Keywords

1 Introduction

This paper presents a systematic approach to the study of models of complex systems.

1.1 What Is a Model?

A simple though not necessarily comprehensive definition of a model is that:

A model is a representation of a postulated relationship among inputs and outputs of a system usually informed by observation and a hypothesis that best explains them.

The definition captures two of the most important characteristics

  • models depend on a hypothesis, and,

  • models use the data from observation to validate and refine the hypothesis.

Errors and uncertainty in the data and limitations in the hypothesis (usually a tractable and computable mathematical construct articulating beliefs like proportionality, linearity, etc.) are immediate challenges that must be overcome to construct useful and credible models.

1.2 Who Needs Them and Why Are There so Many of Them?

A model is most useful in predicting the behavior of a system for unobserved inputs and interpretability or explainability of the system’s behavior. Since, models require a hypotheses implies that the model is a formulation of a belief about the data. The immediate consequence of this that the model may be very poor about such prediction even when sufficient care is taken to use all the available data and information since the subjectivity of the belief can never be completely eliminated. Secondly, the data at hand may not provide enough information about the system to characterize its behavior at the desired prediction. What makes this problem even more acute is that we are often interested in modeling outcomes that are not observed and perhaps sometimes not observable.

The consequence of this lack of knowledge and limited data is the multiplicity of beliefs about the complex system being modeled and a profusion of models based on different modeling assumptions and data use. These competing models lead to much debate among scientists. Principles like “Occam’s razor” and Bayesian statistics [2] provide some guidance but simple robust approaches that allow the testing of models for fitness need to be developed. We present in this paper a simple data driven approach to discriminate among models and the modeling assumptions implicit in each model, given a range of phenomena to be studied. We illustrate the approach by work on granular flow models of large mass flows.

1.3 Models and Assumptions

An assumption is a simple intuitive concept. An assumption is any atomic postulate about relationships among quantities under study, e.g., a linear stress strain relationship \(\sigma =E\epsilon \) or neglecting some quantities in comparison to larger quantities \(\theta \approx \sin (\theta ) \) for small \(\theta \). Models are compositions of many such assumptions. The study of models is thus implicitly a study of these assumptions and their composability and applicability in a particular context. Sometimes a good model contains a useless assumption that may be removed, sometimes a good assumption should be implemented inside a different model - these are usually subjective choices, not data driven. Moreover, the correct assumptions may change through time, making model choice more difficult.

The rest of the paper will define our approach and a simple illustration using 3 models for large scale mass flows incorporated in our large scale mass flow simulation framework TITAN2D [5]. The availability of 3 distinct models for similar phenomena in the same tool provides us the ability to directly compare inputs, outputs and internal variables in all the 3 models.

1.4 Analysis of Modeling Assumptions and Models

Let us define \(\left( M(A), P_{M(A)}\right) \), where A is a set of assumptions, M(A) is the model which combines those assumptions, and \(P_M\) is a probability distribution in the parameter space of M. For the sake of simplicity we assume \(P_M\) to be uniformly distributed on selected parameter ranges. While the support of \(P_M\) can be restricted to a single value by solving an inverse problem for the optimal reconstruction of a particular flow, this is not possible if we are interested in the general predictive capabilities of the model, where we are interested in the outcomes over a whole range.

Stage 1: Parameter Ranges. In this study, we always assume \(P_M\sim \bigotimes _{i=1}^{N_M} Unif(a_{i,M},b_{i,M})\), where \(N_M\) is the number of parameters of M. These parameter ranges will be chosen using information gathered from the literature about the physical meaning of those values together with a preliminary testing for physical consistency of model outcomes and range of inputs/outcomes of interest.

Stage 2: Simulations and Data Gathering. The simulation algorithms can be represented as (Fig. 1):

Fig. 1.
figure 1

Models and variables

The model inputs are the parameters of M, The latent variables include quantities in the model evaluation that are ascribable to specific assumptions \(A_i\). These are usually not observed as outputs from the model. For example in momentum balances of complex flow calculations these could be values of different source terms, dissipation terms and inertia terms. Finally, the model outputs include explicit outcomes, e.g., for flow calculations these could be flow height, lateral extent, area, velocity, acceleration, and derived quantities such as Froude number Fr. In general, for each quantity of interest (QoI), we use a Monte Carlo simulation, sampling the input variables and obtaining a family of graphs plotting their expectation, and their 5\(^{\mathrm {th}}\) and 95\(^{\mathrm {th}}\) percentiles. Our sampling technique of the input variables is based on the Latin Hypercube Sampling (LHS) idea, and in particular, on the improved space-filling properties of the orthogonal array-based Latin Hypercubes.

Stage 3: Results Analysis. These and other statistics can now be compared to determine the need for different modeling assumptions and the relative merits of different models. Thus, analysis of the data gathered over the entire range of flows for the state variables and outcomes leads to a quantitative basis for accepting or rejecting particular assumptions or models for specific outcomes.

2 Modeling of Mass Flows

Dense large scale granular avalanches are a complex class of flows with physics that has often been poorly captured by models that are computationally tractable. Sparsity of actual flow data (usually only a posteriori deposit information is available), and large uncertainty in the mechanisms of initiation and flow propagation, make the modeling task challenging, and a subject of much continuing interest. Models that appear to represent the physics well in certain flows, may turn out to be poorly behaved in others, due to intrinsic mathematical or numerical issues. Nevertheless, given the large implications on life and property, many models with different modeling assumptions have been proposed.

2.1 Three Models

Modeling in this case proceeds by first assuming that the laws of mass and momentum conservation hold for properly defined system boundaries. The scale of these flows, very long and wide with small depth led to the first most generally accepted assumption, shallowness [13]. This allows an integration through the depth to obtain simpler and more computationally tractable equations. This is the next of many assumptions that have to be made. Both of these are fundamental assumptions which can be tested in the procedure we established above. Since, there is a general consensus and much evidence in the literature of the validity of these assumptions we defer analysis of these to future work. The depth-averaged Saint-Venant equations that result are:

$$\begin{aligned} \frac{\partial {h}}{\partial t} + \frac{\partial }{\partial {x}}(h \bar{u}) + \frac{\partial }{\partial {y}}(h\bar{v})= & {} 0 \nonumber \\ \frac{\partial }{\partial {t}} (h\bar{u}) + \frac{\partial }{\partial x}\left( h\bar{u}^2 + \frac{1}{2}k g_{z}h^2\right) + \frac{\partial }{\partial y}(h\bar{u}\bar{v})= & {} S_{x}\\ \frac{\partial }{\partial t} (h\bar{v}) + \frac{\partial }{\partial x}(h\bar{u}\bar{v}) + \frac{\partial }{\partial y}\left( h\bar{v}^2 + \frac{1}{2}k g_{z}h^2\right)= & {} S_{y} \nonumber \end{aligned}$$
(1)

Here the Cartesian coordinate system is aligned such that z is normal to the surface; h is the flow height in the z direction; \(h\bar{u}\) and \(h\bar{v}\) are respectively the components of momentum in the x and y directions; and k is the coefficient which relates the lateral stress components, \(\bar{\sigma }_{xx}\) and \(\bar{\sigma }_{yy}\), to the normal stress component, \(\bar{\sigma }_{zz}\). The definition of this coefficient depends on the constitutive model of the flowing material we choose. Note that \(\frac{1}{2} k g_z h^2\) is the contribution of hydrostatic pressure to the momentum fluxes. \(S_x\) and \(S_y\) are the sum local stresses: they include the gravitational driving forces, the basal friction force resisting to the motion of the material, and additional forces specific of rheology assumptions.

The final class of assumptions are the assumptions on the rheology of the flows – in particular in this context assumptions used to model different dissipation mechanisms embedded in \(S_x, S_y\) that lead to a plethora of models with much controversy on the most suitable model.

Mohr-Coulomb (MC). Based on the long history of studies in soil mechanics [7], the Mohr-Coulomb (MC) rheology model was developed and used to represent the behavior of geophysical mass flows [13].

Shear and normal stress are assumed to obey Coulomb friction equation, both within the flow and at its boundaries. In other words,

$$\begin{aligned} \tau = \sigma \tan \phi , \end{aligned}$$
(2)

where \(\tau \) and \(\sigma \) are respectively the shear and normal stresses on failure surfaces, and \(\phi \) is a friction angle. This relationship does not depend on the flow speed.

We can summarize the MC rheology assumptions as:

  • Basal Friction based on a constant friction angle.

  • Internal Friction based on a constant friction angle.

  • Earth pressure coefficient formula depends on the Mohr circle.

  • The velocity based curvature effects are included into the equations.

Under the assumption of symmetry of the stress tensor with respect to the z axis, the earth pressure coefficient \(k=k_{ap}\) can take on only one of three values \(\{ 0, \pm 1\}\). The material yield criterion is represented by the two straight lines at angles \(\pm \phi \) (the internal friction angle) relative to horizontal direction. Similarly, the normal and shear stress at the bed are represented by the line \(\tau =-\sigma \tan (\delta )\) where \(\delta \) is the bed friction angle.

MC Equations. As a result, we can write down the source terms of the Eq. (1):

$$\begin{aligned} S_x =&g_x h - \frac{\bar{u}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert } \left[ h\left( g_z+\frac{\bar{u}^2}{r_x}\right) \tan (\phi _{bed})\right] - h k_{ap} \ \mathrm{sgn}\left( \frac{\partial \bar{u}}{\partial y}\right) \frac{\partial (g_z h)}{\partial {y}} \sin (\phi _{int}) \nonumber \\ S_y =&g_y h - \frac{\bar{v}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert } \left[ h\left( g_z +\frac{\bar{v}^2}{r_y}\right) \tan (\phi _{bed})\right] - h k_{ap} \ \mathrm{sgn}\left( {\frac{\partial \bar{v}}{\partial x}}\right) \frac{\partial (g_z h)}{\partial {x}} \sin (\phi _{int}) \end{aligned}$$
(3)

Where, \(\underset{^\sim }{\bar{\mathbf{u }}} = (\bar{u} , \bar{v})\), is the depth-averaged velocity vector, \(r_x\) and \(r_y\) denote the radii of curvature of the local basal surface. The inverse of the radii of curvature is usually approximated with the partial derivatives of the basal slope, e.g., \(1/r_x = \partial \theta _x/\partial x\), where \(\theta _x\) is the local bed slope.

Pouliquen-Forterre (PF). The scaling properties for granular flows down rough inclined planes led to a new formulation of the basal friction stress as a function of the flow depth and velocity [6]. PF rheology assumptions can be summarized as:

  • Basal Friction is based on an interpolation of two different friction angles, based on the flow regime and depth.

  • Internal Friction is neglected.

  • Earth pressure coefficient is equal to one.

  • Normal stress is modified by a hydrostatic pressure force related to the flow height gradient.

  • Velocity based curvature effects are included into the equations.

Two critical slope inclination angles are defined as functions of the flow thickness, namely \(\phi _{start}(h)\) and \(\phi _{stop}(h)\). The function \(\phi _{stop}(h)\) gives the slope angle at which a steady uniform flow leaves a deposit of thickness h, while \(\phi _{start}(h)\) is the angle at which a layer of thickness h is mobilized. They define two different basal friction coefficients.

$$\begin{aligned} \mu _{start}(h)=\tan (\phi _{start}(h))\end{aligned}$$
(4)
$$\begin{aligned} \mu _{stop}(h)=\tan (\phi _{stop}(h)) \end{aligned}$$
(5)

An empirical friction law \(\mu _{b}(\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert , h)\) is then defined in the whole range of velocity and thickness.

PF Equations. The depth-averaged Eq. (1) source terms thus take the following form:

$$\begin{aligned} S_{x}= & {} g_{x} h - \frac{\bar{u}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }\left[ h \left( g_z+\frac{\bar{u}^2}{r_x}\right) \ \mu _{b}(\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert , h)\right] \ + g_{z}h\frac{\partial h}{\partial x} \nonumber \\ S_{y}= & {} g_{y} h - \frac{\bar{v}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }\left[ h \left( g_z +\frac{\bar{v}^2}{r_y}\right) \ \mu _{b}(\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert , h)\right] \ + g_{z}h\frac{\partial h}{\partial y} \end{aligned}$$
(6)

Voellmy-Salm (VS). The theoretical analysis of dense snow avalanches led to the VS rheology model [9, 15]. The following relation between shear and normal stresses holds:

$$\begin{aligned} \tau = \mu \sigma + \frac{\rho \Vert \underline{\mathbf{g }} \Vert }{\xi } \Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert ^2, \end{aligned}$$
(7)

where, \(\sigma \) denotes the normal stress at the bottom of the fluid layer and \(\underline{\mathbf{g }} = (g_{x} , g_{y} , g_{z})\) represents the gravity vector. The VS rheology adds a velocity dependent turbulent friction to the traditional velocity independent basal friction term which is proportional to the normal stress at the flow bottom. The two parameters of the model are the bed friction coefficient \(\mu \) and the turbulent friction coefficient \(\xi \).

We can summarize VS rheology assumptions as:

  • Basal Friction is based on a constant coefficient, similarly to the MC rheology.

  • Internal Friction is neglected.

  • Earth pressure coefficient is equal to one.

  • Additional turbulent friction is based on the local velocity by a quadratic expression.

  • Velocity based curvature effects are included into the equations, following an alternative formulation.

The effect of the topographic local curvatures is again taken into account by adding the terms containing the local radii of curvature \(r_x\) and \(r_y\). In this case the formula is considering the modulus of velocity instead than the scalar component [3].

VS Equations. Therefore, the final source terms take the following form:

$$\begin{aligned} S_{x}= & {} g_{x} h - \frac{\bar{u}}{\Vert \underset{^\sim }{\bar{\mathbf {u}}}\Vert } \ \left[ h \left( g_{z} + \frac{\Vert \underset{^\sim }{\bar{\mathbf {u}}} \Vert ^2}{r_{x}} \right) \mu + \frac{\Vert \underset{^\sim }{\mathbf {g}} \Vert }{\xi }\Vert \underset{^\sim }{\bar{\mathbf {u}}} \Vert ^2\right] , \nonumber \\ S_{y}= & {} g_{y} h - \frac{\bar{v}}{\Vert \underset{^\sim }{\bar{\mathbf {u}}}\Vert } \ \left[ h \left( g_{z} + \frac{\Vert \underset{^\sim }{\bar{\mathbf {u}}} \Vert ^2}{r_{y}} \right) \mu + \frac{\Vert \underset{^\sim }{\mathbf {g}} \Vert }{\xi }\Vert \underset{^\sim }{\bar{\mathbf {u}}} \Vert ^2\right] . \end{aligned}$$
(8)

Latent Variables. For analysis of modeling assumptions we need to record and classify the results of different modeling assumptions. These terms are explored in detail in the next sections.

$$\begin{aligned} \varvec{RHS_1} = [g_x h,g_y h], \end{aligned}$$
(9)

it is the gravitational force term, it has the same formulation in all models.

The formula of basal friction force \(\varvec{RHS_2}\) depends on the model:

$$\begin{aligned} \varvec{RHS_2} =&-h g_z\tan (\phi _{bed})\left[ \frac{\bar{u}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }, \frac{\bar{v}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert } \right] ,\text { in MC model.}\nonumber \\ \varvec{RHS_2} =&- h g_z \ \mu _{b}(\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert , h)\left[ \frac{\bar{u}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }, \frac{\bar{v}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }\right] ,\text { in PF model.}\\ \varvec{RHS_2} =&-h g_{z} \mu \left[ \frac{\bar{u}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}}\Vert } , \frac{\bar{v}}{\Vert \underset{^\sim }{\bar{\mathbf{u }}}\Vert }\right] ,\text { in VS model.}\nonumber \end{aligned}$$
(10)

The formula of the force related to the topography curvature, \(\varvec{RHS_3}\), also depends on the model:

$$\begin{aligned} \varvec{RHS_3} =&-h \tan (\phi _{bed})\left[ \frac{\bar{u}^3}{r_x\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }, \frac{\bar{v}^3}{r_y\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }\right] ,\text { in MC model.}\nonumber \\ \varvec{RHS_3} =&-h\ \mu _{b}(\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert ,h)\left[ \frac{\bar{u}^3}{r_x\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }, \frac{\bar{v}^3}{r_y\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }\right] ,\text { in PF model.}\\ \varvec{RHS_3} =&-h \mu \left[ \frac{\bar{u}\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }{r_{x}},\frac{\bar{v}\Vert \underset{^\sim }{\bar{\mathbf{u }}} \Vert }{r_{y}}\right] ,\text { in VS model.}\nonumber \end{aligned}$$
(11)

All the three models have an additional force term, having a different formula and meaning in the three models:

$$\begin{aligned} \varvec{RHS_4} =&- h k_{ap}\sin (\phi _{int})\left[ \ \mathrm{sgn}(\frac{\partial \bar{u}}{\partial y}) \frac{\partial (g_z h)}{\partial y},\ \mathrm{sgn}({\frac{\partial \bar{v}}{\partial x}}) \frac{\partial (g_z h)}{\partial x}\right] ,\text { in MC model.}\nonumber \\ \varvec{RHS_4} =&\, g_{z}h\left[ \frac{\partial h}{\partial x}, \frac{\partial h}{\partial y}\right] ,\text { in PF model.}\\ \varvec{RHS_4} =&-\frac{\Vert \underset{^\sim }{\mathbf {g}} \Vert }{\xi }\Vert \underset{^\sim }{\bar{\mathbf {u}}} \Vert ^2\left[ \frac{\bar{u}}{\Vert \underset{^\sim }{\bar{\mathbf {u}}}\Vert },\frac{\bar{v}}{\Vert \underset{^\sim }{\bar{\mathbf {u}}}\Vert }\right] ,\text { in VS model.}\nonumber \end{aligned}$$
(12)

These latent variables can be analyzed locally and globally for discriminating among the different modeling assumption.

2.2 Monte Carlo Process and Statistical Analysis

For our study, the flow range is defined by establishing boundaries for inputs like flow volume and rheology coefficients. Optionally, these could include also flow initiation site and geometry, and the digital elevation map. The Latin Hypercube Sampling is performed over \([0,1]^3\) for the MC and VS input parameters, and \([0,1]^4\) for PF input parameters. Those dimensionless samples are linearly mapped to fill the required intervals.

Following the simulations, we generate data for each sample run and each outcome and latent variable \(f(\underline{\mathbf{x }},t)\) calculated as a function of time on the elements of the computational grid. This analysis generates tremendous volume of data which must then be analyzed using statistical methods for summative impact. The latent variables in this case are the mass and force terms in the conservation laws defined above.

We devise many statistical measures for analyzing the data. For instance, let \((F_i(x,t))_{i=1,\dots , 4}\) be an array of force components, where \(x\in \mathbf R^2\) is a spatial location, and \(t\in T\) is a time instant. The degree of contribution of those force terms can be significantly variable in space and time, and we define the dominance factors \((p_j)_{j=1,\dots , k}\), i.e., the probability of each \(F_j\) to be the dominant force at (xt). Those probabilities provide insight into the dominance of a particular source or dissipation (identified with a particular modeling assumption) term on the model dynamics.

2.3 Overview of the Case Studies

The first case study assumes very simple boundary conditions, and corresponds to an experiment fully described in [16]. It is a classical flow down an inclined plane set-up, including a change in slope to an horizontal plane (Fig. 2 Left). Four locations are selected among the center line of the flow to accomplish local testing. These are: the initial pile location \(L_1=(-0.7,0)\) m, the middle of the inclined plane \(L_2=(-0.35,0)\) m, the change in slope \(L_3=(0,0)\) m, the middle of the flat plane \(L_4=(0.15,0)\) m.

Fig. 2.
figure 2

[Left] inclined plane description, including local samples sites (red stars). Pile location is marked by a blue dot. [Right] (a) Volcán de Colima (México) overview, including 51 numbered local sample sites (stars) and four labeled major ravines channeling the flow. Pile location is marked by a blue dot. Reported coordinates are in UTM zone 13N. Background is a satellite photo. Six points that are adopted as preferred locations are highlighted in yellow. (Color figure online)

The second case study is a block and ash flow down the slope of Volcán de Colima (MX) - an andesitic stratovolcano that rises to 3,860 m above sea level, situated in the western portion of the Trans-Mexican Volcanic Belt (Fig. 2 Right). The modeling of pyroclastic flows generated by explosive eruptions and lava dome collapses of Volcán de Colima is a well studied problem [4, 10,11,12, 14]. The volcano has been already used as a case study in several studies involving the Titan2D code [8]. We select 51 locations along the flow inundated area to observe model outputs with six of them as preferred locations being representative of different flow regimes.

3 Sample Results

Figure 3 shows the flow height, h(Lt), at the points \((L_i)_{i=1,\dots ,4}\), for the three rheology models. Parameter ranges – outcome of Stage 1 analysis – come from literature and past work in our laboratory. Plot 3 clearly shows the differences in the statistics of the flow outcomes induced by the different choices of rheology at different locations in the plane. Availability of data allows us to subject the data to tests of reasonability both for the means and extremal values. Given a particular type of flow and collected data we can clearly distinguish model skill in capturing not only that flow but also possible flows. Past work [16] allows us to conclude that MC rheology is adequate for modeling simple dry granular flows.

Fig. 3.
figure 3

Records of flow height at four spatial locations of interest. Bold line is mean value, dashed/dotted lines are 5\(^{\mathrm {th}}\) and 95\(^{\mathrm {th}}\) percentile bounds. Different rheology models are displayed with different colors. Plots are at different scale, for simplification. (Color figure online)

While, the above analysis is interesting in helping us accept or reject particular models a lot of insight can be obtained by examining the behavior of latent variables. Figure 4 shows the spatial average of speed and Froude Number, for the three rheology models for flows at Volcan Colima. Ranges of parameters etc. are obtained from our past work at this site [1]. It also shows the inundated area of flow, as a function of time. Similar analysis of model suitability can be conducted here given recorded deposits. In past work [5], we have tuned MC rheology to match deposits for known block and ash flows but a priori predictive ability was limited by inability to tune without knowledge of flow character.

Fig. 4.
figure 4

Comparison between spatial averages of (a) flow speed, and (b) Froude Number, in addition to the (c) inundated area, as a function of time.

The plots 5a, b, c and 5d, e, f are related to point \(L_8\) and \(L_{10}\), respectively. They are significantly similar. \(\varvec{RHS_1}\) related to the gravitational force is the dominant force with a very high chance, \(P_1>90\%\). In MC and PF there is a small probability, i.e., \(P_3\) = 5%–30% at most, of \(\varvec{RHS_3}\) related to topographic curvature effects being the dominant force for a short amount of time, i.e. \({\sim }5\,\text {s}\). This occurs in the middle of the time interval in which the flow is almost surely inundating the points being observed. In VS it is observed a \(P_4=5\%\) chance of \(\varvec{RHS_4}\) related to the turbulent dissipation being dominant, for a few seconds, anticipating the minimum of no-flow probability. Plots 5g, h, i, are related to point \(L_{17}\), and the plots are split in two sub-frames, following different temporal scales. In all the models, \(\varvec{RHS_2}\) is the most probable dominant force, and its dominance factor has a bell-shaped profile, similar to the complementary of no-flow probability. In all the models, \(\varvec{RHS_1}\) has a small chance of being the dominant force. In MC, this is more significant, at most \(P_1=30\%\), for \({\sim }20\,\text {s}\) after the flow arrival, and has again about \(P_1=\%2\) chance to be dominant in \([100, 7200]\,\text {s}\). In PF, the chance is \(P_1=15\%\) at most, and has two maxima, one short lasting at about \(55\,\text {s}\), and the second in \([100,500]\,\text {s}\). Also in VS, the chance is at most \(P_1=15\%\), reached at \([300, 500]\,\text {s}\), but its profile is unimodal in time, and becomes lower than \(P_1=2\%\) after \(2000\,\text {s}\). In MC and PF, \(\varvec{RHS_3}\) has a chance of \(P_3=10\%\) of being the dominant force, for a short amount of time \([30, 50]\,\text {s}\) and \([40, 50]\,\text {s}\), respectively. Figure 5 show the Dominance Factors \((P_i)_{i=1,\dots ,4}\), for the three rheology models and focusing on the RHS terms moduli, at the three selected points \(L_{8}\), \(L_{10}\), and \(L_{17}\), closer than 1 Km to the initial pile (in horizontal projection).

Fig. 5.
figure 5

Records of dominance probabilities of RHS force moduli, at three spatial locations of interest, in the first km of runout. Bold line is mean value, dashed/dotted lines are 5\(^{\mathrm {th}}\) and 95\(^{\mathrm {th}}\) percentile bounds. No-flow probability is also displayed. (Color figure online)

4 Conclusions

In this study, we have introduced a simple, robust statistically driven method for analyzing complex models. We have used 3 different models arising from different rheology assumptions. The data shows unambiguously the performance of the models across a wide range of possible flow regimes and topographies. We analyze local and global quantities and latent variables. The analysis of latent variables is particularly illustrative of the impact of modeling assumption. Knowledge of which assumptions dominate, and, by how much, at the level of assumptions will allow us to construct efficient models for desired inputs. Such model composition is the subject of ongoing and future work.