1. Introduction

Among debris flows, the most devastating phenomena are volcanic flows known as lahars. The name comes from Indonesia and describes flows whose solid phase mostly consist of material of volcanic origin (^{Tilling, 1996}). Lahars are regarded as the second largest destructive volcanic hazard (^{Baxter, 1983}; ^{Blong, 1984}; ^{Tilling, 1989}). During the past century, tens of thousands of people have been killed by volcanic flows and hundreds of thousands forced from their homes (^{Tilling, 1996}; ^{National Research Council, 1991}, ^{1994}). These two-phase mass flows containing water and solid particles are common in volcanic regions. They can be initiated by several mechanisms. For example, a volcanic explosion can be accompanied by large plumes and pyroclastic flows consisting of rock and gas that race along the surface of the mountain at speeds as high as 100 m/s (^{Sheridan, 1979}). The hot ash can melt snow, creating a muddy mixture that knocks down trees and entrains rocks and boulders into the flow. Cotopaxi volcano in Ecuador is an example of a volcano that has produced many large lahars by this process (^{Pistolesi et al., 2013}). Crater lakes on volcanoes can be another source of lahars; a recent example is the 2007 lahar of Ruapehu in New Zealand (^{Procter et al., 2010}). An additional mechanism for initiating lahars is intense rainfall on hillsides that are devoid of vegetation with material like clay soils or volcanic ash exposed. An example of this type of lahar is the 1998 mudflow at Casita Volcano in Nicaragua that occurred during Hurricane Mitch and caused hundreds of deaths (^{Sheridan et al., 1999}). Lahars can carry constituent particles that typically range from clay to boulder size and can propagate tens of kilometers before coming to rest (^{Procter et al., 2010}). As the solid particle component becomes deposited, the resulting deposits can be up to 100 m thick (^{Legros, 2002}). However, the typical deposits left after a debris flow passes are on the scale of meters.

Models that assume one phase, like pure fluid (^{Chow, 1959}) or pure frictional models (^{Savage and Hutter, 1989}), cannot represent the complex behavior of these kinds of flows (^{Iverson et al., 2010}; ^{Iverson, 2014}). In order to develop a complete mathematical model of debris flows, two principal challenges must be overcome: rheology and scale. First, relations among granular material must be developed to describe granular material including clays, sands, pebbles and rocks, with interstitial water. Second, a computational method must be developed that extends over six orders of magnitude, as clay diameters are of the order of *O*(10^{-6}) m and boulders are *O*(1) m. Neither of these challenges can be fully met at this time. This paper tries to strike a balance between fidelity to the physics of mass flows and computational feasibility. We describe a modeling effort that draws on innovation from engineering and geoscience, which postulate constitutive theory and fluid-solid interaction effects, and, through a depth averaging process, results in a system of equations that is computationally tractable.

The modeling effort here has its origins in the pioneering work of ^{Savage and Hutter (1989)}. They began with mass and momentum balance laws based on a Coulomb constitutive description of dry granular material. By scaling and depth averaging, they developed a “thin layer” model for granular flows down inclines. Flow over general topography was addressed in ^{Gray et al. (1999)}, ^{Patra et al. (2005)}, and ^{Pudasaini and Hutter (2003)}. Comparison of thin-layer model results to historic flows was presented for example in ^{Sheridan et al. (2005)}, ^{Charbonnier and Gertisser (2009)}, and ^{Sulpizio et al. (2010)}. In ^{Hutter et al. (2005)}, the appropriateness of these thin layer models was considered for several different types of geophysical flows. Much of the modeling effort in this direction was summarized in ^{Pudasaini and Hutter (2007)}.

^{Iverson (1997)} and ^{Iverson and Denlinger (2001)} argued that the presence of interstitial fluid fundamentally alters the behavior of geophysical flows, and fluid effects should be included in the constitutive behavior of the flowing material. Starting with equations of mixture theory (^{Bedford and Drumheller, 1983}) and through a careful examination of experiments, past studies developed a system of mass and momentum balance laws for the mixture. Unfortunately, in this development a balance equation for the motion of pore fluid was not readily available. Instead, based on experimental data, a transport equation for the fluid was postulated.

A different approach, based on a fully three-dimensional model of two-phase flows, can be found in ^{Dartevelle (2004)} and ^{Meruane et al. (2010)}. Another approach to modeling mud flows employs visco-plastic constitutive assumptions (^{Coussot, 1997}; ^{Balmforth and Craster, 1999}; ^{Mei et al., 2001}; ^{Ancey, 2007}). The choice of a visco-plastic flow model drives the subsequent derivation, as well as the parameter fitting necessary for the constitutive relationship. The process of depth averaging a visco-plastic flow is always difficult. The interface between yielding and non-yielding material is itself a free surface that must be determined. This attribute requires the use of multiple layers in the model system, with all the resulting complexity.

^{Pitman and Le (2005)} developed a two-phase thin-layer model of fluid and granular material. They began with a fully three-dimensional model of two-phase flows, based on model equations in engineering (^{Jackson, 2000}). The model equations are scaled and depth averaged. The resulting system of equations is not complete, and closure assumptions are required. With these assumptions, the mathematical system is shown to be hyperbolic under common conditions, and thus well posed (^{Pelanti et al., 2008}). The model of ^{Pitman and Le (2005)} includes a drag term, which is the only term describing the interaction of the two phases; this coefficient must then be fitted to experiments. That model assumes the fluid is not viscous, and that there is no frictional dissipation in the fluid phase at the basal surface. Both of these assumptions, which are reasonable in bench-scale fluidized bed experiments, are difficult to determine for large mass flows. The ^{Pitman and Le (2005)} model over-simplifies the physics of the fluid phase. When applied to real-scale topographies it only works for constant solid-particle concentrations well above the dilute flow limit, where the physics of two phase flows are mainly governed by pure fluid dynamics. Particle-particle interactions are unimportant for volumetric particle concentrations less than 10% (^{Bagnold, 1962}; ^{Winterwerp, 1999}; ^{Dartevelle, 2003}). In order to address some of the shortcomings of the ^{Pitman and Le (2005)} model, this paper reconsiders the model equations proposed by them and proposes a revision of the model equations that better represents two-phase geophysical flows. For example, the proposed model accounts for the friction at the wall of the fluid phase and no longer assumes a constant volumetric fraction of solids, as in ^{Valentine and Wohletz (1989)} and ^{Dobran (1991)}.

2. Model derivation

The Titan2F model uses a similar framework to that developed in ^{Pitman and Le (2005)}. However, a complete set of model equations for a granular phase and for a fluid phase is written. Phase interaction terms are modeled, and scaling of all terms suggests simplifications that can be made. Depth averaging and closure assumptions complete the derivation.

A note on sign convention: in soil mechanics it is common to consider compressive stresses as positive; by contrast, in fluid mechanics, as an increase in pressure results in a reduction of the volume, compression is negative. We caution the reader to observe the sign convention in the equations below.

2.1 Fundamental assumptions

All symbols used in this paper are described in Table 1, Appendix A. The fundamental theory of two-phase flows used here can be found in ^{Dobran (1991)} and ^{Jackson (2000)}. In two-space dimensions we consider a thin layer of granular material (s) and interstitial fluid (f), each of constant specific density *p*
^{
s
} and *p*
^{
f
} , respectively, flowing over a smooth basal surface, *b*. Erosion and deposition are neglected. Along the basal surface, we define a Cartesian coordinate system *Oxyz*, with origin *O* defined so the *Oxy* is tangent to the basal surface, with *x* in the downstream direction, and *Oz* in the normal direction. Writing *φ*
^{
f
} for the fluid volume fraction. We assume the mass is fully saturated, so the sum of the solid and fluid volume fractions adds to one (*φ*
^{
f
}
*= 1 - φ*). When writing equations in component form, we use subscripts to denote the component of the vectors, and superscripts the phase of the flow (either solids or fluid).

Mass conservation for the two constituent phases may be written as in ^{Anderson and Jackson (1967)}:

The momentum equations are:

Here *T*
^{
s
} and *T*
^{
f
} are the stress tensors for the solid and fluid, respectively, and *f*
^{
s
} and *f*
^{
f
} are the interaction forces between the solid and fluid phases, respectively. We must postulate constitutive relations and an equation for the interphase force to close the system. ^{Jackson (2000)} presented an argument for separating buoyancy from other interphase force terms (such as drag or virtual mass), and for properly accounting for buoyancy in a field with fluid pressure variations. Similar modeling can be found in ^{Valentine and Wohletz (1989)}, ^{Dobran et al. (1993)}, ^{Dobran (2001)}, and ^{Neri et al. (2003}, ^{2007)}.

^{Pitman and Le (2005)} accounted only for the drag in evaluating the interaction force, unlike them, from ^{Dobran (1991)}, (neglecting capillarity, virtual mass, and lift) we account for the total fluid stress as well:

Here, the total fluid stress is *T*
^{
f
}
*= -P*
^{
f
}
*+τ*
^{
f f
} , where *P*
^{
f
} is the fluid pressure and *τ*
^{
f
} is the viscous contribution to the fluid stress. The drag term exchanges momentum between the phases, with a coefficient *D* that is phenomenological. ^{Ergun (1952)}, ^{Wallis (1969)}, ^{Gidaspow (1994)}, ^{Fan and Zhu (1998)}, ^{Dobran (2001)}, ^{Mazzei and Lettieri (2007)} and ^{Panneerselvam et al. (2007)}, among other sources, suggested values. When ∀ *0*, the drag vanishes. Following ^{Mazzei and Lettieri (2007)}, we set

where *d* is the mean particle diameter and β is a constant related to the constant *n* in the Richardson-Zaki equation (^{Khan and Richardson, 1989}). According to ^{Mazzei and Lettieri (2007)}, this constant equals 2.80 either when R∀0 or ∀∞, thus we use β = 2.80 in Equation 6. Finally, by assuming smooth spherical particles on the inertial regime, the drag coefficient is approached as constant C_{d} = 1, which holds for particle Reynolds number up to 500 (^{Sparks et al., 1997}).

2.2 Scaling

The characteristic thickness of the flowing granular material is *H* and the characteristic length is *L*. We scale *x* and *y* by *L*, and *z* by *H*, the time by the free fall time √*Lg*. And we scale the *x*, *y* and *z* velocities by √*Lg* and (*H/L*)√(*Lg*), respectively. The stresses are scaled by *p*
^{
s
}
*gH* for the solids phase and *p*
^{
f
}
*gH* for the fluid phase. After scaling, the mass balance equations are unchanged. Several terms in the momentum equations contain the factor ε = *H/L*, which is small; values of ε from 0.01 to 0.001 are not uncommon (^{Iverson and Denlinger, 2001}). Writing *x*, *y*, *z* for x_{1}, x_{2}, x_{3}, the solid momentum balance equations become (showing here just the x momentum),

Note that components of gravity have been scaled by the magnitude |*g*|, so (g_{x}, g_{y}, g_{z}) is a unit vector.

With the same scaling, the fluid momentum balance equations become

Equation 8 shows an important advance in modeling the physics of the flow related to ^{Pitman and Le (2005)} as we do not neglect the effect of the volumetric fraction of fluids.

In summary, the proposed equation system consists of the solid volume fraction φ, the three solid velocities

2.3 Constitutive assumptions and boundary conditions

The upper surface of the flowing mass at F_{h}(x, y, t) = 0 is assumed to be a material surface and stress free. At the base of the mass, material is assumed to flow tangent to the basal surface, and to satisfy a sliding friction law. For the solid constituent, this friction relation specifies that the shear traction and the normal stress are proportional: _{
bed
} is the basal friction angle and the *−sgn(v)* specifies that the shear traction opposes motion.

We now will discuss constitutive relations. A Coulomb constitutive relation (^{Coulomb, 1776}) is postulated for the material. The Coulomb law is a nonlinear relation among the components of **T**
^{s}, and stipulated that material deforms when the total stress reaches yield, which means *dev*(**T**
^{s}) = **T**
^{s}− ½ *tr* (**T**
^{s})**I** is the stress deviator, *tr* (**T**
^{s}) is the trace of the stress (the sum of the diagonal components), **I** is the identity tensor, and κ is a material parameter, and that as deformation occurs, the stress and strain-rate tensors are aligned. That is, the alignment condition specifies dev(T ^{s}) = λ dev(V), where the strain-rate V = −1/2 *tr* denotes the transpose. To avoid a switching between plastic and non-plastic behavior, we assume the solid material is everywhere in plastic yield.

The full Coulomb relations are too complex to be used here. Two simplifications are proposed. First, at the basal surface, the boundary condition ensures proportionality and alignment of the tangential and normal forces. We assume the same proportionality and alignment holds throughout the thin flowing layer of material. Written in components *ij*, this implies *Tij*
^{
s
} = *ν*
_{
ij
}
*T*
_{
zz
}
^{
s
} , where the proportionality constant ν is a function of φ_{
bed
} . Second, following ^{Rankine (1857)} and ^{Terzaghi (1936)}, an earth pressure relation is assumed for the diagonal stress components, ν_{ii} = k_{ap}, or T_{xx}
^{s} = k_{ap}T_{zz}
^{s}, where

Here **
φ
**

_{int}is the internal friction angle and the choice of the plus or minus sign depends on whether flow is locally contracting (the passive state, with

For the fluid, the stress terms in Equation 8 should be such that in case of **
φ
**

_{s}→ 0 the equations agree with the pure fluid depth-averaged equations.

For pure water, the shear stress at wall can be approached by (^{Chow, 1959})

where *S*
_{
f
} is the slope friction and *R*
_{
h
} is the hydraulic ratio. Note that for shallow-water problems, this assumption could be considered a rough approximation in case the shallow-water theory’s condition is not met. There are several approaches for approximating both the slope friction and the shear stress at wall. For example, ^{Pan et al. (2006)} use the empirical Manning approach, whereas ^{Liu and Leendertse (1978)} and ^{Zhou and Stansby (1999)} use the Chezy equation. Both the Manning and Chezy approaches pose numerical problems when *h* ∀ 0. Thus we use the Darcy equation (^{Zhou, 1995}; ^{Xu, 2006}):

where

is the Chezy coefficient. The friction coefficient depends on the Reynolds number and the roughness of the channel (k_{s}).

3. Depth averaging

The final step in the derivation is a depth averaging of the mass and momentum balance equations. In this and the following sections we will use the same notation used by ^{Savage and Hutter (1989)}. If *h*(*x, y, t*) is the unsteady surface of the flow and *b*(*x, y*) is the terrain surface, for some function *f*, we compute

Where *h − b* is the flow depth at a point (*x, y*) and time *t*. As the function *f* contains partial derivatives, repeated use of the Leibniz rule is made to interchange integration and differentiation, and boundary conditions are employed to evaluate terms at *b* and *h*. In addition, several approximations must be made during the depth averaging process. In what follows, we only briefly sketch the depth averaging process, noting as appropriate those places where approximations are made. ^{Pitman and Le (2005)} provided an estimation of the errors typically made by these assumptions.

The terms of order ε are assumed small and we hope to drop all such terms from the model. However, ^{Savage and Hutter (1989)} argued that diagonal contributions to the solid stress must be retained. Because there is no preferential downslope direction in our applications, and the flow direction may change during a flow, we retain the stress terms in both the *x* and *y* directions, dropping only *O*(ε) terms in the *z* direction. See the discussion in ^{Iverson and Denlinger (2001)}.

3.1 Mass balance equation

As ρ^{s} and ρ^{f} are constants, Equation 1 can be reduced to:

Which says that the volume-weighted mixture flow is divergence free. This equation is integrated from *z* = *b* to *z* = *h*:

The upper free surface *F*
_{
h
} = 0 is a material surface for the mixture. So, following ^{Savage and Hutter (1989)} and ^{Pitman and Le (2005)}, at *z* = *h*(*x, y, t*),

At the basal surface *F*
_{
b
} = *0*, the flow is tangent to the fixed bed, and the bed is fixed in time. Thus, we can drop in Equation 12 the terms in ∂_{t}, taking into account that at the surface *z* = *b*(*x, y*) (^{Pitman and Le, 2005}):

Like in ^{Pitman and Le (2005)}, at this point, we have ignored sedimentation, entrainment, and infiltration of fluid into the bed.

Using Equation 13 and after algebraic manipulation, the depth averaged equation for the total mass of the solid and fluid can be written

In writing this equation, the depth averaged velocities are

with a similar expression for the volume fraction of solids ^{Savage and Hutter (1989)}, *h* ̀= *h - b*.

An additional equation is needed to solve for **
φ
** . After depth averaging Equation 1, we arrive at the mass balance equations for the solid and fluid phases:

Where

**3.2 z momentum**

Note that as *z* momentum equation, the fluid tends to be hydrostatic:

Depth averaging Equation 17:

In the same manner, for the solid *z* momentum we find the equation for an effective stress:

Substituting Equation 17 in 18 and solving for ∂_{z}T ^{s}
_{zz},

Thus the normal solid stress in the *z* direction at any height is equal to the reduced gravity times the volumetric fraction of solids.

In scaling these equations, the *z* velocities have been dropped. Of course neglecting motion in the *z* direction is a central component of a thin layer theory. Furthermore, any contribution to the *z* momentum from fluid shearing -terms such as T ^{f}
_{xz}, T ^{f}
_{yz}- are dropped due to scaling. Thus, only pressure contributions to the fluid stress are important, which is an assumption we will make below, albeit with a modification at the basal surface.

**3.3 x and y momentum**

The nonlinearity of these equations presents difficulties in formulating a depth averaged theory, and complicates the derivation, and in several places it is necessary to take the depth average of products of terms. When necessary, we approximate the required closure relation; for example, as ^{Savage and Hutter, 1989}; ^{Pitman and Le, 2005}).

Considering first the equation for the motion of the solid phase, the left-hand side of the solids *x* momentum Equation 7 can be written

Depth averaging and using boundary conditions, we find

Now the depth average of the right hand side of Equation 7 becomes:

In order to proceed, the following assumptions are made:

This equation governs the motion of the solid phase and we assume the upper free surface for the mixture is a free surface for both of the individual phases.

The drag term

*D*is highly nonlinear and a correct depth average is all but impossible to calculate. We postulate that experiments could fit an averaged phenomenological drag of a similar form. As seen in Equation 6, we assume a drag coefficient*C*_{ d }= 1. In addition, as the range of particle sizes of lahars is huge, at present, tracking all of the particles or even a representative sample of them, is intractable. Nevertheless, the lahar mean class size is about coarse sand (^{Pareschi, 1996}), thus we use a typical mean particle diameter for lahars*d*= 1 mm (^{Schmid, 1981}).The earth pressure relation for the solid phase is employed. That is, the basal shear stresses are assumed to be proportional to the normal stress

Where *i* can be either *x* or *y*, and the velocity ratio enforces that friction opposes motion in the designated direction (^{Savage and Hutter, 1989}; ^{Patra et al., 2005}). The α notation will provide a convenient shorthand that we use in other places. Likewise the diagonal stresses are taken to be proportional to the normal solid stress

Finally, following ^{Iverson and Denlinger (2001)}, *xy* shear stresses are determined by a Coulomb relation

Where the sgn function ensures that friction opposes straining in the (*x, y*) plane.

• For the fluid phase, the basal shear stresses are assumed to be proportional to the square of the depth averaged velocities (^{Zhuo, 1995}; ^{Xu, 2006})

Where *C*
_{
f
} is the Chezy coefficient, which depends on the friction coefficient (see Equation 9). A physical approach for this coefficient is the Colebrook-White equation (^{Colebrook and White, 1937}), which for rough channels can be approximated as

Where *k*
_{
s
} is the roughness of the channel and *R*
_{
h
} is the hydraulic radius, which for shallow-water problems can be approached as the depth of the flow (R_{h}→*h*). Equation 22 is logarithmic, thus large uncertainties in *k*
_{
s
} result only in small variations in *C*
_{
f
} (^{Swaffield and Bridge, 1983}). In the ^{Transport & Road Research Laboratory (1976)} guide, values of *k*
_{
s
} are proposed for different materials and channel types. From that guide, we choose *k*
_{
s
} =1 mm for channels in volcanic environments (as Titan2F is open source software, the user can modify this value). Therefore, here *C*
_{
f
} will depend on the flow depth *ρ*
^{
f
} , and the fluid velocity *u*
_{
i
}
*.*Note that this is a physical approach for T ^{f}
_{iz}, which does not depend on empirical approaches.

From Leibniz’s rule and the stress computations above, we find

Now using the fluid and solid stress relation

term (i) is approximated as

The upper free surface is stress free, so all terms involving

Note that the factor ( ‒ *g*
_{
z
} ) originates in the evaluation and depth averaging of the fluid stress; in typical flows, this factor is positive.

Combining all terms yields a solids *x* momentum equation

The *y* solid momentum equation can be derived in a similar fashion.

For pure fluids, the diagonal stresses and shear stress are zero. Thus, depth averaging the equation for the fluid motion presents fewer difficulties.

The depth averaged *x* momentum equation takes the form

Where *y* momentum equation has a similar form. Note that unlike ^{Pitman and Le (2005)}, if ^{Chow, 1959}). ^{Kowalski (2008)} describes how debris flows become reduced to a shallow-water flow as solid volume fraction vanishes.

As noted above, we solve for volume fractions (**
φ
** ), thus the bulk density can be calculated from

Then, we obtain the dynamic pressure *p* from

Where ^{Fan and Zhu, 1998})

where

The use of the impacting dynamic pressure information on structures and living beings allows us to estimate levels of damage, as in ^{Valentine (1998)}, ^{Valentine et al. (2011)}, and ^{Jones (2013)}, which is useful for vulnerability analysis.

4. Numerical solution

This system of equations (14, 15, 16, 26, 27) is then solved using the finite volume method, whose solution provides results of the velocity field, flow depth, and the volumetric concentration of solids at the center of each finite volume computational grid.

To solve the balance laws, we use the Godunov solver developed by ^{Davis (1988)} already implemented in ^{Patra et al. (2005)} and ^{Pitman and Le (2005)}. The adaptive meshing is used as well, which allows us to have very fine grids where indicators show high gradients, and coarser grids where low gradients are detected. The time step is adjusted from the Courant condition (^{Courant et al., 1967)}. The complexity of the equation system results in typical time steps of the order of 10^{-4} s. As consequence of this small time step, Titan2F becomes a computationally expensive tool.

The numerical solution of the above set of equations presents strong numerical sensitivity to Digital Elevation Model (DEM) errors and the quality of those maps. The DEMscan have regions where elevations are not well defined; they can have crossed contour levels or even infinite holes (^{Wechsler, 2007}; ^{Dalbey et al., 2008}). Abrupt terrain changes, both actual or DEM artifacts, cause computations of gradients and curvatures to become unstable. In order to avoid such numerical problems, patching and intelligent smoothing of the DEMs are needed.

Based on the hyperbolicity analysis done by ^{Pitman and Le (2005)}, we try to ensure it , imposing a minimum _{
max
}
*= φ*
_{
pc
} corresponding to a maximum packing concentration of 0.65, and a minimum φ_{
min
} = 10^{-8} that ensured stability. This constraint makes the program stable if a smooth DEM is used.

The needed initial conditions are the location of the pile of material, its geometry, the volumetric solid concentration, and initial velocity. There is no inflow condition implemented yet. Inflow hydrograph can be approached using several piles, each one with different initial depths. The bed and internal friction are set internally to the fixed values of *φ*
_{
bed
} = 40° and *φ*
_{int}= 42° respectively. Those values correspond to the bed and internal friction of dry smooth spheres moving on the assumed bed roughness size *k*
_{
s
} (see ^{Kirchner et al., 1990}; ^{Miller and Byne, 1966}; ^{Webb, 2004}). Nevertheless, both of these parameters can affect *k*
_{
ap
} . The doted line in Figure 1 shows how *k*
_{
ap
} changes with _{
bed
} for any fixed value of *φ*
_{
s
} (Note that ^{Williams et al. (2008)} shows that Titan2D results are not strongly affected by *k*
_{
ap
} , see Equation 9). _{
bed
} is a function of *h*, which is directly related to φ_{
s
} , as shown by ^{Kowalski (2008)}. The black line in Figure 1 shows how *k*
_{
ap
} evolves when the fixed value of _{
bed
} , is multiplied by 0 < *φ* < 0.65. Thus, in the evaluation of *k*
_{
ap
} we use the result of the multiplication *φ* × _{
bed
} and

instead of the user defined _{
bed
} and

used by ^{Pitman and Le (2005)}. In this way, the range of resulting values of *k*
_{
ap
} are within the ranges used for modeling lahars with Titan2D (^{Sheridan et al., 2005}; ^{Williams et al., 2008}; ^{Procter et al., 2010}, ^{2012}).

5. Differences with Pitman and Le model

There are six major differences between the present paper and ^{Pitman and Le (2005)}:

In

^{Pitman and Le (2005)}, mass and momentum conservation laws are derived for the solid material and for the phase averaged mixture of solid and fluid, whereas here, the Titan2F model uses mass and momentum equations for both individual phases.Any two-phase model system must incorporate several phenomenological functions, such as an interphase drag coefficient. In the present derivation these functions are better adapted to geophysical flows whose fluid phase corresponds to water and the solid phase are rounded solid particles. The drag is calculated from an expression valid for the whole range of Reynolds numbers. It only needs the mean particle diameter of the flow as a parameter. We use a typical mean particle diameter of lahars. As Titan2F is open source, the user can modify the assumed value of the mean particle diameter.

The volumetric particle concentration is no longer a fixed parameter, which in our approach is calculated for every time step and grid point. This means that instead of having constant 𐐐 like in Equations (3.2) in

^{Pitman and Le (2005)}or simply not accounted for, like in the depth averaged fluid mass and momentum balance equations (3.27) and (3.28) in^{Pitman and Le (2005)}, we include the variable 𐐐 within all derivatives and depth averaged equations. Thus, in regions where the particle concentration vanishes, the solid phase role in the equation system vanishes as well. In that way, the equation system becomes the typical hydraulic shallow-water approach, which does not happen in^{Pitman and Le (2005)}.We account for the fluid stresses at wall. We use a physical approach that only needs the roughness of the channel and the flow depth. The last term is calculated by the program at every time step, whereas the former is set as a fixed value.

The interaction forces between the phases now depend on both the drag and the fluid stress. The pressure is no longer the only term in the fluid stress model. Instead, we use the physical Darcy-Weisbach hydraulic approach that accounts for the full fluid stress. This approach allows us to have a more realistic model for the fluid phase than in

^{Pitman and Le (2005)}.The only input parameters needed by the program are the location of the pile of material, its volume, and the volumetric solids concentration. The friction coefficients are no longer needed as input as they are automatically adjusted according to the evolution of the volumetric fraction of solids across the grid and time. The bed and internal friction are set in such a way that when the volumetric fraction of solids tends to an assumed maximum packing concentration (

*φ*^{Sheridan et al. (2005)}and^{Williams et al. (2008)}for ranges of values.

6. Validation and verification

Validation of the accuracy of the code was done with analytical solutions of the Dam Break problem and with several experimental results (*e.g.*^{Dressler, 1954}; ^{Ancey et al., 2008}). Among them, we check the deposited pattern predicted by the program with the results shown by ^{Liu (1996)}. The prediction of the arrival time and the flow-depth profile was compared with the experimental results shown by ^{Iverson et al. (2010)} from his recent work done on his large-channel facility.

Analytical solutions for shallow-water problems are scarce. Only one-dimensional analytical approaches are available in the literature, especially for the well known Dam Break problem (*e.g.*, ^{Dressler, 1954}; ^{Mangeney et al., 2000}; ^{Fernandez-Feria, 2006}; ^{Ancey et al., 2008}). Unfortunately, analytical solutions for geo-mass flows are almost impossible to find due to the complexity of the nonlinear partial differential equations (^{Pudasaini and Hutter, 2007}). Such solutions can be obtained only in special cases like the similarity based solutions proposed by ^{Savage and Hutter (1989)} for dry avalanches. In our test we use the solution proposed by ^{Fernandez-Feria (2006)} for the Dam Break problem on an incline for pure water. In our program we assume *t*-test shows that there is no statistical difference between the analytical solutions and the prediction of Titan2F using the 1D version of the equations. That test shows that about 70% of the predictions of Titan2F can explain the analytical solution. At least for the one dimensional case, the program successfully reproduces analytical solutions for different initial conditions down to very low particle concentrations (less than 1%). However, as one dimensional approaches neglect lateral spreads, and at the initial parts of the flow evolution *H ⁄ L* is not small, that version of the program tends to over-estimate the front advance of the flow, as can be seen in Figure 2.

^{Liu (1996)} performed several experiments for geo-mass flows in an inclined channel. He modified the initial volume, the channel slope, and the particle concentration to find the final size of the debris flows measured by their resulting maximum width and length. We reproduced the experimental final width and length after the simulation reached the same time corresponding to the duration reported by ^{Liu (1996)}. Figure 3 shows the correspondence of the model with the experiments for (a) the width of the deposit and (b) the length. A Pearson correlation shows that 90% of the experimental data for the deposit length can be explained with the predictions of Titan2F, whereas 80% of the data for the deposit width can be explained by the predictions of the program. This illustrates the high accuracy of the program in predicting the deposit characteristics for different initial volumes (ranging from 2.7 m^{3} to 16 m^{3}) and high initial solid concentrations (**
φ
** = 0.53 ‒ 0.65). See

^{Liu (1996)}for details about initial conditions of his experiments.

The experiments performed by ^{Iverson et al. (2010)} done on a 95 m long artificial channel were used to verify the accuracy of the predictions of the flow front arrival time and the temporal evolution of the flow depth. These flows were unsteady and nonuniform. ^{Iverson et al. (2010)} show time-series data for several measured properties: flow thickness, pore pressure, basal normal stress, and arrival time of the front. Raw data provided to us was used to test the Titan2F prediction concerning time evolution of the flow depth and arrival times at the check points located at 32 and 66 m distance from the lock. As shown in Figures 4 and 5, the arrival time and the temporal evolution of the predicted depth fits well within the confidence interval of the experiments. In both of the cases, Titan2F tends to over-estimate the flow depth just after the arrival of the front. This is probably an effect of the slight difference in the shape of the initial pile, as the free surface of the numerical pile follows the same slope of the channel, whereas the actual free surface within the lock is horizontal. Nevertheless, a Pearson correlation shows that more than 90% of this experimental data can be explained with the predictions of Titan2F.

The range of concentrations that the program can cope with are from = 10^{-8} (almost pure water) to = 0.65 (maximum packing concentration). Finally, as expected, the program predicts high particle concentrations at the front of the flow and low particle concentrations at the tail of the flow (in some cases, even near pure water concentrations, or φ → 0), as can be seen in Figure 6 where a longitudinal solids particle distribution predicted by Titan2F is shown. The predictions fit with the qualitative observation of ^{Iverson et al. (2010)} that the tail of the flow remains very watery. Using a predicted concentration of solids, the density is assessed (Equation 28), and, together with the speed of the flow, the dynamic pressure distribution is calculated as well (Equation 29). Despite that ^{Iverson et al. (2010)} do not show information about the dynamic pressure field, Titan2F does. Figure 7 shows a longitudinal section of the dynamic pressure after a 14 s simulation. Knowledge of the dynamic pressure information is of vital importance in risk analysis as structural damage and risk for human life can be assessed from it.

Verification with actual mud flows has been done as well, showing very good fit with field data. For example, ^{Sheridan et al. (2011)} shows that the Titan2F predictions are within 10% of the data shown by ^{Procter et al. (2010)} for the highly channeled mud flow at Ruapehu, New Zealand. In addition, the theory was tested against field data assessed by ^{Williams et al. (2008)} for the 2006 Vazcun creek lahar at Tungurahua volcano, Ecuador, as shown by ^{Córdoba et al. (2015)}, where Titan2F predictions about maximum velocity are within 10% and within 15% for the measured super-elevation. In addition, ^{Rodríguez-Espinosa et al. (2017)} verified the accuracy of Titan2F predictions as well. They used the velocity measurements of the 2001 lahar in the Huiloac creek at Popocatapetl volcano, Mexico, done by ^{Muñoz-Salinas et al. (2007)}, to compare Titan2F predictions. Using the non-parametric ^{Kolmogorov (1933)} confidence test, ^{Rodríguez-Espinosa et al. (2017)} show a significance level of 0.01.

7. Conclusions

We present a new computational two-phase tool for lahar hazard assessment that has no constraints on initial volumetric particle concentration other than the maximum packing concentration (≤ 0.65). The program computes space-time evolution of the particle concentration, flow depth, velocity field, and dynamic pressure at each point of the computational grid.

The model is valid for two-phase flows whose phases consist of solids and water. However, the phenomenological approach used for the interphase drag model assumes an average diameter of the solids, which means individual boulders or particles cannot be tracked. In addition, the model is depth averaged, assuming thin layer and shallow-water approaches. Thus, our model correctly predicts the dynamics of gravity-driven flows providing the depth averaged values for the particle concentration, flow and phases velocities, and flow depth in a three dimensional topography. In order to model other kinds of geophysical mass flows, adjustments to the code must be done. For example, pyroclastic flows could be modeled if the flow density of the fluid phase is appropriately addressed (*e.g.*, the effect of temperature on air density using ideal gas law and an additional equation for temperature).

The proposed mathematical approach allows for the simulation of a range of flow behaviors. Regions with almost pure fluid to regions of friction-dominated flows are correctly described by the equations. Using this information, dynamic pressure is deduced, which becomes a useful tool for risk assessment.

The highly nonlinear coupled-equation system makes the time step very small. The use of this new tool on natural terrains or detailed DEMs requires high computational power. The use of a high performance work station with multiple cores is advised.

Important processes that are not addressed by this tool include the effect of turbulence, incorporation of solid material from the bed of the channel, and incorporation of water into the flow from existing water bodies. Nevertheless, this two-phase flow is an important step forward in forming an acceptable computational model for simulating hazardous natural phenomena.