Modeling of unsaturated granular flows by a two-layer approach

Abstract

Flows of partially saturated grain-fluid mixtures over complex curved topography are commonly observed in nature. However, comprehensive understanding of the physics behind them is to date out of reach. To investigate their dynamic process, a two-layer approach is proposed, in which the fluid-saturated granular layer is overlaid by the pure granular material. More specifically, the lower layer is described by a two-phase mixture theory of density preserving solid and fluid constituents. For the upper layer, the single-phase granular mass is treated as a frictional Coulomb-like continuum, and the dilation effect and the influence of the interstitial air are ignored. The capillarity effects and grains-size segregation are not considered in both the layers. The lower and upper layers interact at an interface which is a material surface for the fluid phase, but across which the mass exchange for the granular phase may take place. The granular mass exchange across the layer interface is parameterized by an entrainment type postulate. In addition, the classical jump conditions are employed to connect both layers at the interface dynamically. Furthermore, we perform the depth-averaged technique for the saturated grain-fluid mixture lower layer and the pure granular upper layer, respectively, to simplify the governing equations established. It is demonstrated that the resulting model equations can be reduced to most of the existing single-layer pure granular flow models and saturated two-phase single-layer debris flow models. Numerical solutions demonstrate that the present two-layer model can describe flows of partially saturated grain-fluid mixtures and the transition process of a saturated grain-fluid mixture into an under-saturated state.

Appendices

Appendix 1: Treatment of boundary conditions

• 1. The kinematic boundary condition (9) at the top free surface can be expanded as

  $$\frac{\partial s}{\partial t} + u_g^s\frac{\partial s}{\partial x} + v_g^s\frac{\partial s}{\partial y} - w_g ^s= 0.$$

  (68)

  where downslope velocity is represented as $u_g^s$, cross-slope velocity as $v_g^s$, and normal velocity as $w_g^s$.

  The down-slope, cross-slope and normal components of the dynamic boundary condition, (10), are respectively

  $$T_{g(xx)}^s\frac{\partial s}{\partial x}+T_{g(xy)}^s\frac{\partial s}{\partial y}-T_{g(xz)}^s =0,$$

  (69)
  $$T_{g(yx)}^s\frac{\partial s}{\partial x}+T_{g(yy)}^s\frac{\partial s}{\partial y}- T_{g(yz)}^s=0,$$

  (70)
  $${T}_{g(zx)}^s\frac{\partial s}{\partial x}+{T}_{g(zy)}^s\frac{\partial s}{\partial y}-{T}_{g(zz)}^s=0.$$

  (71)

• 2. Similarly, at the bottom, the kinematic boundary condition (11) for each constituent of the lower layer can be expanded as follows:

  $$w_s^b = 0, \qquad\qquad w_f^b = 0.$$

  (72)

  For the solid constituent of the lower-layer, the down-slope, and cross-slope components of the bed friction boundary condition (12) are

  $${T}_{e(xz)}^b =-\mu _s^b\frac{u^b_s}{\mid {{\varvec{u}}}^b_s\mid }{T}_{e(zz)}^b -(k_s^b\mu _s^b{\phi _{s}}^b)u_s^b,$$

  (73)
  $${T}_{e(yz)}^b=-\mu _s^b\frac{v^b_s}{\mid {{\varvec{u}}}^b_s\mid }{T}_{e(zz)}^b -(k_s^b\mu _s^b{\phi _{s}}^b)v_s^b,$$

  (74)

  respectively. Similarly, for the fluid constituent of the lower-layer, the down-slope, and cross-slope components of the bed friction condition (13) are

  $$\tau _{f(xz)}^b=k_f^b\,u_f^b,$$

  (75)
  $$\tau _{f(yz)}^b =k_f^b\,v_f^b,$$

  (76)

• 2. At the interface, the kinematic boundary condition (14) is expanded as

  $$\frac{\partial h_m}{\partial t}+u_f^i\frac{\partial h_m}{\partial x}+v_f^i\frac{\partial h_m}{\partial y}-w_f^i=0,$$

  (77)

  and the traction-free dynamic boundary condition (17) in the down-slope, cross-slope and normal directions are

  $$(-p^i+\tau _{f(xx)}^i)\frac{\partial s}{\partial x}+\tau _{f(xy)}^i\frac{\partial s}{\partial y}-\tau _{f(xz)}^i =0,$$

  (78)
  $$\tau _{f(yx)}^i\frac{\partial s}{\partial x}+(-p^i+\tau _{f(yy)}^i)\frac{\partial s}{\partial y}- \tau _{f(yz)}^i=0,$$

  (79)
  $$\tau _{f(zx)}^i\frac{\partial s}{\partial x}+\tau _{f(zy)}^i\frac{\partial s}{\partial y}-(-p^i+\tau _{f(zz)}^i)=0,$$

  (80)

  respectively. Relations (78)–(80) indicate that $p^i=0$ and $\tau _{f(jk)}^i=0, (j, k\in \{x, y, z\})$.

  The mass jump condition (15) (or (16)) is written as

  $$\begin{aligned} J&= -{\widetilde{\rho }}_s\phi _g\left[ \frac{\partial h_m}{\partial t}+u_s^i\frac{\partial h_m}{\partial x}+v_s^i\frac{\partial h_m}{\partial y}-w_s^i\right] \nonumber \\ &= -{\widetilde{\rho }}_s\phi _g\left[ \frac{\partial h_m}{\partial t}+u_g^i\frac{\partial h_m}{\partial x}+v_g^i\frac{\partial h_m}{\partial y} -w_g^i\right] . \end{aligned}$$

  (81)

  Similarly, the components of the momentum jump condition (18) in the down-slope, cross-slope and normal directions are

  $$-J(u_s^i-u_g^i)+T_{s(xx)}^i\frac{\partial h_m}{\partial x}+{T}_{s(xy)}^i\frac{\partial h_m}{\partial y}-{T}_{s(xz)}^i ={T}_{g(xx)}^i\frac{\partial h_m}{\partial x}+{T}_{g(xy)}^i\frac{\partial h_m}{\partial y}-{T}_{g(xz)}^i,$$

  (82)
  $$-J(v_s^i-v_g^i)+{T}_{s(yx)}^i\frac{\partial h_m}{\partial x}+{T}_{s(yy)}^i\frac{\partial h_m}{\partial y}-{T}_{s(yz)}^i ={T}_{g(yx)}^i\frac{\partial h_m}{\partial x}+{T}_{g(yy)}^i\frac{\partial h_m}{\partial y}-{T}_{g(yz)}^i,$$

  (83)
  $$-J(w_s^i-w_g^i)+{T}_{s(zx)}^i\frac{\partial h_m}{\partial x}+{T}_{s(zy)}^i\frac{\partial h_m}{\partial y}-{T}_{s(zz)}^i ={T}_{g(zx)}^i\frac{\partial h_m}{\partial x}+{T}_{g(zy)}^i\frac{\partial h_m}{\partial y}-{T}_{g(zz)}^i,$$

  (84)

  respectively, where ${T}_{s(jk)}^i={T}_{e(jk)}^i$ ($j, k\in (x, y, z)$) due to $p^i=0$ indicated by (80).

  To rewrite (19) in a component form, the tangential velocity ${{\varvec{u}}}_k^{\tau }$ ($k\in \{s, g\}$) is rewritten as ${{\varvec{u}}}_k^{\tau }={{\varvec{u}}}_k^i -({{\varvec{u}}}_k^i\cdot {{\varvec{n}}}^i){{\varvec{n}}}^i$ that yields ${{\varvec{u}}}_s^{\tau }-{{\varvec{u}}}_g^{\tau }={{\varvec{u}}}_s^i-{{\varvec{u}}}_g^i$ due to ${{\varvec{u}}}_s^i\cdot {{\varvec{n}}}^i={{\varvec{u}}}_g^i\cdot {{\varvec{n}}}^i$ from the jump condition (16) together with the continuity assumption of the solid volume fraction across the interface. It follows that the components of relation (19) in the down-slope, cross-slope and normal directions are,

  $$\begin{aligned}&{T}_{e(xx)}^i\frac{\partial h_m}{\partial x}+{T}_{e(xy)}^i\frac{\partial h_m}{\partial y}-{T}_{e(xz)}^i \nonumber \\ &\quad =({{\varvec{n}}}^i\cdot {{{\varvec{T}}}}_e^i{{\varvec{n}}}^i) \left( \frac{\partial h_m}{\partial x}-\frac{u^i_s-u^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i\right) -k_s^i\mu _s^i\phi _g(u^i_s-u^i_g) ,\end{aligned}$$

  (85)
  $$\begin{aligned}&{T}_{e(yx)}^i\frac{\partial h_m}{\partial x}+{T}_{e(yy)}^i\frac{\partial h_m}{\partial y}-{T}_{e(yz)}^i \nonumber \\ &\quad = ({{\varvec{n}}}^i\cdot {{{\varvec{T}}}}_e^i{{\varvec{n}}}^i)\left( \frac{\partial h_m}{\partial y}-\frac{v^i_s-v^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i\right) -k_s^i\mu _s^i\phi _g(v^i_s-v^i_g) , \end{aligned}$$

  (86)
  $$\begin{aligned}&{T}_{e(zx)}^i\frac{\partial h_m}{\partial x}+{T}_{e(zy)}^i\frac{\partial h_m}{\partial y}-{T}_{e(zz)}^i\nonumber \\ &\quad = ({{\varvec{n}}}^i\cdot {{{\varvec{T}}}}_e^i{{\varvec{n}}}^i)\left( -1- \frac{w^i_s-w^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i\right) - k_s^i\mu _s^i\phi _g(w^i_s-w^i_g), \end{aligned}$$

  (87)

  respectively.

Appendix 2: Depth-averaged procedure

1.1 Conservation equations for mass

Integrating the mass equation (1) of the upper-layer granular phase by using Leibniz rule to swap the orders of integration and differentiation, one can obtain

$$\begin{aligned}&\frac{\partial }{\partial t}(h_g{\widetilde{\rho }}_s\phi _g)+\frac{\partial }{\partial x}( h_g{\widetilde{\rho }}_s\phi _g{\overline{u}}_g)+\frac{\partial }{\partial y}( h_g{\widetilde{\rho }}_s\phi _g{{\overline{v}}}_g)\nonumber \\ &\quad -\left[ {\widetilde{\rho }}_s\phi _g\left( \frac{\partial z}{\partial t}+u_g\psi \frac{\partial z}{\partial x}+v_g\frac{\partial z}{\partial y}-w_g\right) \right] _i^s=0, \end{aligned}$$

(88)

where the square bracket denotes the difference of the top and interfacial values of a function q, i.e., $[q]_i^s=q^s-q^i$. The term within $[\cdot ]$ can be further handled by the kinematic boundary condition (68) at the top surface, and jump condition (81), which leads to

$$\frac{\partial }{\partial t}(h_g{\widetilde{\rho }}_s\phi _g)+\frac{\partial }{\partial x}( h_g{\widetilde{\rho }}_s\phi _g{\overline{u}}_g)+\frac{\partial }{\partial y}( h_g{\widetilde{\rho }}_s\phi _g{{\overline{v}}}_g) =J.$$

(89)

Similarly, integrating (3) and (4) can yield the depth-averaged mass-conservation equations for the solid and fluid constituents of the lower layer. They are given by

$$\frac{\partial }{\partial t}({h}_m{\overline{\rho }}_s)+\frac{\partial }{\partial x}({h}_m\overline{\rho _su_s})+\frac{\partial }{\partial y}({h}_m\overline{\rho _sv_s}) =-J,$$

(90)

$$\frac{\partial }{\partial t}({h}_m{\overline{\rho }}_f) + \frac{\partial }{\partial x}({h}_m\overline{\rho _fu_f})+\frac{\partial }{\partial y}({h}_m\overline{\rho _fv_f})=0,$$

(91)

where ${\overline{\rho }}_s={\widetilde{\rho }}_s{\overline{\phi }}_s$ and ${\overline{\rho }}_f={\widetilde{\rho }}_f{\overline{\phi }}_f$. We follow Pitman and Le [30] and Iverson and Denlinger [20] to assume that volume fractions are uniformly distributed along the depth, such that equations (90) and (91) can be transformed into

$$\frac{\partial }{\partial t}({h}_m{\overline{\rho }}_s)+\frac{\partial }{\partial x}({h}_m{\overline{\rho }}_s{\overline{u}}_s)+\frac{\partial }{\partial y}({h}_m{\overline{\rho }}_s{{\overline{v}}}_s) =-J,$$

(92)

$$\frac{\partial }{\partial t}({h}_m{\overline{\rho }}_f) + \frac{\partial }{\partial x}({h}_m{\overline{\rho }}_f{\overline{u}}_f)+\frac{\partial }{\partial y}({h}_m{\overline{\rho }}_f{{\overline{v}}}_f)=0.$$

(93)

1.2 Normal component of momentum equations

Since the Mohr–Coulomb friction criterion, $\mid {\varvec{S}}\mid =N\tan \varphi$, does not explicitly show connections between stress and strain rate tensor, we therefore use shallow-water approximation to simplify the conservation equations. The shallow-water approximation indicates that the typical depth ${\mathcal {H}}$ is much smaller than the typical length ${\mathcal {L}}$, i.e., the geometrical aspect ratio $\epsilon ={\mathcal {H}}/{\mathcal {L}}\ll 1$ is a small quantity. In classical pure granular flow model (see [34]), the normal momentum balance usually reduces to hydrostatic force balance, when shallow-water approximation is employed. Consequently, the upper-layer granular normal momentum equation reduces to

$$\frac{\partial T_{g(zz)}}{\partial z}={\widetilde{\rho }}_s\phi _gg_z,$$

(94)

where $g_z=g\cos \zeta$. Relation (94) can be vertically integrated from any position z to the top free surface, which yields

$$T_{g(zz)}={\widetilde{\rho }}_s\phi _gg_z(s-z)$$

(95)

by using traction-free condition at the free surface.

Similarly, in most of the previous two-phase grain-fluid mixture models, the fluid normal momentum balance implies that the pore pressure is hydrostatic, and the solid normal momentum balance indicates a balance of gravity, internal stress, and buoyancy force, when shallow-water approximation is employed to streamline the normal momentum balances. Consequently, the lower-layer fluid and solid normal momentum balances reduce to

$$\phi _f\frac{\partial p}{\partial z} + {\widetilde{\rho }}_f\phi _fg_z= 0,$$

(96)

$${\widetilde{\rho }}_s{\phi _{s}}g_z + \frac{\partial T_{e(zz)}}{\partial z} +{\phi _{s}}\frac{\partial p}{\partial z}= 0,$$

(97)

respectively. Integrating relation (96) from any position z in the lower-layer mixture to the interface can yield

$$p(z) = {\widetilde{\rho }}_fg_z (h_m-z),$$

(98)

by using the traction-free condition $p^i=0$ at the interface. Similarly, vertical integration of (97) leads to

$$T_{e(zz)} = ({\widetilde{\rho }}_s-{\widetilde{\rho }}_f){\overline{\phi }}_sg_z (h_m-z)+T_{e(zz)}^i,$$

(99)

where the stress $T_{e(zz)}^i$ at the interface remains to be determined. To evaluate $T_{e(zz)}^i$, we appeal to relation (84). In shallow granular flow models, the derivatives $\partial h_m/\partial x$ and $\partial h_m/\partial y$ are small quantities, and the magnitudes of the normal velocities, $w_\eta$ ($\eta \in \{g, s, f\}$), are smaller than the magnitudes of the horizontal velocities, $u_\eta$ and $v_\eta$. With these descriptions, relation (84) reduces to

$$T_{e(zz)}^i=T_{s(zz)}^i = T_{g(zz)}^i = {\widetilde{\rho }}_s\phi _gh_gg_z.$$

(100)

Substitution of relation (100) into (99) leads to

$$T_{e(zz)} = ({\widetilde{\rho }}_s-{\widetilde{\rho }}_f){\overline{\phi }}_sg_z (h_m-z)+ {\widetilde{\rho }}_s\phi _gh_gg_z.$$

(101)

1.3 Horizontal momentum-conservation equations

• 1. For the pure granular phase in the upper layer, integrating downslope component of equation (2) throughout the upper layer by using Leibniz rule to interchange the orders of integration and differentiation can yield the depth-averaged form. The left-hand side terms turn into

  $$\begin{aligned} \int _{h_m}^s\hbox{LHS}\,\hbox{d}z &=\frac{\partial }{\partial t}(h_g{\widetilde{\rho }}_s\phi _g{\overline{u}}_g)+\frac{\partial }{\partial x}(h_g {\widetilde{\rho }}_s\phi _g{\overline{u}}_g^2)+\frac{\partial }{\partial y}(h_g{\widetilde{\rho }}_s\phi _g{\overline{u}}_g{{\overline{v}}}_g)-Ju_g^i , \end{aligned}$$

  (102)

  where we used plug-flow assumption as previous studies did (see [30]). The right-hand side terms turn into

  $$\begin{aligned} \int _{h_m}^s \hbox{RHS}\,\hbox{d}z&= {\widetilde{\rho }}_g\phi _gh_gg_x -\frac{\partial (h_g{\overline{T}}_{g(xy)})}{\partial y}-\frac{\partial ( h_g{\overline{T}}_{g(xx)})}{\partial x}\nonumber \\ &+\left[ {T}_{g(xx)}\frac{\partial z}{\partial x} +T_{g(xy)}\frac{\partial z}{\partial y}-{T}_{g(xz)}\right] _{h_m}^s, \end{aligned}$$

  (103)

  where $g_x=g\sin \zeta$, and the term $-\partial ( h_g{\overline{T}}_{g(xy)})/{\partial y}$ is usually ignored, since it is a high-order term with respect to aspect ratio $\epsilon$ (see [8]). The gradient of lateral stress, $-\partial (h_g{\overline{T}}_{g(xx)})/{\partial x}$, can be determined through the following [34]. We postulate that $T_{g(xx)}$ is proportional to $T_{g(zz)}$ throughout the upper-layer depth, i.e., $T_{g(xx)}=K_x^gT_{g(zz)}$, where $K_x^g$ is an earth pressure coefficient and defined in (28). Then, ${\overline{T}}_{g(xx)}$ can be derived by vertically integrating relation $T_{g(xx)}=K_x^gT_{g(zz)}$. They result in

  $$\begin{aligned} \frac{\partial (h_g{\overline{T}}_{g(xx)})}{\partial x}=\frac{\partial }{\partial x}\left( \frac{1}{2}K^g_x h_g^2{\widetilde{\rho }}_s\phi _gg_z\right) . \end{aligned}$$

  (104)

  Moreover, the terms in $[\cdot ]$ of relation (103) can reduce to

  $$\begin{aligned} \left[ {T}_{g(xx)}\frac{\partial z}{\partial x} +{T}_{g(xy)}\frac{\partial z}{\partial y}-{T}_{g(xz)}\right] _{h_m}^s&=-({{\varvec{n}}}^i\cdot {{{\varvec{T}}}}_e^i{{\varvec{n}}}^i)\left( \frac{\partial h_m}{\partial x}-\frac{u^i_s-u^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i\right) +J(u_s^i-u_g^i)\nonumber \\ &\quad +k_s^i\mu _s^i\phi _g(u^i_s-u^i_g), \end{aligned}$$

  (105)

  by using boundary conditions (69), (82), and (85), where the normal stress $({{\varvec{n}}}^i\cdot {{{\varvec{T}}}}_e^i{{\varvec{n}}}^i)$ at the interface remains to be determined. In shallow flows, the normal stress $({{\varvec{n}}}^i\cdot {{{\varvec{T}}}}_e^i{{\varvec{n}}}^i)$ can be approximated as normal traction, $T_{e(zz)}^i$, plus a centrifugal force term caused by curvature, see [8]. We assume the same curvature between the interface and the bottom that indicates that the centrifugal force is $h_g{\widetilde{\rho }}_s\phi _g\kappa {\overline{u}}_g^2$. These descriptions yield

  $$\begin{aligned} ({{\varvec{n}}}^i\cdot {{{\varvec{T}}}}_e^i{{\varvec{n}}}^i)= h_g{\widetilde{\rho }}_s\phi _gg_z+ h_g{\widetilde{\rho }}_s\phi _g\kappa {\overline{u}}_g^2. \end{aligned}$$

  (106)

  Similarly, integrating the cross-slope component of equation (2) throughout the upper layer can yield the depth-averaged form. The left-hand side terms turn into

  $$\begin{aligned} \int _{h_m}^s\hbox{LHS}\,\hbox{d}z &=\frac{\partial }{\partial t}(h_g{\widetilde{\rho }}_s\phi _g{{\overline{v}}}_g)+\frac{\partial }{\partial x}(h_g{\widetilde{\rho }}_s \phi _g{\overline{u}}_g{{\overline{v}}}_g)+\frac{\partial }{\partial y}(h_g{\widetilde{\rho }}_s\phi _g{{\overline{v}}}_g^2)-Jv_g^i , \end{aligned}$$

  (107)

  and the right-hand side terms are transformed into

  $$\begin{aligned} \int _{h_m}^s\hbox{RHS}\,\hbox{d}z &= -\frac{\partial }{\partial y}\left( \frac{1}{2} K_y^g{\widetilde{\rho }}_sh_g^2\phi _gg_z\right) -{\widetilde{\rho }}_s\phi _gh_gg_z\frac{\partial h_m}{\partial y}+J(v_s^i-v_g^i) \nonumber \\ & \quad +k_s^i\mu _s^i\phi _g(v^i_s-v^i_g) +\frac{v^i_s-v^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i{\widetilde{\rho }}_s h_g\phi _g \left( g_z+\kappa {\overline{u}}_g^2\right) , \end{aligned}$$

  (108)

  where the first term on the right-hand side represents gradient of internal stress; the second term represents the effect of lower-layer height, the third term denote the momentum exchange, and the fourth terms and the fifth term represent viscous and Coulomb friction at the interface.

• 2. For the solid constituent in the lower layer, integrating downslope component of equation (7) can yield the depth-averaged form. The left-hand side terms turn into

  $$\begin{aligned} \int ^{h_m}_0 \hbox{LHS}\,\hbox{d}z&=\frac{\partial }{\partial t}(h_m{\overline{\rho }}_s{\overline{u}}_s)+\frac{\partial }{\partial x}(h_m {\overline{\rho }}_s{\overline{u}}_s^2)+\frac{\partial }{\partial y}(h_m{\overline{\rho }}_s{\overline{u}}_s{{\overline{v}}}_s)+Ju_s^i, \end{aligned}$$

  (109)

  and the right-hand side terms are transformed into

  $$\begin{aligned} \int ^{h_m}_0 \hbox{RHS}\,\hbox{d}z&= -\frac{\partial (h_m\overline{{T}}_{e(xy)})}{\partial y} -\frac{\partial (h_m{\overline{T}}_{e(xx)})}{\partial x}-\int _0^{h_m}{\phi _{s}}\frac{\partial p}{\partial x}\hbox{d}z \nonumber \\ &+\left[ {T}_{e(xx)}\frac{\partial z}{\partial x} +T_{e(xy)}\frac{\partial z}{\partial y}-T_{e(xz)}\right] _0^{h_m}\nonumber \\ &+C_d h_m{\overline{\phi }}_s{\overline{\phi }}_f({\overline{u}}_f-{\overline{u}}_s) +h_m{\overline{\rho }}_sg_x , \end{aligned}$$

  (110)

  where the lateral stress gradient, $\partial (h_m\overline{{T}}_{e(xy)})/\partial y$, can be ignored as we mentioned above. The other stress gradient, $\partial (h_m{\overline{T}}_{e(xx)})/\partial x$, can be determined by introducing earth pressure coefficient. Referring to Savage and Hutter [34] and its modifications (see, e.g., Hutter et al. [17]; Gray et al. [8]), the following relation

  $$\begin{aligned} {T}_{e(xx)}^b=K^s_x{T}_{e(zz)}^b \end{aligned}$$

  (111)

  holds at the base, where the bed earth pressure coefficient $K^s_x$ is defined in (36).

  Moreover, relation ${T}_{e(xx)}^i={T}_{g(xx)}^i$ is postulated, which does not violate the momentum jump condition (82). Since relation ${T}_{g(xx)}^i=K^g_x{T}_{g(zz)}^i$ holds at the interface, relation ${T}_{e(xx)}^i=K^g_x{T}_{e(zz)}^i$ can be obtained. To proceed, we postulate that horizontal stresses vary linearly with normal stress throughout the depth of the lower layer. However, the factors are linear with z-coordinate instead of constant value used previously (see, e.g., Savage and Hutter [34]). This more flexible approach leads to

  $$\begin{aligned}&{T}_{e(xx)}=\left( \frac{K^g_x-K^s_x}{{h}_m}z+K^s_x\right) {T}_{e(zz)}, \end{aligned}$$

  (112)

  which is a generalization of Savage–Hutter theory. When $K^g_x=K^s_x$, relation (112) reproduces the form used in the Savage–Hutter theory. Additionally, relation (112) will produce ${T}_{e(xx)}^i=K^g_x{T}_{e(zz)}^i$ at the interface, and (111) at the base.

  Exploiting relations (101) and (112), one can yield

  $$\begin{aligned} \frac{\partial }{\partial x}(h_m{\overline{T}}_{e(xx)})&= \frac{\partial }{\partial x}\left[ \frac{1}{2} h_gh_m{\widetilde{\rho }}_s\phi _gg_z(K^s_x+K^g_x)\right. \nonumber \\ &\quad \left. +\frac{1}{6}({\widetilde{\rho }}_s-{\widetilde{\rho }}_f) g_zh_m^2{\overline{\phi }}_s(2K^s_x+K^g_x)\right] . \end{aligned}$$

  (113)

  Moreover, the term involving fluid pressure gradient, arising in relation (110), can be formulated as

  $$\begin{aligned} -\int _0^{h_m}{\phi _{s}}\frac{\partial p}{\partial x}\hbox{d}z=-{\widetilde{\rho }}_f{\overline{\phi }}_s\frac{\partial }{\partial x}\left( \frac{1}{2}h_m^2g_z\right) , \end{aligned}$$

  (114)

  by using relation (98). The terms in $[\cdot ]$, arising in relation (110), can be simplified as

  $$\begin{aligned}&\left[ {T}_{e(xx)}\frac{\partial z}{\partial x} +T_{e(xy)}\frac{\partial z}{\partial y}-T_{e(xz)}\right] _0^{h_m} \nonumber \\ &\quad =({{\varvec{n}}}^i\cdot {{{\varvec{T}}}}_e^i{{\varvec{n}}}^i) \left( \frac{\partial h_m}{\partial x}-\frac{u^i_s-u^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i\right) -k_s^i\mu _s^i\phi _g(u^i_s-u^i_g)\nonumber \\ &\qquad -\mu _s^b\frac{u_s^b}{\mid {{\varvec{u}}}_s^b\mid }T_{e(zz)}^b -(k_s^b\mu _s^b{\phi _{s}}^b)u_s^b, \nonumber \\ &\quad =-\frac{u^i_s-u^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^ih_g{\widetilde{\rho }}_s\phi _g (g_z+\kappa {\overline{u}}_g^2)+h_g{\widetilde{\rho }}_s\phi _gg_z\frac{\partial h_m}{\partial x}\nonumber \\ &\qquad -k_s^i\mu _s^i\phi _g(u^i_s-u^i_g)-(k_s^b\mu _s^b{\phi _{s}}^b)u_s^b \nonumber \\ &\quad \quad -\frac{v_s^b}{\mid {{\varvec{u}}}_s^b\mid }\mu _s^b \left[ {\widetilde{\rho }}_sg_z\left( h_m{\overline{\phi }}_s (1-\gamma )+h_g\phi _g\right) \right. \nonumber \\ &\qquad \left. +\kappa \left( h_m{\overline{\phi }}_s({\widetilde{\rho }}_s {\overline{u}}_s^2 -{\widetilde{\rho }}_f{\overline{u}}_f^2 )+{\widetilde{\rho }}_sh_g\phi _g{\overline{u}}_g^2\right) \right] , \end{aligned}$$

  (115)

  by substituting boundary conditions (73) and (85). In (115), the overburden pressure $T_{e(zz)}^b$ at the bottom is written as a normal traction plus a curvature-dependency term, i.e., $T_{e(zz)}^b={\widetilde{\rho }}_sg_z\left( h_m{\overline{\phi }}_s (1-\gamma )+h_g\phi _g\right) +\kappa \left( h_m{\overline{\phi }}_s({\widetilde{\rho }}_s {\overline{u}}_s^2 -{\widetilde{\rho }}_f{\overline{u}}_f^2 )+{\widetilde{\rho }}_sh_g\phi _g{\overline{u}}_g^2\right)$, since the bed curvature also contributes a part of pressure and furthermore affects the friction at the bottom, see [8].

  Similarly, integrating the cross-slope component of equation (7) yields that the left-hand side terms turn into

  $$\begin{aligned} \int _0^{h_m}\hbox{LHS}\,\hbox{d}z =\frac{\partial }{\partial t}(h_m{\overline{\rho }}_s{{\overline{v}}}_s)+\frac{\partial }{\partial x}(h_m {\overline{\rho }}_s{\overline{u}}_s{{\overline{v}}}_s)+\frac{\partial }{\partial y}(h_m{\overline{\rho }}_s{{\overline{v}}}_s^2)+Jv_s^i, \end{aligned}$$

  (116)

  and the right-hand side terms turn into

  $$\begin{aligned}&\int _0^{h_m}\hbox{RHS}\,\hbox{d}z \nonumber \\ &\quad =-\frac{\partial }{\partial y}\left[ \frac{1}{2} h_gh_m{\widetilde{\rho }}_s\phi _gg_z(K^s_y+K^g_y)+\frac{1}{6}({\widetilde{\rho }}_s -{\widetilde{\rho }}_f) g_zh_m^2{\overline{\phi }}_s(2K^s_y+K^g_y)\right] \nonumber \\ &\qquad -\frac{v_s^b}{\mid {{\varvec{u}}}_s^b\mid }\mu _s^b \left[ {\widetilde{\rho }}_s\cos \zeta \left( h_m{\overline{\phi }}_s (1-\gamma )+h_g\phi _g\right) \right. \nonumber \\ &\qquad \left. +\kappa \left( h_m{\overline{\phi }}_s({\widetilde{\rho }}_s {\overline{u}}_s^2 -{\widetilde{\rho }}_f{\overline{u}}_f^2 )+{\widetilde{\rho }}_sh_g\phi _g{\overline{u}}_g^2\right) \right] \nonumber \\ &\qquad -\frac{v^i_s-v^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^ih_g{\widetilde{\rho }}_s\phi _g (g_z+\kappa {\overline{u}}_g^2) -k_s^i\mu _s^i\phi _g(v^i_s-v^i_g)\nonumber \\ &\qquad +C_d h_m{\overline{\phi }}_s{\overline{\phi }}_f({{\overline{v}}}_f-{{\overline{v}}}_s), \end{aligned}$$

  (117)

  where the first term of the right-hand side represents the cross-slope gradient of effective stress; the second term accounts for the bed Coulomb friction; and the rest terms describe interfacial Coulomb friction, interfacial viscous friction, and viscous drag between the two phases of the lower mixture consecutively.

• 3. For the fluid phase of the lower layer, vertically integrating the downslope component of equation (6) allows to derive the depth-averaged form. The left-hand side terms turn into

  $$\begin{aligned} \int ^{h_m}_0\hbox{LHS}\,\hbox{d}z =\frac{\partial }{\partial t}({h}_m{\overline{\rho }}_f{\overline{u}}_f)+\frac{\partial }{\partial x}(h_m {\overline{\rho }}_f{\overline{u}}_f^2)+\frac{\partial }{\partial y}(h_m{\overline{\rho }}_f{\overline{u}}_f{{\overline{v}}}_f). \end{aligned}$$

  (118)

  The right-hand side terms turn into

  $$\begin{aligned} \int ^{h_m}_0\hbox{RHS}\,\hbox{d}z &= -{\overline{\rho }}_f\frac{\partial }{\partial x}\left( \frac{1}{2}h_m^2g_z\right) + \frac{\partial }{\partial x}(h_m{\overline{\phi }}_f{\overline{\tau }}_{f(xx)}) +\frac{\partial }{\partial y}(h_m{\overline{\phi }}_f{\overline{\tau }}_{f(xy)}) \nonumber \\ &-\left[ -\phi _fp\frac{\partial z}{\partial x}+\phi _f\tau _{f(xx)}\frac{\partial z}{\partial x} +\phi _f\tau _{f(xy)}\frac{\partial z}{\partial y}-\phi _f\tau _{f(xz)}\right] _0^{h_m} \nonumber \\ &+{\overline{\rho }}_fg_x\sin \zeta - C_dh_m{\overline{\phi }}_s{\overline{\phi }}_f({\overline{u}}_f-{\overline{u}}_s) \end{aligned}$$

  (119)

  where the first term on the right-hand side, $-{\overline{\rho }}_f\partial (h_m^2g_z)/2\partial x$, accounts for the contribution of fluid pressure. The second and third terms represent the contribution of diffusion terms; however, they are often ignored, since they have a negligible effect on dynamics compared to the bed friction (terms in $[\cdot ]$), see [19]. The terms within $[\cdot ]$ appear due to mathematical operation (Leibniz rule), and they can be simplified as follows

  $$\begin{aligned} \left[ -\phi _fp\frac{\partial z}{\partial x}+\phi _f\tau _{f(xx)}\frac{\partial z}{\partial x} +\phi _f\tau _{f(xy)}\frac{\partial z}{\partial y}-\phi _f\tau _{f(xz)}\right] _0^{h_m} =k_f^bu_f^b \end{aligned}$$

  (120)

  by using boundary conditions (75) and (78).

  Similarly, vertically integrating the cross-slope component of equation (6) can yield the depth-averaged form. The left-hand side terms turn into

  $$\begin{aligned} \int ^{h_m}_0\hbox{LHS}\,\hbox{d}z =\frac{\partial }{\partial t}({h}_m{\overline{\rho }}_f{{\overline{v}}}_f)+\frac{\partial }{\partial x}(h_m {\overline{\rho }}_f{\overline{u}}_f{{\overline{v}}}_f)+\frac{\partial }{\partial y}(h_m{\overline{\rho }}_f{{\overline{v}}}_f^2). \end{aligned}$$

  (121)

  The right-hand side terms turn into

  $$\begin{aligned} \int ^{h_m}_0 \hbox{RHS}\,\hbox{d}z&=-{\overline{\rho }}_f\frac{\partial }{\partial y}\left( \frac{1}{2}h_m^2g_z\right) - k_f^bv_f^b -C_dh_m{\overline{\phi }}_s{\overline{\phi }}_f({{\overline{v}}}_f-{{\overline{v}}}_s), \end{aligned}$$

  (122)

  where the first term on the right-hand side represents the contribution of fluid pressure; the second term stems from the bed friction; and the last term stands for viscous drag between the fluid and the solid phases.

1.4 Depth-averaged dimensional model equations

Combining relation (102) with (103) and using (104)–(106) can yield the downslope component of upper-layer granular momentum equations. Similarly, combining relation (107) with (108) can yield the cross-slope component. They are given by

$$\begin{aligned}&\frac{\partial }{\partial t}( h_g{\widetilde{\rho }}_s\phi _g)+\frac{\partial }{\partial x}( h_g{\widetilde{\rho }}_s\phi _g{\overline{u}}_g)+\frac{\partial }{\partial y}( h_g{\widetilde{\rho }}_s\phi _g{{\overline{v}}}_g) =J, \end{aligned}$$

(123)

$$\begin{aligned}&\frac{\partial }{\partial t}(h_g{\widetilde{\rho }}_s\phi _g{\overline{u}}_g)+\frac{\partial }{\partial x}(h_g {\widetilde{\rho }}_s\phi _g{\overline{u}}_g^2)+\frac{\partial }{\partial y}(h_g{\widetilde{\rho }}_s\phi _g{\overline{u}}_g{{\overline{v}}}_g)+\frac{\partial }{\partial x}\left( \frac{1}{2}K_x^gh_g^2{\widetilde{\rho }}_s\phi _gg_z\right) \nonumber \\ &\quad =Ju_g^i -{\widetilde{\rho }}_sh_g\phi _gg_z\frac{\partial h_m}{\partial x} + s_{x(g)}, \end{aligned}$$

(124)

$$\begin{aligned}&\frac{\partial }{\partial t}(h_g{\widetilde{\rho }}_s\phi _g{{\overline{v}}}_g)+\frac{\partial }{\partial x}(h_g {\widetilde{\rho }}_s\phi _g{\overline{u}}_g{{\overline{v}}}_g)+\frac{\partial }{\partial y}(h_g{\widetilde{\rho }}_s\phi _g{{\overline{v}}}_g^2)+\frac{\partial}{\partial y}\left( \frac{1}{2}K_y^gh_g^2{\widetilde{\rho }}_s\phi _gg_z\right) \nonumber \\ &\quad =Jv_g^i -{\widetilde{\rho }}_sh_g\phi _gg_z\frac{\partial h_m}{\partial y} + s_{y(g)} \end{aligned}$$

(125)

together with mass-conservation equation, where the source terms $s_{x(g)}$ and $s_{y(g)}$ are

$$\begin{aligned}&s_{x(g)} = {\widetilde{\rho }}_sh_g\phi _gg_x + k_s^i\mu _s^i\phi _g(u^i_s-u^i_g) +\frac{u^i_s-u^i_g}{\mid {{\varvec{u}}}^ i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i {\widetilde{\rho }}_sh_g\phi _g\left( g_z+\kappa {\overline{u}}_g^2\right) , \end{aligned}$$

(126)

$$\begin{aligned}&s_{y(g)} = k_s^i\mu _s^i\phi _g(v^i_s-v^i_g) +\frac{v^i_s-v^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i {\widetilde{\rho }}_sh_g\phi _g\left( g_z+\kappa {\overline{u}}_g^2\right) . \end{aligned}$$

(127)

Similarly, combining relation (109) with (110) and using relations (113)–(115) can yield the downslope component of the lower-layer solid momentum equations. Combining relation (116) with (117) can yield the cross-slope component of the lower-layer solid momentum equations. They are given by

$$\begin{aligned}&\frac{\partial }{\partial t}( h_m{\overline{\rho }}_s)+\frac{\partial }{\partial x}( h_m{\overline{\rho }}_s{\overline{u}}_s)+\frac{\partial }{\partial y}( h_m{\overline{\rho }}_s{{\overline{v}}}_s) =-J, \end{aligned}$$

(128)

$$\begin{aligned}&\frac{\partial }{\partial t}(h_m{\overline{\rho }}_s{\overline{u}}_s)+\frac{\partial }{\partial x}(h_m{\overline{\rho }}_s{\overline{u}}_s^2)+\frac{\partial }{\partial y}(h_m{\overline{\rho }}_s{\overline{u}}_s{{\overline{v}}}_s) \nonumber \\ &\qquad \qquad +\frac{\partial }{\partial x}\left( \frac{1}{2}h_mh_g{\overline{\rho }}_sg_z(K_y^s+K_y^g)\right) +\frac{\partial }{\partial x}\left( \frac{1}{6}({\widetilde{\rho }}_s-{\widetilde{\rho }}_f)g_zh_m^2{\overline{\phi }}_s (2K_y^s+K_y^g)\right) \nonumber \\ &\qquad \qquad \qquad \qquad \qquad =-Ju_s^i-\frac{\partial }{\partial x}\left( \frac{1}{2}{\widetilde{\rho }}_f{\overline{\phi }}_sh_m^2g_z \right) +\frac{1}{2}{\widetilde{\rho }}_fh_m^2g_z\frac{\partial {\overline{\phi }}_s}{\partial x} +{\widetilde{\rho }}_sh_g\phi _gg_z\frac{\partial h_m}{\partial x} + s_{x(s)}, \end{aligned}$$

(129)

$$\begin{aligned}&\frac{\partial }{\partial t}(h_m{\overline{\rho }}_s{{\overline{v}}}_s)+\frac{\partial }{\partial x}(h_m {\overline{\rho }}_s{\overline{u}}_s{{\overline{v}}}_s)+\frac{\partial }{\partial y}(h_m{\overline{\rho }}_s{{\overline{v}}}_s^2) \nonumber \\ &\qquad \qquad +\frac{\partial }{\partial y}\left( \frac{1}{2}h_mh_g{\overline{\rho }}_sg_z(K_y^s+K_y^g)\right) +\frac{\partial }{\partial y}\left( \frac{1}{6}({\widetilde{\rho }}_s-{\widetilde{\rho }}_f)g_zh_m^2{\overline{\phi }}_s (2K_y^s+K_y^g)\right) \nonumber \\ &\qquad \qquad \qquad \qquad \qquad =-Jv_s^i-\frac{\partial }{\partial y}\left( \frac{1}{2}{\widetilde{\rho }}_f{\overline{\phi }}_sh_m^2g_z \right) +\frac{1}{2}{\widetilde{\rho }}_fh_m^2g_z\frac{\partial {\overline{\phi }}_s}{\partial y} +{\widetilde{\rho }}_sh_g\phi _gg_z\frac{\partial h_m}{\partial y} + s_{y(s)}, \end{aligned}$$

(130)

together with mass equation, where the source terms $s_{x(s)}$ and $s_{y(s)}$ are given by

$$\begin{aligned} s_{x(s)} &={h}_m{\overline{\rho }}_sg_x +C_d{h}_m{\overline{\phi }}_s{\overline{\phi }}_f({\overline{u}}_f-{\overline{u}}_s) -k_s^b\mu _s^b{\phi _{s}}^b\,u_s^b\nonumber \\ &-\frac{u^b_s}{\mid {{\varvec{u}}}^b_s\mid }\mu _s^b\left[ g_z\left( {h}_m{\overline{\phi }}_s({\widetilde{\rho }}_s -{\widetilde{\rho }}_f) +{\widetilde{\rho }}_s h_g\phi _g\right) +\kappa \left( {h}_m{\overline{\phi }}_s({\widetilde{\rho }}_s{\overline{u}}_s^2 -{\widetilde{\rho }}_f{\overline{u}}_f^2) + {\widetilde{\rho }}_sh_g\phi _g{\overline{u}}_g^2\right) \right] \nonumber \\ &-\frac{u^i_s-u^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^i h_g\phi _g{\widetilde{\rho }}_s\left( g_z+ \kappa {\overline{u}}_g^2\right) -k_s^i\mu _s^i\phi _g(u_s^i-u_g^i) , \end{aligned}$$

(131)

$$\begin{aligned} s_{y(s)} &=C_d{h}_m{\overline{\phi }}_s{\overline{\phi }}_f({{\overline{v}}}_f-{{\overline{v}}}_s) -k_s^b\mu _s^b{\phi _{s}}^b\,v_s^b \nonumber \\ &-\frac{v^b_s}{\mid {{\varvec{u}}}^b_s\mid }\mu _s^b\left[ g_z\left( {h}_m{\overline{\phi }}_s({\widetilde{\rho }}_s -{\widetilde{\rho }}_f) +{\widetilde{\rho }}_s h_g\phi _g\right) +\kappa \left( {h}_m{\overline{\phi }}_s({\widetilde{\rho }}_s{\overline{u}}_s^2 -{\widetilde{\rho }}_f{\overline{u}}_f^2) + {\widetilde{\rho }}_sh_g\phi _g{\overline{u}}_g^2\right) \right] \nonumber \\ &-\frac{v^i_s-v^i_g}{\mid {{\varvec{u}}}^i_s-{{\varvec{u}}}^i_g\mid }\mu _s^ih_g\phi _g{\widetilde{\rho }}_s\left( g_z+ \kappa {\overline{u}}_g^2\right) -k_s^i\mu _s^i\phi _g(v_s^i-v_g^i). \end{aligned}$$

(132)

Combining relation (118) with (119) and using relation (120) leads to the downslope component of the fluid momentum equations. Similarly, combining relation (121) with (122) can yield the cross-slope component. The mass-conservation and momentum-conservation equations are given by

$$\begin{aligned}&\frac{\partial }{\partial t}( h_m{\overline{\rho }}_f)+\frac{\partial }{\partial x}( h_m{\overline{\rho }}_f{\overline{u}}_f)+\frac{\partial }{\partial y}( h_m{\overline{\rho }}_f{{\overline{v}}}_f) =0, \end{aligned}$$

(133)

$$\begin{aligned}&\frac{\partial }{\partial t}(h_m{\overline{\rho }}_f{\overline{u}}_f)+\frac{\partial }{\partial x}(h_m{\overline{\rho }}_f{\overline{u}}_f^2)+\frac{\partial }{\partial y}(h_m{\overline{\rho }}_f{\overline{u}}_f{{\overline{v}}}_f) \nonumber \\ &\quad +\frac{\partial }{\partial x}\left( \frac{1}{2}{\overline{\rho }}_fg_zh_m^2\right) =\frac{1}{2}h_m^2g_z\frac{\partial {\overline{\rho }}_f}{\partial x} + s_{x(f)}, \end{aligned}$$

(134)

$$\begin{aligned}&\frac{\partial }{\partial t}(h_m{\overline{\rho }}_f{\overline{u}}_f)+\frac{\partial }{\partial x}(h_m{\overline{\rho }}_f{\overline{u}}_f{{\overline{v}}}_f)+\frac{\partial }{\partial y}(h_m{\overline{\rho }}_f{{\overline{v}}}_f^2) \nonumber \\ &\quad + \frac{\partial }{\partial y}\left( \frac{1}{2}{\overline{\rho }}_fg_zh_m^2\right) =\frac{1}{2}h_m^2g_z\frac{\partial {\overline{\rho }}_f}{\partial y} + s_{y(f)}. \end{aligned}$$

(135)

Appendix 3: Fluxes and non-conservative products

To derive the fluxes, we need to explicitly express the mass-exchange rate. To this end, we eliminate the temporal derivatives from (56) by using mass equations (22), (29) and (37), so that we obtain

$$\begin{aligned} J&= -c_m\left[ h_m\left( \frac{\partial u_g}{\partial x}+\frac{\partial v_g}{\partial y}\right) -\frac{\partial }{\partial x}(h_m{\phi _{s}}u_s+h_m\phi _fu_f)-\frac{\partial }{\partial y} (h_m{\phi _{s}}v_s+h_m\phi _fv_f)\right] \nonumber \\ &- \frac{c(h_m+h_g)\phi _g}{h_g(1-\phi _g)}(\phi _m-{\phi _{s}}), \end{aligned}$$

(136)

where $c_m=\phi _g/(1-\phi _g)$. Furthermore, the vector system (57) is a compact form of the model equations. For this system, the fluxes can be determined directly from the left-hand sides of Eqs. (22), (23), (29), (30), (37) and (38), which yields

$$\begin{aligned}&\displaystyle {\varvec{F}}= \left( m_g^x-c_m(m^x_{s}+m_f^x), \,\,\,\, m^x_{s}+c_m(m^x_{s}+m_f^x), \,\,\,\, m_f^x, \,\,\,\, (m_g^x)^2/ {\hat{h}}_g+\frac{1}{2}\beta ^x_g {\hat{h}}_g^2, \right. \nonumber \\ &\quad \quad \, \left. (m^x_{s})^2/ {\hat{h}}_s +\frac{1}{2}\beta ^{x_1}_{s} {\hat{h}}_g({\hat{h}}_s+{\hat{h}}_f)+\frac{1}{2}\beta ^{x_2}_{s}{\hat{h}}_s({\hat{h}}_s+{\hat{h}}_f), \,\,\,\, (m_f^x)^2/{\hat{h}}_f+\frac{1}{2}\beta ^x_f{\hat{h}}_f({\hat{h}}_s+{\hat{h}}_f) \right) ^T. \end{aligned}$$

(137)

Eq. (137) indicates a flux Jacobian

$$\begin{aligned} \frac{\partial {\varvec{F}}}{\partial \varvec{\Omega }} = \left[ {\begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 1 &{} -c_m &{} -c_m\\ 0 &{} 0 &{} 0 &{} 0 &{} 1+c_m &{} c_m\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1\\ \displaystyle -\frac{(m_g^x)^2}{{\hat{h}}_g^2}+\beta _g^x{\hat{h}}_g &{} 0 &{} 0 &{} \displaystyle 2\frac{m_g^x}{{\hat{h}}_g} &{} 0 &{} 0 \\ \displaystyle \frac{\beta ^{x_1}_s({\hat{h}}_s+{\hat{h}}_f)}{2} &{} B^*&{} C^*&{} 0 &{} \displaystyle 2\frac{m_s^x}{{\hat{h}}_s} &{} 0\\ 0 &{} \displaystyle \frac{\beta _f^x{\hat{h}}_f}{2} &{} \displaystyle -\frac{(m_f^x)^2}{{\hat{h}}_f^2}+\frac{\beta _f^x(2{\hat{h}}_f+{\hat{h}}_s)}{2}&{} 0 &{} 0 &{} \displaystyle 2\frac{m_f^x}{{\hat{h}}_f} \end{array} } \right] , \, \end{aligned}$$

(138)

where $B^*= {-(m_s^x)^2}/{{\hat{h}}_s^2}+{\beta _s^{x_1}{\hat{h}}_g/2+\beta _s^{x_2}(2{\hat{h}}_s+{\hat{h}}_f)}/2$ and $C^*= \beta _s^{x_1}{\hat{h}}_g/2 + \beta _s^{x_2}{\hat{h}}_s/2$.

Inspect the non-conservative terms of Eqs. (22), (23), (29), (30), (37) and (38) can yield the form of ${\varvec{W}}_x$. It should be noted that in the process of formulating non-conservative terms of Eq.(30), we combine the terms $-\partial (\epsilon \gamma {\overline{\phi }}_sh_m^2\cos \zeta /2)/\partial x$ and $(\epsilon \gamma h_m^2\cos \zeta /2)\partial {\overline{\phi }}_s/\partial x$ on the right-hand side of (30) as $-\epsilon \gamma {\overline{\phi }}_s\partial (h_m^2\cos \zeta /2)/\partial x$, which can be further approximated as $-(\epsilon \gamma {\overline{\phi }}_sh_m\cos \zeta )\partial h_m/\partial x$ due to $\partial \zeta /\partial x \approx 0$ justified by $\lambda<< 1$. With the above instructions, the matrix-valued function ${\varvec{W}}_x$ can be directly formulated as

$$\begin{aligned} {{\varvec{W}}}_x= \left[ {\begin{array}{cccccc} \displaystyle -\frac{c_mh_mu_g}{{\hat{h}}_g} &{} 0 &{} 0 &{} \displaystyle \frac{c_mh_m}{{\hat{h}}_g} &{} 0 &{} 0 \\ \displaystyle \frac{c_mh_mu_g}{{\hat{h}}_g} &{} 0 &{} 0 &{} \displaystyle -\frac{c_mh_m}{{\hat{h}}_g} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ \displaystyle -\frac{c_mh_mu_su_g}{{\hat{h}}_g} &{} \epsilon {\hat{h}}_g\cos \zeta &{} \epsilon {\hat{h}}_g\cos \zeta &{} \displaystyle \frac{c_m \displaystyle h_mu_s}{{\hat{h}}_g} &{} \displaystyle -c_mu_s &{} \displaystyle -c_mu_s\\ \displaystyle \frac{c_mh_mu_su_g}{{\hat{h}}_g}&{} -\epsilon ({\hat{h}}_g-\gamma {\hat{h}}_s)\cos \zeta &{} -\epsilon ({\hat{h}}_g-\gamma {\hat{h}}_s)\cos \zeta &{} \displaystyle -\frac{c_mh_mu_s}{{\hat{h}}_g} &{} c_mu_s &{} \displaystyle c_mu_s\\ \displaystyle 0 &{} \displaystyle \frac{\epsilon {\hat{h}}_f\cos \zeta }{2} &{} \displaystyle -\frac{\epsilon {\hat{h}}_s\cos \zeta }{2} &{} 0 &{} 0 &{} 0 \\ \end{array} } \right] , \, \end{aligned}$$

(139)

where $h_m={\hat{h}}_s+{\hat{h}}_f$, and $u_g=m_g^x/{\hat{h}}_g$ and $u_s=m_s^x/{\hat{h}}_s$.

To formulate the source terms we extract the non-derivative term ${\mathcal {G}}=-c\phi _g(h_m+h_g)(\phi _m-{\phi _{s}})/(h_g(1-\phi _g))$ from mass-exchange rate (136), so that the vector of source terms can be written as

$$\begin{aligned} {\varvec{S}} &=\left( {\mathcal {G}},\,\,\, -{\mathcal {G}},\,\,\, 0, \,\,\, {\mathcal {G}}m_s^x/{\hat{h}}_s + s_{x(g)}, \,\,\, -{\mathcal {G}}m_s^x/{\hat{h}}_s + s_{x(s)}, \,\,\, s_{x(f)} \right) ^T. \end{aligned}$$

(140)