Entrainment of bed material by Earth‐surface mass flows: Review and reformulation of depth‐integrated theory

Earth‐surface mass flows such as debris flows, rock avalanches, and dam‐break floods can grow greatly in size and destructive potential by entraining bed material they encounter. Increasing use of depth‐integrated mass and momentum conservation equations to model these erosive flows motivates a review of the underlying theory. Our review indicates that many existing models apply depth‐integrated conservation principles incorrectly, leading to spurious inferences about the role of mass and momentum exchanges at flow‐bed boundaries. Model discrepancies can be rectified by analyzing conservation of mass and momentum in a two‐layer system consisting of a moving upper layer and static lower layer. Our analysis shows that erosion or deposition rates at the interface between layers must, in general, satisfy three jump conditions. These conditions impose constraints on valid erosion formulas, and they help determine the correct forms of depth‐integrated conservation equations. Two of the three jump conditions are closely analogous to Rankine‐Hugoniot conditions that describe the behavior of shocks in compressible gasses, and the third jump condition describes shear traction discontinuities that necessarily exist across eroding boundaries. Grain‐fluid mixtures commonly behave as compressible materials as they undergo entrainment, because changes in bulk density occur as the mixtures mobilize and merge with an overriding flow. If no bulk density change occurs, then only the shear traction jump condition applies. Even for this special case, however, accurate formulation of depth‐integrated momentum equations requires a clear distinction between boundary shear tractions that exist in the presence or absence of bed erosion.


Introduction
Dense, gravity-driven flows of solid-fluid mixtures that entrain bed material while descending steep slopes and channels have broad importance in geomorphology, volcanology, hydrology, and civil engineering. Examples include dam-break floods, debris flows, rock avalanches, lahars, ground-hugging pyroclastic flows, and wet snow avalanches (Figure 1) [e.g., Plafker and Ericksen, 1978;Sparks et al., 1997;Scott et al., 2005;Sovilla et al., 2006;Rogers et al., 2010;McCoy et al., 2012]. Although these subaerial mass flows may be regarded as gravity currents [Simpson, 1987], their bulk densities typically exceed that of the surrounding air by a factor of 500 or more, fundamentally distinguishing them from more dilute gravity currents in which the surrounding fluid exerts a significant buoyancy force [cf. Bonnecaze et al., 1993;García and Parker, 1993;Ancey, 2004]. Many subaerial mass flows occur abruptly and attain peak speeds that range from 10 to 100 m/s, and the largest have volumes >1 km 3 [e.g., Voight, 1978;Vallance and Scott, 1997;Jakob and Hungr, 2005;Griswold and Iverson, 2008]. In addition to being powerful agents of landscape change, they can pose grave hazards to people and property. The severity of these hazards was evidenced by debris flow disasters in Armero, Colombia, in 1985 andVargas state, Venezuela, in 1999, which each involved more than 20,000 fatalities ( Figure 2).
The size, speed, and destructive potential of geophysical mass flows can be strongly influenced by entrainment of bed material, and accounting for the causes and effects of entrainment is one of the greatest challenges in modeling flow dynamics Crosta et al., 2009]. Field evidence shows that a flow's volume (V) can increase many fold as a consequence of entrainment [e.g., Pierson et al., 1990;Wang et al., 2003;Hungr and Evans, 2004], and data compilations show that areas inundated downslope or downstream from mass flow source areas tend to increase in proportion to V 2/3 [Vallance and Scott, 1997;Legros, 2002;Griswold and Iverson, 2008]. Experimental evidence indicates that entrainment can also be accompanied by feedbacks that may decrease friction and increase flow momentum dramatically ( Figure 3 and Movie S1 in the supporting information). However, progress in understanding the detailed mechanics of entrainment has been hindered by the violent nature of full-scale geophysical mass flows and by the difficulty of making measurements at their evolving basal boundaries. Theory shows that potential exists for constraining entrainment mechanics by comparing measurements of evolving flow velocities and depths to model predictions, but only if the models themselves are soundly grounded . Thus, it is crucial for models of erosive mass flows to be founded on physical conservation laws that are properly formulated.
Our review focuses on mass flow models that employ depth-integrated conservation equations that are generalized from those of classical shallow-water theory [cf. Stoker, 1958;Savage and Hutter, 1989;Vreugdenhil, 1994;Iverson and Denlinger, 2001]. In recent decades depth-integrated mass flow models have become widely used tools in hazard evaluation [e.g., Pudasaini and Hutter, 2007;Medina et al., 2008;Pastor et al., 2014], and they offer three advantages over more detailed, three-dimensional models. First, they generate model output with a level of detail comparable to that of field measurements and large-scale experiments, enhancing the prospects for conclusive model testing [e.g., Iverson et al., 2010Iverson et al., , 2011McCoy et al., 2012;George and Iverson, 2014]. Second, they embed the effects of the evolving positions of the upper and basal flow boundaries directly into the governing conservation equations, eliminating the need to separately resolve domain boundaries [cf. Biscarini et al., 2010;Iverson and George, 2014]. Third, they reduce the degrees of freedom in the conservation equations themselves (commonly by neglecting bed-normal momentum components) and thereby reduce computation time. As a result, use of depth-integrated  [Rogers et al., 2010], (USGS photo by J. D. Spooner); (d) debris flow that incised preexisting landslide deposit, with red arrows identifying remnants of a severed, 5 m high sediment retention dam, Wenjia Gully, Sichuan, China, 2010 [Xu, 2010;Yu et al., 2013], (photo by Xu Qiang, used by permission of the Chinese Journal of Engineering Geology). equations in combination with shock-capturing numerical methods and adaptive mesh refinement enables high-resolution simulations of large-scale geophysical flows to execute quickly on ordinary desktop computers George, 2011;LeVeque et al., 2011]. Fast simulations can be crucial if hazard assessments are performed with short lead times or if probabilistic modeling strategies are used to examine the effects of poorly constrained initial conditions or parameter values [e.g., Ancey, 2005;Dalbey et al., 2008]. Although depth-integrated models of erosive mass flows have been widely used, inclusion of boundary flux terms in their mass and momentum conservation equations has been a source of controversy [Iverson, 2012]. Disagreement exists, for example, on whether boundary flux terms related to mass evolution should appear in  Oblique aerial view of the town of Armero, Colombia, which was largely destroyed by a lahar that originated on Nevado del Ruiz volcano on 13 November 1985. Pierson et al. [1990] reported that the lahar roughly quadrupled in volume between the source area and depositional area shown here. (USGS photo by R. J. Janda.) momentum conservation equations. Fundamental disparities in models' governing equations can lead to large differences in numerical predictions that may involve millions of calculations of mass and momentum fluxes. Therefore, systematic comparison of existing equations and establishment of the correct forms of applicable equations are priorities.

Historical Perspective
The qualitative importance of bed material entrainment by evolving mass flows has been recognized for more than a century. For example, Rickmers [1913] provided a lyrical description of entrainment by debris flows, for which he coined the term "mudspates": When a… slope of grit and shingle has been soaked like a sponge by rain or melting snows there may come a time when it… slides off…Slipping into channels and gullies this mass…attains a higher speed and carries away soft material as well as rocks which it finds on its way. It is during this descent that the mudspate generally acquires its characteristic composition, for only by movement can an even mixture of liquid and solids be maintained.
And two early books in German, Die Muren ("Debris Flows") [Stiny, 1910] and Bergsturz und Menschenleben ("Landslides and Human Lives") [Heim, 1932], included thorough accounts of the scouring action of debris flows and rock avalanches, respectively. Development of quantitative mass flow models that consider coevolution of flow mass and momentum during entrainment of bed material did not begin until the 1960's. Much of the earliest work focused on snow avalanches and appeared in the Russian literature, as summarized by Eglit and Demidov [2005]. In the English language literature, Briukhanov et al. [1967] presented a one-dimensional (1-D) depth-integrated continuum model that addressed entrainment at avalanche fronts by using a Lagrangian reference frame that translated with the front speed. Brugnot and Pochat [1981] analyzed a similar problem by using a fixed (i.e., Eulerian) reference frame and 1-D mass and momentum conservation equations resembling those of classical shallow-water theory. Neither of these pioneering snow avalanche models was accompanied by mathematical derivations that revealed the physical basis of relationships between boundary mass fluxes and evolution of conserved variables, however. A similar lack of detailed derivations has continued to hamper interpretation of many recent models. Takahashi et al. [1987] were perhaps the first to use 1-D Eulerian conservation equations similar to those of shallow-water theory to model landslides and debris flows that entrain mass from erodible basal boundaries. Takahashi et al. [1987] included a flow momentum loss term proportional to the rate of mass entrainment but did not explain the term's physical basis [see also Takahashi, 1991, p. 76]. Over the next several years, Lagrangian point-mass models were introduced with the goal of clarifying the behavior of landslides and debris flows that exchange mass with their beds [Cannon and Savage, 1988;Van Gassen and Cruden, 1989]. The equations used in these models implied that evolution of flow mass could influence flow momentum in a manner analogous to mass ejection by a rocket-an implication that spawned ensuing criticism [e.g., Hungr, 1990aHungr, , 1990b. Importantly, Erlichson [1991] noted that a fundamental shortcoming of these models was a failure to distinguish clearly between the action of internal and external forces that enable mass change [cf. Iverson, 2012].
Multidimensional, depth-integrated continuum models have now almost entirely supplanted 1-D and point-mass models as the standard tools for analyzing the behavior of geophysical mass flows, but the crux issue identified by Erlichson [1991] remains a source of discrepancies and misunderstandings. In many depth-integrated models, distinctions between internal and external forces are ambiguous, because they are linked to nuanced (and commonly unstated) variations in authors' definitions of boundary tractions, boundary fluxes, and conserved variables. Our assessment of depth-integrated models shows that the forms of momentum change terms associated with entrainment are contingent on these definitions as well as on the frame of reference.

Current Objective
We critically evaluate depth-integrated models of erosive mass flows by identifying constraints derived from a more comprehensive, two-layer, depth-integrated theory. The theory considers the interaction of a static layer of bed material with an overriding, mobile layer that may undergo changes in bulk density as it exchanges mass and momentum with the lower layer. In debris flows and rock avalanches the upper layer may span the entire flow thickness. In other cases the upper layer may represent a zone in which bed load motion is driven by both gravitational forces and boundary tractions exerted from above ( Figure 4). In either case the two-layer theory addresses two key questions: how do the dynamics of the moving upper layer respond to entrainment or deposition? And how does conservation of mass and momentum constrain the downward propagation of an erosion front (or upward propagation of a deposition front) where moving material contacts static bed material? We show that precise definitions of boundary tractions and boundary fluxes of conserved variables, as well as precise evaluations of mass and momentum jump conditions at eroding boundaries, play essential roles in answering these questions. These definitions and evaluations also enable reconciliation of some seemingly disparate models, as we describe in section 5 of this paper.

Qualitative Summary of Existing Models
A rapidly expanding body of geoscience, engineering, and applied mathematics literature presents depth-integrated mathematical models of evolving Earth-surface flows that entrain bed material. Table 1 summarizes some key qualitative features of the 31 published models we have examined. Each model considers simultaneous evolution of flow mass and momentum as well as flow-bed momentum exchange. The momentum balances and frames of reference used in the models have important differences, however. For example, some momentum equations are referenced to Lagrangian frames that translate downslope with the local, depth-averaged flow velocity. If entrained mass has a velocity that differs from that of the depth-averaged velocity as it enters this reference frame, then the left-hand sides of these equations are not fully conservative. As a result, the right-hand sides of Lagrangian equations must subtract fluxes of relative momentum at flow-bed boundaries in order to express momentum conservation during entrainment of static bed material [e.g., McDougall and Hungr, 2005]. We quantify this effect in sections 4 and 5 of this paper.
A wider range of possibilities exists in momentum equations formulated in Eulerian reference frames ( Table 1). The left-hand sides of these equations can express momentum conservation precisely, eliminating one complication in Lagrangian models. Eulerian formulations can also distinguish clearly between depth-integrated flow velocities and velocities of material crossing flow-bed boundaries. Nevertheless, existing Eulerian models differ widely. Some models assume that material entrained into the base of a flow automatically attains the depth-averaged downstream velocity of the flow, and an entrainment term derived by using this assumption implies that bed material entering a flow adds momentum to the flow [e.g., Gray, 2001]. A larger number of Eulerian models include an entrainment term that appears to subtract momentum from the flow. Modelers generally justify use of such a term by reasoning that flow momentum must be transferred to the bed to accelerate static bed material during entrainment, thereby diminishing the momentum of the flow itself [Cao et al., 2004;Benkhaldoun et al., 2011;Pirulli and Pastor, 2012]. Other Eulerian models of erosive flows use depth-integrated momentum balance equations that contain no explicit term related to entrainment [e.g., Capart and Young, 1998;Sovilla et al., 2006]. Few models formulated in Figure 4. (a and b) Schematic depictions of a moving sediment-water layer of thickness h 1 and bulk density ρ 1 that can exchange mass and momentum with a static bed layer of thickness h 2 and bulk density ρ 2 . The total thickness of the two-layer system is H = h 1 + h 2 . The position of the evolving boundary where the moving layer contacts static bed material is shown in red and is located at z = z b (x, y, t), where y is normal to the x-z plane. As shown in Figure 4b, an externally imposed shear traction τ may act at z = H. Here z denotes the direction of depth integration, x denotes the direction of momentum flux (normal to z), ρ 1 denotes the depth-averaged flow bulk density, and LHS and RHS, respectively, denote the left-hand slide and right-hand side of the x momentum equation. either Eulerian or Lagrangian frames make clear-and requisite-distinctions between boundary shear tractions that act during entrainment and flow over a fixed bed. In Table 1 we denote the interrelated issues concerning boundary momentum fluxes and boundary tractions as "Ia" and "Ib," and we quantify these issues in sections 3 and 4 of this paper.

Reviews of Geophysics
Other issues we identify in Table 1 and later quantify are (II) neglect of the effect of differences between depth-averaged and boundary values of velocities on interlayer momentum exchange, (III) neglect of the effect of bed-normal velocities that are necessarily associated with changes in bulk density that may accompany entrainment, and (IV) neglect of vertical accelerations when motion on steep slopes is modeled using an Earth-centered coordinate system in which z is vertical. We quantify this last issue in Appendix A. We reserve for Appendix B a detailed mathematical summary of some key models listed in Table 1, because interpretation of the mathematics is facilitated by first establishing a common set of physical principles, definitions, and notation.

Quantitative Model Evaluation
We establish pertinent principles, definitions, and notation by deriving a two-layer model that demonstrates how fluxes of mass and momentum across layer boundaries are constrained by three types of jump conditions. Two of the three conditions are closely analogous to Rankine-Hugoniot conditions for the onedimensional Euler equations that are used to model the dynamics of ideal gases [e.g., Godlewski and Raviart, 1996]. The third jump condition lacks a gas-dynamics analog because grain-fluid mixtures exhibit no clear relationship between pressure and bulk density (i.e., there is no equation of state). Instead, the third jump condition results from the existence of finite shear tractions that transfer tangential momentum between layers. All applicable jump conditions are derived from physical conservation laws, and all must be satisfied by entrainment rate formulas that are used in conjunction with depth-integrated conservation equations.
Our strategy is first to derive depth-integrated mass and momentum conservation equations for an individual layer that can exchange mass and momentum with adjacent layers. The derivation assumes that the evolving bulk density of the layer is not a function of z (Figure 4), although density stratification can be addressed by splitting any single layer into multiple layers with differing bulk densities. The next step in the derivation entails adding conservation equations for adjacent layers to obtain conservation equations that apply to a two-layer system as a whole. To ensure conservation of mass and momentum within the two-layer system, specific terms must cancel from the two-layer equations, and the sum of these canceled terms yields jump equations that must be satisfied at the interface between the layers. From a mathematical perspective, this method of obtaining jump equations is analogous to a more general, vector-calculus procedure used to obtain jump conditions by integrating around a perimeter that encloses a singular surface where two continua contact one another [e.g., Chadwick, 1999]. From a geophysical perspective, our depth-integrated method of obtaining jump conditions has the advantage of linking the conditions directly to the conservation equations commonly used to model Earth-surface flows.

Boundary Definition
All depth-integrated models of erosive flows are based on the premise that the evolving boundary between a moving flow and static bed material can be idealized as a sharp interface. It is worthwhile to examine the conceptual basis of this idealization, because the position of an evolving boundary where a flowing grain-fluid mixture interacts with a static grain-fluid mixture can appear indistinct [Capart and Fraccarollo, 2011]. Definition of the boundary position can also be complicated by shifting and jostling of grains in the bed before they are fully entrained [e.g., Armanini et al., 2005]. Nevertheless, the evolving boundary position must be defined precisely in depth-integrated mathematical models, including those that explicitly consider materials with two-phase compositions [e.g., Spinewine, 2005;Perng and Capart, 2008]. In subsequent sections of this paper, we analyze a scenario in which an arbitrarily shaped boundary forms the interface between flow and bed layers with arbitrary compositions. Here as a prelude to detailed analysis, we illustrate some underlying concepts by considering an example that involves a relatively simple material passing through a planar boundary.
Our example shows how an eroding boundary in a two-phase solid-fluid mixture viewed at the grain scale can be idealized by a sharp interface viewed at the continuum scale. Moreover, it shows that a two-phase material passing through the interface can be idealized as a single-phase material in which the bulk density evolves as entrainment proceeds. The example considers finite rectangular segments of a flow layer (denoted by subscript 1) and underlying bed layer (denoted by subscript 2) that each have a length 3δ, height 2δ, width δ, and volume 6δ 3 before any bed erosion occurs ( Figure 5). The bed layer initially contains six identical spherical grains, each with constant density ρ s , diameter δ, volume (π/6)δ 3 , and mass ρ s (π/6)δ 3 . The space between the grains is filled with fluid of constant density ρ f . Therefore, the bulk density of the bed layer is given by ρ 2 = ρ f + (π/6)[ρ s À ρ f ], where π/6 is the solid volume fraction of the bed. Before erosion occurs, the flow layer contains fluid alone, such that its initial bulk density is given by ρ 1 = ρ f . As erosion occurs and grains are entrained, the value of ρ 1 evolves.

Reviews of Geophysics
Our example assumes that erosion occurs steadily and entails migration of one grain from the bed layer to the flow layer during each of a series of discrete, equal time steps. In some instances, such as at t = 1 and t = 2 in Figure 5, this migration may involve parts of multiple grains that together represent one grain volume, (π/6)δ 3 . A volume of fluid [1 À (π/6)]δ 3 accompanies each grain as it migrates from the bed layer to the flow layer so that the flow layer volume increases by δ 3 and the bed layer volume decreases by δ 3 during each time step. (This proportionate flux of fluid and solid from the bed layer into the flow layer is necessary for the static bed to maintain a fixed bulk density ρ 2 = ρ f + (π/6)(ρ s À ρ f ) as it loses mass during erosion.) The total volume flux across the flow-bed boundary indicates that the boundary position changes by (1/3)δ during each time step, irrespective of whether the boundary passes through individual grains ( Figure 5).
The volume flux described above indicates that migration of material from the bed to the flow causes the flow layer mass m 1 to increase by ρ 2 δ 3 during each time step. More generally, m 1 can be expressed as a continuous function of time as where 6δ 3 ρ f is the initial flow layer mass, h 1 (t) is the flow layer thickness at time t, 3δ 2 h 1 (t) is the flow layer volume at time t, and 3δ 2 ρ 2 [h 1 (t) À 2δ] is the cumulative mass of material entrained into the flow layer up to time t. Some simple algebra shows that these relationships imply that the evolving bulk density of the flow layer ρ 1 (t) = m 1 (t)/[3δ 2 h 1 (t)] is given by Differentiation of (2) then yields an equation showing how mass conservation links the rate of boundary evolution, which is described by dh 1 /dt, to the rate of flow layer bulk density evolution, dρ 1 /dt: Coevolution of the flow bulk density and eroding boundary position accounts for the evolving two-phase composition of material that enters the flow layer from the bed layer. We now consider the same mass conservation and boundary evolution example analyzed above, but we adopt a continuum-scale perspective analogous to that employed in subsequent sections of this paper. During each time increment shown in Figure 5, the change in mass of the flow layer is given by (ρ 1 + Δρ 1 )(h 1 + Δh 1 ) À ρ 1 h 1 , and the change in mass of the bed layer is given by À ρ 2 Δh 1 , where Δ denotes an incremental change in a variable during the time increment Δt.
The mass changes that occur in the flow and bed layers during Δt must total zero, indicating that the difference equation ρ 1 Δh 1 + Δρ 1 h 1 + Δρ 1 Δh 1 À ρ 2 Δh 1 = 0 applies. The difference equation can be converted to a differential equation by dividing each term by Δt, making the substitutions Δρ 1 = ρ 1 (t + Δt) À ρ 1 (t) and Δh 1 = h 1 (t + Δt) À h 1 (t), and then taking the limit as Δt → 0. The result can be written as This equation can be reduced to (3) by using the substitution ρ 2 À ρ 1 = [2δ(ρ 2 À ρ f )]/h 1 , which is obtained from (2). The equivalence of (4) and (3) demonstrates that a continuum-scale perspective that omits consideration of separate solid and fluid phases produces no error in reckoning the effects of mass conservation on the boundary migration rate.
Despite their mathematical equivalence, equations (3) and (4) reflect some important conceptual differences. Unlike (3), (4) results from an assessment of mass conservation alone, without any reference to volume conservation. Indeed, (4) applies to mass conservation across an evolving planar boundary separating continuous materials with arbitrary compositions. Moreover, because it arises from a mathematical limit in which Δt → 0, (4) applies irrespective of whether boundary erosion is steady or unsteady. Equation (4) can also be viewed as a simple jump condition, because it shows how mass conservation relates motion of the boundary separating the flow and bed layers to the evolving contrast (or jump) in the layers' bulk densities, ρ 2 À ρ 1 . More general jump conditions, which we derive below, apply in most geophysical contexts.

Boundary Kinematics
In order to derive jump conditions and general, depth-integrated conservation equations, it is necessary to specify the kinematics of an evolving boundary that has an arbitrary shape. As part of this specification, it is critical to distinguish between the basal erosion rate and the bed material entrainment rate ( Figure 6). The basal erosion rate is simply the negative of the local rate of bottom boundary elevation change, À ∂z bot /∂t. By contrast, the bed material entrainment rate E is the volumetric flux of material per unit area that enters a layer at its boundary. By definition, E is positive upward, normal to the local boundary, such that E > 0 indicates volume gain by a layer that entrains material at its bottom boundary, and also represents volume loss by a layer that erodes at its top boundary. Volume is not necessarily conserved as material passes through a boundary, however, and, in general, E ≠ À ∂z bot /∂t applies.
Kinematic boundary conditions that are somewhat more general than those commonly used in fluid mechanics express relationships between rates of boundary elevation change, rates of material flux through each boundary, and the velocity components of material adjacent to each boundary. If the material consists Figure 6. Schematic depiction of an eroding boundary (red) located at z = z b = z 2top = z 1bot , where subscripts 1 and 2 denote layers 1 and 2. The boundary elevation changes at rate ∂z b /∂t, and the x and z velocity components of material entrained into layer 1 at the boundary are u 1 (z 1bot ) and w 1 (z 1bot ), respectively. The mass entrainment rate per unit area, ρ 1 E 1bot ¼ ρ 2 E 2top , is the bed-normal mass flux of material that leaves the top of layer 2 and enters the bottom of layer 1.
of a two-phase solid-fluid mixture, then its velocity → u is defined by a mass-weighted average of the solid phase velocity → u s and fluid-phase velocity → u f [Iverson, 1997]: where n is the fluid volume fraction, and ρ is the mixture bulk density. Henceforth, we employ the Cartesian components of → u, which are u, v, and w in the x, y, and z directions, respectively.
Appendix C shows how universal forms of the kinematic boundary conditions can be derived from a three-dimensional analysis of boundary evolution. It also shows that specialized forms of the conditions, appropriate for application in depth-averaged models that use Cartesian velocity components, can be expressed as where z top and z bot denote the positions of the top and bottom boundaries of a layer, and ξ top = [1 + (∂z top /∂x) 2 + (∂z top /∂y) 2 ] 1/2 and ξ bot = [1 + (∂z bot /∂x) 2 + (∂z bot /∂y) 2 ] 1/2 are coefficients that account for the fact that E top and E bot describe volumetric fluxes per unit boundary area normal to the local boundary (Appendix C). In simplest terms, ξ top and ξ bot can be interpreted as geometric correction factors, and if (∂z bot /∂x) 2 , (∂z bot /∂y) 2 , (∂z top /∂x) 2 and (∂z top /∂y) 2 are much smaller than 1, then ξ top ≈ ξ bot ≈ 1. (6) and (7) are evaluated at the boundaries z top and z bot , but physically, the velocities apply to material immediately adjacent to boundaries, because no material can occupy a vanishingly thin boundary surface. Mathematically, a discontinuity in at least one of these velocity components must exist if a boundary between static and moving material is present. As shown in section 3.1, a mathematical idealization that uses a sharp boundary to represent the transition from a static bed layer to a moving flow layer is valid at a continuum scale, irrespective of whether the boundary appears somewhat ambiguous when viewed at the grain scale.

Conservation of Mass and Accompanying Jump Condition
We use overstrikes to denote depth-averaged values of the dependent variables u, v, w, and ρ. Although ρ is not a function of z within any individual layer, our conservation equations account for evolution of ρ as a function of x, y, and t. Moreover, ρ can change abruptly as a function of z when mass passes between layers. Under these conditions, mass conservation within an individual layer of thickness h(x, y, t) bounded by moving surfaces at z top and z bot can be expressed by The last line of (8) is obtained from the preceding lines by using the kinematic boundary conditions (6) and (7) as well as the substitution ∂h/∂t = ∂(z top À z bot )/∂t. The terms ρξ top E top and Àρξ bot E bot in the last line of (8) describe fluxes of mass into or out of the layer.

Reviews of Geophysics
Next we consider forms of (8) that apply when mass exchange occurs across a single boundary where an upper layer (designated by subscript 1) contacts an adjacent underlying layer (designated by subscript 2). For the upper layer the last line of (8) reduces to and for the lower layer it reduces to Because mass is conserved in the two-layer system as a whole, the source terms on the right-hand sides of (9) and (10) must sum to zero. This summation yields an interface jump equation, which must be satisfied in addition to (9) and (10). Because the bottom boundary of layer 1 is coincident with the top boundary of layer 2, ξ 1bot = ξ 2top applies. Thus, (11) shows that E 1 bot = E 2 top is not satisfied unless ρ 1 ¼ ρ 2 .
Equation (11) can be expanded to obtain an explicit relationship between velocity jumps, bulk density jumps, and evolution of the interface between layers 1 and 2. By substituting (6) and (7) into (11) and setting z 1 bot = z 2 top = z b , the mass jump condition can be expressed as A simpler form of (12) results from assuming that a planar interface with ∂z b /∂x = ∂z b /∂y = 0 exists: A still-simpler form of this jump condition applies if layer 2 consists of static bed material, such that w 2 (z b ) = 0: If layer 2 is static but the condition ∂z b /∂x = ∂z b /∂y = 0 is not satisfied, then w 1 (z b ) is replaced by /∂y] on the right-hand side of (14).
Equations (13) and (14) show how the bed erosion rate À ∂z b /∂t is related to differences in mass density ρ and the normal component of material velocity w(z b ) across the boundary. They are exactly analogous to one of the Rankine-Hugoniot conditions for the one-dimensional Euler equations [e.g., Godlewski and Raviart, 1996]. Therefore, in (13) and (14) ∂z b /∂t is analogous to the propagation speed of a one-dimensional shock wave in an ideal gas. Like a shock wave, an erosion front represents a mathematical discontinuity, which must move at the speed given by (12)-(14) in order satisfy conservation of mass. However, in the presence of a density contrast ρ 1 < ρ 2 and boundary material velocity w 1 (z bot ) > 0, an erosion front described by (14) moves downward into a static bed layer-opposite to the propagation direction of a shock wave in a gas.

Conservation of x Momentum and Accompanying Jump Condition
Evolution of the downstream component of momentum in a layer with depth-invariant bulk density ρ is described by where ∫ ztop z bot ΣF x dz represents the sum of x direction forces acting on the layer. The boundary momentum exchange terms ρu z top À Á ξ top E top and À ρu z bot ð Þξ bot E bot arise in the last line of (15) as a result of using the kinematic boundary conditions (6) and (7) to combine several terms in the second and third lines of (15). Importantly, these momentum exchange terms involve the boundary velocities u(z top ) and u(z bot ) of material that enters the layer at z top and z bot , not the depth-averaged velocity of the layer, u. Thus, if material with no x direction velocity enters the layer during entrainment, it imparts no x direction momentum to the layer.
The integrals in the fourth line of (15) result from use of the identities which show how the depth integrals of products are related to the products of depth integrals. Physically, the integrands (u À u) 2 and u À u ð Þ v À v ð Þdescribe the effects of differential advection of momentum due to variations of u and v with z [cf. Vreugdenhil, 1994]. Many depth-averaged flow theories neglect these variations, but this neglect can lead to oversights when evaluating the effects of momentum exchange between adjacent layers during entrainment. Therefore, we adopt another approach commonly used in depth-integrated flow theories, which accounts for variations of u and v with z by using Boussinesq momentum distribution coefficients defined as With these substitutions (15) can be rewritten as which is the form we employ below.
The forcing term ∫ ztop zbot ΣF x dz in (17) represents the sum of the downslope gravitational driving force per unit volume ρg x h and the resisting stress gradients that develop in reaction to this forcing. By using Leibniz' rule to evaluate depth integrals, this sum can be expressed as where the normal-stress component σ xx is defined as positive in compression, as is customary in geophysics, but the shear stress components τ yx and τ zx are defined using a continuum mechanics sign (18) Reviews of Geophysics convention [cf. Malvern, 1969;Iverson, 2012]. Hence, terms involving normal-stress and shear stress components have opposite signs in the first line of (18).
Detailed evaluation of the stress components in (18) requires the use of a constitutive model, and in Appendix D we describe one such model. Here, however, to minimize complications, we avoid constitutive assumptions and instead invoke a standard thin-layer scaling argument to infer that a valid approximation of (18) is In this approximation, neglected terms are smaller than retained terms by a factor ε, where ε is the ratio of the characteristic thickness to the characteristic length or breadth of a layer [cf. Savage and Hutter, 1989;Iverson, 2013a]. Values ε ≪ 1 commonly apply in many types of Earth-surface flows [Pudasaini and Hutter, 2007]. Our use of the thin-layer approximation helps highlight the essential features of the results that follow, but in Appendix D we describe effects of terms neglected in (19).
Leading-order (in ε) x momentum equations for adjacent layers that interact at a mutual boundary (where z 1 bot = z 2 top ) are obtained by combining (19) with (17) and using subscripts 1 and 2 to denote the upper and lower layers, respectively: The right-hand side of (20) shows that layer 1 gains momentum at the rate ρ 1 u 1 z 1bot ð Þξ 1bot E 1 bot as a result of bed material entrainment [cf. Gray, 2001], while (21) shows that layer 2 loses momentum at the rate À ρ 2 u 2 z 2top À Á ξ 2top E 2 top . These gains and losses are somewhat illusory, however; the physical implications of the entrainment process become clearer when (20) and (21) are added to obtain an x momentum equation for the two-layer system as a whole: Because the only external forces acting on the two-layer system described by (22) are those due to gravity and the boundary shear tractions τ 1 zx top and τ 2 zx bot , other terms that appear on the right-hand sides of (20) and (21) cancel one another. Thus, (22) shows that no net x momentum is gained or lost as a consequence of entrainment.
The jump condition for x momentum at the interface between layers 1 and 2 is the sum of the terms that were canceled from (22) when it was obtained by addition of (20) and (21): Substituting (6) and (7) into (23) and using the abbreviated notation z 2 top = z 1 bot = z b , τ 1 zx bot = τ 1bot , and τ 2 zx top = τ 2top yields a more explicit form of the equation: Reviews of Geophysics This jump condition can be recast as a relatively simple erosion rate equation if the eroding interface satisfies ∂z b /∂x = ∂z b /∂y = 0: If, in addition, layer 2 is static, then (25) simplifies further to If layer 2 is static, but the condition ∂z b /∂x = ∂z b /∂y = 0 is not satisfied, then /∂y] is added to the right-hand side of (26).
In (26) À ∂z b /∂t describes the downward erosion rate in a fixed, Eulerian reference frame, whereas À ∂z b /∂t + w 1 (z b ) is the bed material entrainment velocity-a quantity that depends on the upward material velocity w 1 (z b ) ( Figure 6). Indeed, by neglecting w 1 (z b ) and using the substitution E = À ∂z b /∂t, (26) can be written in the alternative form E ¼ τ 1 bot À τ 2 top À Á =ρ 1 u 1 z b ð Þ, which is equivalent to entrainment rate equations obtained previously by considering x momentum jumps and neglecting z momentum [e.g., Fraccarollo and Capart, 2002;Iverson, 2012].

Conservation of y Momentum and Accompanying Jump Condition
The depth-integrated equations for the y momentum component and attendant jump condition are precisely analogous to those derived above for the x component. Thus, interchanging x and y as well as u and v in the preceding section yields the y component equations. These equations provide no additional insight, however, and we do not address them further.

Conservation of z Momentum and Accompanying Jump Condition
The equation describing conservation of the bed-normal component of momentum in a layer with depth-invariant bulk density ρ is The forcing term ∫ ztop zbot ΣF z dz in (27) can be expressed as where the normal-stress component σ zz is defined as positive in compression. Note that ∂σ zz /∂z is negative if σ zz counteracts the effects of gravity, because g z < 0.
Reviews of Geophysics We next adopt a rationale like that used to obtain the x momentum jump condition (23). This procedure yields leading-order z momentum equations for layers 1 and 2, and addition of these two equations yields a merged, two-layer, z momentum equation, Identification of terms canceled from (29) yields an equation describing the z momentum jump between layers 1 and 2: Expanding (30) by using (6) and (7) yields a more explicit form of the equation, For a planar interface with ∂z b /∂x = ∂z b /∂y = 0, (31) reduces to Like the mass jump equation (14), (32) is closely analogous to one of the Rankine-Hugoniot jump conditions for the one-dimensional Euler equations. However, consistent with (14), the erosion front speed À ∂z b /∂t in (32) describes motion in a direction opposite to that of a shock wave in an ideal gas. If layer 2 is static, then the jump condition (32) reduces further to If layer 2 is static but the condition ∂z b /∂x = ∂z b /∂y = 0 is not satisfied, then /∂y] is added to the right-hand side of (33).
Equation (33) shows that if w 1 (z b ) > 0, then bed erosion occurs only if σ 2top À σ 1bot > ρ 1 w 1 2 (z b ). On the other hand, if w 1 (z b ) < 0, then bed erosion occurs only if σ 2top À σ 1bot < ρ 1 w 1 2 (z b ). These possibilities illustrate the crucial role of z momentum in constraining the behavior of an erosional or depositional interface.

Summary of Jump Conditions and Their Implications
The foregoing analysis identifies three jump conditions that must be satisfied simultaneously at an evolving interface between adjacent layers of material with ρ 1 ≠ρ 2 . To examine the implications of these conditions for bed sediment entrainment, we focus on forms that apply if the bottom layer is static. Table 2 summarizes dimensional and normalized forms of the conditions for cases in which ∂z b /∂x = ∂z b /∂y = 0 is satisfied and also lists normalized forms that apply if ∂z b /∂x = ∂z b /∂y = 0 is not satisfied and the additional term Moving upper layer is denoted by subscript 1, stationary lower layer is denoted by subscript 2, and dimensionless quantities denoted by asterisks are defined in the text

Reviews of Geophysics
All normalizations employ the velocity scale u 1 (z b ), because u 1 (z b ) must be finite for the notion of a flow-bed boundary to have any physical meaning.
The normalized jump conditions in Table 2 contain six independent, dimensionless quantities, defined as Generally, these quantities have straightforward physical interpretations that require no elaboration. However, w * has a special interpretation in the context of granular mechanics, because it can be regarded as a dilatancy angle that relates the bed-normal and tangential velocity components of material entrained at the interface between layers 1 and 2. More specifically, w * > 0 indicates dilative behavior, w * < 0 indicates contractive behavior, and w * = 0 indicates behavior without any dilation or contraction. Table 2 can be summarized graphically for cases in which A = ∂z b /∂x = ∂z b /∂y = 0 is satisfied. Figures 7-9 depict the signs of interface discontinuities required by the jumps as wells as hypothetical profiles of dimensional variables on either side of the jumps for three distinct cases in which bed erosion occurs. (For simplicity the figures assume that the profiles of all variables are linear, but this need not be the case; the jump conditions dictate only the changes in values of variables at the interface, where z = z b .) Figure 7 illustrates a baseline case in which the upper layer has the same bulk density as the underlying layer, such that entrainment involves jumps in only the boundary shear traction and tangential velocity. Figure 8 illustrates a case in which the lower layer is denser than the upper layer, such that material dilates during entrainment and jumps in all variables occur at the boundary. Figure 9 illustrates the opposite case in which the lower layer is loose and contracts as it is entrained. Owing to the effects of w 1 (z b ) < 0, the contractive case necessarily involves smaller jumps in shear traction τ and normal traction σ than in a dilative case with an identical erosion rate. Sediment contraction during entrainment also implies the existence of a higher normal traction σ on the upper side than on the lower side of the interface. This counterintuitive condition reflects the lack of balance between gravitational forcing and normal stress that is associated with finite downward momentum, ρ 1 w 1 z b ð Þ < 0.

Some basic implications of the jump conditions in
Under ideal conditions, observable properties of eroding interfaces include the position of the boundary between moving and static material, velocity components of material adjacent to the boundary, and bulk density contrasts [Capart and Fraccarollo, 2011].  . Graphical depiction of jump conditions for a baseline case in which shearing layer 1 maintains the same bulk density as underlying static layer 2. As a consequence, erosion at the boundary between layers 1 and 2 (red) involves jumps in only shear stress and downstream velocity. By algebraically manipulating the normalized jump conditions for ∂z b /∂x = ∂z b /∂y = 0 listed in Table 2, these properties can be expressed entirely in terms of the jumps in normal and shear tractions acting at the eroding boundary:

Reviews of Geophysics
Equation (34) shows that the normalized erosion rate Àż* depends on the normalized excess shear traction acting at the boundary τ * minus the resisting effects of dilatancy expressed by (35), σ */τ * = w *. Indeed, (35) can be viewed as a stress-dilatancy relationship applicable at an eroding interface [cf. Rowe, 1962]. If the effects of dilatancy are large enough that σ * = τ * 2 , then (34) indicates that no erosion occurs. However, (36) shows that the condition σ * = τ * 2 cannot be realized except in the limit ρ * → 0. Therefore, for any material with finite ρ 1 , the conditions (34)-(36) not only show that erosion will inevitably occur if τ * > 0 but also show that the rate of erosion will be restricted if it is accompanied by positive dilatancy-irrespective of whether sediment dilation during entrainment involves pore pressure feedback of the type considered by van Rhee [2010] or Iverson [2012].
In cases in which sediment entrainment involves negative dilatancy (i.e., contraction), (34) and (35) show that the rate of bed erosion is enhanced because w * < 0. Contractive motion during entrainment can occur, for example, when loose sediment beds are overrun by debris flows or avalanches [Wang et al., 2003;Iverson et al., 2011]. In this case σ */τ * < 0 applies, and (36) therefore requires that ρ * < À 1. This restriction indicates that ρ 1 > ρ 2 must apply-as is necessary if bed material contracts as it is entrained and merged with a denser overriding flow.

Implications for Depth-Integrated Conservation Equations
Substitution of the jump conditions into the depth-integrated mass and momentum conservation equations demonstrates that the effects of entrainment can be expressed in several different ways. For example, substitution of (14) into (9) yields the mass conservation equation for layer 1: Three forms of the source term on the right-hand side of this equation are shown in order to highlight relationships between E, w 1 (z b ), and the bulk density contrast ρ 2 ≠ρ 1 . The first form shows that mass entering layer 1 during bed erosion depends on the bulk density ρ 2 of layer 2, whereas the second form shows that if ρ 2 ≠ρ 1 , then the mass transfer between layers necessarily involves dilatancy (i.e., w 1 (z b ) ≠ 0). The third form shows that if ρ 2 ≠ρ 1 , then mass conservation couched in terms of the volumetric entrainment rate E depends on the bulk density ρ 1 rather than on the bulk density ρ 2 . The physical reasons for these differing representations of entrainment terms in (37) are quite transparent, but reasons for analogous differences in momentum conservation equations are subtler and a source of considerable disagreement (Table 1). Figure 9. Graphical depiction of jump conditions for a case in which layer 1 is denser than layer 2, indicating that bed material contracts as it is entrained. Jumps in all quantities therefore exist at the eroding boundary (red) between layer 1 and layer 2. The signs of three of the five jumps are opposite of those illustrated in Figure 8. The x momentum conservation equation for layer 1 is obtained by substituting (26) into (20). As in (37), three forms of the terms associated with entrainment are shown on the right-hand side of the resulting equation:

Reviews of Geophysics
The first form of the right-hand side of (38) includes no momentum source term. Instead, all effects of momentum exchange with the static bed are summarized by the basal shear traction, τ 2 zx top , which expresses the basal flow resistance engaged during the interaction of layer 1 with the eroding upper surface of layer 2. By contrast, the second and third lines of (38) each contain momentum gain terms proportional to ρ 1 u 1 z b ð Þ. These apparent gains arise only because the second and third lines of (38) also contain the basal shear traction term, τ 1 zx bot , which misrepresents the basal flow resistance engaged during entrainment. (The term τ 1 zx bot would be correct in the absence of entrainment, but if entrainment occurs, then the resisting traction is τ 2 zx top , and the jump conditions derived above show that τ 2 zx top ≠ τ 1 zx bot is mandatory.) Thus, although each of the forms of the right-hand side of (38) is mathematically correct, any change in momentum associated with entrainment results physically from attendant reduction of the resisting basal shear traction from τ 1 zx bot to τ 2 zx top . In no case does this reduction imply a loss of x momentum by layer 1 as a result of momentum transfer to the static bed, as is erroneously indicated by some of the models listed in Table 1.
The z momentum conservation equation for layer 1 is obtained in an analogous manner, yielding As in (38), the first form of the right-hand side of (39) shows that layer 1 gains or loses no momentum as a result of entrainment of static bed material, because the external force per unit area exerted on layer 1 by layer 2 is expressed entirely by the normal traction σ 2 zz top . The second and third forms of the right-hand side of (39) have implications like those of the second and third lines of (38): if the normal traction acting on layer 1 is expressed in terms of σ 1 zz bot instead of σ 2 zz top , then illusory momentum gain terms proportional to ρ 1 w 1 z b ð Þ appear.
The conservation equations (37)-(39) can be written in alternative forms that facilitate comparisons with the diverse models listed in Table 1. If the equations' left-hand sides are written in forms like those used in standard shallow-water equations (which assume that ρ 1 is constant), then additional source terms must appear on the equations' right-hand sides to account for the effects of any evolution of ρ 1 . Table 3 summarizes both the modified left-hand sides and the additional source terms that arise when the conservation equations are written in this way.
Equations (37)-(39) may also be written in Lagrangian forms (or closely analogous primitive Eulerian forms)provided that velocity profiles are assumed uniform, such that β = 1 applies everywhere. In this case some The notation d=dt ¼ ∂=∂t þ u 1 ∂=∂x þ v 1 ∂=∂y is used to denote a depth-averaged material time derivative.

10.1002/2013RG000447
terms cancel from the left-hand sides of (38) and (39) through use of (37), but the cancelation leads to additional momentum source terms on the right-hand sides of the equations (Table 4). Importantly, the additional source terms À ρ 1 u 1 E and À ρ 1 w 1 E shown in Table 4 cancel other terms on the right-hand sides of (38) and (39) if all motion is concentrated as basal slip (i.e., u 1 = u(z b ) and w 1 ¼ w z b ð Þ, such that β = 1) and the basal tractions τ 1 zx bot and σ 1 zz bot are used. On the other hand, use of the correct basal tractions τ 2 zx top and σ 2 zz top in the Lagrangian forms of (38) and (39) shows that the effects of the source terms À ρ 1 u 1 E and À ρ 1 w 1 E are physically meaningful: they represent the inertia of static bed material that must be accelerated to attain the velocity of the Lagrangian reference frame, as inferred by Hungr [1995]. No analogous inertial effect exists in an Eulerian reference frame-a contrast that has led to some confusion in the literature [Iverson, 2013b].

Implications for Bed Erosion Formulas
The foregoing results show that conservation of mass and momentum requires that the erosion rate À ∂z b /∂t of a planar bed with ∂z b /∂x = ∂z b /∂y = 0 satisfies (34). In dimensional form, (34) can be expressed as Any erosion rate formula that does not consider the influence of bed topography must be consistent with (40), and in Table 5 we use this inference to categorize various published erosion formulas. The second term on the right-hand side of (40) represents the effects of dilatancy in suppressing the erosion rate, and if this term is zero, then (40) reduces to an erosion rate formula À ∂z b =∂t ¼ τ 1 bot À τ 2 top À Á =ρ 1 u 1 z b ð Þ deduced by Fraccarollo and Capart [2002].
The presence of u 1 (z b ) in the denominator of the first term on the right-hand side of (40) might seem puzzling in two respects. First, it indicates that erosion rates decrease as the basal flow velocity u 1 (z b ) increases, provided that τ 1 bot À τ 2 top remains constant. This effect is a necessary consequence of x direction momentum conservation. However, it is important to recognize that τ 1 bot À τ 2 top might increase faster than u 1 (z b ) increases and that the relative rates of change of τ 1 bot À τ 2 top and u 1 (z b ) cannot be ascertained from conservation laws alone. Another implication of the first term on the right-hand side of (40) is that À ∂z b /∂t → ∞ as u 1 (z b ) → 0. Physically, this singular behavior implies that the flow layer merges with the static bed layer as the flow layer ceases to move. In this case the notions of a flow-bed boundary or erosion rate lack any conventional meaning.
Equation (40) can be written in an alternative form that may be more illuminating: Here u 1 (z b )[τ 1bot À τ 2top ] is the excess power per unit bed area that causes erosion (similar to the "excess stream power" of Harsine and Rose [1992]), and τ 1 ref ≡ ρ 1 [u 1 (z b )] 2 is a reference value of τ 1 bot . Physically, τ 1 ref can be interpreted as the value of τ 1 bot that would enable erosion to occur at a steady rate given by À ∂z b /∂t = u 1 (z b ) À w 1 (z b ) if entrainment were resisted only by the inertia of the entrained sediment, and no other forces were acting (an inference drawn from the second line of (38)). Thus, erosion rates described by (41) increase with flow power, but they are restricted by noninertial forces that influence τ 1 bot , and also by the effects of dilatancy.
The notation d=dt = ∂=∂t þ u 1 ∂=∂x þ v 1 ∂=∂y is used to denote a depth-averaged material time derivative. Erosion formulas that are less abstract and more explicit than (40) and (41) can be obtained only if constitutive assumptions are invoked to evaluate the magnitude of the traction jumps τ 1 bot À τ 2 top and σ 2 top À σ 1 bot . Here we do not attempt to find definitive constitutive relationships, but we describe implications of some logical possibilities. The simplest possibility is that there is no dilatancy or normal traction jump (i.e., the second term on the right-hand side of (40) is zero) and that τ 1 bot and τ 2 top each obey Coulomb friction equations,

Reviews of Geophysics
where ϕ 1 and ϕ 2 are friction angles for layers 1 and 2, respectively, and p 1 bot and p 2 top are pore-fluid pressures that may counteract the effects of σ 1 bot and σ 2 top . In this case an erosion rate formula like that of Iverson [2012] results from postulating lithostatically balanced total normal tractions (i.e., σ 1 bot ¼ σ 2top ¼ ρ 1 g Àz H À z b ð Þ , where z = H denotes the top of the two-layer system), and then substituting (42) and (43) into (40) to obtain This equation is well posed only if the conditions ϕ 1 ≤ ϕ 2 and p 1 bot tan ϕ 1 ≤ p 2 top tan ϕ 2 are both satisfied. (If they are not satisfied, then (44) implies either that z b grows exponentially with time when quantities other than z b (t) are constant or that z b > H when ∂z b /∂t = 0. Neither behavior is physically plausible.) We focus on the simplest well-posed case, in which ϕ 1 = ϕ 2 applies and (44) reduces to Formula is based on contrast between flow sediment concentration, c, and an equilibrium concentration, c e . E lacks explicit dependence on boundary tractions.
Formula is based on contrast between equilibrium bed slope θ e and ambient slope θ. E lacks explicit dependence on boundary tractions. Fraccarollo and Capart [2002] À Formula is compatible with (40) if bulk density contrast and dilatancy at eroding boundary are zero.
where α is a fitting parameter Empirical formula. E lacks explicit dependence on boundary tractions.

Cao et al. [2004]
where θ is Shields parameter and θ c is its critical value Empirical formula postulated to describe bed erosion in response to dam-break floods. E lacks explicit dependence on boundary tractions.

McDougall and Hungr
Empirical formula postulated to describe bed erosion by rapid landslides. E lacks explicit dependence on boundary tractions Formula postulated to describe bed erosion by snow avalanches. Similar to z momentum jump condition in Table 2.
where τ 1zxbot is obtained from Bingham; Herschel-Bulkey; or Voellmy model; Formula is compatible with (40) if bulk density contrast and dilatancy at eroding boundary are zero.
This formula indicates that erosion occurs if p 2top exceeds p 1bot , a condition that may develop as a result of undrained bed deformation caused by rapid loading by an overriding flow [e.g., Wang et al., 2003;Iverson et al., 2011;Reid et al., 2011]. Normalization of (44) using a length scale H and pore pressure scale ρ 1 g Àz H shows that the timescale of this erosion is u 1 (z b )/g À z . If u 1 (z b ) < 10 m/s, then values of u 1 (z b )/g À z are smaller than about 1 s, implying that erosion described by (44) responds almost instantly to changes in u 1 (z b ). Thus, the principal factor limiting erosion described by (44) is the capacity of overridden bed sediment to develop and sustain high pore pressures.
An alternative model of erosion results from retaining the assumption of zero dilatancy but assuming that the basal shear traction τ 1 bot in (40) increases as u 1 (z b ) increases. Several studies summarized by Fraccarollo and Capart [2002] have inferred that the traction obeys τ 1 bot ¼ Cρ 1 u 1 z b ð Þ ½ 2 , where C is a dimensionless proportionality coefficient that is typically smaller than 0.1. This quadratic dependence of shear traction on slip velocity (or shear rate) might be a consequence of grain collisions [e.g., Bagnold, 1954], but irrespective of this interpretation, τ 1 bot ¼ Cρ 1 u 1 z b ð Þ ½ 2 is a rational postulate based on dimensional reasoning. Use of this postulate in conjunction with (40), (43), and the assumptions of a hydrostatic pore-pressure distribution (i.e., p(z b ) = ρ w g À z (H À z b ) + p| z = H , where ρ w is the pore-fluid density and p| z = H denotes the fluid pressure at z = H) and lithostatically balanced total normal tractions (i.e., σ 1 bot ¼ σ This formula is equivalent to that derived by Fraccarollo and Capart [2002, p. 194], who used it successfully to predict key features of erosion patterns measured in laboratory experiments involving dam-break water floods. The formula indicates that Coulomb friction of bed material restricts erosion but that the restriction diminishes as u 1 (z b ) increases and the erosion rate becomes dictated mostly by the values of C and u 1 (z b ).
If all quantities except z b on the right-hand side of (45) are treated as constants, then (45) has an exponential-decay solution, z b t ð Þ ¼ z 0 À z steady À Á e Àt=t0 þ z steady , which indicates that an increment of erosion due to an imposed change in u 1 (z b ) causes z b to relax from the initial value z 0 toward the steady equilibrium value z steady ¼ H À C u 1 z b ð Þ ½ 2 = 1 À ρ w =ρ 1 ð Þg Àz tanϕ 2 ð Þ ½ . The timescale of this relaxation is , which differs little from the erosion timescale u 1 (z b )/g À z identified in conjunction with (44). Indeed, t 0 < 1 s generally applies for reasonable ranges of u 1 z b ð Þ; ρ w ; ρ 1 ; g Àz and ϕ 2 in (45). This finding implies that bed erosion described by (45) can largely keep pace with changes in flow dynamics [Fraccarollo and Capart, 2002]. Unlike erosion described by (44) with ϕ 2 = ϕ 1 , however, erosion described by (45) is a self-limiting relaxation process if u 1 (z b ) is constant.
For cases in which entrainment involves sediment dilation or contraction, plausible bed erosion formulas are more complicated than (44) or (45) because they must satisfy the three jump conditions (34)-(36) simultaneously. If the normal traction jump condition satisfies σ 1 bot < σ 2top ¼ ρ 1 g Àz H À z b ð Þþpj z¼H (implying lithostatic bed stresses and positive dilatancy during entrainment), then combination of (35), (40), and (43) with the postulate τ 1 bot ¼ Cρ 1 u 1 z b ð Þ ½ 2 of Fraccarollo and Capart [2002] leads to a generalized version of (45): In this formula w * is used to isolate the influence of bed sediment dilatancy, defined as w * = w 1 (z b )/u 1 (z b ), and the notation pj z¼z b is used to denote the pore pressure at z = z b . Equation (46) implies that w * > 0 causes a reduction in erosion rates, irrespective of pore pressure, but pore pressure feedback might produce additional effects owing to coupling between the dilation rate and pj z¼zb [Iverson, 2012].
Development of erosion rate formulas constrained by jump conditions remains in its earliest stages. However, it is important to emphasize that formulas such as (44)-(46) differ fundamentally from formulas generally used in fluvial erosion studies, which have focused on laboratory settings or river systems in which driving and resisting forces are essentially balanced and bed load transport rates equilibrate as part of this balance [e.g., Gilbert, 1914;García, 2008;Berzi and Fraccarollo, 2013;Revil-Baudard and Chauchat, 2013]. Most investigations of the bed sediment entrainment process have also focused on a narrow range of conditions close to equilibrium [e.g., Pugh and Wilson, 1999;Lamb et al., 2008;Capart and Fraccarollo, 2011]. Early analyses of bed mobilization by mass flows took a similar approach; they implicitly adopted Bagnold's [1956] postulate, which states that basal shear tractions exerted by overriding flows are balanced by resisting shear tractions engaged during the onset of bed sediment motion [e.g., Takahashi, 1978;Bovis and Dagg, 1992]. However, this postulate can be inappropriate even for steady states [Seminara et al., 2002;Parker et al., 2003], and it is clearly inappropriate for unsteady mass flows in which entrainment can occur under conditions far from equilibrium. In such flows momentum and mass evolve interdependently, as reflected by the unbalanced boundary tractions that are explicit in (40) and (41) and implicit in (44)-(46).

Reconciliation of Disparate Models
Equations (37)-(39) differ from the conservation equations employed in many depth-integrated models of Earth-surface flows that entrain bed sediment (Appendix B). The model of Cao et al. [2004], for example, is commonly used to simulate the effects of bed sediment entrainment by dam-break floods, but it employs conservation equations with forms quite different from those presented here. To reconcile these disparate forms, we first note that equation (1) of Cao et al. [2004] is essentially the same as our mass conservation equation (37). Furthermore, manipulation of mass conservation equations (3) and (5) of Cao et al. [2004] using our notation shows that the equations imply that the evolving sediment concentration c 1 in a flow that entrains bed sediment with constant concentration c 2 (i.e., with bed porosity 1 À c 2 ) obeys where d=dt denotes the depth-averaged material time derivative defined in Tables 3 and 4 and (E À D) Cao is the net volumetric sediment influx per unit area due to bed erosion. Note that (E À D) Cao differs from our quantity E, which describes the entrainment rate of a solid-fluid mixture, rather than that of sediment alone. Moreover, Cao et al. [2004] use D to explicitly represent sediment deposition rates, whereas we use E < 0 to implicitly indicate net deposition.
Through use of (47), the last term in momentum conservation equation (6) of Cao et al. [2004], identified by the authors as a "momentum transfer" term, can be recast as Àu 1 h 1 ρ s À ρ w ð Þ =ρ 1 ½ dc 1 =dt À Á . The origin and meaning of this term can be clarified by noting that dρ 1 =dt ¼ ρ s À ρ w ð Þ dc 1 =dt À Á applies and thereby making the identification This result demonstrates that the last term in equation (6) of Cao et al. [2004] is the same as the additional source term listed here in the second row of Table 3. The term accounts for changes in flow bulk density, which are independent of any effects of momentum exchange with the bed.
Another widely used entrainment model is that of McDougall and , whose x momentum conservation equation is cast in a Lagrangian form similar to that shown here in the second row of Table 4. Unlike Cao et al. [2004], McDougall and Hungr [2005] assume that all flow material and entrainable bed material have the same, invariant bulk density (i.e., ρ 1 ¼ ρ 2 ). As a result, in the x momentum equation of McDougall and Hungr [2005], Àρ 1 u 1 E is the only source term analogous to those derived here and listed in the second row of Table 4. This term expresses the apparent negative momentum of entrained, static bed material when it instantaneously enters a moving reference frame that has velocity u 1 . The term Àρ 1 u 1 E is correct in the model of McDougall and Hungr [2005] only if β = 1, implying that the basal slip velocity is the same as the depth-integrated flow velocity. By contrast, use of an Eulerian frame to model the entrainment process obviates the need for this assumption and reveals the importance of differences between depth-averaged and basal velocities.
Finally, it is worthwhile to note that the models used by Cao et al. [2004], McDougall and Hungr [2005], and in the other entrainment studies we have examined make the implicit assumption of a nearly planar erosional interface that satisfies (∂z b /∂x) 2 < < 1 and (∂z b /∂y) 2 < < 1. Under this assumption, the jump conditions listed in the first two rows of Table 2 and in the conservation equations (37)-(39) apply exactly. For irregular erosional interfaces in which (∂z b /∂x) 2 < < 1 and (∂z b /∂y) 2 < < 1 do not apply, the jump conditions have the forms listed in the final row of Table 2, and their implications for conservation equations are more complicated. Despite these complications, (37)-(39) demonstrate that a clear distinction must be made between the boundary shear traction τ 1 zx bot acting at the base of a flow crossing an effectively rigid bed and the smaller boundary shear traction τ 2 zx top engaged during entrainment.

Summary and Outlook
Depth-integrated conservation equations are widely used to model the behavior of Earth-surface mass flows that entrain bed material. A two-layer depth-integrated analysis that accounts for mass and momentum conservation as material moves across an evolving flow-bed boundary yields three jump conditions that constrain the correct forms of single-layer depth-integrated equations. Many published single-layer equations are inconsistent with these jump conditions. Sources of disparities in published equations are diverse, but many stem from imprecise accounting for the relationships between boundary tractions and flow-bed momentum exchange during entrainment.
Identification of mass and momentum jump conditions applicable at an evolving boundary also enables any continuous medium to be represented as a stack of layers that interact at multiple boundaries. As the number of layers increases, the mathematical description provided by a set of depth-integrated conservation laws and boundary jump conditions can be used to approximate a full, three-dimensional description of the physical system. The key step in this procedure is identification of all conditions that must be satisfied at an individual, evolving boundary.
Several alternative forms of the single-layer depth-integrated conservation equations are each correct mathematically, but some are better suited than others for clarifying the physics represented by the mathematics. In our view the first lines of (37)-(39) provide the clearest exposition of the underlying physics. The left-hand sides of these equations accurately represent evolution of conserved quantities (mass and momentum), and the right-hand size sides contain only simple source terms and boundary tractions that are engaged during entrainment. Source terms that are more complicated and more difficult to interpret arise if different boundary tractions are used or if the left-hand sides of the equations are expressed in nonconservative forms such as those shown in Tables 3 and 4. Use of a nonconservative form and Lagrangian reference frame can be particularly problematic because it requires the use of a frame-dependent source term to account for the relative negative momentum of static bed sediment that is entrained. In this context the distinction between relative and total momentum is crucial. In any valid model, regardless of its frame of reference, basic physics dictates that the inertia of entrained static bed material can affect the velocity but not the total momentum of a flow that entrains it.
Any depth-integrated model that considers entrainment of bed material must also consider the jump in boundary shear traction that necessarily exists at the interface between moving and static material (i.e., τ 1bot À τ 2top > 0). Such a jump implies that if all other factors are constant, then an actively eroding bed exerts less shear resistance than that exerted in the absence of erosion. Thus, to correctly represent evolution of flow momentum during bed material entrainment, models must clearly distinguish between resisting shear tractions engaged in the presence of static and eroding beds. Otherwise models can confuse the effects of constitutive assumptions (which influence calculated boundary shear tractions) with those of momentum conservation. In depth-integrated models, decreased basal shear resistance can be the only source of flow momentum gains that sometimes accompany entrainment [e.g., Iverson et al., 2011].
The jump conditions used to obtain (37)-(39) invoke no physics beyond that which is contained in mass and momentum conservation laws themselves. Indeed, the jump conditions simply embed into single-layer equations the constraints that two-layer conservation laws impose. In this way, the jump conditions can enhance the scope and rigor of single-layer models without adding the computational burdens of a two-layer formulation.
To our knowledge no previous study has identified a z momentum jump condition that links boundary normal tractions and bed sediment dilatancy during entrainment. We have shown that this condition is

10.1002/2013RG000447
inextricably coupled to a mass jump condition that must be satisfied if moving material and static bed material have differing bulk densities. The z momentum jump can have particularly important ramifications when viewed in the context of bed erosion formulas, owing to its possible association with pore pressure feedback accompanying dilative or contractive motion of water-saturated sediment.
The three jump conditions summarized in Table 2 indicate that erosion formulas for planar beds with ∂z b /∂x = ∂z b /∂y = 0 must be consistent with the form of equation (40). If a change in bulk density and z momentum jump during entrainment are absent, then (40) reduces to a relatively simple form derived by Fraccarollo and Capart [2002]. On the other hand, some commonly used erosion formulas are difficult to reconcile with (40), regardless of whether changes in bulk density occur during entrainment (Table 5). To provide testable erosion rate predictions, (40) must be accompanied by constitutive information that describes shear and normal tractions at eroding interfaces. High-resolution experiments with idealized, well-sorted sediment beds have identified some relevant constitutive information [Capart and Fraccarollo, 2011], but similarly detailed experiments with poorly sorted natural debris are needed. A critical issue concerns identification of the magnitude and location of basal slip in complicated sediment beds that contain natural grains with diverse shapes and sizes.
In our opinion, further progress in using depth-integrated models to predict the behavior of evolving, erosive Earth-surface flows will hinge on meeting four objectives. First, consistent use of rigorously derived conservation equations that embed the effects of basal jump conditions will minimize the potential for confusion about the roles of momentum exchange and boundary tractions during bed sediment entrainment. Second, use of high-resolution, shock-capturing numerical methods to solve the conservation equations will minimize the potential for spurious computational results. Third, and perhaps most challenging, refinement of erosion rate formulas that satisfy the constraints imposed by basal jump conditions will ensure that predictions are well founded mechanistically. Finally, and perhaps most importantly, predictions must be subjected to unequivocal tests that employ data from controlled experiments with geophysically relevant materials. Such work is in its earliest stages.

Appendix A: Effects of Vertical Accelerations in Depth-Integrated Models
Many models of dam-break floods, debris flows, and rock avalanches use depth-integrated conservation equations formulated in terms of Earth-centered coordinates with z vertical and x horizontal, whereas our equations employ rotated coordinates with z normal to the spatially averaged basal slope and x parallel to the slope (Figures 4 and 6). This difference in coordinate systems can complicate comparisons, but it also helps reveal a fundamental difficulty that can arise when using an Earth-centered formulation to model flows on steep slopes.
To appreciate the source of the problem, it is useful first to consider a scenario in which both formulations give correct results: hydrostatic equilibrium of a liquid of uniform density ρ in a pond with sloping bottom topography ( Figures A1a and A1b). In a typical formulation with z vertical [e.g., Cao et al., 2004], the local x direction static force balance in the pond is expressed by the equation ρgh[∂(h + z b )/∂x] = ρgh[tan θ À tan θ] = 0, where θ is the local angle of bed inclination. In our rotated coordinates the same balance is expressed by the equation ρghsin θ À ρghcos θ[∂(h + z b )/∂x] = ρgh[sin θ À cos θ tan θ] = 0. Thus, although the quantities h and z b have different meanings in the two formulations, the formulations give equivalent results.
Now consider a scenario in which the two formulations give differing results: a tabular mass of uniform density ρ and thickness h accelerates as it descends a uniform frictionless slope inclined at the angle θ ( Figures A1c and  A1d). For this scenario the source terms in the x momentum equations are the same as in the static balances for the ponds. The formulation with z vertical describes evolution of the horizontal velocity, u hor , given by By contrast, the rotated-coordinate formulation describes evolution of the slope-parallel velocity, u par , which obeys the downslope momentum equation Equation (A1) is not equivalent to (A2), because (A2) implies that the horizontal component of u par obeys ρh(du hor /dt) = ρghsin θ cos θ. The difference between this result and (A1) arises from the fact that (A1) disregards the existence of a vertical acceleration component that obeys ρh (du ver /dt) = À ρghsin 2 θ and necessarily exists when a frictionless mass descends a slope.
Corrections to (A1) (and to classical shallow-flow equations) to account for the effects of different coordinate systems can be considered in a general way [e.g., Keller, 2003], but in the context of the simple problem illustrated here, the value of g can be adjusted to account for the effects of vertical acceleration, such that in (A1) g is changed to g(1 À sin 2 θ).
The consequent reduction of g is 25% for slopes of 30°, which are routinely encountered by debris flows and avalanches and occasionally by water floods. Without this adjustment, the error in equation (A1) is substantial.

Appendix B: Mathematical Summary of Some Key Models
See Table B1 on the following page.

Appendix C: Derivation of Kinematic Boundary Conditions
In a fixed Cartesian reference frame, consider a bed or flow layer with a bounding interface that has a geometry that can evolve in response to both motion of material along the interface and passage of material through the interface. The interface position can be described mathematically by F x; y; z; t ð Þ¼z À z b x; y; t ð Þ¼0; where z = 0 defines a reference plane, x and y are coordinates normal to z, and z b is the height of the interface above z = 0 ( Figure 6). The function F is defined everywhere, but it satisfies (C1) only on the interface. Vectors defined by ∇F are necessarily normal to surfaces of constant F, and unit vectors normal to F are defined by where |∇F| denotes the magnitude of ∇F. Together, (C1) and (C2) imply that positive unit vectors ( → n > 0) are directed outward from the interface of interest, which in this case forms the top boundary of a layer that occupies a region where z < z b .
The instantaneous velocity of the interface → i depends on the velocity of material at the interface → u and on the rate at which material passes outward through the interface during erosion (or inward during deposition), E → n. Here E is the outward volumetric flux of material per unit area normal to the interface. Therefore, the interface velocity can be expressed as where → i, → u, E and → n are each functions of x, y, and t on the interface, where F = 0 and z = z b (x, y, t).
A differential equation describing evolution of the interface can be obtained by noting that the interface definition F = 0 implies that the time derivative of F must vanish in a reference frame Figure A1. Diagrams illustrating relationships between bed elevation z b and layer thickness h evaluated using coordinate systems with z vertical and z rotated at the angle of the local bed slope, θ. (a and b) A statically ponded layer. (c and d) A tabular layer descending a frictionless slope.

Reviews of Geophysics
Evolution of flow bulk density is considered, but attendant jump in z momentum is neglected. Use of horizontal x coordinate omits effects of vertical acceleration. Uniform velocity profile normal to x is assumed.

McDougall and Hungr [2005]
x direction parallel to bed owing to use of Lagrangian reference frame and assumptions of constant bulk density, uniform velocity profile normal to x direction, and no z momentum.

Pitman et al. [2003]
∂ h x direction parallel to bed in direction of maximum local slope = u r j j ð Þ tanδ; All variables are nondimensional. The physical variables, which are originally expressed in Cartesian coordinates, are converted into local quantities based on varying topography. The momentum equation is consistent with that derived here.
moving with the interface velocity → i. Expressing this time derivative as D=Dt ¼ ∂=∂t þ → iÁ ∇, we stipulate that applies at F = 0. This equation is similar to a standard kinematic boundary condition commonly used in fluid mechanics, except that the relevant velocity is the interface velocity → i rather than the material velocity → u. By using (C3), (C4) can be expressed in terms of → u as A more useful from of (C5) results from employing (C2) and the vector identify ∇F Á ∇F = |∇F| 2 to find that → nÁ∇F ¼ ∇F j j. Substitution of this relationship in (C5) reduces the equation to This equation describes three-dimensional motion of an evolving interface at F = 0 with a flux of material E passing normally through it, and it differs from a standard kinematic boundary condition owing to inclusion of the interface flux term, À E|∇F|.
A more specific form of (C6) that employs only Cartesian scalar quantities is needed for use in depth-integrated models. This form can be obtained by first differentiating (C1) to find that the x, y, and z components of ∇F on the interface at F = 0 are given by ∂F=∂x ¼ À∂z b =∂x; ∂F=∂y ¼ À∂z b =∂y; ∂F=∂z ¼ 1 ; and also to find that the time derivative of F is given by Next, the magnitude of ∇F at the interface can be expressed by using the magnitudes of the vector components in (C7) and the Pythagorean theorem to find that Substitution of (C7), (C8), (C9) in (C6) then reduces the kinematic boundary condition to where u(z b ), v(z b ), and w(z b ) are the x, y, and z components of the material velocity → u along the interface. Equation (C10) it is the kinematic boundary condition we express as equations (6) and (7) of the main text.
The terms involving (∂z b /∂x) 2 and (∂z b /∂y) 2 in (C10) account for the fact that the material volume flux directed normally across the interface, E, can have components in the x and y directions, and not merely in the z direction. These terms can be interpreted further by considering the relationship between → n and the x, y, and z direction unit vectors, → e x , → e y , and → e z at the interface where F = 0 [cf. Gray, 2001]. Using these unit vectors in conjunction with (C7), the vector ∇F at the interface can be expressed as The magnitude of ∇F at the interface can also be expressed in terms on these unit vectors: This equation simply recapitulates (C9) because → e x ¼ → e y ¼ → e z ¼ 1. Thus, substitution of (C11) and (C12) in (C2) yields This result shows how the terms in brackets in (C10) account for differences in the orientation of the unit normal vector → n and the orientations of the unit vectors → e x , → e y , and → e z .

Appendix D: Evaluation of Stress Components Affecting x Momentum
Explicit evaluation of the stress components in (18) reveals their potential influences on bed sediment entrainment and also facilitates comparisons of our x momentum conservation equation with equations used in other models. Evaluation of stress components requires constitutive assumptions, and we restrict the scope of our evaluation by assuming that stresses within any individual layer are gravitationally induced and are linearly dependent on depth H À z, where z = H denotes the top of the multilayer system. This assumption is routinely used in geomechanics, but it does not necessarily require that z direction forces are statically balanced (i.e., "lithostatic"). However, it implies that for layer 1, the stress components σ xx and τ yx in (18) can be expressed as where κ 1N and κ 1S are proportionality coefficients and g À z = À g z is used to clarify the fact that σ xx and τ yx are positive (because z is positive upward and, consequently, g À z > 0). If z direction forces are unbalanced, upward or downward acceleration modifies the influence of g À z , and in (D1) and (D2) this influence is absorbed into κ 1N and κ 1S . The simplest special case of (D1) and (D2) is the hydrostatic case assumed in many shallow-water models. In this case κ 1N = 1 and κ 1S = 0 apply. In layer 2, σ xx and τ yx depend on tractions imposed by the base of layer 1, which has thickness h 1 , and on the weight of material within layer 2. Consequently, the stress components in layer 2 can be expressed as σ 2xx ¼ κ 1N ρ 1 g Àz h 1 þ κ 2N ρ 2 g Àz H À h 1 À z ð Þ ¼ κ 1N ρ 1 À κ 2N ρ 2 ð Þ g Àz h 1 þ κ 2N ρ 2 g Àz H À z ð Þ ; (D3) τ 2yx ¼ κ 1S ρ 1 g Àz h 1 þ κ 2S ρ 2 g Àz H À h 1 À z ð Þ ¼ κ 1S ρ 1 À κ 2S ρ 2 ð Þ g Àz h 1 þ κ 2S ρ 2 g Àz H À z ð Þ : The second lines of (D3) and (D4) show the stresses in layer 2 involve a term proportional to H À z, analogous to (D1) and (D2), but additionally involve a term proportional to the differences κ 1N ρ 1 À κ 2N ρ 2 and κ 1S ρ 1 À κ 2S ρ 2 at the interface between layers 1 and 2. Nonzero values of these differences imply that jumps in σ xx and τ yx exist at the interface. If no jump exists, then (D3) and (D4) reduce to forms like (D1) and (D2).
For the special case without a jump in σ xx and/or τ yx , we evaluate the effects of stress components by substituting (D1) and (D2) into (18), and then generalizing the algebraic procedure of Iverson [2012, p. 15] to account for nonzero gradients of ρ. The resulting equation may be written as where L x ¼ κ N g À z h ρ ∂H ∂x þ H À z mid ð Þ ∂ρ ∂x ! and (D6a) summarize the effects of stress gradients in the x and y directions, respectively, and z mid = (z top + z bot )/2 denotes the height of the midpoint of a layer bounded by the surfaces at z top and z bot . Subscripts denoting layer 1 or layer 2 are omitted from (D5) and (D6a) and (D6b) because the equations apply equally to either layer, provided that κ 1N ρ 1 À κ 2N ρ 2 ¼ κ 1S ρ 1 À κ 2S ρ 2 ¼ 0 applies. The terms H À z mid ð Þ∂ρ=∂x ð Þand H À z mid ð Þ ∂ρ=∂y ð Þin (D6a) and (D6b) account for the effects of bulk density gradients on stress gradients, and for a single-layer model with z bot = 0, (D6a) and (D6b) reduce to If κ N = 1 and κ S = 0 apply, then (D7a) and (D7b) can be manipulated to obtain the concentration gradient source term in momentum equation (6) in the model of Cao et al. [2004]. This term arises irrespective of whether bulk density gradients are caused by entrainment.
For more complicated cases that involve jumps in σ xx and τ yx at the interface between layers 1 and 2, additional terms must be included in (D5) to account for derivatives of κ 1N ρ 1 À κ 2N ρ 2 ð Þ g Àz h 1 and κ 1S ρ 1 À κ 2S ρ 2 ð Þ g Àz h 1 that arise when (D3) and (D4) are substituted into (18). These additional terms are relevant only for layer 2, and after some algebraic manipulation they can be reduced to the form Use of these expressions generalizes (D5) to ∫ ztop zbot Σ F x dz ¼ ρ 2 g x h 2 À L 2x þ L 2y À J 2x þ J 2y þ τ 2zx top À τ 2zx bot ; where subscripts are used to distinguish quantities in layers 1 and 2. Physical interpretation of J 2x and J 2y is straightforward if κ values are the same in each layer. In this case (D8a) reduces to J 2x ¼ g Àz h 2 κ N ∂ ρ 1 À ρ 2 ð Þ h 1 ½ =∂x, demonstrating that longitudinal variations in ρ 1 À ρ 2 associated with the mass jump condition (14) can produce a longitudinal thrust that affects the net x direction force acting on layer 2. This force, in turn, could influence entrainment by modifying the x momentum jump condition (26). However, the effect of this force is higher order (i.e., of order ε) than are the effects already represented in (38), and it is likely to be negligible in many cases. Similar higher-order effects are neglected in the z momentum equation (39) and jump condition (33).