## Abstract

To simulate debris-flow behaviour from initiation to deposition, we derive a depth-averaged, two-phase model that combines concepts of critical-state soil mechanics, grain-flow mechanics and fluid mechanics. The model's balance equations describe coupled evolution of the solid volume fraction, *m*, basal pore-fluid pressure, flow thickness and two components of flow velocity. Basal friction is evaluated using a generalized Coulomb rule, and fluid motion is evaluated in a frame of reference that translates with the velocity of the granular phase, *v*_{s}. Source terms in each of the depth-averaged balance equations account for the influence of the granular dilation rate, defined as the depth integral of ∇⋅*v*_{s}. Calculation of the dilation rate involves the effects of an elastic compressibility and an inelastic dilatancy angle proportional to *m*−*m*_{eq}, where *m*_{eq} is the value of *m* in equilibrium with the ambient stress state and flow rate. Normalization of the model equations shows that predicted debris-flow behaviour depends principally on the initial value of *m*−*m*_{eq} and on the ratio of two fundamental timescales. One of these timescales governs downslope debris-flow motion, and the other governs pore-pressure relaxation that modifies Coulomb friction and regulates evolution of *m*. A companion paper presents a suite of model predictions and tests.

## 1. Introduction

Debris flows are water-saturated masses of soil and fragmented rock that rush down mountainsides, funnel into stream channels and form lobate deposits when they spill onto valley floors. Their geological nature and mechanical behaviour make them intermediate in character between rock avalanches and flash floods. Because debris flows have solid grain concentrations that exceed 0.4, peak speeds that commonly surpass 10 m s^{−1}, and volumes that range up to approximately 10^{9} m^{3}, they can denude slopes, bury floodplains and devastate people and property [1,2]. Notable recent debris-flow disasters involved more than 20 000 fatalities in Armero, Columbia, in 1985 and in Vargas state, Venezuela, in 1999 (figure 1).

Alternative terms such as mudflow, mudslide, debris torrent and lahar are sometimes used to describe debris flows, but the terms ‘debris’ and ‘flow’ have precise geological meanings. ‘Debris’ implies that grains with greatly differing sizes and irregular shapes are present. This trait—and the consequent lack of a characteristic grain size—distinguishes debris-flow mixtures from most man-made granular mixtures. In debris flows, the largest grains can have linear dimensions exceeding 10 m, but the presence of at least several weight per cent mud-sized particles smaller than about 62 *μ*m is generally more critical [3]. Hydrodynamic suspension of these small particles increases the effective viscosity of the muddy water that fills pore spaces between larger grains [4], and this enhanced viscosity can promote development and persistence of high pore-fluid pressure that facilitates debris-flow motion by reducing grain contact forces [5–7]. The term ‘flow’ implies that rearrangement of grain contacts is pervasive during debris-flow motion. Indeed, granular debris that is liquefied by high pore-fluid pressure can appear to flow almost as fluidly as water [8].

Spatial and temporal changes in macroscopic material behaviour that result from local rearrangements of grains and attendant evolution of pore-fluid pressure pose fundamental challenges for continuum mechanical modelling of debris flows. The most conspicuous transitions occur as debris mobilizes from static material on slopes, liquefies and flows rapidly, and later regains rigidity during consolidation of deposits [9,10]. Most debris-flow models neglect these transitions, and instead treat rheology or depth-averaged flow resistance as inherent properties of debris [11]. With this approach, flow dynamics simulations typically employ basal flow resistance coefficients less than half as large as those necessary to statically balance forces at flow initiation sites [12–14]. Unbalanced initial states are held in check by assuming that an imaginary dam restrains the debris until the modeller issues a command, but use of this type of initial condition compromises physical relevance.

By contrast, most natural debris flows commence when balanced forces are infinitesimally perturbed. Flow onset commonly results from rainfall or snowmelt that triggers failure of debris-mantled slopes or mobilization of scree in steep rills and gullies. As masses of water-saturated grains begin to move, however, the governing force imbalance can evolve dramatically owing to pore-pressure feedbacks that modify the apparent rheology of the debris [15–17]. Differences between this type of behaviour and behaviour that arises from a stipulated rheology and force imbalance have great practical as well as theoretical significance, because pore-pressure feedbacks can determine whether a rapid debris flow develops at all—as opposed to a creeping landslide that moves imperceptibly or intermittently downslope [18,19].

## 2. Objectives

Our chief goal is seamless simulation of debris-flow motion from initiation to deposition without any redefinition of governing equations, re-evaluation of parameters or restructuring of numerical methods. An ancillary goal is efficient machine computation of solutions for use in practical applications. We pursue these goals by formulating a depth-averaged model that allows feedbacks to develop during coupled evolution of solid and fluid volume fractions, pore-fluid pressure and debris-flow velocity and thickness. The feedbacks involve dilatancy – the state-dependent propensity of granular materials to undergo changes in solid volume fraction as they shear. Well-known since the time of Reynolds [20,21], variable dilatancy underpins the critical-state theory of soil mechanics [22–24], and it plays a key role in determining of the continuum-scale rheology of dense granular flows [25].

Previous depth-averaged debris-flow models have included effects of solid–fluid interactions, and a few have accounted for evolution of solid and fluid volume fractions [5,26–35]. However, no previous model has explicitly considered dilatancy coupled to pore-pressure feedbacks that mediate transitions between static and dynamic states. Dry granular materials can undergo analogous state transitions [36], which can be characterized by using a continuum mechanical model that employs an evolving order parameter to account for changes in mobility due to changes in solid volume fraction [37]. The order parameter is not physically measurable, however. By contrast, direct measurements of pore-pressure evolution due to variable dilatancy have led to development of a model that successfully simulates a transition from failure to flow that occurs when water-saturated granular avalanches are manually triggered by increasing the steepness of submerged slopes [38]. Additional experiments with concentrated grain–fluid mixtures have shown that a state-dependent Coulomb friction rule that embeds the effects of variable dilatancy can successfully describe bulk flow resistance over a broad range of normalized shear rates, including quasi-static rates [39]. Our model is founded partly on these findings.

All physical principles and empiricisms used in our model are implemented within a depth-averaged framework. Use of depth-averaged equations in conjunction with shock-capturing numerical methods and adaptive mesh refinement yields a combination of computational speed and accuracy well suited for assessment of diverse flow hazards [40,41]. Like all depth-averaged models, however, our model neglects some details of three-dimensional behaviour. It nevertheless provides predictions that can be rigorously tested because it computes flow depths, velocities and basal pore pressures with a resolution similar to that of the most-detailed debris-flow data collected to date [10,16,42–45]. Use of the solid volume fraction as an additional prognostic variable expands the possibilities for future model tests.

Section 3 describes physical principles and observations that guide formulation of our model. Subsequent sections describe derivation of the depth-averaged model equations and analytical results that help clarify their physical implications. A companion paper [46] examines mathematical properties of the model, describes our numerical method for solving the model equations, and tests computational predictions of initiation of debris-flow motion as well as subsequent flow and deposition.

## 3. Physical principles and definitions

To represent debris-flow behaviour, our model combines principles from soil mechanics, grain-flow mechanics and fluid mechanics, which we quantify below. We also quantify our rationale for treating two-phase debris flows as granular flows with evolving solid volume fractions and pore-fluid pressures. In this approach the fluid-phase velocity is evaluated in a frame of reference that translates with the granular phase, and variations of fluid velocity result in variations of fluid pressure. Our model excludes the effects of grain-size segregation and bed-sediment entrainment, which are important in many debris flows [44,47]. It is designed to incorporate the effects of segregation and entrainment in a manner consistent with the depth-integrated analyses of Gray & Kokelaar [48] and Iverson [49], however.

### (a) Continuum conservation laws

Within each phase of a debris-flow mixture, and also within the mixture as a whole, conservation of mass is expressed by
*ρ* is the mass density and ** v** is the velocity vector. Similarly, conservation of linear momentum is expressed by

**is the acceleration due to gravity and**

*g**T*is the stress tensor, which is assumed to be symmetric [50]. We use the standard sign convention of continuum mechanics to define shear components of

*T*[51,52]. However, because granular debris can sustain negligible tension, we define normal components of

*T*using a soil mechanics sign convention in which compressive stresses are positive. Thus, in rectangular Cartesian coordinates we express the

*x*component of ∇⋅

*T*as −∂

*τ*

_{xx}/∂

*x*+∂

*τ*

_{yx}/∂

*y*+∂

*τ*

_{zx}/∂

*z*. Use of this convention simplifies our definition of constitutive equations and evaluation of boundary tractions.

Our model includes no energy-conservation equation because we treat debris flows as isothermal systems. By contrast, Bartelt *et al.* [53] have developed a depth-averaged model of granular snow avalanches that includes an energy equation roughly analogous to those used in kinetic grain-flow theories. In the model of Bartelt *et al.* [53], the effective frictional resistance of an avalanche evolves because it is mediated by evolution of kinetic energy associated with grain velocity fluctuations. The effective frictional resistance evolves in our model as well, but it does so as a result of evolution of dilatancy and pore-fluid pressure. The difference between our formulation and that of Bartelt *et al.* [53] is logical, because debris flows contain muddy pore water that is roughly 10^{3} times more dense and viscous, as well as 10^{4} times less compressible, than pore air at Earth's surface. Therefore, as demonstrated by dimensional analysis, pore fluid is likely to play a pivotal role in debris flows but not in most dry avalanches [54].

### (b) Definitions of volume fractions, densities and velocities

Our model assumes that the solid and fluid constituents of debris are incompressible, with mass densities *ρ*_{s} and *ρ*_{f}, respectively. Therefore, the mixture bulk density *ρ* varies exclusively as a result of grain rearrangements and elastic shear distortion localized at grain contacts (cf. [55]). Our definition of the mixture bulk density,
*m* and fluid volume fraction 1−*m* can be regarded as continuous variables. Commonly *m* varies from about 0.4 in the most dilute debris flows (or parts of debris flows) to 0.8 in the densest (cf. [56]). This variation implies that *ρ* ranges from about 1700 to 2400 kg m^{−3}, because the solid constituents of debris are largely rock fragments with *ρ*_{s}≈2700 kg m^{−3}, and the fluid constituents are mostly muddy water with *ρ*_{f} ranging from about 1000 to 1200 kg m^{−3} [5,57].

Like the mixture density, the mixture momentum is weighted by the mass of solid grains *ρ*_{s}*m* and mass of pore fluid *ρ*_{f}(1−*m*) per unit volume, but it also depends on the velocities of each phase. Thus, the linear momentum of the mixture is defined as *ρ*** v**=

*v*_{s}

*ρ*

_{s}

*m*+

*v*_{f}

*ρ*

_{f}(1−

*m*), where

*v*_{s}is the velocity of the solid grains, and

*v*_{f}is the velocity of the pore fluid. This definition of mixture momentum implies that the mixture velocity is defined as

*v*_{f}−

*v*_{s}. To an observer moving with the local solid velocity

*v*_{s}, the apparent fluid velocity is the volumetric flux of pore fluid per unit area of mixture,

**:**

*q***is termed the specific discharge [58].**

*q*To simplify our model, we approximate the mixture velocity as ** v**≈

*v*_{s}. The basis for this approximation can be established by using (3.3)–(3.5) to find that

**≈**

*v*

*v*_{s}, because values

*ρ*

_{f}/

*ρ*≈1/2 and ∥

**∥/∥**

*q*

*v*_{s}∥≪1 are typical [5]. For the initial and final stages of debris-flow motion, with ∥

*v*_{s}∥≈0, (3.6) instead reduces to

**≈(**

*v**ρ*

_{f}/

*ρ*)

**. Nevertheless,**

*q***≈**

*v*

*v*_{s}remains a good approximation because ∥

**∥≪0.1**

*q**ms*

^{−1}is typical, implying that the mixture momentum balance (3.2) reduces to the quasi-static form

*ρ*

**+∇⋅**

*g**T*≈0 [59].

### (c) Definitions of total stress, effective stress and pore-fluid pressure

A formulation that approximates the mixture momentum as *ρ**v*_{s} requires that the total mixture stress *T* must account not only for solid- and fluid-phase stresses but also for all effects of solid–fluid interactions. These effects can be diverse, but in mixtures of fluid and macroscopic grains (subject to negligible Brownian forces), the effects are commonly dominated by buoyancy and viscous drag [33,60]. Other interaction forces include hydrodynamic-lift, added-mass and Basset forces, but these forces are probably muted in debris flows owing to restriction of fluid flow caused by high concentrations of grains with differing shapes and sizes [5,6]. Indeed, a linear Darcian drag law can explain relaxation of non-hydrostatic pore-pressure gradients measured in freshly agitated, liquefied debris-flow slurries [61]. On this basis, our model assumes that buoyancy and viscous drag wholly account for the effects of solid–fluid interactions on the mixture stress, *T*. As shown below, the net effect of these interactions is manifested on a bulk continuum scale by gradients of pore-fluid pressure, *p*.

Our model isolates the effects of *p* by using a principle of soil mechanics, which states that *T* can be expressed as the sum of *p* and an effective stress *T*_{e} [55,62]. By convention *T*_{e} and *p* are treated mathematically as if they act throughout the entire solid–fluid mixture, although this practice is largely a matter of convenience and is not mandated by mechanics [63,64]. By combining *T*_{e} and *p* with a viscous deviatoric stress component *T*_{vis} (which acts only in the fluid volume fraction 1−*m* as a result of shearing of the fluid phase), we define the total mixture stress as
*I* is the identity tensor. Equation (3.7) indicates that, in the presence of constant *T*, increases in *p* cause commensurate reductions in the mean effective normal stress *σ*_{e}, defined as
*σ*=(1/3)*tr*(*T*) is the mean total normal stress.

To distinguish the effects of buoyancy from those of solid–fluid drag, we subdivide *p* into a hydrostatic component *p*_{h} and a non-hydrostatic ‘excess’ component *p*_{e}
*p*_{h} depend only on the evolving geometry of a mass of fluid-saturated debris, whereas values of *p*_{e} are affected by ** q**. In turn,

**is a consequence of debris dilation or contraction.**

*q*### (d) Definition of dilation rate

We define the local debris dilation rate as the divergence of *v*_{s}, which is related to evolution of the solid volume fraction *m* by
*d*/*dt*=∂/∂*t*+*v*_{s}⋅∇ is a material time derivative in a frame of reference that moves with the granular velocity *v*_{s} [58]. Thus, *dm*/*dt*<0 indicates a positive dilation rate, whereas *dm*/*dt*>0 indicates a negative dilation rate (i.e. a rate of contraction). Although ∇⋅*v*_{s} does not appear explicitly in our final model equations, *dm*/*dt* plays a prominent role.

Because we assume that incompressible fluid fills all pore spaces in the debris (except where grains may protrude from the free upper surface of a flow, as detailed in §4), the debris dilation rate is locally balanced by the divergence of ** q**. This fact can be demonstrated formally by manipulating the mixture mass-conservation equation (3.1) while embedding the definitions (3.3)–(3.5) to obtain

A central postulate of our model is that the dilation rate depends on two mechanical phenomena that modify *m*: changes in *σ*_{e} that cause reduction or enlargement of the debris pore space due to bulk compressibility of the mixture, and changes in mixture shear strain that cause pore space contraction or expansion due to dilatancy. We assume that these two phenomena are independent of one another, except insofar as they each involve interactions with other variables such as pore-fluid pressure. We therefore define the total rate of change of *m* as
*γ* is a scalar measure of natural shear strain, which is referenced to evolving granular-phase material coordinates that move with velocity *v*_{s} (cf. [52, p. 151]), and *dm*/*dt* by considering idealized deformations consisting of isotropic volume strain and homogeneous simple shearing in *x*−*z* and/or *y*−*z* planes (figure 2). This idealization yields a parsimonious model that can be tested using existing datasets [46].

### (e) Definition of compressibility

Application of (3.12) requires a constitutive relationship to determine ∂*m*/∂*σ*_{e}. At effective confining stresses typical of debris flows less than 50 m thick (*σ*_{e}≤10^{6} Pa), bulk compression or expansion of granular debris caused by changes in *σ*_{e} results almost entirely from changes in pore volume rather than from changes in the densities of the debris’ solid and fluid constituents [55,62]. Therefore, the debris compressibility *α* can be defined by using a relationship commonly employed in soil mechanics [58]
*α*=[*m*(*ρ*_{s}−*ρ*_{f})]^{−1}(∂*ρ*/∂*σ*_{e}), a form that reveals an analogy with ideal-gas compressibility. Just as compression of a gas can be a nearly reversible process, *α* can be regarded as a nearly reversible, elastic compressibility, because the pore volume change described by (3.13) arises largely from localized elastic shear straining of grain asperities that exist where irregular grains contact one another [55,67]. On the other hand, (3.13) implies that *α* cannot be a constant. Indeed, if *α* were constant and *m*/d*t*=*mα*(d*σ*_{e}/d*t*), indicating that exponential growth of *m* would occur in response to any growth of *σ*_{e}.

Existing data provide guidance for evaluating variation of *α* as a function of *m* and *σ*_{e}. Gravity-driven consolidation of unconfined, nearly liquefied debris-flow slurries with *m*≥0.4 generally exhibits behaviour consistent with *α*≈10^{−5} Pa^{−1} [61,68]. Compressibilities as large as *α*∼10^{−2} Pa^{−1} have been measured for relatively dilute, mud-rich slurries and dredged sludges with *σ*_{e}<10^{3} Pa and *m*<0.4 [69–71], whereas values *α*∼10^{−7} Pa^{−1} are more typical of loosely packed sand and sand–gravel mixtures with *σ*_{e}>10^{4} Pa and *m*>0.55 [72,73]. For virtually all soils and debris slurries, values of *α* increase monotonically as values of *σ*_{e} and *m* decrease (cf. [5,74]), although *α* must approach some upper limit when *σ*_{e} is sufficiently small. To reconcile these observations with (3.13), we postulate the constitutive relationship
*a* is a proportionality coefficient and *σ*_{0} is a reference stress that establishes the maximum debris compressibility, *a*/*mσ*_{0}, which applies when *σ*_{e}=0.

The implications of our constitutive relationship for *α* can be seen most clearly by combining (3.14) with (3.13) and (3.12), and then integrating the resulting equation for a special case with *m*/∂*γ*=0 (i.e. zero shear rate or zero dilatancy). This operation yields the solution
*m*=*m*_{min} when *σ*_{e}=0. Graphs of (3.15) can be compared with data obtained in laboratory compression and slurry-consolidation tests of diverse debris-like materials (figure 3). Because none of the tests involved gravity-free conditions in which an ideal limiting state with *m*=*m*_{min} and *σ*_{e}=0 was attained, we extrapolated the data for *m* as a function of *σ*_{e} in order to estimate the values of *m*_{min} and *σ*_{0} necessary for normalizing the data in figure 3 (table 1). Use of alternative values of *m*_{min} and *σ*_{0} would not affect the slopes of the data trends depicted in the figure, however. These trends show that the most-plausible values of *a* in (3.15) range from about 0.01 to 0.05, irrespective of the debris composition. The largest values of *a* as well as the smallest values of *m*_{min}and *σ*_{0} apply to mixtures with the highest mud contents. Owing to the ability of (3.15) to mimic the data trends in figure 3, we adopt (3.14) as our definition of debris compressibility, and we infer that values 0.01≤*a*≤0.05 and 10 *Pa*≤*σ*_{0}≤1000 Pa in (3.14) are commonly suitable. Model predictions presented in our companion paper [46] help test the validity of this inference.

### (f) Definition of dilatancy

The dilatancy of granular materials is traditionally expressed as an angle, as discussed in detail by Rowe [77], Nemat-Nasser [78], Bolton [23], Houlsby [24] and Rao & Nott [79]. In our model, evolving values of the dilatancy angle *ψ* are not specified. Instead they are the by-product of coupled evolution of the solid volume fraction, flow velocity, flow thickness and pore-fluid pressure. Nevertheless, *ψ* can also be interpreted as a state-dependent material property defined by
*V*/*V* is the part of the natural volume strain caused by an increment of shear strain ∂*γ* (figure 2). This definition of *ψ* is similar to that used in theories of quasi-static plastic deformation of rocks and soils [55], but our model avoids explicit evaluation of strains and instead combines (3.16) with (3.12) and (3.13) to recast the definition of *ψ* in terms of rates
*ψ* is related to the total dilation rate (−1/*m*)d*m*/d*t* and the portion of the total dilation rate caused by changes in the mean effective stress, *α*(d*σ*_{e}/d*t*).

Some key implications of (3.17) can be illustrated by considering special cases. For example, if *α*=0 or d*σ*_{e}/d*t*=0, then (3.17) reduces to *m* would occur in the presence of constant *ψ* and *ψ* must evolve as values of *σ*_{e} and *m* evolve during deformation. For experiments in which steady values of *m* are imposed, (3.17) implies that *ψ* evolves during a start-up phase in which *σ*_{e} evolves until d*σ*_{e}/d*t*=0 and *ψ*=0 are satisfied. This equilibration of *σ*_{e} to an imposed *m* is tantamount to development of steady dispersive normal stress in the enclosed shear cell experiments of Bagnold [80] (cf. [81,82]). On the other hand, if a steadily shearing grain–fluid mixture freely dilates or contracts when a step change in *σ*_{e} is externally imposed, then (3.17) implies that *m*/d*t*→0 and *ψ*→0. Boyer *et al.* [39] observed this type of equilibration in experiments more relevant than those of Bagnold [80]—at least with respect to gravity-driven geophysical flows—because Boyer *et al.* [39] employed a stress-controlled shear cell that allowed unrestricted evolution of *m* and attendant pore-fluid flow.

Although (3.17) is consistent conceptually with broad patterns of behaviour observed in experiments, it gives an incomplete description of dilatancy. It provides no relationship between instantaneous values of *ψ*, *m* and *σ*_{e}, and it contains no information about behaviour in the static limit *ψ* can be viewed as a material property that depends on the evolving difference between the ambient solid volume fraction *m* and a volume fraction *m*_{eq} that is equilibrated to the ambient confining pressure and shear rate. Dilatancy can thereby exhibit memory effects as well as dependence on the instantaneous state of the material (cf. [84]). This approach effectively extends concepts that were originally developed in critical-state models of quasi-static soils [22], and we employ it in our model.

To evaluate the instantaneous relationship between *m*−*m*_{eq} in disequilibrium states, we adopt a linear formula used by Roux & Radjai [85] and Pailha & Pouliquen [38], who inferred that tan *ψ*=*K*(*m*−*m*_{eq}), where *K* is a calibration factor. Moreover, to limit the number of adjustable parameters in our model, we tentatively set *K*=1, and thereby define
*ψ* generally falls in the range −0.2<*ψ*<0.2 (radians) in our model, because debris-flow materials typically have 0.4≤*m*≤0.8 and *m*_{eq}≈0.6.

To gauge the combined effects of the stress state and shear rate on *m*_{eq}, we use a dimensionless parameter similar to the ‘viscous number’ *et al.* [39] (cf. [86,87]). However, because our model allows for the development of fully liquefied states with *σ*_{e}=0, we use a generalized dimensionless parameter, *N*, which has a value that remains finite in such states
*μ* is the effective shear viscosity of muddy pore fluid, and *δ* is a length scale associated with generation of normal stresses by grain collisions during shearing (cf. [80,88]). For flows of geological debris containing a great diversity of grains, the value of *δ* is not well constrained. Therefore, in our model we simply set *δ*=0.001 m to match the size of typical sand grains.

Physical interpretation of *N* is facilitated by noting that (3.19) can be recast as
*S* is a Stokes number defined as *ρ*_{s}=2700 kg m^{−3}, *δ*=0.001 m, 0.001≤*μ*≤0.1 Pa-s and 10≤*σ*_{e}≤10^{6} Pa apply, the values 0≤*I*_{v}≤0.1 and 0<*S*<30 apply in (3.20). These values indicate that the full range of typical *N* values is 0≤*N*≤0.1, and that the range of *N* commonly mirrors that of *I*_{v}. The only significant difference in the values of *N* and *I*_{v} arises if *σ*_{e}<10 Pa, in which case *N*→1/*S* as *σ*_{e}→0 and *N*=1/*S* is an upper bound value that applies if grain motion is almost fully damped by fluid viscosity. (The mathematical singularity at *S*=0 is physically irrelevant unless there is perfect damping, which cannot occur in moving debris described by our model.)

The function *m*_{eq}(*N*) can be inferred from experimental data and from physical constraints that imply behavioural limits. In the static limit with *I*_{v}=0 and *N*=0, the value *m*_{eq}=*m*_{crit} applies, where *m*_{crit} is the critical-state volume fraction defined for lithostatically stressed states in normally consolidated soils [22]. Few direct measurements of *m*_{crit} for poorly sorted geological debris exist, but limited data indicate that a reasonable range is 0.5≤*m*_{crit}≤0.7 [16]. The relationship between *m*_{crit}, *m*_{eq} and *I*_{v} for non-static cases has been revealed most clearly by Boyer *et al.* [39], who found that the function
*N*≈*I*_{v} with (3.21) with (3.18) to pose what we believe is the simplest plausible constitutive relationship for state-dependent debris-flow dilatancy
*m* tend to relax towards an equilibrium value *m*_{eq} that depends on *N*. If *N*=0, then *m*_{eq}=*m*_{crit} applies. The relaxation time for *m* can be evaluated precisely for a special case in which *σ*_{e} and *N* are constants [38]. In this case, combination of (3.17) and (3.22) yields
*m*(0)=*βm*_{crit}, where *β*>1 describes initial states denser than the static critical state, and *β*<1 describes initial states looser than the static critical state. This solution
*m* approaches *m*_{eq} in a manner that depends on both *β* and *N*, but in every case the approach occurs with a characteristic relaxation time *m*_{eq} values of roughly 0.6 and ^{−1} [5], (3.24) implies that the relaxation *m*→*m*_{eq} typically would occur in less than 2 s if *N* were constant and no pore fluid were present. However, in our model relaxation of *m* is strongly regulated by evolution of pore-fluid pressure, as described below.

### (g) Definition of pore-fluid viscosity and solid–fluid drag

On the basis of diverse evidence summarized by Iverson [5], we define the debris-flow fluid phase as water plus persistently suspended clay and silt-sized particles (collectively called mud particles) with effective diameters smaller than about 62 *μ*m. In most natural debris flows and in large-scale experiments, the mud suspension has inferred solid volume fractions less than 0.2 [10,57]. These relatively dilute mud suspensions generally have effective viscosities that exceed that of pure water (approx. 0.001 Pa-s) by a factor ranging from about 5 to 100, whereas suspensions with solid volume fractions >0.4 can have viscosities that exceed that of water by a factor more than 1000 [3,89]. Furthermore, as the mud volume fraction increases, nonlinear rheological effects become increasingly prominent, particularly if mud particles consist of colloidal or chemically reactive clay [4,90,91]. Here, for the sake of parsimony, we neglect possible nonlinearities and treat pore fluid as an incompressible, Newtonian viscous material with a constant shear viscosity *μ*. Moreover, we assume that on a bulk scale, the liquid's shear deformation rate can be represented adequately by a deformation-rate tensor *R* that applies to the solid–fluid mixture as a whole. We therefore use the fluid-phase constitutive equation [50]
*k*(*m*) is the hydraulic permeability of the granular aggregate, and *p*_{e} is the excess pore-fluid pressure defined in (3.9). For the function *k*(*m*), we use the empirical formula
*k*_{0} in (3.27) has values that range from about 10^{−10} m^{2} for debris with approx. 2% mud content to 10^{−13} m^{2} for debris with about 50% mud content by dry weight. Large volumes of rapidly shearing debris are probably characterized by larger values of *k*_{0} owing to the effects of agitation that opens transient pathways for fluid flow, and also to the effect of scale on the effective permeabilities of porous media [92]. However, no measurements have been performed on large volumes of agitated, poorly sorted debris. Therefore, we infer that *k*_{0} may have values larger (but not smaller) than those portrayed in figure 5.

Because values of *k* are probably no larger than about 10^{−7} m^{2} in realistic, agitated debris, and because values of *μ* in muddy pore water are about 10^{−2}−10^{−1} Pa-s, we infer that the excess pore-pressure gradients described by (3.25) can readily have magnitudes comparable to those of hydrostatic gradients (∇*p*_{e}∼*ρ*_{f}*g*∼10^{4} Pa m^{−1}). Gradients of at least this size occur if the magnitude of ** q** exceeds 10

^{−1}m s

^{−1}in debris-flow materials with the largest

*k*/

*μ*values and if it exceeds 10

^{−8}m s

^{−1}in debris-flow materials with the smallest

*k*/

*μ*values. Development of excess pore-pressure gradients can inhibit debris dilation and contraction, and these gradients serve as a volume-averaged surrogate for the effects of local lubrication forces where individual grains displace fluid in adjacent pores [90,93,94]. In turn, inhibition of dilation or contraction produces feedbacks that can lead to regulation of ∇

*p*

_{e}by diffusive redistribution of excess pore-fluid pressure.

### (h) Diffusive redistribution of excess pore-fluid pressure

The relationship between the excess pore-pressure gradient described by (3.26) and the dilation rate described by (3.17) implies that pore-pressure evolution is mathematically analogous to forced diffusion (cf. [95]). This analogy can be established by first combining (3.10), (3.11) and (3.26) to find that
*m* can be eliminated from (3.28) by combining it with (3.17) to obtain
*k*/*αμ* plays the role of a pore-pressure diffusivity. Finally, definitions (3.8) and (3.9) can be used to infer that d*σ*_{e}/d*t*=d*σ*/d*t*−d*p*_{e}/d*t*−d*p*_{h}/d*t*, and substitution of this relationship in (3.29) enables the equation to be recast as a forced, advection–diffusion equation for *p*_{e},
*t*=∂/∂*t*+*v*_{s}⋅∇. If all time derivatives are zero and

The forcing effects described by the right-hand side of (3.30) can drive rapid pore-pressure change, but *p*_{e} can also evolve in the absence of forcing owing to diffusion described by the left-hand side of (3.30). A standard diffusion-equation normalization shows that *p*_{e} relaxes diffusively with a time-scale *h*^{2}*αμ*/*k*, where *h* is the debris-flow thickness (i.e. the distance over which pore-pressure diffusion typically occurs). For the range of *k*/*μ* values considered above (10^{−12}−10^{−5} m^{3} kg^{−1} s) and the representative value *α*=10^{−5} Pa^{−1}, *h*=1 m implies that timescales of unforced pore-pressure diffusion can range from roughly 1 s to 30 years. Thus, in an extreme case involving flow of gravel and water with little finer sediment, excess pore pressure may dissipate diffusively on a timescale much shorter than that of a flow duration. In the opposite extreme (e.g. a large, clay-rich volcanic debris flow), pore-pressure dissipation may proceed so slowly that its effects are modest until long after downstream motion has ceased. As a result of the dependence of the diffusion timescale on *h*^{2}, large debris flows can readily maintain higher fluid pressures, lower Coulomb friction, and greater mobility than can small flows with the same composition [5].

### (i) Coulomb friction

An abundance of evidence indicates that Coulomb friction generates most of the resistance to debris-flow motion—particularly if Coulomb behaviour is defined as a proportionality between shear stresses and normal stresses on planes of shearing, irrespective of the shear rate [11]. With this definition, Coulomb friction characterizes the shear-to-normal stress ratio in most quasi-static soils [96], in dense granular flows with enduring grain contacts [25], in more dilute granular flows characterized by brief grain collisions [97], and in shear flows of concentrated solid–liquid mixtures [39,80]. Indeed, on a bulk continuum scale, a proportionality between normal and shear stresses on planes of shearing appears to be a nearly universal property of granular materials.

In our model, the Coulomb friction rule represents the influences not only of grain-boundary friction but also of dilatancy and modification of normal stresses at grain contacts by pore-fluid pressure. On the basis of well-established principles of soil mechanics and empirical evidence [98,99], we use the effective stress *T*_{e} defined in (3.7) as the stress that influences bulk frictional resistance to shearing. Thus, if all else is constant, frictional resistance decreases as pore pressure locally increases. Dilatancy *ψ* influences bulk shear resistance through geometrical effects on grain rearrangement. Following the rationale of Rowe [77], Nemat-Nasser [78], Bolton [23] and Houlsby [24], we infer that the effects of *ψ* contribute additively to those of a constant-volume friction angle *ϕ*, which characterizes friction under conditions with *m*=*m*_{eq}. By embedding this effect of *ψ* and also embedding the effect of pore pressure, we express the Coulomb friction rule as
*τ*_{s} is the granular-phase shear stress vector acting on planes of shearing, *σ*_{n} is the total normal stress on these planes, and ** v**/∥

**∥ indicates that positive shear stresses resist motion. The Coulomb rule can be expressed more completely in three-dimensional tensor form [100], but this complication is unnecessary here because we do not evaluate three-dimensional stress states in our depth-averaged model. A further complication is that Coulomb friction may exhibit rate-dependence, manifested as a gradual increase of**

*v**ϕ*as a function of

*I*

_{v}[39]. Here, for the sake of parsimony, we neglect this effect, although (3.31) implies that granular shear resistance exhibits some implicit rate-dependence owing to the effects of the shear rate on

*ψ*and

*p*.

## 4. Depth-averaged model equations

Our depth-averaged model describes coupled evolution of the flow thickness *h*(*x*,*y*,*t*), solid volume fraction *m*(*x*,*y*,*t*), basal pore-fluid pressure *p*_{b}(*x*,*y*,*t*) and depth-averaged components of the mixture velocity, which are *x*-direction and *y*-direction. The *z* component of the mixture velocity, *w*, is not predicted explicitly, but it is related to *h*, *u*(*h*) and *v*(*h*) through the kinematic boundary conditions
*z*=*h*(*x*,*y*,*t*) and at a rigid basal boundary, where *z*=0.

The necessity of specifying the direction of depth integration *a priori* makes the choice of a coordinate system crucial. Some depth-averaged models use an Earth-centred, orthogonal Cartesian coordinate system with *z* vertical, which has the advantage of being independent of terrain geometry. This approach leads to complicated mechanical considerations when accounting for the effects of steep, irregular topography, however [101]. Other models, including the one we present here, use a *z*-coordinate normal to the local ground surface, such that the *x*-coordinate is directed downslope and the *y*-coordinate is directed cross-slope (figure 6). This approach, also, has limitations, because curvilinear coordinates are needed to adapt it to irregular terrain, limiting the utility of conventional digital elevation models of Earth's topography. Here, to emphasize the main new features of our model, we omit the effects of bed curvature and focus on flows that traverse planar terrain. Our companion paper [46] addresses bed curvature effects.

### (a) Mass balances

Our model includes two mass-balance equations: one for the solid–fluid mixture as a whole and one for the granular solid phase. Because no data exist that document variations of the bulk density *ρ* or solid volume fraction *m* as a function of *z* in realistic debris flows, the equations assume that *ρ*(*x*,*y*,*t*) and *m*(*x*,*y*,*t*) are uniform at all depths below the surface at *z*=*h*(*x*,*y*,*t*). However, the definition of *h* itself involves a potential ambiguity: because rocks may protrude variable distances above the adjacent liquid at the surface of debris flows, solid and liquid volume fractions may be ill defined in the vicinity. To avoid this ambiguity, we define *z*=*h*(*x*,*y*,*t*) as the position of a *virtual* free surface, such that the debris-flow mass per unit basal area Δ*x*Δ*y* is given by *ρh*(*x*,*y*,*t*) (figure 7). In effect, this definition replaces some combination of solid or liquid mass immediately above or below *h*(*x*,*y*,*t*) with a thin, equivalently massive, homogeneous layer with a density *ρ* and upper surface at *h*(*x*,*y*,*t*). In this way, *ρh*(*x*,*y*,*t*) is conserved, despite the fact that solid or liquid constituents may pass through *h*(*x*,*y*,*t*) during dilation or contraction.

The depth-integrated dilation rate *D*(*x*,*y*,*t*) plays a key role in our mass-balance equations. It is defined by integrating (3.10) through the flow thickness to obtain
*Dm*/*h* accounts for the effects of dilation or contraction on *m*.

A depth-averaged mass-conservation equation that describes evolution of *ρh* for the solid–fluid mixture as a whole results from integrating (3.1) using Leibniz’ rule and the kinematic boundary conditions (4.1) to find that
*ρ*=*ρ*_{s}*m*+*ρ*_{f}(1−*m*) in conjunction with (4.3) to obtain
*ρ*_{f}=0 (which implies that *ρ* and *h* describe attributes of only a granular phase) or with *ρ*_{f}=*ρ* (which implies that *ρ* and *h* are uninfluenced by dilation or contraction because the debris effectively behaves as a homogeneous, incompressible material). In the more-general case with 0<*ρ*_{f}<*ρ*, (*ρ*−*ρ*_{f})*D* describes the simultaneous effects of a source *ρD* and sink −*ρ*_{f}*D* as fluid mass replaces solid mass in dilating debris beneath *h*. In this case, the surface *h*(*x*,*y*,*t*) coincides with neither the fluid surface nor the solid surface, but rather with that of the virtual free surface described above (figure 7).

Use of (4.6) also enables (4.3) to be recast in a form that accounts for simultaneous evolution of *m* and *h*. Multiplication of (4.6) by *m* and multiplication of (4.3) by *h*, followed by addition of the resulting equations, yields the depth-averaged equation we use to evaluate evolution of the solid volume fraction

The presence of finite source terms in (3.7) and (3.8), even for the special case with *ρ*_{s}=*ρ*_{f}=*ρ*, reflects the fact that *ρ*_{s}*hm* and *ρ*_{f}*h*(1−*m*) are not conserved beneath the virtual free surface at *z*=*h*(*x*,*y*,*t*). Rather, *ρh* is the quantity conserved beneath this surface, and *ρh*=*ρ*_{s}*hm* or *ρh*=*ρ*_{f}*h*(1−*m*) are valid only if *ρ*_{f}=0 or *ρ*_{s}=0, respectively.

### (b) Momentum balances

Depth-integrated equations for conservation of the *x*, *y* and *z* components of linear momentum are exactly analogous to one another. Therefore, we illustrate integration of only the *x* component, which is accomplished by using (3.2) and Leibniz’ rule to obtain the left-hand side terms
*x*-momentum through the virtual free surface at *h*(*x*,*y*,*t*) in response to dilation or contraction, but it does not imply that momentum enters or leaves the debris flow. Rather, the term reflects the fact that momentum is not conserved beneath *h*(*x*,*y*,*t*) unless *D*=0 or *ρ*_{f}=*ρ*_{s}=*ρ*. Although (4.9) can be written in alternative forms that lack *D*, we use the last line of (4.9) to illustrate the role of *D* and to clarify the relationship between our equations and standard, depth-integrated momentum equations for flows of incompressible material.

Depth integration of the *x* component of the right-hand side of (3.2) using Leibniz’ rule yields
*g*_{x} is the *x* component of ** g**, and

*τ*

_{xx}and other normal stress components as positive in compression). Except for the resisting basal shear traction

*τ*

_{zx}(0), all boundary terms introduced during depth integration vanish from (4.10), because we assume a traction-free upper boundary at

*z*=

*h*and rigid basal boundary at

*z*=0.

Substitution of (3.9) and (4.9) into (3.2) yields the depth-averaged *x*-component momentum-conservation equation
*y*-component momentum equation is obtained by interchanging *x* and *y* in (4.11), yielding
*z*-momentum equation can be written in an Eulerian form analogous to that of (4.11) and (4.12), but it is useful to manipulate its left-hand side by employing (4.6) to obtain the Lagrangian form

If shallow-flow scaling (detailed in §5) is used to justify neglect of *z* component of the effective weight per unit area of the moving debris, *ρg*′_{z}*h* [101]
*z* is reckoned positive upward (such that *g*_{z}<0), a positive acceleration *g*_{z} and thereby increases the effective weight of the debris and the basal normal traction *τ*_{zz}(0). Lack of precise accounting for the effects of *θ* is bed slope angle (figure 6), *g* is the magnitude of ** g** and

*θ*on flow dynamics.

### (c) Basal pore pressure

We use depth integration and rational approximations to obtain a reduced pore-pressure evolution equation that contains the basal pore pressure *p*_{b} rather than *p*_{e}. The first step entails recasting (3.30) in terms of the total pore-fluid pressure, *p*=*p*_{e}+*p*_{h}, where *p*_{h}=*ρ*_{f}*g*_{z}(*h*−*z*) is the hydrostatic pressure that balances the fluid weight. Next, we invoke shallow-flow scaling that applies if *H*/*L*≪1 (see §5), where *H* is the characteristic thickness and *L* is the characteristic length of a debris-flow surge (figure 6). This scaling indicates that ∂^{2}*p*/∂*z*^{2}≫∂^{2}*p*/∂*x*^{2}, ∂^{2}*p*/∂*y*^{2} because ∂^{2}/∂*z*^{2} scales with 1/*H*^{2}, whereas ∂^{2}/∂*x*^{2} and ∂^{2}/∂*y*^{2} scale with 1/*L*^{2}. Analogous reasoning indicates that ∂/∂*x*, ∂/∂*y*∼1/*L* and ∂/∂*z*∼1/*H*. Using these scales to identify and neglect terms of order *H*/*L* or smaller in (3.30) yields the reduced equation
*w*≈(*z*/*h*)d*h*/d*t*, *k*/*μ*)/∂*z*=0 (because we assume that no material property varies as a function of *z*), and thereby reduce (4.16) to
*p*(*h*)=*σ*(*h*)=0, yielding
*k*/*αμ*)*ρ*_{f}*g*_{z} arises in (4.18) from depth integration of the pore-pressure diffusion term in (4.17) and application of a zero-flux basal boundary condition stipulating that the pore-pressure gradient at the bed is hydrostatic: [∂*p*/∂*z*]_{z=0}=−*ρ*_{f}*g*_{z}. The term *σ*−*p*)/∂*z* in (4.17) by parts. This term cancels some other terms and thereby reduces (4.18) to
*p* and *p*_{b}.

To express (4.19) in terms of *p*_{b}, we find substitutions for *p*/∂*z*]_{z=h} that are consistent with our assumption that *m* is uniform at all depths below *z*=*h*. This condition implies that ∇⋅*v*_{s} and ∇⋅** q** are not functions of

*z*, further implying that ∂

^{2}

*p*/∂

*z*

^{2}is not a function of

*z*in (4.16) and (4.17). With this restriction, we solve ∂

^{2}

*p*/∂

*z*

^{2}=const. and employ the zero-flux basal boundary condition ([∂

*p*/∂

*z*]

_{z=0}=−

*ρ*

_{f}

*g*

_{z}) and pressure-free surface boundary condition

*p*(

*h*)=0 to find that

*p*(

*z*,

*t*) satisfies the quadratic equation

*p*(

*z*,

*t*) is represented by evolving values of the basal pressure

*p*

_{b}(

*t*) and flow thickness

*h*(

*t*) (figure 8). Integration and differentiation of (4.20) show that the equation also implies that

*N*in (4.24) employs the depth-averaged shear rate estimate

A direct connection between *p*_{b} and the dilation rate *D* can be established by using (3.28), (4.2), and the pore-pressure distribution described by (4.20) to find that
*p*_{b} influences all of the other dependent variables in the model, because each of the depth-averaged evolution equations (4.6), (4.7), (4.11) and (4.12) contains a source term involving *D*. This influence can be viewed in a different light by combining (4.25) with (4.22) and (3.26) to obtain
*q*_{z}|_{z=h} is the specific discharge (positive upward) of fluid moving through the surface at *z*=*h* during dilation or contraction (figure 7). Thus, in our model the dilation rate, excess basal pore pressure and fluid flux through the virtual free surface are different manifestations of the same phenomenon.

The relationships expressed by (4.25) and (4.26) are contingent on our assumption that all debris, except that passing through the virtual free surface, remains fully saturated with liquid. This assumption would be violated if debris were to dilate so rapidly that the liquid phase cavitates. Energy requirements for creation of liquid–gas interfaces in porous media tend to mitigate against this behaviour [55], but they do not preclude it. Our model does not attempt to simulate air entry and instead enforces the restriction that *p*_{b}≥0, which, according to (4.25), effectively limits *D* to the range *D*≤(2*k*/*μ*)*ρ*_{f}*g*_{z}. The model implies that small negative fluid pressures (i.e. tension-saturated zones) can develop internally as *p*_{b}→0, however (figure 8).

### (d) Stress estimation

Our estimates of stress components in (4.11) and (4.12) assume that Coulomb friction in the granular fraction of the debris is fully engaged if driving forces per unit area (*ρg*_{x}*h* and *ρg*_{y}*h* in (4.11) and (4.12)) equal or exceed the sum of the other right-hand side terms in (4.11) and (4.12). If *ρg*_{x}*h* and *ρg*_{y}*h* do not suffice to fully engage Coulomb friction, then the debris is static and the stress state is indeterminate—as it is in classical limit-equilibrium methods of slope-stability analysis. To obtain a formula for the fully engaged basal shear traction, we combine (3.31) with (4.16) to model the Coulomb component of the traction, and then add a term derived from (3.25) to account for depth-integrated viscous shear resistance of fluid in the volume fraction 1−*m*. The resulting formulae for the basal shear tractions in the *x*- and *y*-directions are
*θ*_{x} and *θ*_{y} are the *x* and *y* components of the slope angle *θ*, and *x* and *y* components of the fluid shear rate. In (4.27) and (4.28), the leading factors *τ*_{zx}(0) and *τ*_{zy}(0) resist motion in the direction given by the vector sum of

Typically, the viscous terms in (4.27) and (4.28) account for little of the total basal shear traction. For example, in a scenario with hydrostatic basal pore pressure and the values *h*=1 m and *p*_{b}=*ρg*_{z}*h*). Dissipation associated with inelastic grain collisions may contribute additional flow resistance. For the sake of parsimony, we do not include this effect in our model, although experiments indicate that the effect can be mimicked by treating *ϕ* as a rate-dependent quantity [39].

We represent the solid- and fluid-phase contributions to the lateral stress gradient terms *τ*_{xx} and *τ*_{yy} are each proportional to *τ*_{zz}, leading to the expressions
*ρg*_{z}*h*−*p*_{b} account for lateral effective normal stresses acting in the granular phase, whereas the terms *h*(∂*p*_{b}/∂*x*) and *h*(∂*p*_{b}/∂*y*) account for fluid pressure contributions to the lateral stress. Additional terms can be added to (4.29) and (4.30) to account for depth-averaged deviatoric fluid normal stresses (cf. [27]), but for plausible values of the fluid viscosity *μ* these terms are negligible, and we exclude them here.

The granular components of the lateral normal stress terms in (4.29) and (4.30) involve a lateral pressure coefficient *κ*, which has been discussed at length previously [27,102–104]. For an ideal fluid, *κ*=1 applies, and in this case (4.29) and (4.30) reduce to conventional forms used in shallow-water theory. Similarly, if a fully liquefied state exists (i.e. *p*_{b}=*ρg*_{z}*h*), then (4.29) and (4.30) reduce to the conventional shallow-water form. On the other and, for a uniform slab of unliquefied Coulomb material that undergoes elongation or compression as well as basal sliding, theory indicates that *κ* has one of two values given by (see [5], for a derivation)
*ϕ* is the basal friction angle and *ϕ*_{int} is the friction angle that describes Coulomb resistance to internal deformation. These two friction angles may or may not be the same, but *ϕ*≤*ϕ*_{int} is necessary to ensure that basal sliding accompanies internal deformation. The larger of the two values of *κ* described by (4.31) is associated with compressive longitudinal deformation and always exceeds 1, whereas the smaller *κ* value is associated with extensional deformation and is commonly smaller than 1. Our model includes the option of using these variable *κ* values or, as suggested by the findings of Gray *et al.* [103,105], setting *κ*=1 as a simplifying approximation.

## 5. Normalization and dimensionless parameters

We normalize the depth-averaged model equations by using scales appropriate for a debris flow of finite length *L* and thickness *H* (figure 6). Because downslope debris-flow motion is driven by gravitational potential, the scale for the velocity components in the *x*- and *y*-directions is (*gL*)^{1/2}, whereas the scale for the *z*-direction velocity component is (*gH*)^{1/2} [102,106]. The scale (*gH*)^{1/2} also applies for *D*, because *D*≠0 indicates relative motion of solid and liquid phases in the *z*-direction. The length-scale *L* divided by the downslope velocity scale yields the timescale for downslope debris-flow motion, (*L*/*g*)^{1/2}. The scale for *ρ* is the initial static debris bulk density *ρ*_{0} associated with the initial solid volume fraction *m*_{0}:*ρ*_{0}=*m*_{0}*ρ*_{s}+(1−*m*_{0})*ρ*_{f}. The value *m*_{0} also serves as the scale for *m*. Finally, the scale for all depth-averaged stress components *p*_{b} is *ρ*_{0}*gH*. Use of these scales leads to definition of the following dimensionless quantities, denoted by asterisks:

The scaled mass-conservation equations (4.6) and (4.7) are
*ε*≪1 provides a basis for comparing the relative magnitudes of terms in these scaled equations. The stress gradient terms in (5.2)–(5.4) are of order *ε*^{1}, and consequently they may be negligible in many circumstances. The bed-normal acceleration term involving *ε*^{1/2}. All terms that contain *D** are effectively of order *ε*^{0}, because the factors *ε*^{1/2} and *ε*^{−1/2} in these terms cancel one another if (5.9) is substituted into (5.2), (5.3) and (5.6)–(5.8). Thus, no term containing *D** can be neglected.

Although shallow-flow scaling indicates that stress gradient terms are small, modelling surge-like motion of debris flows requires retention of the normal stress gradient terms *ε*, however—thereby distinguishing our approach from that used in some previous debris-flow and avalanche models [27,29,101]. In our model, the terms containing *D** play a more essential role, and we aim to evaluate this role without juxtaposing unnecessary complications.

The scaled model equations contain relatively few dimensionless parameters. Three of these parameters, *g*_{x}/*g*, *g*_{y}/*g* and *g*_{z}/*g*, reflect the extrinsic influence of the local slope angle and orientation, which are independent of the properties of the debris flow. The remaining six parameters express the influence of intrinsic debris-flow properties:
*H*=0.1 m and *L*=100 m) and large-scale flows (with *H*=10 m and *L*=10^{4} m).

The parameter *k*(*L*/*g*)^{1/2}/*αμH*^{2} plays the most important role in affecting model predictions of flow behaviour, because it has values that can vary by many orders of magnitude, even in flows of fixed size (table 2). (Values of the parameter *μ*/*ρ*_{0}*H*(*gL*)^{1/2} also can range widely, but they are universally small, indicating that viscous shearing plays a relatively minor role in resisting motion.) Model predictions also depend strongly on *m*_{crit}, because the value of this parameter establishes whether the initial state of debris with *m*=*m*_{0} is dilative (*ψ*>0 if *m*_{0}>*m*_{crit}) or contractive (*ψ*<0 if *m*_{0}<*m*_{crit}). The sign of *m*_{0}−*m*_{crit} consequently determines the sign of pore-pressure changes that occur during the first stages of downslope debris motion. Debris-flow size also has an important influence, because as flow size increases, values of *k*(*L*/*g*)^{1/2}/*αμH*^{2} decrease, whereas those of *αρ*_{0}*gH* increase. In summary, the behaviour of model predictions depends largely on topography, initial conditions and the value of *k*(*L*/*g*)^{1/2}/*αμH*^{2}.

The parameter *k*(*L*/*g*)^{1/2}/*αμH*^{2} can be viewed as a time-scale ratio, in which (*L*/*g*)^{1/2} is the timescale for downslope debris-flow motion and *αμH*^{2}/*k* is the timescale for pore-pressure relaxation to a hydrostatic state (cf. [54]). Within these timescales, the greatest source of variability arises from the pore-pressure diffusivity, *k*/*αμ*, which depends chiefly on the grain-size distribution and degree of debris dilation. A very broad range of diffusivities (10^{−7}–10^{2} m^{2} s^{−1}) results from considering the full ranges of plausible values of *k* (10^{−13}–10^{−7} m^{2}), *α* (10^{−7}–10^{−4} Pa) and *μ* (10^{−3}–10^{−1} Pa-s) for debris ranging from mud-rich and gravel-poor to gravel-rich and mud-poor. The typical range of *k*/*αμ* values may be considerably smaller, however. Values of *k*/*αμ* determined experimentally for quasi-static debris-flow materials with compositions dominated by mud or by sand and gravel have ranged only from about 10^{−7} to 10^{−4} m^{2} s^{−1} [5,61,68,72].

## 6. Discussion

Our model comprises a set of five simultaneous partial differential equations describing coupled evolution of *h*(*x*,*y*,*t*), *m*(*x*,*y*,*t*) and *p*_{b}(*x*,*y*,*t*). A concise recapitulation of these equations and ancillary constitutive formulae is provided in our companion paper [46].

Application of our model typically entails simulating debris-flow motion that begins from a statically balanced initial state with *p*_{b}(*x*,*y*,*t*) are by definition too small to satisfy (4.27) or (4.28) together with (4.29) or (4.30). Motion can then be triggered either by gradually increasing *p*_{b} (simulating the effect of rainfall or snowmelt infiltration), gradually reducing the basal friction angle, *ϕ* (simulating the effects of rock weathering or decay of roots that help bind soil), gradually changing the slope geometry (simulating erosion or human intervention) or rapidly changing ** g** (simulating earthquakes). In all scenarios, motion begins when and where an infinitesimal force imbalance first develops. In this way, the model simulates the onset of slope instability and a transition to dynamic behaviour that may or may not entail rapid flow—contingent on dilatancy and pore-pressure feedback.

The model's representation of flow dynamics can be compared with that of some better-established models. For example, the depth-averaged evolution equations reduce to a form like that of the Savage–Hutter [102] granular avalanche model if the dilatancy, compressibility and pore-fluid viscosity and density are zero, implying that the restrictions *D*=0, *m*=*m*_{0}=*m*_{crit} and *p*_{b}=0 are enforced. On the other hand, if the restrictions *D*=0, *m*=*m*_{0}=*m*_{crit}, *ρ*=*ρ*_{f} and *ϕ*=0 apply, the equations provide a conventional description of shallow-water flow. In this case, the debris behaves as an incompressible liquid, and the pore-pressure evolution equation (4.24) reduces to the hydrostatic balance *p*_{b}=*ρg*_{z}*h*. Moreover, the longitudinal stress gradients given by (4.29) and (4.30) reduce to *ρg*_{z}*h*(∂*h*/∂*x*) and *ρg*_{z}*h*(∂*h*/∂*y*), and the basal shear tractions given by (4.27) and (4.28) contain only viscous components that represent fluid effects. In shallow-water models, these basal tractions can be approximated in a variety of other ways in an effort to account for the effects of turbulence.

Behaviour that temporarily or locally mimics shallow-water flow can occur as *D*, *m* and *ρ* evolve in our model, provided that *k*(*L*/*g*)^{1/2}/*αμH*^{2}≪1 and *m*<*m*_{eq} apply. In this case, contractive shear deformation produces a persistent state of near-liquefaction because debris-flow motion occurs much more rapidly than pore-pressure relaxation. A caveat, however, is that elongation and thinning of a moving debris flow can reduce the effective value of *H*, thereby increasing the speed of pore-pressure relaxation in proportion to the change in *H*^{2}. This effect causes frictional resistance to increase where the flow thickness becomes small. As thickness approaches zero at debris-flow snouts, for example, pore pressure dissipates readily and increased friction develops.

Reduction of pore pressure and consequent growth of flow resistance also can result from increased shear rates, because increased shear rates reduce the value of *m*_{eq} through the action of dilatancy. Indeed, the order-*ε*^{−1} term on the right-hand side of (5.8) indicates that a lowest order, steady-state approximation of the physical processes represented in our model consists of a balance between the effects of dilatancy and generation of excess basal pore pressure. This balance can be expressed in dimensional terms as
*m*<*m*_{eq} applies, then the steady state defined by (6.1) implies that *p*_{b}−*ρ*_{f}*g*_{z}*h*<0, indicating that the effects of positive dilatancy reduce pore-fluid pressure and help stabilize motion. Faster motion promotes pore-pressure reduction not only by increasing the magnitude of *m*_{eq}—as demonstrated in idealized experiments [82]. According to (6.1), the propensity for evolution of *m* to stabilize motion can also manifest itself during contractive behaviour with *m*<*m*_{eq} and *p*_{b}−*ρ*_{f}*g*_{z}*h*>0, because faster motion reduces *m*_{eq} and thereby reduces *p*_{b}.

In another special case, the model equations provide an approximate description of consolidation of quasi-static debris-flow deposits. In this case, all velocities vanish, and the pore-pressure evolution equation (4.24) reduces to
*D* given by (4.26). Equation (6.2) has a particularly simple solution if *h* and *ρ* are assumed to be constants, as is typically assumed in soil consolidation modelling
*ρ*_{f}*g*_{z}*h*, and it satisfies the initial condition *p*_{b}=*p*_{b}(0). The relaxation coefficient in (6.3) has a complicated form, but in general (*αρ*_{f}*g*_{z}*h*/2)[(*ρ*−*ρ*_{f})/*ρ*]≪3 applies for all but the largest debris flows. Therefore, a useful approximation of (6.3) is
*αμh*^{2}/*k*. In a depth-averaged sense, this exponential decay approximates pore-pressure relaxation predicted by classical soil consolidation theory (figure 9). As excess basal pore pressures relax, quasi-static debris becomes progressively more rigid.

The pore-pressure relaxation timescale has another important implication, which arises from its relationship to the relaxation rate of *m*. The relaxation *m*→*m*_{eq} has a characteristic time-scale *αμh*^{2}/*k* that typically ranges from a few minutes to many months. Owing to these disparate timescales and the strong coupling of *m* and *p*_{b}, the relaxation *m*→*m*_{eq} can proceed only as rapidly as the relaxation *p*_{b}→*ρ*_{f}*g*_{z}*h* permits. In this way, the model equations mirror experimental findings demonstrating that coupled evolution of pore space and pore-fluid pressure plays a dominant role in regulating debris-flow dynamics [16,65].

Complications that are not addressed by the model include several related to evolution of *m*. For example, d*m*/d*t*<0 might occur near the surface of a debris flow while d*m*/d*t*>0 occurs near the bed, yielding a commensurately variable pore-pressure response. Moreover, the value of *m*_{crit} might evolve in response to a variety of factors identified in quasi-static soil tests [10,23,108,109]. Our aim is to omit rather than include complicating factors wherever possible, however, and we therefore assume that *m*_{crit} is constant. In our companion paper, we show that such a parsimonious model can make useful predictions [46].

## 7. Conclusion

We have derived a depth-averaged debris-flow model aimed at seamlessly simulating all stages of flow behaviour, from initiation to post-depositional debris consolidation. The model formalizes the hypothesis that the evolving debris dilation rate, coupled to evolution of pore-fluid pressure, plays a primary role in regulating debris-flow dynamics. The model's representation of this role involves three key postulates. One postulate is that changes in the solid volume fraction *m* result from the interaction of the depth-averaged shear rate, dilatancy and effective stress. In turn, the evolving dilatancy angle *ψ* obeys tan *ψ*=*m*−*m*_{eq}, where the equilibrium solid volume fraction *m*_{eq} depends on the ambient stress state and shear rate. The second key postulate, interrelated with the first, is that a Darcy drag formula describes the effect of solid–fluid interactions on the relaxation of *m* towards *m*_{eq}. The third postulate is that flow resistance is dominated by basal Coulomb friction, which is affected by dilatancy and by pore-fluid pressure mediated by Darcy drag. The Darcy and Coulomb postulates are consistent with behaviour observed previously in replicable experiments [16,61,65]. Therefore, relationships involving dilatancy represent the primary new hypothesis embodied by the model. The fact that dilatancy produces a leading-order effect in the normalized model equations enhances the prospects for conclusive model tests.

The normalized model equations contain only one dimensionless parameter that varies significantly among debris flows with differing compositions and sizes, *k*(*L*/*g*)^{1/2}/*αμH*^{2}. Therefore, model predictions require little calibration. A typical debris-flow simulation might employ values of other dimensionless model parameters estimated as *ϕ*≈0.7, *ρ*_{f}/*ρ*_{0}≈0.6, *m*_{crit}≈0.6, *αρ*_{0}*gH*≈0.1 and *μ*/*ρ*_{0}*H*(*gL*)^{1/2}≈10^{−6}. If these values are fixed, then variations in predicted debris-flow behaviour depend only on variations in the value of *k*(*L*/*g*)^{1/2}/*αμH*^{2} and on the initial value of *m*−*m*_{crit} (in addition to extrinsic factors such as the geometry of the debris-flow source area and path). In stringent model tests, such as those performed in flume experiments, values of these parameters can be constrained by independent measurements of material properties [16,65]. In scenarios that lack such constraints, the initial value of *m*−*m*_{crit} and the value of *k*(*L*/*g*)^{1/2}/*αμH*^{2} may be used as calibration parameters or as stochastic variables in probabilistic forecasting.

The key role played by the time-scale ratio *k*(*L*/*g*)^{1/2}/*αμH*^{2} implies that debris-flow mobility predicted by our model involves inherent scale-dependence. For flows of varying size but fixed aspect ratio *H*/*L*, the motion timescale (*L*/*g*)^{1/2} grows as *L*^{1/2}, whereas the pore-pressure relaxation timescale grows as *H*^{2}. Thus, large flows can retain high pore pressures longer than can small flows of similar composition. Consequent reduction of effective basal friction can help explain the extraordinarily high mobility exhibited by many large debris flows [5].

The structure of our model implies that if contractive initial conditions with *m*<*m*_{crit} exist, then slope failure leads to positive pore-pressure feedback, making partial liquefaction and runaway debris-flow motion almost inevitable. On the other hand, dilative initial conditions with *m*>*m*_{crit} lead to negative pore-pressure feedback. This feedback may lead to slow or intermittent landslide motion [19], but it does not preclude debris-flow initiation. According to our model, however, debris-flow initiation that commences with dilative deformation is a relatively gradual process that does not involve abrupt liquefaction and runaway behaviour. Progressive destruction of soil aggregates during the early stages of shearing can promote this type of flow initiation by causing a transition from dilative to contractive behaviour [10]. This complication, as well as complications due to grain-size segregation and entrainment of boundary material, remains to be incorporated in our model.

- Received December 10, 2013.
- Accepted June 30, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.