A Depth‐Averaged Description of Submarine Avalanche Flows and Induced Surface Waves

This paper develops a depth‐averaged theory to investigate submarine landslides and resulting water waves. The problems here consist of a pure fluid regime and a mixture regime of grains and fluid. Both regimes separate from one another by an interface, which is a material surface for grains. While the downslope velocities of the both phases are assumed to be identical in the mixture regime, the velocity shear causes a rearrangement of grains, which induces a vertical relative motion between the phases. The established theory consists of five coupled conservation equations, which describe the evolution of the pure fluid thickness, the mixture thickness, the solids volume fraction, and depth‐averaged velocities. To handle nonconservative products emerging in the equations, a new coordinate system is introduced to rewrite the equation system in an equivalent form, so that numerical solutions are insensitive to the choice of discretization of nonconservative products, which enables us to accurately characterize the dynamic behaviors of particles in the collapse experiments of underwater particles and describe free‐surface wave profiles. It is shown that the computed results are in good agreement with the experiments reported in previous literatures.

distance. These poorly understood phenomena are often observed in nature and laboratory settings (see Bryn et al., 2005;Locat & Lee, 2002;Viroulet et al., 2013). This paper focuses on the dense layer of submarine landslides and generated water waves, as shown in Figure 1. This problem has became a scientific hot subject due to the catastrophic hazard posed by tsunami waves generated by submarine landslides (see Lange et al., 2020;Rauter et al., 2022;Rzadkiewicz et al., 1997;Si et al., 2018). Early studies usually treated dense-layer particles as a rigid box sliding on an incline to simulate landslide (see Grilli & Watts, 2005;Walder et al., 2003). Given that the rigid box model is incapable of accounting for the shear-induced granular deformation and the fluid-solid interaction, deformable materials such as fine gravels were used to mimic the motion of submarine landslides in the experiments (Mohammed and Fritz (2012) and Takabatake et al. (2020)). Ma et al. (2013) numerically simulated submarine landslide as a dense water-sediment mixture, whose motion is driven by the baroclinic pressure forcing due to spatial density variation. E. Fernández-Nieto et al. (2014) considered landslide and ambient water as a two-layer medium. The deformable granular materials behave like a rigid body as the shear stress acting on them is below a critical stress and above that like a viscoplastic fluid. In order to model the complex behavior of granular materials, numerical analyst employed Bingham fluids (Bingham, 1922), Herschel-Bulkley rheology (Coussot, 1994), Bagnold fluids (Bagnold, 1954), and Coulomb plasticity and μ(I)-rheology (Jop et al., 2006), respectively, to describe constitutive response of the granular flows to shear.
Large-flume experiments conducted by the U.S. Geological Survey (USGS), see Iverson et al. (2000), showed that loosely packed wet granular material accelerates rapidly when the pore fluid pressure rises to a threshold value, whereas densely packed wet granular material deforms slowly and moves intermittently. Experiments of underwater granular flows, see for example, Pailha et al. (2008), Pailha and Pouliquen (2009), Rondon et al. (2011), and C. , have also evidenced certain fast-moving and some slow-moving underwater granular flows that differ only in compactness. This implies that in addition to a reliable constitutive description to characterize the response of the granular material to shear, additional physics is also needed in order to rigorously model the motion of submarine landslide and the induced water waves.
Progress in understanding granular flows shows that the combination of the granular dilatancy law and mixture theory (or two-phase averaging theory) is a promising approach to account for submarine landslides, which show distinct behaviors in response to a small change in the initial solids volume fraction. In this context, recent representative works include sediment-fluid mixture theory model of Baumgarten and Kamrin (2019), Rauter's (2021) OpenFOAM implementation of compressible μ(J)-and ϕ(J)-rheology of Boyer et al. (2011), and two-phase flow model of Montellà et al. (2021). These two-phase models based on mixture theory incorporate both the granular dilatancy and the viscoplastic constitutive relation, and hence they might be subject to tremendous computational burden when applied to geophysical submarine granular flows and generated water waves. This provides The geometries of submarine granular flows are analogous to subaerial debris flows for which depth-averaged theories are usually developed. Hence, this paper aims to develop a depth-averaged theory that eliminates the dependence of the normal coordinate from the original mixture theory equations to describe underwater granular flows and generated waves. Depth-averaged theories make full use of the flow characteristic that the typical flow thickness H is much smaller than the typical flow length L, which allows an asymptotic expansion to be performed with respect to aspect ratio ɛ = H/L ≪ 1 for mass and momentum conservation equations. Depth-averaged theories have a long tradition to be applied to flood flows (see Abril & Knight, 2004;Mignot et al., 2006), snow avalanche flows (see Gray & Edwards, 2014;Savage & Hutter, 1989), debris flows (see Iverson, 1997), and submarine landslides (see Bouchut et al., 2016;E. D. Fernández-Nieto et al., 2008), as a few examples. Lynett and Liu (2002) derived a depth-averaged model to describe the generation and propagation of water waves induced by a submarine landslide, in which the motion of landslide is not explicitly described and instead modeled only as a forcing function. The resulting model is capable of describing wave propagation from relatively deep water to shallow water. Iverson and George (2014) particularly focused on the motion of fluid-saturated sediment, in which the tangential percolation of fluid through the skeleton of grains is not considered but the vertical relative motion between grains and fluid, caused by the granular dilatancy, is modeled. Since the upper pure fluid regime is not taken into account in the model of Iverson and George (2014), their model is more appropriate for the motion of subaerial debris flows. Recently, Bouchut et al. (2016) have combined the averaging theory of Anderson and Jackson (1967) and the dilatancy law of Roux and Radjai (2001) to construct a two-phase two-layer model, which captures distinct transient behaviors between uniform-thickness experiments of loose and dense underwater granular flows (see Pailha & Pouliquen, 2009). The theory of Bouchut et al. (2016) allows for tangential and normal percolation of fluid through the granular matrix, and the fluid exchange between two layers is a function depending on both the tangential and normal velocity differences between grains and fluid. The resulting model is relatively more complex compared to the model of Iverson and George (2014). The system of their model equations consists of seven coupled partial differential equations (PDE) for two-dimensional flow problems, which are challenging to perform numerical simulation. To simulate the experiments presented in Rondon et al. (2011) on the nonuniform collapse of underwater granular flows, Bouchut et al. (2017) assumed that upper pure fluid does not move, so that the mass and momentum balance equations are trivially satisfied in the pure fluid regime, and hence a numerical implementation is achievable in this case. Moreover, Grilli et al. (2017) modeled the upper fluid layer flow using a nonhydrostatic model and the submarine landslide as either a dense Newtonian fluid or a Coulomb friction medium to assess the influence of the rheology of the submarine slide layer on the kinematics of the tsunami waves. More relevant depth-averaged models are documented in Yavari-Ramshe and Ataie-Ashtiani (2016). This paper derives a depth-averaged theory toward the goal to improve the prediction of existing models for submarine landslides and generated water waves. The problems of interest consist of a pure fluid regime and a mixture regime of grains and fluid, separated by an interface, which is a material surface for the granular phase. The fluid can penetrate this interface. Analogous to the model of Iverson and George (2014), the present theory does not explicitly take into account the tangential percolation of fluid through the grain skeleton in the mixture regime, but the vertical relative motion between grains and fluid is characterized using the granular dilatancy law (see Meng et al., 2020;Pailha & Pouliquen, 2009;Roux & Radjai, 1997). In this way, the form of the current depth-averaged equations is manipulable, though nonconservative products emerge. By using a new coordinate system, the derived equations are rewritten into an equivalent form, in which nonconservative products are now proportional to the free-surface fluctuation of fluid, which is usually a small quantity. In this equivalent form, numerical solutions are insensitive to the choice of discretization of nonconservative products. It is the combination of this relatively simple mathematical form and proper numerical treatment that allows the model to provide a prediction that is in better agreement with the experiments of Rondon et al. (2011) for the collapse of underwater particles, compared to the depth-averaging theory of Bouchut et al. (2017). Additionally, the rigorous model of particle-particle interactions in this paper enables a significant improvement in the prediction of generated wave profiles compared with the theory of Ma et al. (2013)

Field Formulation
As shown in Figure 1, submarine landslides involve a pure fluid regime and a mixture regime. The pure fluid regime overlies the mixture consisting of grains and fluid, which is described by mixture theory presented below.

Mixture Theory
Mixture theory of continuum mechanics has a long tradition in describing multiphase interacting media (see de Boer & Ehlers, 1990;Truesdell, 1984). It is assumed that each spatial point is occupied simultaneously by all the phases with individual volume fraction ϕ ν , where ν = g, f for the granular and the fluid phase, respectively. This makes it possible to define overlapping partial density ϱ ν , partial velocity u ν , and partial stress σ ν . The individual mass and momentum conservation equations for each constituent in mixture theory are as follows: where β ν represents the interaction force acting on the ν phase. In mixture theory, the partial density ϱ ν and the partial velocity u ν are related to the intrinsic density ϱ ν⋆ and the intrinsic velocity u ν⋆ as follows: The partial stress tensor for the fluid phase, σ f , is defined as follows: where p f⋆ is the pore fluid pressure and the fluid deviatoric stress, τ f , satisfies Newtonian fluid law and it is in mixture theory (see Nunziato et al., 1986), where η f is the dynamic viscosity and the superscript T is the transpose. The form Equation 5 is a potential source of discrepancy between mixture theory and averaging approaches, see Joseph and Lungdren (1990). The form Equation 5 implicitly assumes that the fluid deviatoric stress is proportional to its local concentration. In mixture theory of Baumgarten and Kamrin (2019) for sediment-fluid flows, the effective viscosity of Einstein's (1906), η f (1 + 5ϕ g /2), is applied. Nevertheless, in the ensemble averaging theory of Jackson (2000) and the phase-averaging theory of Drew (1983) and Ishii and Hibiki (2010), the forms of the effective viscosity are more complex and differ from one another. So far, it is by no means clear how to express this deviatoric stress.
Following de Boer and Ehlers (1990), the partial stress for grains is as follows: where σ e is the solid effective stress defined by Terzaghi (1943), which has opposite sign to that of the Cauchy stress.
The interaction drag usually takes the following form: where the first right-hand side term p f⋆ ∇ϕ g in combination with −∇(ϕ g p f⋆ ) results in −ϕ g ∇p f⋆ , which is usually called the buoyancy force in classical fluid mechanics. The second right-hand side term is the Darcy interaction drag acting on particles, where C d is the Darcy drag coefficient and d is the grain diameter, see Pailha and Pouliquen (2009).
The bulk density, ϱ m , and the bulk velocity, u m , are defined as follows: for the center of mass. The bulk stress, σ m , is as follows: Summation of the grain and fluid mass conservation law (Equation 1) gives the mass conservation law for the bulk. The combination of the grain and fluid momentum balance equations to eliminate the Darcy interaction drag gives the momentum balance equation for the bulk: where σ′ is Geophysical statistical data imply that the percolation speed of the fluid through the grain skeleton is relatively small compared to the granular convective velocity, and hence u g ≈ u m . Iverson and his collaborators have frequently employed this postulation to create a model with a relatively simple mathematical form; for details, see Iverson and Denlinger (2001) and Iverson and George (2014). This paper will use this approximation. Volume fractions sum to unity, that is, ϕ g + ϕ f = 1. The momentum Equation 12 in this regime reduces to

Single Phase Regime
Equations 11 and 14 from the two-phase mixture theory are capable of reducing to those describing a single phase regime. In the pure fluid regime, the solids volume fraction vanishes and the fluid volume fraction is unity, that is, ϕ g = 0 and ϕ f = 1. The partial density and the bulk density are therefore identical to the fluid intrinsic density, ϱ f = ϱ m = ϱ f⋆ , and u m = u f in this regime. The solid effective stress vanishes in this regime, and hence Equations 11 and 14 naturally reduce to the incompressible Navier-Stokes equations in this case, In the pure particle regime, ϕ f = 0 and the pore fluid pressure vanishes. In this case, Equations 11 and 14 from the mixture theory reduce to those employed to model granular flows by Savage and Hutter (1989). The pure particle regime does not appear in this work, and hence the forms of the governing equations will not be shown here.

Dilatancy Model
When particles are subject to shear, particles will rearrange themselves in respond to the shear, leading to a change in the solids volume fraction. This behavior was found by Reynolds (1986) and termed dilatancy thereinafter. Recent experiments, see Iverson et al. (2000) and Pailha et al. (2008), have shown that the granular dilatancy has a crucial influence on the generation of the excess pore fluid pressure, which in turn mitigates or enhances the friction of particles. The dilatancy effect is incorporated by introducing a dilatancy angle ψ to characterize the change of volume in a granular material sheared at a shear rate and describes the evolution of the solids volume fraction in terms of shear rate and dilatancy angle, where d/dt is the substantial derivative. For details of the derivation of Equation 17, see Roux and Radjai (1997) and Roux and Radjai (2001).
With mass balance (Equation 1), the substantial derivative on the left-hand side of Equation 17 is eliminated, yielding the following equation: which defines the sign of dilatancy angle. When the bulk volume increases under shear, the dilatancy angle ψ > 0. In contrast, the dilatancy angle ψ < 0 as the bulk volume shrinks under shear. Combination of the grain and fluid mass balance equations implies that shear-induced granular dilatancy is also related to the relative motion between grains and fluid through the following relation: The dilatancy angle ψ governs the rate of Reynolds dilation and a linear expression, given in Pailha and Pouliquen (2009), is employed here, where K 3 is a nondimensional parameter and ϕ eq is the equilibrium packing fraction corresponding to steady-state shearing (see Baumgarten & Kamrin, 2019;Bouchut et al., 2016). In viscous regime, ϕ eq is a decreasing function of viscous number =̇∕ , where p p is the pressure acting on grains (see Boyer et al., 2011). In the transient regime toward the inertial regime, however, recent studies, see Trulsson et al. (2012) and Amarsid et al. (2017), show ϕ eq depends on both the viscous number J and the inertial number =̇∕ √ ∕ . Accurately modeling equilibrium packing fraction ϕ eq is still an ongoing subject in the field of granular flows, and it has not reached an agreement how to express ϕ eq . Pailha and Pouliquen (2009) shows that a linear decreasing function of the viscous number is capable of capturing leading-order behavior of wet uniform-thickness granular flows in low shear rates, which lends weight to the use of Equation 21 in this paper. The parameter ϕ c determines if initial packing is dense or loose. K 2 is a parameter that controls compactness of particles under steady shear, and its value will be calibrated below.

Boundary Conditions
A terrain-following curvilinear coordinate system oxz is used with the x-axis following the change of the topography and the z-axis pointing to normal direction of the basal surface, as shown in Figure 1. For simplicity, we assume the flow is uniform in the lateral direction, but it is straightforward to add the lateral coordinate onto the model if needed. The current problem involves three surfaces: (a) the free-surface , which separates the pure fluid regime from the mixture regime of grains and fluid, and (c) the basal surface F b (x, z) = z − b(x) = 0. In each surface, the upward pointing unit normal is On the free-surface, kinematic boundary condition is as follows: The interface between layers is a material surface for grains is used to express the conservation of the fluid mass flux across this surface. To simplify the derivation, we introduce the notation M to represent the fluid mass transfer, M = ϱ f⋆ (u f+ − u g ) ⋅n m , where the superscript "+" denotes the quantity on the upper side of the interface.
The nonpenetration condition for the bulk consisting of grains and fluid on the basal surface is as follows: Traction-free condition holds on the free-surface On the interface, the upper pure fluid experiences tangential friction from the underlying mixture consisting of grains and fluid. A pragmatic approach is to assume the tangential friction acting on the upper pure fluid is proportional to the difference of velocities between layers, that is, which has been proposed by Beavers and Joseph (1967), who designed an experiment to investigate the boundary condition corresponding to a fluid flow above a porous medium, where the superscript "-" represents quantities on the lower side of the interface and the value of the parameter ϑ depends on permeability of the grain skeleton. For simplicity, this paper assumes ϑ = 0, since the interlayer friction is not important for the problems of interest here, as we shall see later.
Recent progress in suspension rheology, for example, Boyer et al. (2011), shows that the friction acting on the bulk consisting of grains and fluid takes a form similar to the dry granular μ(I)-rheology. The bulk shear stress at the basal surface is equal to the particle pressure multiplying by a rate-dependent coefficient μ b , that is, where the friction coefficient μ b carries rheology information and its expression will be discussed in Section 2.8.
On the interface, the momentum-jump condition between layers is as follows: see Hutter (2008), where the stress 〚σn m 〛 = σ f n m − σ m n m . Using the notation of mass transfer M, the momentum-jump condition (29) is rewritten as follows:

Two-Layer Model Framework
Equations 11 and 14, alongside with the aforementioned boundary conditions, constitute a two-phase system, which could be numerically solved as long as the solid effective stress σ e is prescribed. Accurately modeling solid effective stress, especially in the presence of the viscous liquid, is still a not well solved scientific problem (Guazzelli & Pouliquen, 2018). Alternatively, flow geometries where typical flow thickness H is much smaller than typical flow length L referred to as shallow flow assumption imply that the integration of Equations 11 and 14 through the mixture thickness and Equations 15 and 16 through the thickness of pure fluid allows to derive a set of tractable thickness-averaged equations. The shallow flow assumption made here requires a slowly varying basal topography (Chiou et al., 2005;Gray et al., 1999), that is, the characteristic radius of curvature of the basal topography R is much larger than the typical flow length L. It should be noted that it is not common in the community of granular flows to refer to such models that are averaged over the flow thickness, that is, normal to the basal topography, as thickness-averaged models but as depth-averaged models. For this reason, we refer to such an averaged model as a depth-averaged model also in this paper, 8 of 32 even though the averaging is performed over the flow thickness. The system of depth-averaged equations delineate the thickness h pf of the pure fluid regime, the thickness h m of the mixture of grains and fluid, the depth-averaged tangential velocity of the pure fluid regime, and the depth-averaged tangential velocity of the mixture of grains and fluid as well as the solids volume fraction ϕ g that is assumed to be constant in the normal direction of the topography within the mixture regime but varies with downslope coordinate and time. The thicknesses h pf and h m are respectively, and the depth-averaged tangential velocities are defined as follows: The model equations, derived by integrating conservation equations in the normal direction of the topography, corresponding to the regimes of the pure fluid and the mixture of grains and fluid, for mass are as follows: and for momenta are as follows: The derivation processes of Equations 33-37 are shown in Appendices (A-C) for the sake of brevity. The key assumption made to derive Equations 33-37 is that the solids volume fraction ϕ g is independent of the z-coordinate within the mixture regime. The shape factor χ ν , ν = pf, m for the pure fluid and the mixture of grains and fluid, respectively, in the convective terms of the momentum balance equations characterizes the ratio of the depth-averaged square of the velocity, ( ) 2 , to the depth-averaged velocity squared, (̄) 2 . For simplicity, it is postulated that the shape factor χ ν = 1, which is a common assumption in depth-averaged models (see Bouchut et al., 2016;Gray & Edwards, 2014;Meng et al., 2022;Pitman & Le, 2005). Equation 33 delineates the mass balance corresponding to the pure fluid regime, in which the right-hand side term represents the fluid mass transfer due to change of voids between grains under shear. As grains under shear compact, an amount of interstitial fluid is expelled from voids into the pure fluid regime, thereby resulting in a positive mass product M on the right-hand side of Equation 33 and a negative fluid mass product −M on the right-hand side of Equation 34, which describes mass balance of the bulk consisting of grains and fluid. The right-hand side terms of the momentum balances, Equations 36 and 37, are as follows: where the term Mu f+ is the interfacial momentum transfer induced by the fluid mass transfer across the interface. The source terms in Equations 38 and 39 require expressions of the fluid mass transfer M, the particle basal pressure , and the depth-averaged fluid viscous stress and , which are presented below.

Fluid Mass Transfer
As particles dilate under shear, the ambient fluid is sucked into the voids, resulting in a fluid mass flux from the pure fluid regime into the mixture regime of grains and fluid. Contrarily, as particles compact under shear, the pore fluid is expelled, resulting in a fluid efflux from the mixture regime. Equation 23 is rewritten as u g ⋅ n m = −(1/|∇F m |)∂F m /∂t. Using the mass-jump condition (24), the fluid mass transfer M is as follows: For shallow flows, the derivative, ∂s m /∂x, is a small quantity and is of the order O(ɛ) with respect to aspect ratio ɛ = H/L ≪ 1. It gives |∇F m | = 1 + O(ɛ 2 ), in which |∇F m | = 1 to leading order in the small parameter ɛ (see Gray & Edwards, 2014;Meng et al., 2022).
The assumption of vanishing difference of the tangential velocity implies that Equation 19 reduces to Integrating Equation 41 from the base to the mixture surface of grain and fluid gives the following equation: As mentioned in Pailha and Pouliquen (2009), rigorous calculation of the integration on the left-hand side of Equation 42 is not possible. An assumption has to be made to proceed. We then follow Pailha and Pouliquen (2009) to assume that the dilatancy at the bottom gives the right order of magnitude of dilatancy inside the mixture regime, which implies the following equation: where ψ b is the basal dilatancy angle of grains. By substituting Equations 42 and 43 into Equation 40, the fluid mass transfer can be expressed as: where the basal dilatancy angle ψ b involves the velocity profile to express the solid basal shear rate . To this end, steady uniform suspension flows, see Pailha and Pouliquen (2009) and Boyer et al. (2011), have already shown that a parabolic velocity profile prevails, and hence the basal shear rate is = 3̄∕ℎ . Moreover, substitution of into Equations 20 and 21 shows the basal dilatancy as follows:

Solid Basal Pressure and Depth-Averaged Viscous Stress
Terzaghi's effective stress principle from soil mechanics expresses that the solid effective stress is equal to the pore fluid pressure extracted from the bulk pressure, that is, where the basal pore fluid pressure, , includes hydrostatic and excess components. The excess pore fluid pressure is induced due to the vertical relative motion between grains and fluid, see Pailha and Pouliquen (2009), Iverson and George (2014), Bouchut et al. (2016), and Meng and Wang (2018). The basal pore fluid pressure is as follows: In depth-averaged models (see Bouchut et al., 2017;Pitman & Le, 2005;Savage & Hutter, 1989), the depth-averaged viscous stress is usually ignored, since it is a small quantity in most flows. Nevertheless, viscous stress plays a crucial role in maintaining well-posedness of the model equations (see Baker et al., 2016), in explaining cut-off frequency of roll wave instability (see Gray & Edwards, 2014), and in properly capturing levee formation (see Rocha et al., 2019) and lateral velocity profiles in channel flows (see Iverson & Denlinger, 2001;Meng & Wang, 2018). The fluid in the problem later is much more viscous than normal water, and hence this paper takes account of fluid viscosity effect.
In the pure fluid regime, the depth-averaged in-plane deviatoric stress is as follows: Inclusion of all the terms on the right-hand side of Equation 48 is very challenging and it would lead to a very complex theory. A pragmatic approach is to assume that the bottom boundary layer does not have an influence on the pure fluid regime, that is, =̄ for s m ≤ z ≤ s f . Thus, Equation 48 becomes Analogously, the depth-averaged fluid deviatoric stress is in the mixture regime of grains and fluid.

Friction Coefficient μ b
Wet granular μ(J)-rheology revealed that the basal friction coefficient μ b is a function of the basal viscous number =̇∕ at low shear rates, see Boyer et al. (2011). A linear relation is employed here. It is originally proposed by Pailha and Pouliquen (2009) and used by Bouchut et al. (2016) to model underwater granular flows. The incorporation of the basal dilatancy angle, ψ b , into the friction coefficient expresses that friction is enhanced, as particles are dilated. Conversely, the friction of particles weakens when particles are compressed. In Equation 51, K 1 is a constant and its value will be discussed in Section 4.1.

Final Model Equations in the Terrain-Following Coordinate System
For the description of numerical method in the following section, the system composed of Equations 33-37 has to be put in a more standard mathematical form. The intrinsic densities, ϱ f⋆ and ϱ g⋆ , are canceled out to further 11 of 32 simplify mass balance Equations 33 and 35. By repeatedly using the grain mass balance Equation 35, one can eliminate the dependency of the bulk mass balance Equation 34 on the bulk density ϱ m . Finally, the system of mass conservation laws becomes Using mass balance Equations 33 and 34 to eliminate the dependence of the left-hand side terms of the momentum Equations 36 and 37 on the fluid and grain intrinsic densities, respectively, one obtains the following equations: where the right-hand side terms S pf⋆ and S m⋆ are as follows: where ν f = η f /ϱ f⋆ is the kinematic viscosity of the fluid. The bulk basal friction  is as follows: In momentum Equations 55 and 56, the gradients of the basal topography are omitted because the cases discussed later do not involve elevation of the basal topography. In Equations 57 and 58, the terms −(h pf g cos ζ)∂h m /∂x, − 1 2 (ℎ ) 2 (cos )∕ ∕ , and −ϱ f⋆ /ϱ m h m ∂(h pf g cos ζ)/∂x represent the pressure coupling between layers (see the "layered interplay" terms in Equations 38 and 39), and they cannot be rewritten in a conserved form so that the system of model equations is nonconservative. The existence of nonconservative products implies that the jump conditions are not uniquely defined without involving additional physics. One has to construct more complex scheme to solve nonconservative PDE (see Abgrall & Karni, 2009;Castro Diaz et al., 2007;Diaz et al., 2019;Kim & LeVeque, 2008). Due to the fact that it is challenging to derive a unique weak solution in the presence of nonconservative products, different choices of discretization for the nonconservative products tend to give qualitatively different computed solutions, and hence no agreement has yet been reached on what a robust scheme for discretizing the nonconservative products should hold. Additionally, it has been assumed that all interfacial quantities are equal to the corresponding depth-averaged quantities. Moreover, the source terms s pf and s m , emerging in Equations 57 and 58, are as follows: 12 of 32

Reformulated Model Equations in Global Coordinate System
In this paper, the proposed two-layer system is numerically solved using a high-resolution central-upwind scheme with a total variation diminishing slope limiter to prevent spurious numerical oscillation (Kurganov & Tadmor, 2000), which has proved its ability to solve similar systems of equations, including the models of dry granular flows (Baker et al., 2016;Edwards & Gray, 2015;Y. Wang et al., 2004) and the models of debris flows (Meng et al., , 2022Tai et al., 2019). Moreover, we introduce a new coordinate system, depicted in Figure 2, to rewrite Equations 52-56 in an equivalent form in which the nonconservative products are small compared to leading-order terms. In this case, the numerical solutions are not sensitive to the choice of numerical methods for discretizing nonconservative terms. Such a coordinate transformation has shown promise in solving similar systems (see Kurganov & Miller, 2014;Kurganov & Petrova, 2009).
In the coordinate system OXZ, the X-axis of the coordinate system points along the horizontal direction and the Z-axis points upward, as shown in Figure 2. The central idea of Kurganov and Miller (2014) lies in the fact that the free-surface fluctuation of fluid ξ(X, t) is usually much smaller than the thickness h pf , where B(X) is the elevation of the basal topography in relation to the horizontal plane, see Figure 2. In the framework of this redefined coordinate, for any given function f(x(X)), the transformation hold.
Then, the new notations are introduced, where w pf , w m , and w mg are the projected thicknesses of the pure fluid, the mixture and the granular phase of mixture to the OZ-axis, q pf , and q m are the momenta of the pure fluid and the mixture along the OX-axis, for which it follows and the final "lake-at-rest" solution is given by the following equation: Using the conservative variables = ( ) T , the system composed of Equations 52-56 can then be rewritten as follows: where γ = ϱ f⋆ /ϱ m . The horizontal gradient of density ϱ m in Equation 72, ∂ϱ m /∂X, is small compared with other terms and is hence neglected in George and Iverson (2014) and the following numerical implementation. By using this transformation, all the nonconservative products emerging in the model equations are now proportional to the free-surface fluctuation of fluid ξ. When numerically solving the model equations, the derived forms Equations 71 and 72 are advantageous, since the free-surface fluctuation of fluid ξ is usually a small quantity, hence the influence of the present nonconservative products is very small compared to the original nonconservative products emerging in Equations 57 and 58. This is remarkable for the numerical implementation because the influence of the nonconservative products is so small that the numerical solutions are insensitive to the choice of discretization of the nonconservative products in the current coordinate system.

Semidiscrete Central-Upwind Scheme
The central-upwind scheme requires that the system of Equations 68-72 is written in vector form as follows: and source terms: The nonconservative products are as follows: the friction source terms: and the viscous terms; In the central-upwind scheme, the computation domain is discretized using uniform cells = [ , where − 1 2 = − Δ ∕2 and + 1 2 = + Δ ∕2 with ΔX = X j+1 − X j . In each cell I j , the cell-averaged conservative variable ( ) , defined as is updated through the semidiscrete ordinary differential equation, where the cell-averaged source, viscous, nonconservative, and friction terms are as follows: The numerical fluxes − 1 2 ( ) and + 1 2 ( ) in Equation 80 maintain nonoscillatory property by using a flux limiter and they are explicitly given in Kurganov et al. (2001) (see (3.15) in Kurganov et al. (2001)).
In this work, the ordinary differential Equation 80 is numerically solved using second-order Runge-Kutta method documented in Shu and Osher (1988). The source terms ( ) , the nonconservative products ( ) , and the viscous terms ( ) were first applied to iterate for an intermediate solution. Then, the friction source terms ( ) were applied to update the solution, see Kurganov and Miller (2014) for details. Moreover, local wave speeds, a 1 , a 2 , …, a 5 , are calculated by solving eigenvalue problems of the Jacobian matrix A = ∂F/∂U, The characteristic Equation 82 is a fifth-order polynomial equation. The nature of the solutions of this polynomial equation determines whether the system of Equation 73 is well-posed or ill-posed. As the velocity difference between the layers becomes large, the system of equations in a multilayer model usually tends to loss its hyperbolic property (Pelanti et al., 2008). This is also the case for the present system Equation 73. However, recent studies show that the introduction of mass-and momentum-exchange across the interface between the layers as well as viscous terms could make the multilayer model exempt from ill-posedness (see Baker et al., 2016;Sarno et al., 2021). This lends weight to conduct a numerical simulation. In this paper, the Lagrange theorem was used to establish bounds on the roots of the polynomial equation and then a real value in this interval was chosen (see Mignotte & Stefanescu, 2002). Finally, since the central-upwind scheme is an explicit scheme, the Courant-Friedrichs-Levy (CFL) condition needs to be satisfied to assure numerical stability, where Cr is the Courant number and is usually chosen less than one. Given the fact that a small time step is needed due to constrains of the source term, in the following numerical simulations the Courant number Cr = 0.1 is chosen. Rondon et al. (2011) designed small-scale experiments of underwater granular collapse, which exhibits many features that are closely akin to what is commonly observed in large-scale experiments and geophysical flows. These experiments showed that small variations in the solids volume fraction influence particle behavior profoundly, which was also evidenced in large-scale experiments (Iverson et al., 2000), and hence they have been widely used to test theories, including the mixture theory of fluid-sediment flow (see Baumgarten & Kamrin, 2019;, the dense suspension μ(J) and ϕ(J) theory (see Boyer et al., 2011;Rauter, 2021) and the depth-averaged theory of submarine granular avalanche flows (see Bouchut et al., 2016Bouchut et al., , 2017. This provides a strong motivation to use the experiments of Rondon et al. (2011) to test the present theory.

Granular Collapse in a Viscous Fluid
In the experiments of Rondon et al. (2011), glass beads of density ϱ g⋆ = 2,500 kg m −3 with an average diameter of d = 225 μm are immersed in a fluid of density ϱ f⋆ ≈ 1,000 kg m −3 with viscosity η f = 12 cP to create a loose packing column (initially more porous than critical) and a dense packing column (initially less porous than critical). The fluid, used in the experiments, is composed of 87% water and 13% Ucon oil, making it much more viscous than tap water at room temperature. These glass beads are initially deposited behind a vertically removable plate and the initial solids volume fraction is controlled between 0.55 and 0.62 in the experiments. The SUN ET AL.

10.1029/2022JF006893
16 of 32 column packed with ϕ g = 0.55 is looser than that in the critical state, while the column packed with ϕ g = 0.62 is denser than that in the critical state. As the plate is removed, the column is allowed to collapse, and the collapse process is recorded using a high-speed camera, thereby capturing run-out profiles. The experimental setup is shown in Figure 3.
The computation domain follows the experimental setup depicted in Figure 3, and the values of physical parameters, required in the computation, strictly follow the experimental measurements. In the experiments, the values of ϕ c , K 1 , K 2 , K 3 , and K 4 were not measured. Nevertheless, the experimental findings in Rondon et al. (2011) show that the transition between dense and loose regimes seems to occur around a critical volume fraction of 0.58. Thus, ϕ c = 0.582 is selected in the computation. In the μ(J) and ϕ(J) suspension rheology (Boyer et al., 2011), ϕ eq = 0.57 when the viscous number is J = 5.4 × 10 −3 , which implies K 2 = 25.0 by relation (21). Pailha and Pouliquen (2009) investigated uniform-thickness flows of glass beads submerged in a viscous fluid, in which the properties of particles and the roughness of the bed are similar to those used in this paper. The calculation of Pailha and Pouliquen (2009) clearly shows the choice of K 1 = 90.5, K 3 = 4.09, and K 4 = 1.8 provided a good prediction for the particle velocity and the basal pore fluid pressure. Thus, the choice of these parameters is maintained in this work, and all the parameters are summarized in Table 1. In particular, none of them is a fitting parameter.

Temporal Evolution of the Thickness Profiles
In this section, an overall impression is first given to highlight that the present theory works remarkably well to predict the dynamic process of both dense and loose columns from release to final deposit, prior to a detailed discussion of the influence of aspect ratio and fluid mass transfer. Experiments with both dense and loose packing are presented, in which the dense column initially has a thickness of h ini = 4.2 cm, a length of l ini = 6 cm and a solids volume fraction of 0 = 0.6 , and the loose column has a thickness of h ini = 4.8 cm, a length of l ini = 6 cm and a solids volume fraction of 0 = 0.55 . Figure 4 shows the comparison between theory and experiments with respect to the temporal evolution of the particle-surface profiles. After the mass is released, the loose column spreads much faster than the dense column, and the run-out distance in the loose case is approximately twice that in the dense case. This behavior has been captured by the present theory, as shown in Figures 4a and 4c, which reveals the significant influence of a small change in the initial solids volume fraction on the dynamic behaviors of the underwater avalanche collapse. Both theory and experiments show that loosely packed column spreads rapidly immediately after the mass is released, and the spreading process lasts less than 3 s. Conversely, the spreading process in the dense column takes more than 10 s and the final deposit exhibits a trapezoidal shape.
It is quite interesting to compare the simulation results with those obtained using mixture theory (Baumgarten & Kamrin, 2019) and depth-averaged theory (Bouchut et al., 2017). Figure 5 shows thickness profiles at the beginning and the end of collapse from the present theory (black solid lines), Baumgarten and Kamrin (2019) (blue dash-dotted lines), Bouchut  3. Experimental setup used in Rondon et al. (2011). A column of glass beads with thickness h ini and length l ini is held in place by a plate. By rapidly moving the plate upward, the column is released and allowed to collapse. Using a high-speed camera, the run-out profiles are recorded (courtesy of Baumgarten and Kamrin (2019) (2019), the current depth-averaged theory works remarkably well in predicting the collapse process of the loose column, though the depth-averaged theories do not take account of the vertical variation of the physical quantities. For the collapse of the dense column, the experiments show that the mass hardly deforms as the plate is removed, and only the particles in the top right-hand corner yield and subsequently move slowly. The present theory does not predict this behavior (see Figure 5a) as well as the full dimensional theory of Baumgarten and Kamrin (2019), mainly because this behavior might be linked to higher order terms of the suspension rheology (see Rauter, 2021) and the equilibrium packing fraction (21), used in this work, only captures leading-order behavior. Nevertheless, the dilatancy law (21), incorporated into the mixture theory, predicts the deposit morphology of the dense case quite well, as evidenced in Figure 5b. Compared with the depth-averaged theory of Bouchut et al. (2017); Figure 5 shows that the present theory provides a better prediction to most cases except for the deposit morphology of the dense column at t = 18 s. The theory of Bouchut et al. (2016), used in Bouchut et al. (2017), consists of seven coupled PDE for two-dimensional flows, which allows the fluid to easily percolate through the grain skeleton downslope. Bouchut et al. (2017) assumed = 0 , so the governing equations reduce to five  tractable PDEs. The postulation = 0 , however, implies that particles are exempt from the drag of the upper pure fluid during their motion, so the particles collapse more rapidly and the run-out distance is farther in Bouchut et al. (2017). Moreover, the temporal evolution of the front position is shown in Figure 6. The simulation results from the present theory show a nontrivial improvement compared to those from Baumgarten and Kamrin (2019) and Bouchut et al. (2017). The front position of the granular collapse predicted by the present theory is closer to the experimental measurements in both the dense and loose cases.

Depth-Averaged Particle Velocity and Fluid Mass Transfer
To better understand the dynamics of the initial dense and loose packing as well as the behavior of dilation or compaction during the flow, the velocity of the depth-averaged velocity of the mixture of grains and fluid and the fluid mass transfer M are investigated. Figures 7a and 7b show longitudinal profiles of of dense and loose cases at different times, respectively. At the beginning of the collapse, the velocity remains at zero for the initially dense packing at t = 3.0 s near the back wall, while for the initial loose packing at t = 0.66 s, almost the entire mass is in motion. The largest velocity appears near the front and decreases with increasing time in both cases. However, as time increases, the motion for the initial dense packing spreads to the interior of the flow, while for the initial loose packing the motion of the mixture of grains and fluid near the rear of the flow has already stopped.
A key feature of the present theory, differing from the depth-averaged dilatancy model of Iverson and George (2014), lies in the fact that in this paper the upper pure fluid is taken into account which percolates downward and upward through the grain skeleton with resuspension and sedimentation of the particles. The vertical relative motion between grains and fluid induces fluid mass transfer across the particle surface in addition to the excess pore fluid pressure. For the dense cases, most particles dilate after the release, with an exception at the flow front where particles compress, as shown in Figure 7c. Thus, at the flow front there is a fluid efflux from the mixture regime into the pure fluid regime, whereas in most parts of the flow there is a net fluid mass flux from the pure fluid regime into the mixture regime. For the loose case, as indicated in Figure 7d, the particles at the flow front move sufficiently fast and the flow thickness at the front is sufficiently thin that the basal shear rate is strong, thereby leading to a dilatant behavior at the flow front at the beginning of the collapse. As a result, for the loose case there is a net fluid mass flux into the mixture regime at the flow front. In most of the flow, however, particles contract for the loose case, thereby resulting in a mass efflux from the mixture regime into the pure fluid regime.

Influence of the Initial Aspect Ratio
Depth-averaged theories are derived based on the assumption of the aspect ratio ɛ = H/L ≪ 1. One of the main aims of this case is to check how far the present depth-averaged theory can be pushed to accurately describe problems in which the characteristic flow thickness H is not sufficiently small compared to the characteristic flow length L. Rondon et al. (2011) have conducted experiments with different initial aspect ratios and investigated how the deposit morphology varies in response to the change of the initial aspect ratio. This provides reliable data to validate to what extent the present theory is capable of predicting underwater granular-flow problems. Figure 8 shows the comparison of the deposit morphology between the theory and the measurement. A quantitatively good agreement with the experiments is found for the simulation results in the prediction of the front position, no matter whether the aspect ratio is relatively small (ɛ = 0.5) or large (ɛ = 2). For the dense cases, the theory accurately predicts a quasi-trapezoidal deposit morphology for the small aspect-ratio case and a shape similar to triangle for the deposit morphology of the high aspect-ratio case, both of which were well observed in experiments. For the loose cases, the simulated deposit profiles are qualitatively consistent with experimental measurements, though some discrepancies exist. In the loose cases, the particles move rapidly immediately after the plate is removed, which is dramatically different from the slow movement in the dense cases. In the present theory, the dilatancy law (21) plays a crucial role in the determination of different response of densely and loosely packed particles to shear, and the linear dependence of the viscous number J does not include second-and high-order terms with respect to J in the limit of a small viscous number J → 0, which is a potential source of discrepancy. There is therefore potential to improve the model (at this point) in future. Nevertheless, it should be noted that the comparison of Figures 8a-8d implies that the present theory is not necessarily restricted to cases with small aspect ratios. It looks promising to apply the present theory to nonshallow flows in the future.

Waves Generated by Underwater Granular Flows
The present theory is also capable of simulating water waves generated by underwater granular flows in addition to the collapse of underwater particles. In this section, the experimental measurements of Rzadkiewicz et al. (1997) on the temporal and spatial evolution of free-surface waves, generated by underwater granular flows, and the corresponding simulation results using the present theory are presented and compared. The experiments of Rzadkiewicz et al. (1997) are among the most widely used ones to validate theories or numerical methods, including the incompressible-smoothed particle hydrodynamics (I-SPH) formulation (see Ataie-Ashtiani & Shobeyri, 2008) and the nonhydrostatic wave model (see Ma et al., 2013). In the experiments of Rzadkiewicz et al. (1997), water waves are generated by a sand mass sliding down a slope. The slope consists of an inclined plane at an angle of 45° to the horizontal, which extends 2.26 m downstream and sharply connects to a horizontal plane. In the experiments, the initial vertical profile of the particles is triangular with two identical sides perpendicular to each other with a length of 0.65 m. Particles are held in this initial position by a vertical water gate, which crosses the water free-surface and is lifted up very quickly to release the mass. The water depth is 1.6 m and the upper surface of the triangular mass is initially 10 cm below the quiescent water surface, see Figure 9a.
The computational domain is set −1 m ⩽ x ⩽ 4 m and is sufficiently large that the particles will not touch the boundaries of the domain. The computational domain is discretized into 500 grid points. On the left-and right-boundaries of the domain, the free outflow condition is set, that is, the gradients of the unknown quantities (flow depth and depth-averaged velocities) vanish, which apparently hold for the particle flow and successfully prevents the reflection of water waves on the boundaries. The experiments provided the value of the bulk density consisting of grains and water ϱ m = 1,950 kg m −3 , instead of the individual densities of particles and fluid. For dense flows like here, ϕ g = 0.58 is a good approximation (see Ma et al., 2013). Since ϱ f⋆ = 1,000 kg m −3 , it follows ϱ g⋆ = 2,638 kg m −3 from ϱ g⋆ = (ϱ m − ϱ f⋆ (1 − ϕ g ))/ϕ g . The remaining parameters follow those presented in Table 1. Figure 9 shows the simulated results of underwater avalanche flows and the generated water waves. In the first instance after the gate is lifted up, the sand mass mainly collapses and slightly spreads out. At t = 0.4 s, the deformation of the sand mass at the tail induces a wave trough, whose shape is very close to the experimental measurement. Figure 10 shows that the predicted free-surface wave profile agrees very well with the experiment overall, but the wave amplitudes predicted by the present theory and the theory of Ma et al. (2013) are larger than the measurement at t = 0.4 s. In the theory of Ma et al. (2013), the shear strength of the particles was not considered. As a result, the deformation of the granular material during the release of the grains is more pronounced than in the experimental flows, and hence the predicted wave amplitudes induced by the deformation of the granular material are accordingly larger. In the present numerical simulation, the particle-water mixture is initially somewhat loosely packed because the initial solid volume fraction ϕ g = 0.58 is slightly smaller than the critical volume fraction ϕ c = 0.582. From the numerical results of the submarine granular column collapse experiment, the initially loosely packed granular-fluid mixture moves faster at the beginning of the collapse than in the experiment (see Figure 5c), while at the end of the collapse, the deformation of the mixture agrees well with the experiment (see Figure 5d). Analogously, in this numerical simulation, at the early phase t = 0.4 s, the motion of the particle-water mixture could be faster than in the experiment, so that the generated water waves also propagate faster than in the experiment. At t = 0.8 s, the predicted shape of the avalanche closely matches what was observed in the experiments (see Figure 7b in Rzadkiewicz et al. (1997)), and hence the discrepancy for the free-surface wave profile between the present simulation and the measurement becomes smaller. This behavior, well observed in experiments, was not accurately predicted by Ma et al. (2013), mainly because of the inaccurate prediction of the avalanche shape in their model. When the particles finally settle in the horizontal plane, the water wave slowly evolves into a quiescent free-surface, see Figure 9d. Since the numerical scheme, used in this paper, possesses a well-balanced property, the present results accurately captures the final "lake-at-rest" state.

Conclusions
In this paper, a two-layer model based on the depth-averaging in the normal direction of the topography is developed to study submarine avalanche flows and generated water waves. The current framework has two regimes consisting of (a) a pure fluid and (b) a mixture regime of fluid and particles. To maintain simplicity and capture key physics, this paper assumes that the tangential percolation of the fluid through the grain skeleton is sufficiently small , allowing a detailed study of the vertical relative motion between grains and fluid that remain. Shear causes a rearrangement of particles, which induces a relative motion between grains and fluid in the normal direction of the topography. On the one hand, the vertical relative motion is accompanied by the fluid mass transfer across an interface which is material surface for grains (Bouchut et al., 2016;Meng & Wang, 2018). Moreover, the vertical relative motion between grains and fluid leads to the generation of an excess basal pore fluid pressure that reduces or consolidates the bed friction of the particles (Iverson, 2005;Pailha et al., 2008). The model equations, derived in this paper, describe the temporal and spatial evolutions of the pure fluid thickness h pw , the grain-fluid mixture thickness h m , the solids volume fraction ϕ g , and the associated depth-averaged tangential velocities and .
The derived equations contain nonconservative products, which usually appear in the multiphase or multilayer models (see Abgrall & Karni, 2009;Bouchut et al., 2016;E. D. Fernández-Nieto et al., 2008;Meng et al., 2017Meng et al., , 2022Pitman & Le, 2005). To minimize the influence of the nonconservative products on the accuracy of numerical solutions, a new coordinate is introduced, which makes possible a favorable handling of the nonconservative products. The core of this coordinate transformation lies in the fact that the free-surface fluctuation of fluid ξ is much smaller than the fluid thickness (Kurganov & Miller, 2014). It turns out that the nonconservative products in this new coordinate system are proportional to the variable ξ. In this case, nonconservative products are insensitive to the choice of discretization of the nonconservative products.
To illustrate the performance of this model and the robustness of the numerical method, a comparison is made to the collapse experiments of underwater particles performed by Rondon et al. (2011). It is shown that the present theory predicts the influence of a small change in the initial solids volume fraction on the dynamic 22 of 32 behavior of the collapse remarkably well. The present theory predicts the temporal evolution of the front position and the deposit morphology as well as the results obtained by the full dimensional model of Baumgarten and Kamrin (2019), without the need of a tremendous computation burden. Moreover, direct comparison of simulation results with experimental measurements for cases whose aspect ratios violate the shallow approximation (ɛ ≪ 1) clearly shows that the present theory is not necessarily confined to shallow flows, which makes promising to apply the present theory to investigate more complex nonshallow flows encountered in industrial and in nature. The present theory is also applied to revisit the experiments of Rzadkiewicz et al. (1997), performed with a mixture of monodisperse sand and water. The theory captures the key qualitative features of the avalanche flows as well as provides a good quantitative match to Rzadkiewicz's free-surface wave profile data without any fitting parameters.
The theory, derived in this paper, provides a framework to incorporate the two-phase mixture theory and the granular dilatancy law, which is an ongoing subject in the community of granular flows. In this paper, the equilibrium packing fraction ϕ eq is a monotonically decreasing linear function of the viscous number J, which captures the leading-order behavior in the limit of a small viscous number J → 0. Higher order terms of the viscous number can be incorporated into ϕ eq by referring to the suspension rheology of Boyer et al. (2011), which works mainly at low shear rates. The simulation data in Trulsson et al. (2012) imply that the equilibrium packing fraction ϕ eq is a combined function of viscous number J and inertial number I in the regime between viscous and inertial regimes. The two-layer framework proposed in this paper is sufficiently general to incorporate those sophisticated dilatancy laws, which is a very important point to improve the present model and will be carried out in the future.
where the hats represent nondimensional variables. The typical flow thickness H is much smaller than the typical flow length L, and hence the aspect ratio ɛ = H/L ≪ 1. Lithostatic balance implies that the scaling of the pore fluid pressure and the grain normal stresses is ϱgH, where ϱ is the characteristic density. Coulomb-type friction implies that the scaling for the grain shear stress is μϱgH with characteristic friction coefficient μ. The curvature of the basal reference surface is defined as κ = −dζ/dx, and hence the scaling of the curvature is 1/R, where R is the characteristic radius of curvature of the basal topography with the assumption L/R ≪ 1. Note that the velocity scaling differs from that of Gray and Edwards (2014) and Meng et al. (2022), which used the slower scale U = (gH) 1/2 .
The nondimensional form of the grain mass balance Equation 1 is as follows: where the nondimensional parameters λ and Ψ associated with the curvilinear coordinate system emerge, defined by the following equation: The nondimensional form of the mass balance equation for the bulk, Equation 11, is as follows: The nondimensional form of the mass balance Equation 15 corresponding to the pure fluid regime is as follows: and the nondimensional downslope and normal components of the momentum balance Equation 16 are as follows:̂(̂) The nondimensional forms of the kinematic boundary conditions on the free-surface (22), the interface (23), and the bottom (25) are as follows: The nondimensional form of the mass jump condition (24) is as follows: with the nondimensional fluid mass transfer as follows: SUN ET AL.

10.1029/2022JF006893
24 of 32 The nondimensional downslope and normal components of the traction free boundary condition on the free-surface (26) are as follows: On the interface, the nondimensional downslope and normal components of the friction condition (27) are as follows: On the bottom, the nondimensional downslope and normal components of the friction law (28) are The nondimensional downslope and normal components of the momentum jump condition (30) are respectively.

Appendix B: The Pore Fluid Pressure and the Grain Normal Stress
We consider that the basal topography varies gently, and hence the variable λ = L/R, the function |∇̂| and the factor Ψ are, respectively, of magnitude where 0 < α < 1.
In shallow flows, the assumption of hydrostatic pore fluid pressure is usually made. Nevertheless, the vertical relative motion between grains and fluid plays a crucial role in the generation of the excess pore fluid pressure, which has a significant influence on the motion of particles. To leading order in the small parameter ɛ, the normal component of Equation 16 reduces to in the pure fluid regime and Equation 2 to The integration of Equation B2 to the free-surface, alongside the traction-free condition, Equation A10a, shows that the pressure is hydrostatiĉ= in the pure fluid regime. On the interface = ( ) , the pore fluid pressure is as follows: To leading order in ɛ the condition Equation A15b reduces to the following equation: Integration of Equation B3a in the mixture regime of grains and fluid to the interface gives the following equation:̂= where the excess pore fluid pressure is expressed as follows: Then, the basal pore fluid pressure is as follows: Substituting Equation B3a into Equation B3b gives the following equation: It is assumed that the grain normal stress vanishes on the interface, that is, ( ) = 0 . Thus, the integration of Equation B10 gives the following equation: Moreover, provided the earth pressure coefficient is equal to unity (Meng et al., 2022;Savage & Hutter, 1989), the downslope grain normal stress takes the following form: The basal grain normal stress is as follows: Vertical integration of Equation 41 from the base to any position within the mixture regime of grains and fluid, in which the no-slip boundary condition is applied, yields the following equation: where the small parameter ɛ emerges, implying that the dilatancy angle ψ is small. Since permeability is small in geophysical flows, the Darcy drag is strong enough that the product of the drag coefficient ̂ and the difference in vertical velocity − is of the same order as the hydrostatic pressure.

Substitution of Equation B14
into Equations B9 and B13 gives the following equations: where volume fractions are assumed to be constant in the mixture regime of grains and fluid.
Substitution of Equation B17 into Equations B15 and B16, respectively, shows the basal pore fluid pressure and the solid basal stress takes the following forms: As the scaling Equation A1 is reversely used to transform Equation B18 into a dimensional form, the dimensional basal pore fluid pressure, that is, Equation 47 in the main text, is obtained.

C1. Depth-Averaged Mass Balance Equations
The grain mass balance Equation A2 is integrated from the base to the interface using Leibniz's rule to swap the order of differentiation and integration, and then simplified by using conditions Equations A10b, A10c, which gives the following equation: The mass balance equation of the bulk consisting of grains and fluid, Equation A4, is integrated from the base to the interface and simplified by using Equations A10c and A11, showing the following equation: Analogously, integration in the direction normal to the topography of the mass balance equation of the pure fluid Equation A7 from the interface to the free-surface gives the following equation: 27 of 32

C2.1. Mixture Regime of Grains and Fluid
The momentum balance Equation A5 can now be integrated through the mixture thickness by using Leibniz's rule to swap the order of differentiation and integration. The resulting equations can then be simplified by using the surface and basal kinematic conditions, Equations A10b, A10c-A12, the downslope basal and surface traction conditions Equations A15a-A17a, and the ordering (B1) to give the following equation: where the basal bulk normal stress ( ⋅̂) and the basal grain normal stress Integration of the normal momentum balance Equation A6 through the mixture thickness to leading order in ɛ shows the following equation: The shape factor χ m is introduced to characterize the ratio of (̂) 2 to (̂) 2 , that is, = (̂) 2 ∕ (̂) 2 . Substituting the grain downslope normal stress Equation B12, the basal bulk normal stress Equation C6, the basal grain normal stress Equation C7 and the shape factor into Equation C5 results in the depth-averaged momentum equation of the bulk consisting of grains and fluid in the following form: (C8)

C2.2. Pure Fluid Regime
Analogously, the momentum balance Equation A8 is integrated from the interface to the free-surface. Leibniz's integration rule is used to swap the order of differentiation and integration, and the surface and the massjump kinematic conditions, Equations A10a and A12, the downslope surface traction conditions Equations A14a and A15a, and the ordering Equation B1 are used to simplify the resulting equations. It follows: The shape factor χ pf is introduced to characterize the ratio of (̂) 2 to ̂ 2 , that is, Substituting the pore fluid pressure Equation B4 and the shape factor Equation C10 into Equation C9 shows that the depth-averaged momentum equation of the pure fluid is as follows: As the scaling Equation A1 is reversely used to transform the nondimensional mass balance Equations C2-C4 and the momentum balance Equations C8 and C11 into a dimensional form, the dimensional Equations 33-37 in the main text are obtained.