## Abstract

Geophysical mass flows—debris flows, avalanches, landslides—can contain *O*(10^{6}–10^{10}) m^{3} or more of material, often a mixture of soil and rocks with a significant quantity of interstitial fluid. These flows can be tens of meters in depth and hundreds of meters in length. The range of scales and the rheology of this mixture presents significant modelling and computational challenges. This paper describes a depth-averaged ‘thin layer’ model of geophysical mass flows containing a mixture of solid material and fluid. The model is derived from a ‘two-phase’ or ‘two-fluid’ system of equations commonly used in engineering research. Phenomenological modelling and depth averaging combine to yield a tractable set of equations, a hyperbolic system that describes the motion of the two constituent phases. If the fluid inertia is small, a reduced model system that is easier to solve may be derived.

## 1. Introduction

The risk of volcanic eruption is one that public safety officials throughout the world face several times each year (Tiling 1989). The US Geological Survey reports that globally there are about 50 volcanoes that erupt every year. Across the world, in the 1980s, approximately 30 000 people were killed and almost a half million were forced from their homes because of volcanic activity (see National Research Council 1991, 1994).

Consequent to an eruption can be hazards, including passive gas emission or the slow effusion of lava, or violent explosions accompanied by stratospheric plumes and pyroclastic flows of red-hot ash, rock and gas that race along the surface of the mountain at speeds of up to hundreds of meters per second. These flows can melt snow, creating a muddy mixture of ash, rock and water. Other potential hazardous activity includes giant debris flows, coarse block and ash flows, and avalanches.

A different kind of geophysical mass flow is caused by intense rainfall on hillsides that are devoid of vegetation, generating massive mudslides. The 1998 mudflow at Casita Volcano in Nicaragua caused thousands of deaths. In the San Bernardino area of California during December 2003, hills that were scorched by summertime wildfires gave way in massive mudslide destroying homes and businesses.

In these flows, constituent particles are typically centimetre to metre sized, and the flows, sometimes as fast as tens of meters per second, propagate distances of tens of kilometres. As these flows slow, the particle mass sediments out, yielding deposits that can be tens to 100 meters deep and many kilometres in length (Legros 2002). These flows can contain *O*(10^{6}–10^{10}) m^{3} or more of material, often including a significant amount of water.

The range of scales and the complexity of the rheology for geological materials, coupled with the mathematical problem of describing a free surface flow, makes modelling and computing geophysical mass flows a significant challenge. We lack a full understanding of how these flows are initiated, but there is a growing understanding of processes governing flow, once that motion has been initiated. This paper describes one effort in modelling geophysical mass flows. Our efforts try to strike a balance between fidelity to physics on the one hand and mathematical and computational tractability on the other.

The modelling effort here has its origins in two distinct traditions of science. The first is the pioneering modelling of Savage & Hutter (1989). They began with mass and momentum balance laws based on a Coulomb constitutive description of dry granular material. Making use of the small aspect ratio of typical flows—flows are far longer than they are deep—Savage & Hutter developed a ‘thin layer’ model for granular flows down inclines. This work was later extended to two dimensions (Hutter *et al*. 1993), and to flows over more general topography (Gray *et al*. 1999; Pudasaini & Hutter 2003; Patra *et al*. 2005). An effort to relate thin layer model results to historic flows is presented in Sheridan *et al*. (2005). Hutter *et al*. (2005) examine the appropriateness of these thin layer models for several different types of geophysical flows. Iverson (1997) and Iverson & Denlinger (2001) argue that the presence of interstitial fluid alters the behaviour of flows and should be included in the constitutive behaviour of the flowing material. Starting with equations of mixture theory (Bedford & Drumheller 1983) and through a careful examination of experiments, these papers developed a system of mass and momentum balance laws for the mixture. Unfortunately, an equation for the motion of pore fluid was not readily available. Based on experimental data, a transport equation for the fluid was postulated. The two papers, Iverson (1997) and Iverson & Denlinger (2001), propose two different equations for this pore fluid motion. Computed solutions to the model equations have been compared with flume experiment results (Denlinger & Iverson 2001). A different approach to modelling mixtures may be found in Wang & Hutter (1999), in which a thermodynamically consistent model for solid and fluid flows is developed, and the resulting equations are applied to special flows on inclines. Related contributions can be found in Hutter & Kirchner (2003).

The second point of departure is the mathematical and engineering work on fluid–particle flows (Anderson & Jackson 1967; Drew 1983; Anderson *et al*. 1995; Jackson 2000). These so-called ‘two-phase’ or ‘two-fluid’ models are based on an averaging of mass and momentum balance laws for fluid and solid constituents. Interaction terms are modelled in one of two ways: either flow is assumed to have certain characteristics and particle motion is determined accordingly (e.g. particle motion in a potential flow), or experimental results suggest phenomenological interaction effects (Jackson 2000). In the chemical engineering literature, these models have achieved notable success in providing insight into the physics of fluidization, and in particular, into the stability and bifurcation of fluidized beds (Anderson *et al*. 1995; Glasser *et al*. 1996). Mixture theory and the phase-averaged theory are two distinct approaches to the problem of describing the movement of interacting materials. Historically, mixture theory largely grew out of solid mechanics, while phase averaging has its origins in fluid mechanics. Not surprisingly, these approaches to multiphase flow represent the constituent stresses and interaction forces in different ways. Indeed, what is meant by the constituent stresses is different in the two theories, and care must be exercised in describing and contrasting model systems. The interested reader may consult Joseph & Lundgren (1990) for a comparison of the mixture theory and phase-averaged formulations for a model flow composed of particles and fluid.

The modelling in the present paper begins with the two-fluid equations of Anderson & Jackson (1967), which contain phase interactions that are modelled phenomenologically; we follow Anderson *et al*. (1995) and include only the simplest of these interactions. Nevertheless, the model equations provide a rich enough description of solids and fluid flow and interaction while still being amenable to mathematical analysis. In particular, momentum equations for both fluid and solid phases are retained, providing equations for the velocities of both phases. After a brief review of the Savage–Hutter theory and the model of Iverson, §3 describes the two-phase equations. Through a combination of depth averaging and ad hoc modelling, a ‘thin-layer’ model is developed for the solid–fluid mixture. To arrive at a tractable system of equations, several assumptions are made. These modelling assumptions may later be subject to experimental validation. The resulting equations constitute a system of hyperbolic equations, as demonstrated in §4. If one assumes fluid inertia effects are small, a Darcy-like reduced system can be derived, and may be easier to solve computationally. Section 5 briefly discusses this reduced model. The paper concludes with remarks regarding computational solutions to the model equations, and the appropriateness of the models to real flows. We mention that our thin-layer, two-phase model may apply equally well to slush-snow avalanches (Bozhinskiy & Nazarov 1998).

## 2. Modelling flows

We begin by briefly summarizing the derivation of the original Savage & Hutter (1989) equations and the mixture model of Iverson (1997) and Iverson & Denlinger (2001). This summary is contained in two subsections. These subsections address several issues that will be revisited in §3, where the two-fluid model system is developed. Several mathematical issues that require care in their application are merely stated in this section; these issues are addressed more fully in the derivation of the two-fluid model. Savage & Hutter (1989) introduced a depth-averaging paradigm to model geophysical mass flows. Although they considered only dry granular flows, the subsequent modelling by Iverson examined debris flows containing a mixture of fluid and solids.

A note on sign convention: in soil mechanics, it is common to consider compressive stresses as positive, introducing a minus sign into the stress terms. This is the convention we follow in §3, and which Savage & Hutter use. Iverson follows the usual convention of denoting compressive stresses as negative. We caution the reader to observe the different conventions.

### (a) The Savage–Hutter model

Savage & Hutter (1989) considered an incompressible granular mass of constant density *ρ*^{s} in two dimensions flowing down an incline. They assumed the granular material could be described as frictional and modelled as a Coulomb continuum. By noting that typical flows tend to be much longer than deep, a scaling suggested omitting several terms in the governing mass and momentum balance laws. Depth averaging offered further simplifications, yielding a tractable model of hyperbolic equations describing granular mass flows. We briefly highlight this derivation now and refer the interested reader to the original paper for details.

The governing balance laws for mass and momentum may be written(2.1)Here * v* is the velocity of the granular material down a slowly varying incline. The incline is described by

*F*

_{b}(

*)=0, where*

**x***=(*

**x***x*,

*z*) is a two-dimensional coordinate system with

*x*the downslope direction and

*z*normal to the incline, and

*t*is time. T

^{s}is the solids stress tensor.

Several steps are required to complete the model. First, boundary conditions must be specified. The upper surface of the flowing mass at *F*_{h}(* x*,

*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

*F*

_{b}=0, and to satisfy a sliding friction law. This friction relation specifies that the shear traction and the normal stress are proportional: where

*ϕ*

_{bed}is the basal friction angle and the −sgn(

*v*) specifies that the shear traction opposes motion.

Second, a Coulomb constitutive relation is postulated for the material. The Coulomb law is a nonlinear relation among the components of T^{s}, and may be expressed by two relations.

Material deforms when the total stress reaches yield, , where 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.As material deforms, the stress and strain-rate tensors are aligned: , where the strain-rate and † denotes the transpose.

Implicit here is the assumption that the flowing material is everywhere in plastic yield. It is important to recognize that the mass and momentum balance laws (2.1) with a Coulomb constitutive relation are linearly ill-posed (Schaeffer 1987). Nevertheless, the Coulomb theory and variations of it are widely used in engineering practice.

The full Coulomb relations are too complex to be used in the Savage–Hutter derivation. Two simplifications are proposed. First, the proportionality (and alignment) of the tangential and normal forces that is imposed as a basal boundary condition is assumed to hold throughout the thin flowing layer of material. Written in components with superscripts denoting the tensor element, this implies *T*^{sxz}=*νT*^{szz} where the proportionality constant *ν* is a function of *ϕ*_{bed}. Second, following Rankine (1857) and later Terzaghi (1936), an earth pressure relation is assumed for the diagonal stress components, *T*^{sxx}=*k*_{ap}*T*^{szz}, where*ϕ*_{int} is the internal friction angle—a measure of material roughness—and the choice of the plus or minus sign depends on whether flow is locally expanding (the active state, with , and the − sign) or contracting (the passive state, with , and the + sign). (Note: these simplifications partially regularize the governing equations. More specifically, with these assumptions the governing equation (2.1) are no longer ill-posed, but rather are degenerate.)

The third step to specify the Savage–Hutter model is a scaling of the equations. The characteristic thickness of the flowing granular material is *H* and the characteristic length *L*. Scale *x* by *L* and *z* by *H*, time by the scale of free fall potential , and the *x* and *z* velocities by and , respectively. Stresses are scaled by *ρ*^{s}*gH*. The resulting equations contain the factor *ϵ*=*H*/*L*, which is small; *ϵ*≈0.001 is not uncommon (Iverson 1997). The scaled equations may be written in component form as(2.2)In these equations, cos(*θ*), sin(*θ*) are the local components of gravity, where *θ* is the angle that the basal surface *F*_{b}=0 makes with the horizontal. The limit *ϵ*→0 reduces the *z*-momentum equation to an expression of lithostatic equilibrium , and suggests ignoring the normal velocity *v*^{z}.

Taking *ϵ*→0 in the *x*-momentum equation yields an equation too simple to describe mass flow; the term *ϵ*∂_{x}*T*^{sxx} must be retained. We refer the reader to Savage & Hutter (1989) for a more detailed examination of this term.

The final step in deriving the Savage–Hutter model is a depth averaging of the incompressibility and *x*-momentum equations. That is, the equations are integrated in *z*, from the basal surface *z*=*b*(*x*) to the upper free surface *z*=*h*(*x*, *t*). Repeated use of the Leibnitz rule is made to interchange integration and differentiation, and boundary conditions are employed to evaluate terms at *b* and *h*. One additional approximation is required during this depth averaging. Specifically, if and the depth averaged velocity , the inertia term . Nonetheless, equality is assumed, recognizing that deviation from equality is due to non-uniformity of downstream velocity and is probably small.

Combining the four steps yields the thin layer equations(2.3)In writing these equations, we have dropped the overbar and overhat for simplicity of notation.

It is worth examining the various terms especially in the momentum equation, to appreciate the essential physics of model system. On the left-hand side of the equation are inertial terms and a pressure-like term *ϵk*_{ap}*h*^{2} whose factors arise from: (i) scaling; (ii) the earth pressure assumption; (iii) depth averaging *T*^{sxx}, which is proportional to *h* owing to lithostatic equilibrium in the *z*-direction. On the right of the equation is a basal frictional term, an acceleration due to changes in the topography of the basal surface *b*(*x*), and the driving force of gravity. It is also worth noting that, in the normal direction, . Steps similar to those outlined here will be employed in our two-fluid derivation below. Before that derivation, however, a brief summary of the mixture model of Iverson is presented.

### (b) The Iverson model

In an extensive review of debris flows, Iverson (1997) presented a thin layer model for a mixture of granular material and interstitial fluid in one spatial dimension. Later, Iverson & Denlinger (2001) extended the modelling to two dimensions, and included several refinements to the equations. We summarize Iverson's review and comment on the extensions of Iverson & Denlinger (2001).

Neglecting any generation of mass, the mixture theory equations of motion read(2.4)In writing this system, *ϕ* is the volume fraction of the solid phase, * u* is the fluid velocity, T

^{f}the fluid stress, and

*is a force density that includes all phase interactions. (See the two-fluid model derivation in appendix A, where buoyance effects are separated from other interaction terms.)*

**f**By adding together the species mass and momentum equations, Iverson derives conservation laws for the mixture. Guided by experimental results, several approximations are then made. Important for our purposes are the following.

The specific discharge

can be estimated, and it is found thatTo this end, the solid velocity can be substituted for the fluid velocity, simplifying the mixture momentum equation, and providing a simplified continuity equation**q**The mixture stress is taken to be the sum T

^{s}+T^{f}.When comparing with experiments, the viscous components of the fluid stress are often ignored compared with the solid stress contributions.

Through depth averaging, Iverson derives a system of equations similar to equation (2.3)(2.5)Several modifications to the Savage–Hutter equations are present, including an intergranular shearing force, proportional to the fluid viscosity *μ* and a parameter *η* that is a function the flow profile; this term is sometimes dropped when comparing predictions of the theory to experimental results. *P* is the pore fluid pressure, and is proportional to *hg*^{z} which is the hydrostatic pressure. Unfortunately, in using the solid velocity as an estimate of the fluid velocity, the fluid is constrained to move in lockstep with the solid material, and no change of the fluid pressure distribution is possible. That is, the fluid stress contribution is constant in time in Iverson (1997). The two-dimensional theory of Iverson & Denlinger (2001) postulates an advection–diffusion equation for the fluid pressure, allowing for changes in the pressure in response to the movement of solid. This paper also includes an intergranular shearing term that is important for flows over rapidly changing topography. Both papers provide comparisons with experiments and support for the validity of the proposed model.

## 3. The two-fluid model

We now turn to the derivation of a system of equations based on the two-fluid model. In three space dimensions, consider a thin layer of granular material and interstitial fluid, each of constant specific density *ρ*^{s} and *ρ*^{f}, respectively, flowing over a smooth basal surface. Neglect any erosion of the base. At any location, consider a Cartesian coordinate system O*xyz*, with origin O defined so the plane O*xy* is tangent to the basal surface and O*z* the normal direction. As noted in Iverson & Denlinger (2001), in considering flow over variable terrain, there should be no preferential direction in the *xy*-plane. We do assume a globally consistent orientation of this plane. Again, write * v*,

*for the velocities of the solid and fluid constituents (now with three components), respectively, and*

**u***ϕ*for the solid volume fraction. When writing equations in component form, we use superscripts to denote the component.

From Anderson & Jackson (1967), mass conservation for the two constituent phases may be written as(3.1)The momentum equations for the two species explicitly account for buoyance effects, and take the form(3.2)We note the difference between these two-fluid equations and the mixture model (2.4). See appendix A for a discussion of the phase-averaged equations and how the non-standard form of equation (3.2) accounts for buoyance, and Jackson (2000) and Joseph & Lundgren (1990) for a comparison of mixture theory and phase-averaged theory. Here * f* includes all non-buoyancy interaction forces of the fluid on the particle. The reader is referred to appendix A and Jackson (2000) for a discussion. We follow Anderson

*et al*. (1995) and Glasser

*et al*. (1996) and adopt a simple drag interaction (see also the discussion in Iverson 1997)where the leading factor of (1−

*φ*) accounts for the volume of the fluid acting on the entire particle phase,is a phenomenological function based on the experimental results of Richardson & Zaki (1954),

*v*

_{T}is the terminal velocity of a typical solid particle falling in the fluid under gravity,

*g*is the magnitude of the gravitational force and

*m*is related to the Reynold number of the flow. Appendix A provides details regarding these parameters.

In writing these equations, no a priori assumptions are made about the constitutive behaviour of the material, and specification of the stresses is still required. See appendix A for typical forms used in the context of fluidized beds. The solid is assumed to be an incompressible Coulomb granular material. Much of the derivation below could be generalized to a viscous Newtonian fluid, but, as Iverson (1997) argues, removing viscosity sometimes allows for a more direct comparison of theory with experiment and may not require as much fitting of parameters to data. Thus we consider the only fluid stress to be a pressure.

Finally, at several junctures in the derivation below, it is necessary to depth average a product of terms. In particular, for some function *f* we have occasion to calculate a term such as . We will assume this may be approximated as . Certainly, this is not justified mathematically. In some instances, it is possible to quantify the error made in the approximation. For example, in the Savage–Hutter theory sketched above, it is possible to say that and *γ* can be determined for certain flows (e.g. *γ*=6/5 if the flow has a parabolic profile). In other cases, it is not possible to provide such quantification. However, in appendix B, we provide an estimate of the error . This estimate provides a basis for judging whether the interchange of multiplication and depth averaging is justifiable.

### (a) Scaling

To proceed, we scale the equations in a manner similar to the Savage–Hutter development. The characteristic length in the *z*-direction is *H*, and in *x*, *y* it is *L*. The time-scale is taken to be *t*^{2}=*L*/*g*. Typical stresses are on the order of the weight of the material, so scale the solid stresses by *ρ*^{s}*gH* and the fluid stress by *ρ*^{f}*gH*. After clearing the coefficients, the continuity equations for the solid and fluid are unchanged. Several terms in the momentum equations are multiplied by *ϵ*=*H*/*L* which is small. Specifically, the solid momentum balances in the *x*, *y* and *z* directions become(3.3)

(3.4)

(3.5)In writing these equations, we have explicitly scaled the momentum interaction force, leaving a scaled parameter *v*_{T} in the drag coefficient. Also, we note that components of gravity have been scaled by the magnitude *g*, so (*g*^{x}, *g*^{y}, *g*^{z}) is a unit vector. Usually, one would drop all terms of order *ϵ*; again we refer the reader to Savage & Hutter (1989) for a discussion of the need to retain the diagonal stress contributions. Because no preferential downslope direction is prescribed, and the flow direction may change during a flow, we retain all the stress terms in both the *x*- and *y*-directions, dropping only *O*(ϵ) terms in the *z*-direction; see the discussion in Iverson & Denlinger (2001).

The fluid equations are likewise scaled, and the common factor (1−*φ*) removed(3.6)(3.7)(3.8)We note that the different scaling for the solid and fluid momentum equations has produced phase interaction terms with different coefficients. Thus, adding the component momentum equations will not result in the force terms cancelling each other.

### (b) Mass balance

Divide the mass conservation equations for the solid and fluid phase by the constant specific density *ρ*^{s} and *ρ*^{f}, respectively, and add the two equations together to find(Observe for the time derivatives .) That is, the volume-weighted mixture flow is divergence free. This equation is not an assumption, but rather a direct consequence of mass balance. That the mixture is isochoric allows us to depth average,(3.9)Use the Leibniz rule to interchange differentiation and integration. The upper free surface *F*_{h}(* x*,

*t*)=0 is assumed to be a material surface for the mixture, so at

*z*=

*h*(

*x*,

*y*,

*t*)(3.10)Likewise, at the basal surface,

*F*

_{b}(

*)=0, the flow is tangent to the fixed bed, so at*

**x***z*=

*b*(

*x*,

*y*)(3.11)In arriving at these equations, we have ignored sedimentation, entrainment and infiltration of fluid into the bed.

After algebraic manipulation, we find an equation for the total mass of the solid and fluid(3.12)In writing this equation, the depth-averaged velocities , with a similar expression for the components , , , where again , and is the depth-averaged volume fraction.

We need to follow the evolution of the depth-averaged volume fraction , but there is no obvious candidate equation to which we could appeal. We turn to the mass balance law for the solid phase and depth average. The lack of phase-specific basal and free surface conditions would seem to preclude such an integration, but reasonable assumptions allow us to proceed. More specifically, beginning with the integration(3.13)we find(3.14)In order for the boundary terms to vanish here, the solid and fluid phases must separately satisfy material-free surface and basal surface conditions, just as the mixture does ((3.10) and (3.11)). If no erosion or infiltration occurs, it is reasonable to assume the motion of each species is tangent to the basal surface. Imposing the material surface condition on the upper free surface is less justified. If the fluid and solid velocities do not differ by much, consideration of flows for which the *x*- or *y*-velocity is zero suggest special conditions under which each phase does in fact satisfy a material-free surface condition. We make the assumption, recognizing it is subject to error: at *z*=*h*,(3.15)(3.16)With this assumption, we arrive at a solids mass balance equation(3.17)To be sure, the imposition of individual phase mass balance requires experimental validation. The justification offered here breaks down when the difference * v*−

*is not small—such as at the time of final deposition, where interesting phase separation occurs.*

**u**### (c) *z*-Momentum

Observe that, upon setting *ϵ* to zero in the fluid *z*-momentum equation, we find the fluid to be hydrostatic,Integrating and imposing boundary conditions, we find(3.18)and the average(3.19)

In the same manner, for the solid *z*-momentum (3.5) we find an equation for an effective stress,Substituting,Thus, the normal solid stress in the *z*-direction at any height is equal to the (buoyancy) reduced weight of the solid material overburden.

Observe that, owing to scaling, the *z*-velocities have been dropped from the *z*-momentum equations. Of course, neglecting accelerations in the *z*-direction is a central component of the thin layer theory.

### (d) *x*- and *y*-Momentum

We now must depth average the remaining momentum equations. In several steps of the derivation that follows, we will use the approximation , and we do not indicate every place at which this is performed. Again, we refer the reader to appendix B.

First, consider an equation for the motion of the solid phase. The left-hand side of the *x*-momentum equation (3.3) can be writtenDepth averaging and using boundary conditions yields(3.20)Now, depth averaging the right-hand side of equation (3.3) yields(3.21)In order to proceed, several assumptions are made, namely

A free surface condition for the individual species is specified, as above.

The drag term is phenomenological and highly nonlinear, and expressing the depth-averaged function in terms of , , is all but impossible. However, it is reasonable to expect that experiments could fit a different phenomenological-averaged drag of a similar form:1 . This provides term (iii).

The earth pressure relation for the solid phase is employed. To this end, basal shear stresses are assumed to be proportional to the normal stress,where * can be either

*x*or*y*, and the velocity ratio determines the force opposing the motion in the *-direction, to the extent that force is mobilized (Hutter*et al*. 1993; Patra*et al*. 2005) in that direction. Thus, sliding velocity and basal traction are collinear. The*α*notation will provide a convenient shorthand. Likewise, the diagonal stresses are taken to be proportional to the normal solid stresswhere the same index*x*or*y*is used in both *s. Finally, following Iverson & Denlinger (2001),*xy*shear stresses are determined by a Coulomb relationwhere the sgn function ensures that friction opposes straining in the (*x*,*y*)-plane.Observe that

From the assumptions above and the formulae for the fluid and solid normal stresses, we find(3.22)Finally, using the fluid and solid stress relation , term (i) is approximated as(3.23)Because the upper free surface is stress free, all terms *T*^{szz}|_{z=h} vanish. This provides the expression(3.24)The (−*g*^{z}) factor appears here and below, originating in the evaluation and depth averaging of the fluid stress. In typical flows, this factor is positive.

Combining all terms yields the solids *x*-momentum equation(3.25)The *y*-solid momentum equation can be derived in a similar fashion and yields(3.26)

We now proceed to derive an equation for the fluid motion. Depth averaging for the fluid presents fewer difficulties than in the solid equations. The *x*-momentum equation takes the form(3.27)The fluid *y*-momentum equation is similar,(3.28)

### (e) The full model equations

The complete set of model equations are rewritten here. To simplify notation, hats and overbars are dropped.(3.29)The model equation (3.29) constitutes a system of six equations in the variables *h*, *hϕ* and the four momenta and (the ^{∼} denoting the *x*- and *y*- components of the vector only). Because of the buoyancy term in the solid momentum balance, the equations cannot be put in conservation form. The structure of the equations is symmetric upon interchange of the *x*- and *y*-coordinates. Of course, it is possible to derive other forms of the solid or fluid momentum equations by algebraic manipulations of the equations presented here. In §4, we show this system is hyperbolic, under common conditions. Together with initial and boundary conditions, the system constitutes a well-posed set of equations describing a two-fluid thin layer system.

Exact solutions to this system are difficult to find. We refer the reader to Savage & Hutter (1989) and Chugunov *et al*. (2003) for a discussion of solutions to the model equations for dry flows, and to Iverson (1997) and Iverson & Denlinger (2001) for solutions to the mixture equations. For purposes here, note the special case of flow in one space dimension. That is, consider the model equation (3.29) confined so that all motion is in the *x*-direction only, and set all *y*-derivatives, *y*-stress terms and *y*-velocities to zero. Further assume the basal surface has constant derivative in the *x*-direction. Physically, one might imagine flow in a narrow channel with a rectangular cross-section, or flow down an infinite plane with uniform cross-section. Under such assumptions, a constant solution may exist, depending on prescribed quantities such as the fluid–solid density ratio, the components of gravity (which are determined by the basal surface), and the thickness ratio *ϵ*. Examining the equations, the volume fraction satisfiesand the species velocities are related through the drag term. The value of *h* is not specified by the equations. This constant solution ceases to exist if the right-hand side is negative.

## 4. Hyperbolicity

The analysis of hyperbolicity is most familiar in one space dimension, and we show that analysis here. The full two-dimensional equations share this hyperbolicity, but the algebra required to demonstrate it is significantly more complicated than what follows.

In this section, consider the model equation (3.29) confined so that all motion is in the *x*-direction only. Because the motion is in one direction, we drop the superscripts that denote the components of velocity. The one-dimensional equations may be written in non-conservative form aswhereTo analyse the hyperbolicity of the system, we must compute det(*B*−*λA*) for the eigenvalue *λ*. After tedious algebra, we find the characteristic polynomialwhere

Certainly if *h*→0 or if *ϕ*→1 or *ϕ*→0, the determinant vanishes and hyperbolicity is lost. Away from these degeneracies, we claim the system is hyperbolic. One can show that *a*, *b* and *c* are all positive away from *h*=0 (or the unphysical case *g*^{z}=0). Observe that when *ϵ* is zero, there are two pairs of double real roots; the roots remain real, and thus the system is hyperbolic, at least for *ϵ* sufficiently small (one can appeal to Descartes' rule of signs). However, except under special conditions, deriving the roots explicitly is non-trivial. Knowing the roots are real, we will appeal to a perturbation argument to find their approximate values.

### (a) Case *u*=*v*

The case *u*=*v* is special and requires separate consideration. One has a quartic equation, sayThe quadratic formula can be applied,The discriminant isBy considering this as a quadratic formula in *ϕ*, the discriminant can be shown to be positive, at least for physically relevant values of *α*_{xx} and *ρ*^{f}/*ρ*^{s},2 and the square root gives real values. Because *c*>0, this root is less than (*a*+*b*), and another square root can be taken to findThus the system is hyperbolic in this special case.

### (b) Case *u*≠*v*

More generally, we find a perturbation expansion for *λ* when *u*≠*v*. An order of magnitude analysis suggests the expansion of *λ* in terms of , so let .

*O*(*1*)* problem:* , which implies *λ*_{0}=*v* or *λ*_{0}=*u*.

*problem:* . Because *λ*_{0} is either *u* or *v*, this equation gives no information about *λ*_{1}.

*O*(*ϵ*)* problem:* . When *λ*_{0}=*u*, the equation becomes , which implies *λ*_{1}=±*b*. Similarly, if *λ*_{0}=*v*, *λ*_{1}=±*a*.

Further terms may be computed: up to , the eigenvalues are(4.1)

## 5. Reduced model

The two-fluid model, described in §2, assumes the interstitial fluid is inviscid, the only stress being a pressure. The model does include inertial effects in the fluid momentum equations. We now present an alternative model system, one that simplifies the effort required to compute solutions numerically.

If the fluid motion relative to the solid motion is not large (Iverson 1997), and if the pressure is the only fluid stress, one might consider a Darcy-like approximation and ignore the fluid acceleration. The fluid momentum equations in the *x*- and *y*- directions reduce to the pair(5.1)This pair of equations replaces the fluid momentum equation (3.29)_{5,6}, and one should view this reduced system (consisting of (3.29)_{1,2,3,4,7,8} and (5.1)) as a set of evolution equations for *h*, *hϕ* and . Given , (5.1) is an equation for .

Equation (5.1) for contains derivatives of *h*, a remnant of the gradient of the fluid pressure. When substituting back into the total mass conservation equation (3.29)_{1}, a derivative of appears, leading to second derivatives of *h*. Other equations contain only first derivative terms. Calculations indicate that the coefficient matrices of the first-order derivative terms retain a hyperbolic structure. Thus the reduced model is degenerate parabolic—that is, the coefficient matrices for second derivatives have zero eigenvalues. Nonetheless, the reduced model system is a well-posed evolution system once appropriate initial and boundary conditions are assigned.

This Darcy model is clearly a consequence of the modelling assumption that drag is the principal (non-buoyancy) contribution to interphase momentum transfer. As mentioned earlier, and as noted in appendix A, the precise form of the permeability coefficient—borrowed from the Richardson–Zaki theory—presents difficulties in the limiting cases of pure fluid and pure granular material. Other forms of the coefficient may be better suited to these special cases.

Numerical computations with the reduced Darcy model demonstrate how the presence of a fluid increases the length of the flow runout, relative to the no-fluid case. In the model, this occurs primarily because the effective solid stresses are lower, relative to the stress without pore liquid, generating a smaller overburden. Basal friction, proportional to this overburden, is then applied only to the solid volume fraction. Figure 1 illustrates this longer runout, showing the computed solution to the reduced model equations for two cases of initial volume fraction. The mass flows down an incline that is divided into three sections. The upper part of the incline (from *x*=0 to 60 in figure 1) makes an angle of 30° with the horizontal, the middle part an angle of 15° and the third section (from *x*=80 onwards) is horizontal. In the computations, a slight smoothing of the corners was performed. The granular material was relatively free-flowing, with an internal friction angle of 35° and a basal friction angle of 20°.3 The initial pile was a parabola of horizontal extent 10 and maximum height 6 (thus a relatively large pile on the scale of this numerical experiment). The computational domain is sufficiently large that the pile never approaches the edges during the simulation, allowing us to set boundary conditions to zero for these simulations. (Of course, other boundary conditions could be applied, for example, absorbing conditions would mimic the pile flowing out of the computational domain.) The plots show the total pile height at the beginning, intermediate and long times. In the plot on the left, the initial solid fraction was 1.0; for the plot on the right, the initial solids fraction was uniform in *x* with a value of 0.6. We see that a mixture flow travels much further in a fraction of the time, compared with a solids-only flow. Indeed in these simulations, the *ϕ*=1 pile had essentially stopped at *x*≈96 by the time *T*=40; the *ϕ*=0.6 pile had propagated to *x*≈140 by the time *T*=10, and was still travelling at a non-negligible rate (approximately 40% of its maximum velocity) when we stopped the simulation.

## 6. Conclusion

We have presented a ‘thin layer’ system of equations modelling the depth-averaged flow of granular material coupled to interstitial fluid, as the two-fluid mix moves over smooth topography. Several assumptions are made to develop these equations. A detailed comparison with experiments is needed to demonstrate whether the assumptions and the resulting model equations adequately describe the physical flows. Using reasonable parameters, the model equations constitute a hyperbolic system. We have also presented a reduced model that ignores inertial terms in the fluid phase, much like Darcy's law. Computations with the reduced system show that the presence of a pore fluid greatly increases the runout of a pile of fluid–solid material, relative to the runout for a pure granular material. Both the full model equations and the Darcy system can be integrated into numerical software (Denlinger & Iverson 2001; Patra *et al*. 2005), allowing simulation of flows over realistic terrain. This computational framework will allow for careful comparison with experimental flows both in the laboratory and in large flumes, and with geologically important flows.

Relative to earlier models of mixture flow, the two-fluid model equations developed here (as well as the reduced model equations) offer one significant advantage: velocities for both solid and fluid phases are determined from physically proper momentum balance relations. Certainly, approximations are made to derive a tractable system. However, the momentum balance laws in the two-fluid theory can be traced back to balance laws that are well accepted in fluid mechanics.

The reader may ask whether the two-fluid model (3.2) is an appropriate place to begin development of a thin layer model. Appendix A offers a summary of the origins of this form of momentum balance. As that appendix shows, there is an alternative form that replaces the non-conservative product *ϕ*∇.T^{f} in the solid momentum balance with a convective derivative . Further algebraic manipulations can bring this alternative formulation into a more benign looking form. However the additional convective derivative presents computational challenges. Algebraic operations to remove the derivative introduces other complicating terms. Thus, we content ourselves with the non-conservative form in equation (3.2). It is interesting to note that a formal mathematical derivation of two-phase equations is given by Drew (1983) and yields essentially the same non-conservative system (3.2).

There are several forces present in equation (3.2) acting on the constituent phases, including solid and fluid stresses, gravity, buoyance and drag. In the original derivation of the equations by Anderson & Jackson (1967), the fluid stress that appears in equation (3.2) is not the simple Newtonian fluid stress, but rather it is the usual fluid stress at every point weighted by a solid phase indicator function and integrated. The solid stress in equation (3.2) is the sum of point contact forces exerted on a specified particle, again weighted by a solid phase indicator function and averaged. In both cases, the stress T^{f} and T^{s} that appears in equation (3.2) are not the usual Cauchy stresses of the pure phase constituents. Nevertheless, after postulating the equations, it is common to identify the stress terms that appear with their pure phase analogues and to postulate pure phase constitutive relations for those terms. Such an identification is a modelling assumption.

In the exposition, we assume pressure is the only fluid stress. This fact, coupled with the presence of the buoyancy term in the basic equations of motion, introduces a net force akin to Terzaghi's original postulate of an effective stress—that is, the net force is the solid stress minus a pore fluid pressure, . It is the buoyance effect that produces the long runout demonstrated in the reduced model calculations—the normal force of the solid is reduced, causing a reduction in energy dissipation by basal friction (which is active only on the solid phase, whereas the fluid phase does not contribute at all to basal friction in the model here). This buoyancy effect is not present in earlier models of solid–fluid avalanche and debris flows.

## Acknowledgments

This work has been supported by NSF grants ITR 0121254 and EHR 0087665. Computing support from the Center for Computational Research at the University at Buffalo is gratefully acknowledged. For their diligent reading of this manuscript, and for their suggestions which led to improvements in the paper, we gratefully acknowledge the assistance of the referees.

## Footnotes

One contribution of 11 to a Theme ‘Geophysical granular and particle-laden flows’.

↵In studies of fluidized beds in a cylinder using one-dimensional ‘slice analysis’—that is, averaging the model equations over the cross-section of the cylinder—such an averaged drag is often postulated; see Jackson (2000).

↵One can demonstrate that the discriminant is positive for special cases (e.g.

*α*_{xx}=1 and*ρ*^{f}/*ρ*^{s}=1/2, and show there is no sign change as*α*_{xx}changes while remaining*O*(1) and while 1>*ρ*^{f}/*ρ*^{s}>0.↵Note a very thin pile of granular material without pore fluid released on the middle section would not flow. In these computations, momentum generated on the steeper first section carries the pile through the middle section.

↵These constitutive relations for the granular phase are similar to the granular kinetic theory relations of Jenkins & Savage (1983).

- © 2005 The Royal Society