## Abstract

We evaluate a new depth-averaged mathematical model that is designed to simulate all stages of debris-flow motion, from initiation to deposition. A companion paper shows how the model's five governing equations describe simultaneous evolution of flow thickness, solid volume fraction, basal pore-fluid pressure and two components of flow momentum. Each equation contains a source term that represents the influence of state-dependent granular dilatancy. Here, we recapitulate the equations and analyse their eigenstructure to show that they form a hyperbolic system with desirable stability properties. To solve the equations, we use a shock-capturing numerical scheme with adaptive mesh refinement, implemented in an open-source software package we call D-Claw. As tests of D-Claw, we compare model output with results from two sets of large-scale debris-flow experiments. One set focuses on flow initiation from landslides triggered by rising pore-water pressures, and the other focuses on downstream flow dynamics, runout and deposition. D-Claw performs well in predicting evolution of flow speeds, thicknesses and basal pore-fluid pressures measured in each type of experiment. Computational results illustrate the critical role of dilatancy in linking coevolution of the solid volume fraction and pore-fluid pressure, which mediates basal Coulomb friction and thereby regulates debris-flow dynamics.

## 1. Introduction

In a companion paper, we derive a new depth-averaged model aimed at seamlessly simulating debris-flow motion from initiation to deposition [1]. Here, we describe mathematical properties of the model, our numerical solution technique, and tests of numerical predictions against experimental data. Our model can be viewed in a broad context by recognizing that debris flows belong to a large class of gravity-driven mass flows in which a relatively thin layer of material with an evolving free upper surface moves over basal topography. Such phenomena include rock avalanches, snow avalanches, pyroclastic flows, lava flows and various types of water flooding (tsunamis, dam breaches, etc.). To simulate these phenomena, researchers have proposed many mathematical models, most of which are depth averaged. These models have the advantage of computational efficiency, due to their restricted dimensionality as well as inclusion of the free-surface height as a dependent variable in the governing equations. The classical depth-averaged model is the St. Venant or shallow-water model, which describes a pure fluid with a hydrostatic vertical pressure gradient (e.g. [2]). The shallow-water equations are hyperbolic conservation laws with source terms—an extensively studied class of mathematical problems—and they represent a special case of the more-general set of equations we use to model debris flows.

A feature distinguishing debris flows from many other shallow flows is the presence of roughly equal volumes of diverse solid grains and incompressible, viscous interstitial fluid that interacts with the solids. This interaction has a profound effect on the apparent rheological properties of the bulk mixture [3], but accounting for the interaction in depth-averaged mathematical models has proven challenging. Pitman & Le [4] developed a depth-averaged version of Jackson's [5] three-dimensional, two-phase granular–fluid equations, and various extensions of their model have been proposed (e.g. [6–8]). These models are expressed as hyperbolic systems of partial differential equations, yet they possess an unstable elliptic degeneracy associated with loss of hyperbolicity in particular flow regimes [6,9].

An alternative approach towards encapsulating granular–fluid interactions in depth-averaged models involves calculating the bulk mixture velocity (i.e. as if the mixture were a one-phase material) and accounting for the effects of phase interactions in the mixture stresses. Iverson [10] and Iverson & Denlinger [11] pursued this approach by including diffusive dissipation models for the pore-fluid pressure and advecting this pressure with the bulk mixture velocity. Solid stresses are then mediated by fluid pressure through the use of Terzaghi's effective stress principle, such that normal stresses at grain contacts are proportional to bulk normal stresses minus the fluid pressure. Nonetheless, these models do not fully account for the dynamic pressure response to inter-granular rearrangements accompanied by evolution of solid and fluid volume fractions.

Recently, Kowalski & McElwaine [12] presented a mixture model that parallels our present model because it includes feedback due to changes in solid and fluid volume fractions. Beginning with two-phase equations in three dimensions, they derive equations for the evolving and equilibrium volume fractions resulting from suspension or settling of solid grains. The fluid pressure component of the bulk mixture stress is algebraically coupled to the volume fraction, based on whether there is net particle settling or resuspension. The vertical structure of the solid volume fraction is then projected onto a vertical centre of mass during depth-averaging. The resulting equations employ a bulk mixture velocity, but because the centre of mass is modelled, the equations retain some information about the vertical structure of the evolving solid volume fraction. Therefore, as in our model, Kowalski & McElwaine [12] account for the coevolution of changing volume fractions and pore-fluid pressure.

A fundamental distinction between our model and that of Kowalski & McElwaine [12] involves the mechanisms regulating coevolution of volume fractions and fluid pressure; ours is based on granular dilatancy, whereas theirs is based on the vertical centre of mass as a measure of particle settling versus resuspension. Additionally, the models differ in how pore-fluid pressure influences the flow momentum balance. In our model, the primary influence is through the basal frictional resistance, which is reduced in the presence of excess pore-fluid pressure via Terzaghi's effective stress principle. Kowalski and McElwaine's 2013 model could potentially accommodate a similar effect, although currently the basal friction is proportional to the total bulk normal stress.

The equations derived in our companion paper [1] are also distinct from those of previous debris-flow models in other regards. The role of pore-fluid pressure in our model is unlike that in the models mentioned above, because it is one of the evolving dependent variables. In our model evolving granular dilatancy relates the evolution of excess pressure and solid volume fractions, and closure is provided by a Darcy drag law for solid–fluid interactions [1].

Because of these novelties, one of our objectives in this paper is to evaluate the mathematical structure of our model equations in a general context. The equations comprise a non-conservative hyperbolic system with source terms, but unlike the equations of most two-phase models, hyperbolicity is guaranteed without elliptic degeneracy (as is true in Kowalski and McElwaine's recent model). Additionally, the eigenstructure can be determined in a simple explicit form. The influence of dilatancy is administered through the source term, which tightly couples the set of evolution equations.

Other primary objectives of this paper are to provide a description of numerical issues unique to our equations and a comparison of numerical simulations and data from large-scale debris-flow experiments at the US Geological Survey (USGS) debris-flow flume [13–15]. These comparisons test the model's ability to simulate the complex, dynamic properties of debris flows. In particular, they test the ability to seamlessly simulate processes involving slope failure, liquefaction, downslope motion and deposition.

## 2. Summary of the mathematical model

### (a) Governing equations

Our model equations describe a debris flow on an inclined planar basal surface, with flow velocity components parallel to the surface. We use bed-normal coordinates, where *x* and *y* are local orthogonal directions tangential to the bed, and *z* is the bed-normal direction in which quantities are depth averaged (see fig. 6 in [1]). We adopt conventional shallow-flow assumptions and thereby assume that the depth of the flow is small relative to the tangential length scales of the flow (e.g. [16]). In §2*c*, we further discuss these assumptions as well as issues regarding more arbitrary basal surfaces.

Our equations govern five dependent variables, *h*(*x*,*y*,*t*), *hu*(*x*,*y*,*t*), *hv*(*x*,*y*,*t*), *hm*(*x*,*y*,*t*) and *p*_{b}(*x*,*y*,*t*), where *h* is the height normal to the bed of a virtual free surface defined in our companion paper (see fig. 7 in Iverson & George [1]), *u* and *v* are the depth-averaged flow velocities in the *x* and *y* directions, respectively, *m* is the depth-averaged solid volume fraction and *p*_{b} denotes the pore-fluid pressure at the basal boundary. The equations can be written as
*g*_{z} is the gravitational acceleration in the bed-normal direction (i.e. ** g**=(

*g*

_{x},

*g*

_{y},

*g*

_{z})

^{T}). (See appendix A for the exact correspondence between the system (2.1) and the equations presented in [1].) The depth-averaged bulk density,

*ρ*, is a function of the solid volume fraction

*ρ*

_{s}and

*ρ*

_{f}are the intrinsic material densities of the solid and fluid constituencies, respectively, which are assumed to be constant. The variable coefficient

*χ*is introduced in (2.1e) for notational convenience and is given by

*κ*in (2.1b) and (2.1c), is described in detail by Iverson & George [1], and appendix A discusses the implications of a variable versus constant

*κ*.

The source terms on the right-hand side of (2.1) are defined as
*τ*_{s}=(*τ*_{s,x},*τ*_{s,y})^{T} and *τ*_{f}=(*τ*_{f,x},*τ*_{f,y})^{T} are the basal shear traction exerted by the solid phase and fluid phase, respectively, *α* is the debris' elastic compressibility, *ψ* is the granular dilatancy angle and *D* is the depth-averaged granular dilation rate, equal to the depth integral of the divergence of the solid-phase velocity field [1]. The variable coefficient *ζ* is introduced in equation (2.4e) for notational convenience and is given by

As shown by Iverson & George [1], our formulation indicates that *D* is proportional to the difference between basal pore-fluid pressure and hydrostatic pressure
*k* is the hydraulic permeability of the debris, and *μ* is the pore-fluid viscosity. Equation (2.6) implies that the pore-fluid pressure *p*_{b} relaxes towards hydrostatic pressure, *ρ*_{f}*g*_{z}*h*, at the same rate as the relaxation *ζD* contributes to pore-fluid pressure relaxation, with a rate constant dependent on the model parameters and solution state. This relaxation arises naturally as a consequence of pore-pressure diffusion, and the timescale of this process is analysed and discussed in detail in our companion paper [1].

The hydraulic permeability of the debris, *k*, is a variable parameter. Various formulae for the dependence of *k* on properties of sediment mixtures have been proposed in the literature, but most formulae assume relatively idealized distributions of grain sizes and geometries. For poorly sorted geological debris containing a great variety of grains, we use the empirical formula
*k*_{0} can be measured directly (see fig. 5 in [1]) or estimated for a natural debris flow.

The elastic compressibility of the debris, *α*, is also a variable parameter that is constrained by data. Iverson & George [1] show that a suitable empiricism is
*a* and *σ*_{0} are specified constants (e.g. *σ*_{0}=1000 *Pa* and *a*=0.01–0.05 are typical for debris used in flume experiments).

Because the shear stress is depth-integrated, *τ*_{f} and *τ*_{s} account for the total basal resistance to motion. For the fluid component of the resistance, we use a standard estimate
*μ*, accounts for the presence of suspended fine particles (e.g. clay and silt). In most cases, the contribution of *τ*_{f} to flow resistance is minor, however. The dominant contribution comes from the solid component *τ*_{s}. For non-stationary states, we estimate *τ*_{s} by using the Coulomb–Terzaghi equation for granular friction [1]
*σ*_{e} is the basal effective normal stress, given by
*ϕ*, plus the granular dilatancy angle, *ψ*—an evolving quantity that depends on *m*. As explained in detail by Iverson & George [1], we define the dependence of *ψ* on *m* by using the simple relationship
*m*_{crit} is the critical-state solid volume fraction commonly employed in soil mechanics (i.e. the quasi-static equilibrium solid volume fraction such that there is no change in *m* during shearing, with *σ*_{e} constant). The quantity *N* is a dimensionless state parameter, which we define as
*δ* is a length scale associated with grain collisions [1]. In our simulations, *m*_{crit} ranges from 0.56 to 0.64, and we set *δ*=1 mm.

For stationary debris, *τ*_{s} represents static friction, which cannot be defined by equation (2.10). For this case, we bound *τ*_{s} by the shear traction for deforming debris
*τ*_{s} should oppose the net driving force due to gravity and longitudinal stress gradients in the momentum equations (2.1b,*c*). If we denote this net driving force per unit area as ** F**=(

*F*

_{x},

*F*

_{y})

^{T}, we let

**∥ exceeds the maximum shear traction defined by (2.15), the local force balance is no longer in equilibrium and the mass locally becomes mobile. Otherwise, the local force balance is in equilibrium and the mass locally remains immobile. Therefore, in general,**

*F*

*τ*_{s}is defined by equation (2.16) when ∥(

*u*,

*v*)

^{T}∥=0 and by equation (2.10) otherwise.

### (b) Hyperbolic structure of the equations

The system of governing equations (2.1) can be written compactly as
*q*=(*h*,*hu*,*hv*,*hm*,*p*_{b})^{T}, *φ*=(*φ*_{1},*φ*_{2},*φ*_{3},*φ*_{4},*φ*_{5})^{T} and the functions *q* in the *x*- and *y*-directions, respectively. The matrix-valued functions *f*_{x}(*q*), *f*_{y}(*q*),

Desirable stability properties and well-posedness of systems having the form of (2.19) (neglecting the source term *φ*) are established by hyperbolicity in the eigenstructure of *n*_{x},*n*_{y})^{T} and admissible solutions *q*, the matrix
*ϑ*=*χ*+(1−*χ*)*κ*. Equation (2.3) and *ρ*_{f}≤*ρ* imply that 3/4<*χ*≤1, which, with *κ*>0, implies that *ϑ* must remain positive. All of the eigenvalues (2.24) are therefore real, with
*u* has an algebraic multiplicity of three, hyperbolicity requires three independent eigenvectors associated with *λ*^{2,3,4}(*q*). The eigenvectors are
*r*^{2,3,4}(*q*) are linearly independent for all *q*, implying that *q*. Therefore, (2.1) is hyperbolic.

The wave structure of the five characteristic fields associated with the eigenpairs is revealing, although not surprising. The first and fifth characteristic fields associated with *λ*^{1,5}(*q*) and *r*^{1,5}(*q*) are genuinely nonlinear fields associated with propagation of waves at speeds *ϑ*. This modification is due to the parameter *κ*, where *ϑ*=1 if either *κ*=1 or *χ*=1 (such as with a pure fluid with *ρ*=*ρ*_{f}). The remaining fields are linearly degenerate and represent advection of *v*,*m* and *p*_{b} at the flow speed *u*. In Riemann problems discussed in §3, the nonlinear fields carry shocks or rarefactions, while the linearly degenerate fields carry moving contact discontinuities between adjacent material states. This structure is similar to that of many hyperbolic problems for which numerical techniques have been devised. For an overview of these topics, see e.g. LeVeque [17].

The hyperbolicity of system (2.1) is due, in part, to our modelling of only a single velocity field rather than independent fields for solid and fluid phases. Many two-phase depth-averaged systems are known to lose hyperbolicity due to the appearance of complex eigenvalues when phase velocity fields diverge past some threshold (e.g. [4,7,9]). This property can render the governing equations unstable, and the loss of hyperbolicity must be prevented numerically, often through techniques that limit the differences between solid and fluid velocity fields. The unconditional hyperbolicity of our model equations (2.1), and the simple explicit expressions (2.24) and (2.26) for the eigenpairs, facilitates the development of stable numerical schemes, as briefly described in §3. For many two-phase depth-averaged models (e.g. [4,7,9]), the eigenvalues must be estimated numerically. The model recently developed by Kowalski & McElwaine [12] also maintains hyperbolicity and has a simple eigenstructure.

### (c) Variable basal topography

The classical shallow-water equations [22] are derived by vertically depth-averaging three-dimensional Euler equations in a standard Cartesian coordinate system and retaining the leading-order terms implied by the thin-layer assumption. Variable basal surface topography, described by *z*=*b*(*x*,*y*), introduces source terms into the momentum equations. For shallow flows on steep slopes, errors from neglected terms in the standard derivation can grow with the slope magnitude, and depth-averaging in a coordinate system with *z* normal to the basal surface can give a better approximation. This is the approach advanced by Savage & Hutter [16,23] for granular avalanches, and frequently adopted since. Such derivations commonly neglect some terms arising from bed curvature, i.e. from variations in the slope inclination angle, *θ* [24]. More recently, depth-averaged equations in more general surface-following curvilinear coordinate systems have been derived using asymptotic analysis (e.g. [24–28]).

Regardless of the choice of coordinate system, errors related to the flow aspect ratio, the basal slope and basal curvature are introduced during the depth-averaging process [16,24,27,29]. The relative magnitude of these errors depends on the basal surface and flow properties, and on simplifying assumptions used in the model derivation, which may be violated in the model's numerical application. The errors can be placed in two categories: errors related to neglected physics and errors related to neglected metric terms. For instance, the slope inclination angle *θ*(*x*,*y*) can be allowed to vary in order to accommodate a curving basal surface, while neglecting terms needed to accurately express conservation. The resulting errors grow with the magnitude of the curvature, but may be smaller than other terms neglected in the model derivation [16,24,27]. For flow on steep and curving surfaces, an approximation used in nearly all depth-averaged models diminishes their physical fidelity: flow acceleration in the direction of depth-integration is assumed to be negligible. In general, curvilinear coordinates, some bed-normal centripetal acceleration effects can be retained by expressing the basal normal force as a curvature-dependent pseudo-force. For Coulomb materials, this basal pseudo force can modify the basal frictional resistance, which can generate the leading-order effect of curvature [16,26].

Natural debris flows occur in topography with potentially large curvatures and irregularities. However, the suitability requirements of a curvilinear coordinate system impose restrictions on the surface in terms of smoothness and curvature relative to the flow thickness. For a coordinate system following a surface *z*=*b*(*x*,*y*), where *z* is measured normal to *S* from *b* at the curvilinear coordinates (*x*,*y*) (e.g. [26]). This representation introduces source terms in the momentum equations similar to those in the standard shallow-water equations. The momentum source terms in (2.4) then become
*et al.* [26] by including a correction to basal Coulomb resistance due to an apparent centrifugal force by replacing the *x*-component of equation (2.10) with
*x*) direction. Similarly, as in [16,26], we neglect other geometric terms arising from curvature, and simply allow *θ* to vary as a function of only *x*. This variation implies a spatial dependence of *g*_{z}(*x*), which does not alter the conclusions drawn in §2*b*. Deriving improved depth-averaged models that more accurately accommodate the physical effects of complex topography is an active research pursuit (e.g. [24–28]). Because our model is tested in §4 with the simple basal surface of the concrete flume, we postpone this issue for future studies.

## 3. Numerical methods and D-Claw software

General hyperbolic systems such as (2.19) present well-known numerical challenges and require specialized algorithms in order to obtain stable and convergent solutions (e.g. [17,30–33]). Most notably, the equations allow the development of non-unique discontinuous solutions requiring suitable numerical schemes that can converge stably to the correct solution. A variety of such methods have been developed over the past several decades (e.g. LeVeque [17] for an overview), many based on shock-capturing Godunov methods [30]. Our algorithms belong to this class and are implemented in a software package (named D-Claw), which contains features tailored to modelling depth-averaged flows over topography. D-Claw is an extension of the Clawpack software package [34], an open-source project for solving general hyperbolic problems (e.g. [17,35–37]).

In the following sections, we provide a brief overview of our computational algorithms, with emphasis on features particular to our debris-flow equations (2.1). Background material and a detailed description of the algorithms are provided in the references cited.

### (a) Shock-capturing finite volume methods

The numerical schemes in D-Claw are finite volume methods implemented on logically rectangular grids. We denote an arbitrary grid cell *x*_{i−1/2},*x*_{i+1/2})×(*y*_{j−1/2},*y*_{j+1/2}), with area

Our numerical integration schemes are based on wave-propagation algorithms [17,38] developed for systems of form (2.19). These are high-resolution (i.e. formally second-order accurate and shock-capturing) Godunov-type methods, meaning that they use solutions to Riemann problems computed at grid-cell interfaces. For instance, between grid cells *t*^{n} to *t*^{n+1}. (For convenience, in this section we focus on the Riemann problem between *y*-direction, modification of the formulae that follow is straightforward.)

Nonlinear Riemann problems are commonly approximated by using, for example, a local linearization of the matrix

For conservative systems with *et al.* [39] for details). Numerical conservation of quantities transported by fluxes is required for general convergence to discontinuous solutions [17,40,41].

If the system has a source term, *φ*(*q*,*x*,*y*), (3.5) can be modified to include its effect
*Φ*_{i−1/2,j} is a consistent approximation to *φ*(*q*,*x*,*y*) at (*x*_{i−1/2},*y*_{j}). The advantage of (3.6) over treating the source term via operator splitting is due to the numerical preservation of steady states that arise due to a balance of the spatial gradients on the left-hand side of (2.19) and the source term on the right-hand side. Numerically maintaining these equilibria poses a subtle but significant challenge, as the countervailling terms might be large in magnitude compared with their difference (and hence ∂*q*/∂*t*). This well-known difficulty has given rise to ‘well-balanced schemes’, which can accurately resolve balanced steady states, or small perturbations to such (e.g. [39,42–46,47]). Prevalent equilibria of interest depend on the particular application. Examples include a motionless pool of water lying over variable topography, a constant river-like flow, or, in the case of debris flows, a stationary mass of debris poised on a slope. In the former cases, the balance is dominated by pressure gradients and the momentum sources (2.27) due to variable bottom slope. In the latter case of a static debris mass, the dominant balance includes normal-force gradients, the large gravitational driving forces and Coulomb friction. For debris flows, initiation from a static state occurs as a small perturbation to these large counterbalanced forces, leading to a dynamic feedback that influences the ultimate flow mobility.

We distinguish terms in *φ* (equation (2.4)) by the following splitting:
*φ*^{B}=(0,0,0,0,*φ*_{5})^{T} is the pressure source term, and *φ*^{A} contains all other terms in (2.4) (and in (2.27) in the case of variable topography *b*(*x*,*y*)). The dominant terms arising in *φ*^{A} represent gravitational driving forces and frictional resistance—potentially large terms for the balanced steady states described previously. The term *φ*^{B} represents relaxation of basal pore-pressure towards a steady state, a term that is large when the solution is far from a steady state.

Our development of the Riemann solver for system (2.1) was guided by several criteria. First, it is desirable to treat the transport of quantities due to fluxes in a numerically conservative fashion (without introducing spurious mass, momentum, etc.) as well as include terms from the non-conservative products directly in the Riemann solution. Additionally, the solver must possess well-balanced properties with regard to important balanced equilibria, particularly static states. Finally, the solver must generalize to Riemann problems at the flow front, between vacant and inundated cells, without the appearance of spurious results such as negative depths or unphysical oscillations. (There is extensive literature on solving similar ‘wet-dry’ Riemann problems for the shallow-water equations (e.g. [37,47–50]).) Our solver is based on a decomposition of the form
*r*^{p}(*q*), *p*=1,…,5, (equation (2.26)) based on the solution states *f*_{x}(*q*,*x*,*y*) and *g*_{z}, if *θ* is allowed to vary.) The decomposition (3.8) uniquely determines the coefficients *p*=1,…,5. The Riemann solution can then be constructed from this decomposition and from approximations for the eigenvalues,

The solver is simplified somewhat if *κ*=1, which was used for the computations described in §4; here we consider that case exclusively. Note from (2.26) with *κ*=1, that all of the eigenvectors *r*^{1,5}_{i−1/2,j}, denoted by (*r*^{1,5}_{i−1/2,j})_{4,5}. For these final two components, we use *v*, *m*, *ρ* and *χ*, appearing in *r*^{1}(*q*) and *r*^{5}(*q*), are constant through the waves in the nonlinear fields, and vary only across the contact discontinuities carried by *r*^{2,3,4}(*q*). The derivation of *p*=1,…,5, is beyond the scope of this paper but endows the Riemann solver with properties necessary for wet-dry problems, shock capturing and physical admissibility (energy) requirements. George [47] provides further details.

The term *φ*(*q*)^{A} appearing in (3.7) and described previously. The treatment of variable topography, the approximation of (2.27) in the definition of *θ*. Other terms approximated by *Q*^{n}_{i−1/2,j} and *p*_{b}, at each cell centre the system,

The numerical update is constructed from the Riemann solutions using the wave-propagation algorithms of LeVeque [51]. Details are beyond the scope of this paper, but are thoroughly described for general hyperbolic problems by LeVeque [17,51]. Multiple splitting options are available in D-Claw to account for wave propagation and (3.9). For the simulations described in §4, we simply used Godunov splitting (e.g. [17]).

### (b) Adaptive mesh refinement

Our numerical methods use adaptive mesh refinement (AMR), which allows efficient resolution of flow dynamics by providing a discretization that evolves based on the numerical solution. Our algorithms are an extension of patch-based, block-structured AMR originally developed for general hyperbolic problems (e.g. [52,53]), and later specialized for depth-averaged flows across topography, such as tsunamis and overland floods (e.g. [35–37]). With patch-based, block-structured AMR, the domain is discretized on numerous rectangular grids, each associated with a given level, *l*=1,…,*L*, of spatial and temporal resolution. Preceding any time step *l* grid can be dynamically refined into multiple level-(*l*+1) cells, by predefined cell-length ratios, *l*+1) grids. Figure 2 illustrates a single level-*l* grid cell, *l*+1) cells, with

Our simulations of debris-flow flume experiments presented in §4 illustrate the use AMR. However, to appreciate the efficiency enhancements conferred by AMR, one must consider a large-scale natural debris flow in mountainous terrain. Such a debris flow typically begins in a relatively localized region where mass failure of a slope or steep stream bed occurs. The flow then follows a complex route that may be difficult to anticipate in advance. At any given time, the debris flow typically occupies only a small portion of even the smallest rectangular region that bounds the ultimate path of the debris. During a computation, debris is absent in the vast majority of the domain, implying that grid cells there should be as coarse as possible. Simply using static discretizations with finer grid cells in regions where the debris is likely to go, based on the terrain, confers some computational efficiency. However, ideally, grid cell sizes should evolve with the flow because the optimal grid resolution is highly spatially and temporally dependent. For debris flows, our AMR algorithms generate fine grids in regions just prior to the arrival of the flow front. It is also possible to use even finer grids surrounding complex features where spatial gradients in the solution are highest. The use of AMR for floods in natural and highly variable terrain, using the shallow-water equations, is illustrated by George [37].

## 4. Numerical results

To test D-Claw, we present simulations of two types of large-scale experiments conducted at the USGS debris-flow flume (figures 3 and 4). Detailed descriptions and video documentation of these experiments are available elsewhere [14,15,54]. In contrast to natural debris flows, flume debris flows involve tightly controlled initial and boundary conditions as well as material properties, and thereby afford opportunities for stringent model tests. Moreover, experimental flows in the USGS flume are large enough (typically 6–10 m^{3}) to obviate most scaling problems [15].

In the first type of test, we address short-timescale dynamics during onset of debris-flow motion triggered by gradually rising pore-fluid pressures in static slopes. In the second type, we address dynamics over the longer timescales that are relevant during subsequent downslope motion, runout and deposition. Traditionally, these two distinct problems have been treated with dissimilar models, one for onset of slope instability and one for the dynamics of ensuing motion. By contrast, D-Claw simulates both flow onset from an initially balanced equilibrium state and flow dynamics in states that are far from equilibrium.

### (a) Predictions of slope failure and flow mobilization

In this section, we investigate D-Claw's ability to predict debris-flow mobilization from an initially stable sediment mass. In each of the mobilization experiments of Iverson *et al.* [14], a 6 m^{3} rectangular prism of sediment composed mostly of sand, with 11% finer grains by dry weight, was placed behind a retaining wall positioned normal to the bed near the top of the flume (figures 4 and 5). A ramp downslope of the retaining wall provided a continuous basal running surface for flows that overtopped it. Rising pore-water pressures that triggered slope failure were generated by various combinations of surface sprinkling, infiltration from an upslope pond and groundwater injection (cf. [55]). As pore pressures rose, water tables were maintained nearly parallel to the impermeable flume bed by using a valve to adjust the water discharge from a horizontal drain pipe buried in the sediment adjacent to the base of the retaining wall. Evolving pore-water pressures and slope deformation were monitored continuously by a network of piezometers, tensiometers, subsurface tiltmeters and surface extensometers [14].

The chief variable manipulated in these experiments was the initial bulk density of the sediment prisms. In some experiments, the sediment was loosely placed by dumping it from buckets with minimal disturbance, resulting in initial solid volume fractions *m*_{0}<0.5. In other cases, the sediment was densified by placing it in a series 10 cm thick layers and compacting each layer with high-frequency mechanical vibrations, resulting in initial solid volume fractions *m*_{0}≈0.6. Ancillary testing revealed that the critical-state solid volume fraction for this sediment at appropriately low confining stresses was about *m*_{crit}=0.56 [14]. Therefore, initial states of sediment in experiments with *m*_{0}<*m*_{crit} were classified as ‘loose’, while those in experiments with *m*_{0}>*m*_{crit} were classified as ‘dense’.

During slope failure, loosely packed sediment quickly liquefied and transformed into rapid debris flows (see video documentation of [54] for 6 June 1999). This behaviour resulted from sediment contraction upon shearing, evidenced by a dramatic, contemporaneous increase in pore-fluid pressure. With densely packed sediment, on the other hand, slope failures exhibited self-stabilizing behaviour, because pore-fluid pressures decreased as the deforming sediment dilated. In such cases, landslides crept slowly and sometimes intermittently downslope [14,54]. These contrasting styles of behaviour test not only D-Claw's ability to seamlessly simulate the transition of a stable mass into a moving mass, but also the model's ability to regulate evolution of pore-fluid pressure and thereby influence the degree and nature of landslide mobilization.

#### (i) Mobilization experiment 1: loosely packed sediment (*m*_{0}<*m*_{crit})

We first consider model performance in simulating an experiment with loosely packed sediment. To simulate the rising water table prior to slope failure, we gradually and uniformly increased basal pore-water pressures so that they represented an increasing percentage of the total lithostatic pressure acting on the basal slip surface. Tiltmeter data showed that most of the initial failure surface ran along the concrete flume bed, but near the retaining wall the failure surface curved upward to intersect the wall's crest. We used this basal slip surface to define the initial thickness distribution, *h*(*x*,*y*,0), in D-Claw simulations (figure 5).

Measured soil properties and corresponding values of parameters used to simulate the loose-soil experiment are listed in the upper half of table 1. The only complication in using these parameter values involves interpretation of the effective friction angle. Based on triaxial unloading tests of loose sediment specimens, Iverson *et al.* [14] inferred values of 29°±2° for the effective sediment friction angle during failure. This value differs negligibly from the basal friction angle measured as sediment begins to slide along a smoothly finished concrete surface like that underlying the soil prisms in the flume [14]. However, some of our computations used *ϕ*=38° as the basal friction angle in order to reproduce the pore-fluid pressure head ≈0.26 m measured when slope failure occurred. The high effective friction might be due to sediment settling observed during the pre-failure wetting process, with resultant increases in sediment bulk density and frictional strength. Another contributing factor might be non-negligible sidewall friction. In our simulations, the flume sidewalls are modelled as topographic features without other distinctive properties, and because the walls are vertical, no basal friction is engaged there. Fortunately, the use of *ϕ*=29° versus *ϕ*=38° in our simulations has little effect on short-timescale dynamics immediately following slope failure, because the effects of intrinsic friction are overwhelmed by those of evolving pore pressure.

Computational results are illustrated in figures 6 and 7. Upon failure, the value *ψ*=−4.57° applies in order to satisfy equation (2.12) with the measured values *m*_{0}=0.48 and *m*_{crit}=0.56. Physically, however, *ψ* has no relevance until an infinitesimal displacement occurs. As failure begins, we impose the negative dilatancy at locations where there is non-zero motion, and the resultant slight weakening of the sediment qualitatively mimics behaviour observed in laboratory testing of loose, contractive soils (e.g. [56,57]). More significantly in our computations, however, negative dilatancy can lead to pore-pressure feedback that results in profound weakening associated with liquefaction.

Slope failure in the model is detected when non-zero momentum occurs at any location, after which changes in pore-fluid pressure are no longer imposed, but instead are allowed to evolve according to the model equations. Upon failure, sediment in the region −2 m<*x*<0 begins to deform, accompanied by a co-localized rise in pore pressure. The region of deforming sediment then expands upslope, accompanied by rapidly increasing pore pressures that reach lithostatic levels (*p*_{b}=*ρg*_{z}*h*) within a few seconds (see figure 6 and electronic supplementary material).

Displacement and pore-pressure time-series data collected during the experiment are compared to D-Claw output in figure 7. Because instruments moved downslope with the sediment, the D-Claw results were produced by following Lagrangian material points positioned at the same initial locations as the instruments. Determining the position of, and pressures at, the moving material points involved dynamically integrating and interpolating the Eulerian numerical solution. It is impossible to ascertain the exact time of failure in the experiment, owing to some preliminary settling of the soil and the positioning of isolated detection instruments. In figure 7, *t*=0 *s* corresponds to *t*=2781.1 *s* in the experimental dataset presented by Iverson *et al.* [14]. The temporal axes for the model results are chosen to align with the experimental data, based on when the basal pore pressure sensor moved over the downslope retaining wall. Displacements of the extensometers are plotted on these same temporal axes.

The evolution of the pore-fluid pressure at Lagrangian points was found to vary based on the initial value of the elastic compressibility, *α*_{0}. Because this parameter is difficult to constrain experimentally, we used two values for *a* in equation (2.8), which is used to calculate *α*. Although most data from quasi-static soil testing suggest that 0.01≤*a*≤0.05 [1], we found that a larger value reproduced the measured slope-failure dynamics more closely (figure 7). Nevertheless, the initially higher compressibility, *α*_{0}=7.2×10^{−5} *Pa*^{−1} with *a*=0.3, lies well within the range of plausible values for shearing debris [1], and in simulations the evolving value of *α* maintains the same order of magnitude.

Figure 7*a*,*c* compare computed pore-fluid pressures with data collected using a piezometer initially positioned at the base of the sediment (65 cm deep) 1.9 m upslope from the retaining wall, and with data from two piezometers directly above it at vertical depths of 30 and 50 cm [14]. The simulated pressure for these depths is determined by scaling the computed basal pressure according to the vertical pressure profile derived by Iverson & George [1]. The right panels show the position of two ground surface extensometers, initially located 3.07 and 0.71 m upslope from the retaining wall. The upper and lower panels correspond to the different values of *a* and *ϕ* employed in alternative simulations (*a*=0.03, *ϕ*=29° and *a*=0.3, *ϕ*=38°). In both cases, the key features of interest are the elevated values of pore-fluid pressures that rapidly developed as the mobilized sediment contracted and then overtopped the retaining wall. The larger pore-fluid pressures calculated with *a*=0.03 and *ϕ*=29° (lower panels) resulted from persistence of a greater thickness of sediment, and hence higher basal pressure head, as sediment moved towards the retaining wall. The case with *a*=0.3 and *ϕ*=38° resulted in greater variability in velocities at different locations and a corresponding longitudinal extension and thinning of the moving sediment slab. This resulted in reduced depths of sediment flowing over the retaining wall. In both simulations, the debris became fully liquefied, and pore pressures dropped abruptly after sediment moved over the retaining wall.

#### (ii) Mobilization experiment 2: densely packed sediment (*m*_{0}>*m*_{crit})

Next, we consider a dense soil experiment described by Iverson *et al.* [14]. Our main motivation for including this experiment is to explore the degree to which D-Claw can simulate fundamentally disparate dynamics due to different initial states of the sediment alone. In the dense case, the initial porosity was significantly reduced due to systematic compaction during sediment placement. This compaction in turn affected some of the soil properties, such as the permeability and effective basal friction angle. Otherwise, the experimental procedure was very similar to that of the loose soil experiment described above. The experimentally measured and model parameters are shown in the lower panel of table 1. In contrast to the loose soil experiment, the initial solid-volume fraction, *m*_{0}, is higher than the critical solid-volume fraction, *m*_{crit}, indicating positive dilatancy.

In the experiment, deformation first occurred when the basal pore-fluid pressure head exceeded 45 cm (considerably higher than for the loose soil experiment). The resulting landslide exhibited dynamics that differed dramatically from those of the loose soil case [14]. The sediment underwent slow episodic motion over a protracted period of time (15 min). Brief episodes of incremental deformation produced displacements of less than 0.3 m, with a total displacement of less than 1.5 m. This behaviour was due to the dilative response of the denser soil, which reduced the pore-fluid pressure during deformation, acting to retard motion [14]. The pore-fluid pressure slowly rebounded after each precipitous drop, leading to semi-periodic stick–slip cycles.

As in the loose soil case, we simulated the dense soil experiment by gradually raising the pore pressure to the point of slope failure. The basal pressure head triggering failure at *t*=0 was approximately 43 cm, slightly lower than the pressure observed in the experiment. For comparison, we also performed a simulation with an initial basal pressure head of 45 cm, the pressure observed experimentally when failure commenced. The results were similar in each case. A cross-sectional perspective of the simulation is shown in figure 8. Displacements observed in the first 60 s of simulated time are nearly undetectable in figure 8, but are responsible for slight downslope movement.

The D-Claw simulation exhibited the slow, episodic creeping motion similar to that observed experimentally. Figure 9 shows a time series of the pore pressure and displacement of a Lagrangian material element. Note that the sediment moves only centimetres in 120 s—a very slow rate compared with the loose soil case. The time-averaged velocity (the total displacement per duration) of ≈1 mm s^{−1} was similar in the experiment and the simulation [14]. Small stick–slip cycles that are evident in figure 9 have periods more brief than those observed in experiments, probably because motion was continually driven in the experiments by continually adding pore water. In the simulation, such water addition is not considered, and the evolving pore pressure slowly declines. More-detailed analysis of stick–slip behaviour in the absence of water addition is provided by Schaeffer & Iverson [58].

### (b) Predictions of downslope debris-flow dynamics

In this section, we investigate D-Claw's ability to predict debris-flow dynamics, runout and deposition measured in a set of flume experiments reported by Iverson *et al.* [15]. In these experiments, loose sediment was placed in a hopper behind a vertical gate near the top of the flume, saturated with water, and released by suddenly opening the gate (figures 3, 10 and 11). Upon releasing the headgate, the saturated mixtures quickly liquefied and developed into debris flows that descended the flume channel before emerging onto a nearly planar runout pad and forming deposits.

We focus on the mean behaviour measured in an aggregated set of eight experiments referred to as sand–gravel–mud (SGM) experiments due to the composition of the sediment, which contained about 6% mud-sized particles by dry weight [15]. Once motion commenced, persistent suspension of these particles enhanced the effective pore-fluid viscosity. This enhanced viscosity has been shown to increase debris-flow mobility by facilitating maintenance of high pore pressures [10,15].

The simulation parameters and corresponding values of sediment properties measured in conjunction with the SGM experiments are shown in table 2. Values of all parameters used in the D-Claw simulation matched experimentally determined values, except for that of the initial hydraulic permeability *k*_{0}, which was chosen to be about two orders of magnitude larger in the simulation (although still within the typical range for loose soils [59]). Our rationale for using a higher value was partly that *k*_{0} was measured in static debris, whereas agitation of rapidly moving debris transiently opens pores and thereby facilitates fluid flow with respect to adjacent solid grains. We discuss our choice of *k*_{0} further below.

Oblique views of the D-Claw simulation during the initial stages of a gate release experiment are shown in figure 12 and in the accompanying animation. The opening of the gate was simulated by using an evolving function, *b*(*x*,*y*,*t*), for the basal elevation in the model, representing the height of the laterally pivoting doors. Because of our model coordinate system, the doors were assumed to be normal to the bed rather than vertical. Videotape footage from several experiments was analysed in order to approximate the speed of gate opening, which occurs in about 0.85 s. Based on comparisons with dam-break initial conditions (i.e. an instantaneously opened gate), we found that including the effect of the finite gate opening speed had an important influence on the first stages of flow evolution (figure 13).

Aggregated time-series data reported by Iverson *et al.* [15] provide a basis for straightforward comparison with our D-Claw simulation. In figures 13 and 14, grey shading represents the mean ±1 s.d. measured in eight experiments. Figure 13*a* compares the measured flow depth at *x*=2 m with model results obtained with and without the simulated motion of the headgate. Note that the gate slightly impedes the flow, resulting in a delayed and more abrupt rise in thickness as most of the flow reaches *x*=2 m. Figure 13*b* shows the computed basal pore-fluid pressure *p*_{b}, compared to the lithostatic and hydrostatic pressure for the computed thickness *h*, illustrating the high degree of liquefaction at early times. Measurements of the basal pore-fluid pressure in the experiments were not obtained at *x*=2 m.

Figure 14 compares measured and computed time series that depict flow thicknesses and pore-fluid pressures at two downslope locations, *x*=32 and 66 m. Amplitudes compare well in most of the time series, but the speeds of computed flows begin to outpace those of experimental flows by the time they reach *x*=66 m. The resultant error in flow-front arrival-time is only about 1 *s*, however. Another timing discrepancy exists with respect to pore-fluid pressure. In the experimental data, particularly for *x*=66 m, the arrival of the flow front precedes the arrival of positive pore-fluid pressures by roughly 1 s. This delay in pore-pressure response is absent in the computational results but is commonly observed in debris flows. It results from particle-size segregation that creates a highly permeable, dilated, coarse-grained front that cannot sustain high pore-fluid pressure [10,15,60]. High-friction fronts thereby act to impede motion of the more liquefied flow body that trails them, resulting in reduced flow-front propagation speeds. This behaviour partly motivated our use of an enhanced permeability in the model, which globally increased the relaxation rate of fluid pressures. Indeed, in the computational results pore pressures relaxed to nearly hydrostatic levels after 4 s, while they remained at nearly lithostatic levels in regions behind the high-friction fronts of the experimental flows. In essence, the high model permeability partly compensated for D-Claw's inability to simulate high-flow resistance in pressure-depleted flow fronts, but it caused excessive pore-pressure relaxation in flow tails. Remedying this shortcoming by extending the model to simulate particle size segregation is therefore a priority.

The final state of debris-flow motion involves deceleration and formation of static deposits. Figure 15 illustrates our simulation of this process on the run-out pad at the foot of the flume. The flow front reaches its ultimate downslope distance from the headgate (*x*≈112 m) at *t*≈12 s. Subsequently, additional debris continues to advance onto the run-out pad, leading to a gradual accumulation of material and a maximum deposit thickness of *h*≈18 *cm* at *x*≈106 m. This behaviour mimics that observed in experiments [15,54], but direct quantitative comparisons are complicated by the variability of the experimental results, illustrated in figure 16. A consistent pattern is that the experimental deposits underwent less lateral spreading than in the simulation, however. This discrepancy is due to the presence of lateral levees of coarse-grained high-friction material that helped confine and channel runout of the experimental flows. This phenomenon, which results from grain-size segregation, has been observed and studied extensively (e.g. [60]) but is not yet included in D-Claw. Despite the lack of a segregation mechanism, the model does possess some ability to produce levee-like effects owing to higher basal shear resistance in shallow lateral flanks of the decelerating debris. In the model, this behaviour results from enhanced relaxation of pore pressure where low values of velocity and depth exist. As can be seen in figure 15 after ≈12 s, wing-shaped levees form as deposited lateral material remains immobile while most of the mass continues to move forward near the centre of the flow (bottom panels of figure 15).

### (c) Gauging model behaviour with dimensionless parameters

Representation of the dynamics of flowing debris by D-Claw is fundamentally affected by three dimensionless combinations of parameters, as well as by an initial condition characterized by *m*_{0}−*m*_{crit} (i.e. the initial dilatancy) [1]. Values of these parameters for the three types of simulations described above are shown in table 3. From top to bottom, the quantities represent the tangent of the initial dilatancy angle, a ratio of timescales for downslope motion and pore-pressure diffusion, a normalized debris compressibility and a reciprocal Reynolds number. Because the initial stages of motion were of primary concern for the mobilization experiments, values of *H*=0.5 m and *L*=5 m were used in table 3 for depth and length scales based on the initial debris profile for those experiments. For the gate release experiments, downslope motion was of primary concern and hence *H*=0.1 m and *L*=100 m were used, roughly corresponding to the dimensions of a fully developed debris flow in the flume.

Among all of the values of the dimensionless parameters listed in table 3, the greatest contrast exists in the values of the motion-to-diffusion time-scale ratios. For the gate release experiments, a relatively large value of permeability as well as small value of *H*/*L* leads to a much larger motion-to-diffusion time-scale ratio than in the mobilization experiments. A consequence of this large time-scale ratio was relatively rapid depletion of the high pore-fluid pressures that were generated by debris collapse and contraction during opening of the flume headgate. This rapid depletion increased basal friction and regulated flow dynamics in the manner described in the preceding section. As noted there, widespread pore-pressure depletion only partially mimics the effects of grain-size segregation, which are currently lacking in D-Claw.

The values of the final three dimensionless parameters in table 3 are of similar magnitude for the mobilization experiments with dense and loose initial states. The small difference in values is primarily due to the larger baseline permeability for the loose soil case (*k*_{0}=2.6×10^{−11} versus 2.3×10^{−12}), which is a result of higher initial porosity. Therefore, the dramatically disparate results for the loose and dense case can be attributed to the initial dilatancy angles, which are opposite in sign as indicated by the values of *m*_{0}−*m*_{crit} in table 3.

## 5. Concluding remarks

In this paper, we show that the depth-averaged debris-flow model derived in our companion paper [1] constitutes a non-conservative hyperbolic system of partial differential equations with source terms. Additionally, the eigenstructure of the system has a familiar form, with two nonlinear characteristic fields similar to those of the shallow water equations and three contact discontinuities with propagation speeds equal to the flow velocity. The eigenstructure is simpler than that of many two-phase models, such as that derived by Pitman & Le [4], and it avoids elliptic degeneracy.

Notwithstanding the simplicity of its eigenstructure, as a hyperbolic system the model presents well-known challenges for numerical solution schemes. Additional numerical challenges arise owing to the presence of moving flow boundaries and variable topography. Numerous specialized techniques have been developed previously to overcome these challenges for the shallow-water equations, and we have adapted such techniques to our debris-flow model in the form of the software package D-Claw.

D-Claw employs the concept of granular dilatancy, which encapsulates the tendency of the solid volume fraction to evolve towards a value that is equilibrated to the ambient flow rate and effective-stress state. This evolution enables D-Claw to unify models for slope failure and flow initiation with those for debris-flow dynamics. Inclusion of dilatancy also results in novel evolution equations for the solid volume fraction and pore-fluid pressure, which are tightly coupled through source terms. Because of the inherent sensitivity associated with the transition from a stable debris mass to a dynamic debris flow, it is necessary for the numerical schemes in D-Claw to reliably and stably resolve balanced steady states. For debris flows, the prevalent steady state arises from the counterbalance of longitudinal normal-force gradients, gravitational driving forces and basal shear resistance, the latter of which is strongly dependent on the evolving fluid pressure. Flow initiation occurs when this balance is slightly perturbed. In order to resolve this perturbation and the resulting dynamics accurately, D-Claw requires a specialized Riemann solver and a carefully designed splitting of the source term.

Comparisons of D-Claw model output with results from two types of large-scale experiments indicate that the model can successfully simulate key features of debris-flow mobilization as well as debris-flow dynamics, runout and deposition. Model results show that in many cases the effects of persistent high pore pressure strongly suppress the effects of basal Coulomb friction, which is the primary source of flow resistance. This finding is consistent with experimental and field observations indicating that reduction of friction by high pore-fluid pressure is largely responsible for debris flows' remarkable mobility.

Model comparisons with results of natural-release experiments focus on triggering of debris-flow motion by externally imposed, gradually rising pore-water pressures. In the case of loosely packed soil, D-Claw simulations mimic debris-flow behaviour in mobilization experiments in which negative dilatancy and positive pore-pressure feedback allow a local instability to grow rapidly, thereby causing wholesale slope liquefaction and runaway motion within a few seconds. In the case of densely packed soil, by contrast, simulations show that negative pore-pressure feedback due to positive dilatancy stabilizes motion and gives rise to slow, intermittent downslope creep. Model results imply that debris-flow mobilization in such scenarios is, at most, a relatively slow, piecemeal process.

Model comparisons with results of gate release experiments involve initial conditions that closely parallel those in classical dam-break problems. However, D-Claw output shows that model results best match experimental data if the finite opening speed of the gate is included in the simulation. Flow dynamics during the first approximately 1 s of motion, when debris settlement and liquefaction are triggered by opening of the gate, set the stage for ensuing flow behaviour downslope. Model predictions of evolving flow depths and basal pore-fluid pressures closely match those measured 32 m downslope from the gate, but deviations between predictions and measurements increase as flows proceed further downslope. We believe the primary source of these deviations is the model's lack of grain-size segregation effects. In experimental flows, this segregation produces coarse-grained flow fronts that have high frictional resistance owing to high fluid permeability and rapid relaxation of pore-fluid pressure. Our simulations lack the capacity to generate this highly focused flow-front resistance, and further lack the capacity to reproduce its effect on lateral levee development during flow runout and deposition. We anticipate that the mathematical structure of D-claw will enable us to incorporate size-segregation effects in the future, however.

All of our model predictions and comparisons provide evidence of a sensitive dependence on initial conditions. In particular, the initial solid volume fraction and its implications for dilatancy and liquefaction during the onset of motion during natural flow releases are crucial. Indeed, a fundamental bifurcation of behaviour (runway acceleration versus slow, regulated motion) can result from rather minor differences in the initial state. Dependence of gate-release simulations on the details of the gate-opening rate provides further evidence of the importance of initial conditions. We infer that future modelling directed at practical hazard assessment should account for uncertainties in initial states.

Along with their dependence on initial conditions, our D-Claw simulations depend strongly on the values of just a few dimensionless parameters. Paramount among these is a time-scale ratio that involves the characteristic times for downslope motion and pore-pressure diffusion. Within these timescales, two physical properties play crucial roles. One is the Darcian hydraulic permeability of the debris and the other is its compressibility (or reciprocal bulk modulus). Both properties are likely to vary as a consequence of dilation and agitation of moving debris. Virtually, no data or compelling theories exist to firmly constrain these variations, and D-Claw currently treats the variations using simple postulates.

## Appendix A. Correspondence of model equations

In this appendix, we demonstrate the correspondence of the governing equations, comprising system (2.1) and source terms (2.4), with the equations derived in our companion paper [1]. In [1], the governing equations were presented in their most physically revealing form. In this paper, for analysis and description of numerical schemes, the hyperbolic system has been put in more standard mathematical form. Additionally, form (2.1) provides the relevant fluxes for physically based Rankine–Hugoniot conditions. Note that in this paper the overbars used to denote depth-averaged quantities in [1] have been dropped.

**(a) Momentum** **and**

Beginning with equation (4.11) for *ρ* is neglected in the final expression. Inserting the final expression of (A 2) into (A 1), dividing by *ρ* and grouping spatial derivatives on the left-hand side leads to
*κ* is constant (i.e. if *ϕ*_{int}=*ϕ* in eqn (4.32) in [1]) then it can be absorbed into the second spatial derivative and become part of the flux *f*_{x}(*q*), similar to the momentum flux in standard shallow-water models. The correspondence of eqn (4.12) in [1] for

**(b) Basal pore pressure, p_{b}**

Eqn (4.25) in [1] for the basal pressure, *p*_{b}, is
*d*/*dt* indicates a Lagrangian material time derivative. By noting that the second term on the left-hand side of (A 4) is proportional to *D* (equation (2.6)) and by using equation (2.12) for *D* and defining *χ* and *ζ* by equations (2.3) and (2.5) reduces (A 5)
*p*_{b} and using the product differentiation rule for the terms in brackets results in

Finally, eqn (4.6) for *h* in [1] is identical to equations (2.1a) and (2.4a) here, and eqn (4.7) for *hm* in [1] is identical to equations (2.1d) and (2.4d) for *hm* here.

## Appendix B. Fluxes and non-conservative products

System (2.1) is an example of the more compact form (2.17)—a general hyperbolic balance law with distinct flux-gradients and non-conservative products. For system (2.1)
*f*_{x}(*q*)/∂*q* and

If the term *κ* is constant, then the second component of *f*_{x}(*q*) can become *hu*^{2}+1/2*κg*_{z}*h*^{2} by including *κ* inside the derivative in (2.1b). In this case, the sum *κg*_{z}*h* appearing in *f*_{x}(*q*)/∂*q*. Although mathematically equivalent, including the extra term in the flux facilitates the identification of proper weak solutions and the application of conservative numerical treatment and shock-capturing. If instead *κ* is piecewise constant, it contributes discontinuous coefficients to the non-conservative product. This additional complexity may require specialized numerical treatment in order to converge to appropriately defined weak solutions (e.g. [61]).

Owing to symmetry in the momentum equations, the forms of ∂*f*_{y}(*q*)/∂*q* and *u* and *v* and permuting the second and third rows and columns in (B 2) and (B 3).

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

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