2.1. Overview of the microscopic and macroscopic descriptions

As in the two-phase case, three-phase mixtures can also be described as flows composed of regions occupied by single phases separated from each other by mobile surfaces. This description, commonly defined as ‘microscopic’ or discrete, characterizes the flow inside each region and the motion of each interface. The balance principles inside a region are obtained by applying the conservation equations to an infinitesimal control volume that is fully contained in the phase. Conversely, the equations valid across an interface, namely a surface that separates a given phase from another one, are obtained from an infinitesimal control volume that straddles the interface itself (Ishii & Hibiki Reference Ishii and Hibiki2006). These equations are called jump conditions. This approach can be used when the number of interfaces is limited, but in the case of a large number of interfaces (as in the case we are interested in), its application presents overwhelming difficulties.

The approach can be simplified by describing the flow of each phase at the macroscopic level. This description, also called the ‘multiphase continuum approach’, has the ability to filter the details of the flow, thus maintaining only its essential aspects. It can be derived by applying, to the microscopic equations, suitable averages that satisfy specific properties, as specified by Drew (Reference Drew1983). Time, space and ensemble averages are commonly employed by many authors (Anderson & Jackson Reference Anderson and Jackson1967; Drew & Segel Reference Drew and Segel1971; Morland Reference Morland1992; Ishii & Hibiki Reference Ishii and Hibiki2006). As a consequence, the averaging procedure leads to a formal coexistence of all the phases in any point of the field. With respect to the original equations, these macroscopic balances contain some new terms that arise as a contribution of the interfaces to the average. Conversely, with reference to the jump conditions, the averaging procedure leads to relations, called transfer conditions, that are formally equivalent to the original ones.

With regard to the mathematical symbolism used in this work, while in the microscopic description the instantaneous and local variables are denoted by the prime symbol, in the macroscopic description the averaged variables are not denoted by any superscript. In addition, the subscript $k$, with $k=1,2,3$, is used to identify one of the different phases that compose the mixture, while the double subscript $kj$ identifies quantities related to the interfaces between the phases $k$ and $j$. Unless otherwise indicated, the double subscript holds $\forall (k,j)$ with $k,j\in \lbrack 1, 2, 3]$ and $k\neq j$. Finally, the dependence on space and time of the functions will not be explicitly indicated.

2.2. Microscopic description

The equations of the microscopic description of a three-phase mixture are reported hereafter. While the equations valid inside a phase are formulated in the same way as in the two-phase case, more attention must be paid to the jump conditions. In fact, in three-phase solid–liquid mixtures the interfaces are not of one type but of two types, according to the state of the phases that they separate. More precisely, the interfaces can be defined ‘solid–solid’ if they separate two distinct solid phases and ‘solid–liquid’ if they separate a solid phase from the liquid one. This aspect requires a more detailed treatment of the jump conditions with respect to the two-phase solid–liquid case.

2.2.1. Equations valid inside each phase

The mass and momentum PDEs valid inside each region of the $k$th phase (Truesdell & Toupin Reference Truesdell and Toupin1960) are

(2.1) \begin{equation} \left. \begin{gathered} \dfrac{\partial}{\partial t}\rho_{k}^{\prime}+\boldsymbol{\nabla\cdot} (\rho_{k}^{\prime}\boldsymbol{v}_{k}^{\prime})=0 \\ \dfrac{\partial}{\partial t}(\rho_{k}^{\prime}\boldsymbol{v}_{k}^{\prime })+\boldsymbol{\nabla\cdot}(\rho_{k}^{\prime}\boldsymbol{v}_{k}^{\prime }\otimes\boldsymbol{v}_{k}^{\prime})=\boldsymbol{\nabla\cdot}{\boldsymbol{\mathsf{T}}}_{k} ^{\prime}+\rho_{k}^{\prime}\boldsymbol{g} \end{gathered} \right\},\end{equation}

where $\rho _{k}^{\prime },\boldsymbol {v}_{k}^{\prime }$ and ${\boldsymbol{\mathsf{T}}}_{k} ^{\prime }$ are, respectively, the density, velocity and stress tensor related to the $k$th phase, while $\boldsymbol {g}$ is the gravity vector.

2.2.2. Jump conditions

According to Ishii & Hibiki (Reference Ishii and Hibiki2006), the mass jump condition valid between phase $k$ and phase $j$ reads

(2.2)\begin{equation} {-}\rho_{k}^{\prime}(\boldsymbol{v}_{k}^{\prime}-\boldsymbol{v}_{kj}^{\prime })\boldsymbol{\cdot}\boldsymbol{N}_{kj}+\rho_{j}^{\prime}(\boldsymbol{v} _{j}^{\prime}-\boldsymbol{v}_{kj}^{\prime})\boldsymbol{\cdot}\boldsymbol{N} _{kj}=0 ,\end{equation}

where $\boldsymbol {v}_{kj}^{\prime }$ represents the velocity vector of the interface separating phase $k$ from phase $j$, while $\boldsymbol {N}_{kj}$ identifies the normal unit vector in a point of this interface, pointing outside phase $k$ towards phase $j$. Each term expresses a mass flux that is associated with a phase transformation.

The momentum jump condition (Ishii & Hibiki Reference Ishii and Hibiki2006) is defined as follows:

(2.3)

where $\sigma ^{\prime }$ and ${K}$ represent the surface tension and the surface curvature, respectively. There appear two different contributions: one referred to the momentum flux associated with the mass exchange and the other related to the stresses exerted by the two phases on the interface.

2.3. Macroscopic description

The macroscopic equations reported hereafter are derived by considering, as a reference average operator, a spherical volume $\mathcal {V}$ with a radius larger than the solid particles. In this description, the fraction of reference volume occupied by the generic phase $k$ is defined by $\alpha _k$. Since all the space occupied by the mixture is completely filled by the phases, the following equality holds:

(2.4)\begin{equation} \sum_{k=1}^{3}\alpha_{k}=1. \end{equation}

2.3.1. Gradients of the volume fractions

Following Drew (Reference Drew1983), the gradient of a volume fraction represents the oriented surface (per unit volume) that arises by averaging the local contact areas of a phase with the others (see Appendix A). With reference to the liquid phase, here identified by the subscript $k=3$, the quantity $\boldsymbol {\nabla }\alpha _3$ identifies the oriented average interfacial surface that separates the liquid phase from the two solid phases. Thus, it can be split into two terms

(2.5)\begin{equation} \boldsymbol{\nabla}\alpha_3=[\boldsymbol{\nabla}\alpha_3]_1+ [\boldsymbol{\nabla}\alpha_3]_2,\end{equation}

where $[\boldsymbol {\nabla }\alpha _3]_k$ with $k=1, 2$ represents the contribution relevant to each solid phase. With regard to the gradients of the two solid phases, the corresponding oriented average interfacial surfaces need to take account of both the solid–liquid and solid–solid contact areas. In this way, they can be expressed as follows:

(2.6a)$$\begin{gather} \boldsymbol{\nabla}\alpha_1=[\boldsymbol{\nabla}\alpha_1]_3+ [\boldsymbol{\nabla}\alpha_1]_2, \end{gather}$$

(2.6b)$$\begin{gather}\boldsymbol{\nabla}\alpha_2=[\boldsymbol{\nabla}\alpha_2]_3+ [\boldsymbol{\nabla}\alpha_2]_1, \end{gather}$$

where $[\boldsymbol {\nabla }\alpha _k]_3$ with $k=1, 2$ represents the oriented average interfacial surface that separates the solid phase $k$ from the liquid one, while $[\boldsymbol {\nabla }\alpha _1]_2$ and $[\boldsymbol {\nabla }\alpha _2]_1$ correspond to the oriented average interfacial surfaces separating the two solid phases.

Considering that, in a point of an interface, the surface that separates two phases is unique, the oriented surfaces relevant to these phases are vectors equal in absolute value and opposite in direction. This relation holds also in terms of averaged variables (see Appendix A) and, therefore, the following equalities hold:

(2.7) \begin{gather} {[}\boldsymbol{\nabla}\alpha_3]_k={-}[\boldsymbol{\nabla}\alpha_k]_3 \quad\text{with }k=1, 2, \end{gather}

(2.8) \begin{gather}{[}\boldsymbol{\nabla}\alpha_2]_1={-}[\boldsymbol{\nabla}\alpha_1]_2. \end{gather}

The quantity $\boldsymbol {\nabla }\alpha _3$ can be linked to the gradients of the volume fraction of the other phases by applying the gradient operator to (2.4), thus obtaining the following equation:

(2.9)\begin{equation} \boldsymbol{\nabla}\alpha_3={-}\boldsymbol{\nabla}\alpha_1- \boldsymbol{\nabla}\alpha_2.\end{equation}

It is worth noting that adding (2.6a) to (2.6b) and applying (2.7) and (2.8) to the resulting equation, it is possible to derive an equation that coincides with (2.9).

Finally, assuming that the averaged contact area between the two solid phases is much smaller than the averaged surface between a solid phase and the liquid one, the relation (2.6) can be simplified as follows:

(2.10)\begin{equation} \boldsymbol{\nabla}\alpha_k\simeq[\boldsymbol{\nabla}\alpha_k]_3 \quad\text{with } k=1, 2.\end{equation}

It is worthy of note that this assumption may not occur in all situations, such as when the liquid concentrations are low and the mixture behaves like a solid (Takahashi Reference Takahashi2007; Schneider et al. Reference Schneider, Huggel, Haeberli and Kaitna2011; Andreotti, Forterre & Pouliquen Reference Andreotti, Forterre and Pouliquen2013; Pudasaini & Krautblatter Reference Pudasaini and Krautblatter2014). For more details on the mathematical proof of (2.5)–(2.10), we refer the reader to Appendix A.

2.3.2. Mass balance equations

The averaged mass balance equation for the $k$th phase can be derived from the first equation of system (2.1) and reads

(2.11)\begin{equation} \frac{\partial}{\partial t}( \alpha_{k}\rho_{k}) +\boldsymbol{\nabla\cdot}( \alpha_{k}\rho_{k}\boldsymbol{u}_{k}) =\varGamma_{k}, \end{equation}

where $\rho _{k}$ and $\boldsymbol {u}_{k}$ indicate the averaged density and velocity of the $k$th phase, respectively, while $\varGamma _{k}$ expresses the global mass transfer (per unit volume) between the $k$th phase and all the other components of the mixture. It can be subdivided as

(2.12)\begin{equation} \varGamma_{k}=\sum_{j=1}^{3}(1-\delta_{kj})\varGamma_{kj},\end{equation}

where $\delta _{kj}$ identifies the Kronecker delta function and allows us to delete the meaningless term $\varGamma _{kk}$. We define $\varGamma _{kj}$ as positive whether the transfer of mass occurs from phase $j$ to phase $k$. Since it cannot be written automatically in terms of averaged variables, a proper closure relation needs to be provided.

2.3.3. Momentum balance equations

The momentum balance equation for the $k$th phase can be derived from the second equation of system (2.1) and reads

(2.13)\begin{equation} \frac{\partial}{\partial t}( \alpha_{k}\rho_{k}\boldsymbol{u} _{k}) +\boldsymbol{\nabla\cdot}( \alpha_{k}\rho_{k} \boldsymbol{u}_{k}\otimes\boldsymbol{u}_{k}) =\boldsymbol{\nabla\cdot }( \alpha_{k}{\boldsymbol{\mathsf{T}}}_{k}) +\alpha_{k}\rho_{k} \boldsymbol{g} +\boldsymbol{M}_{k}.\end{equation}

In this equation, ${\boldsymbol{\mathsf{T}}}_{k}$ is the average intraphase tensor that is viscous in the case of a fluid phase and frictional/collisional in the case of a solid phase. In both cases, proper closure relations need to be introduced for this term. In addition, in (2.13) $\boldsymbol {M}_{k}$ is the resultant of the interphase stresses exerted on the $k$th phase by the other components of the mixture

(2.14)\begin{equation} \boldsymbol{M}_{k}=\sum_{j=1}^{3}(1-\delta_{kj})\boldsymbol{M}_{kj}.\end{equation}

Due to the presence of solid–solid and solid–liquid interfaces, the interphase stress $\boldsymbol {M}_{kj}$ needs to be specified according to the type of interfaces that separate the phases $k$ and $j$.

Solid–liquid interphase stress. In agreement with Drew (Reference Drew1983), Jackson (Reference Jackson2000) and Ishii & Hibiki (Reference Ishii and Hibiki2006), the solid–liquid interphase stress can be split into a momentum transfer and into static and dynamic stresses exerted on the interfaces,

(2.15)\begin{equation} \boldsymbol{M}_{kj}=\varGamma_{kj}\boldsymbol{u}_{kj}- [\boldsymbol{\nabla}\alpha_{k}]_{j}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{kj}+\boldsymbol{M}_{kj}^{D}.\end{equation}

On the right-hand side of (2.15), the first term, namely $\varGamma _{kj}\boldsymbol {u}_{kj}$, identifies the momentum transfer between phase $k$ and phase $j$. The variable $\boldsymbol {u}_{kj}$ represents the averaged interfacial velocity and needs to be defined by proper (constitutive) assumptions. The second term namely $-[\boldsymbol {\nabla }\alpha _{k}]_{j}\boldsymbol {\cdot }{\boldsymbol{\mathsf{T}}}_{kj}$, represents the static contribution of the averaged stress exerted at the interfaces by phase $j$ on phase $k$. The third term, namely $\boldsymbol {M}_{kj}^{D}$, identifies the dynamic contribution of the averaged stress exerted on the interfaces and it arises from the relative motion of the solid and liquid phases. It takes account of the effects of drag, virtual mass and lift (Ishii & Hibiki Reference Ishii and Hibiki2006). For brevity, in the following we call this term as the drag stress.

Both ${\boldsymbol{\mathsf{T}}}_{kj}$ and $\boldsymbol {M}_{kj}^{D}$ need to be specified by proper closure relations.

Finally, thanks to (2.10) and to (2.7), the solid–liquid interphase stresses can be rewritten as follows:

(2.16) \begin{equation} \boldsymbol{M}_{kj}=\left\{ \begin{array}{@{}ll} \varGamma_{kj}\boldsymbol{u}_{kj}-\boldsymbol{\nabla}\alpha_{k}\boldsymbol{\cdot }{\boldsymbol{\mathsf{T}}}_{kj}+\boldsymbol{M}_{kj}^{D} & \text{if } k \text{ identifies a solid phase},\\ \varGamma_{kj}\boldsymbol{u}_{kj}+\boldsymbol{\nabla}\alpha_{j}\boldsymbol{\cdot }{\boldsymbol{\mathsf{T}}}_{kj}+\boldsymbol{M}_{kj}^{D} & \text{if }k \text{ identifies the liquid phase}.\end{array} \right. \end{equation}

Solid–solid interphase stress. Similarly to the solid–liquid case, the interphase stress can be split into a momentum transfer and an averaged stress exerted on the interfaces. Conversely, while the splitting of the averaged stress into a static and a dynamic effect is quite accepted for dry granular flows (Gray & Thornton Reference Gray and Thornton2005; Gray & Chugunov Reference Gray and Chugunov2006; Tripathi & Khakhar Reference Tripathi and Khakhar2013; Tunuguntla, Bokhove & Thornton Reference Tunuguntla, Bokhove and Thornton2014), it does not seem to be the case for fully saturated flows. To the best of our knowledge, the only model for granular–liquid flows that considers the splitting of the solid–solid averaged stress into static and dynamic effects corresponds to the approach proposed by Pudasaini & Mergili (Reference Pudasaini and Mergili2019). However, we prefer to express this averaged stress in terms of frictional and collisional actions by analogy with the solid intraphase tensor (Jenkins & Savage Reference Jenkins and Savage1983; Armanini et al. Reference Armanini, Larcher, Nucci and Dumbser2014). We do not discuss further the differences between these two approaches, since this is outside the scope of this work, as explained further on in the paper. In this way, for the solid–solid interphase stresses we get

(2.17)\begin{equation} \boldsymbol{M}_{kj}=\varGamma_{kj}\boldsymbol{u}_{kj}+\boldsymbol{M}_{kj} ^{F}+\boldsymbol{M}_{kj}^{C}, \end{equation}

where $\varGamma _{kj}\boldsymbol {u}_{kj}$ represents the momentum transfer; $\boldsymbol {M}_{kj}^{F}$ identifies the action due to the long-term frictional contacts between particles of different phases; $\boldsymbol {M}_{kj}^{C}$ corresponds to the action due to short-term collisions between particles of different phases.

Both $\boldsymbol {M}_{kj}^{F}$ and $\boldsymbol {M}_{kj}^{C}$ need to be specified by proper closure relations.

2.3.4. Mass transfer conditions

The mass transfer conditions are derived averaging (2.2) and read as follows:

(2.18)\begin{equation} \varGamma_{kj}+\varGamma_{jk}=0. \end{equation}

Adding together the mass transfer conditions related to all the phases, it is possible to derive the following expression:

(2.19)\begin{equation} \sum_{k=1}^{3}\sum_{j=1}^{3}(1-\delta_{kj})\varGamma_{kj}=0, \end{equation}

which, thanks to (2.12), becomes

(2.20)\begin{equation} \sum_{k=1}^{3}\varGamma_{k}=0. \end{equation}

2.3.5. Momentum transfer conditions

The momentum transfer conditions can be specified in two ways, whether the surface tension is assumed to be relevant or not.

Presence of surface tension. The momentum transfer conditions can be expressed by

(2.21)\begin{equation} \boldsymbol{M}_{kj}+\boldsymbol{M}_{jk}=[ \sigma\kappa\boldsymbol{n} ]_{kj}, \end{equation}

where $\kappa$ represents the averaged curvature and $\boldsymbol {n}$ the averaged normal unit vector. It is worth noting that the effect of the surface tension $\sigma$ can be considered relevant in the case of mixtures consisting of solid particles submerged in a fluid. In this situation, the surface curvature $\kappa$ can be linked to the particle diameter. These transfer conditions can be expressed by using (2.16) as follows:

(2.22)\begin{equation} (\varGamma_{kj}\boldsymbol{u}_{kj}+\varGamma_{jk}\boldsymbol{u}_{jk} )+\boldsymbol{\nabla}\alpha_{k}\boldsymbol{\cdot}({\boldsymbol{\mathsf{T}}}_{jk} -{\boldsymbol{\mathsf{T}}}_{kj})+(\boldsymbol{M}_{kj}^{D}+\boldsymbol{M}_{jk}^{D} )=[\sigma\kappa\boldsymbol{n}]_{kj}, \end{equation}

where $k$ identifies a solid phase and $j$ the liquid phase.

Absence of surface tension. The momentum transfer conditions are expressed by

(2.23)\begin{equation} \boldsymbol{M}_{kj}+\boldsymbol{M}_{jk}=0. \end{equation}

These relations can be considered as representative of the interphase stresses between two solid phases. By using (2.17), they become

(2.24)\begin{equation} \varGamma_{kj}\boldsymbol{u}_{kj}+\boldsymbol{M}_{kj}^{F}+\boldsymbol{M}_{kj} ^{C}={-}(\varGamma_{jk}\boldsymbol{u}_{jk}+\boldsymbol{M}_{jk}^{F}+\boldsymbol{M} _{jk}^{C}). \end{equation}

3.1. Interphase terms

With regard to the mass transfers, in RIW mixtures they can occur only between ice and water. Therefore, the mass transfers related to these phases differ from zero,

(3.2)\begin{equation} \varGamma_{iw}\neq0\cup\varGamma_{wi}\neq 0, \end{equation}

while all the other terms are null. Thus, considering (2.12) we can write

(3.3a)$$\begin{gather} \varGamma_{r} =\varGamma_{ri}+\varGamma_{rw}\Longrightarrow\varGamma_{r}=0, \end{gather}$$

(3.3b)$$\begin{gather}\varGamma_{i} =\varGamma_{ir}+\varGamma_{iw}\Longrightarrow\varGamma_{i}=\varGamma_{iw}, \end{gather}$$

(3.3c)$$\begin{gather}\varGamma_{w} =\varGamma_{wi}+\varGamma_{wr}\Longrightarrow\varGamma_{w}=\varGamma_{wi}. \end{gather}$$

With regard to the interphase stresses, we can distinguish between: (i) solid–liquid interactions that are described by (2.16) and occur between water and rock and between water and ice; and (ii) solid–solid interactions that are described by (2.17) and occur between rock and ice.

By using (3.3), these terms become

(3.4a)$$\begin{gather} \boldsymbol{M}_{r}=\boldsymbol{M}_{ri}+\boldsymbol{M}_{rw}\quad \text{where } \begin{cases} \boldsymbol{M}_{ri} =\boldsymbol{M}_{ri}^{C}+\boldsymbol{M}_{ri}^{F}\\ \boldsymbol{M}_{rw} =\boldsymbol{M}_{rw}^{D}-\boldsymbol{\nabla}\alpha _{r}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{ri} \end{cases}, \end{gather}$$

(3.4b)$$\begin{gather}\boldsymbol{M}_{i}=\boldsymbol{M}_{ir}+\boldsymbol{M}_{iw}\quad \text{where } \begin{cases} \boldsymbol{M}_{ir} =\boldsymbol{M}_{ir}^{C}+\boldsymbol{M}_{ir}^{F}\\ \boldsymbol{M}_{iw} =\boldsymbol{M}_{iw}^{D}-\boldsymbol{\nabla}\alpha _{i}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{iw}+\varGamma_{i}\boldsymbol{u}_{iw} \end{cases}, \end{gather}$$

(3.4c)$$\begin{gather}\boldsymbol{M}_{w}=\boldsymbol{M}_{wr}+\boldsymbol{M}_{wi}\quad \text{where } \begin{cases} \boldsymbol{M}_{wr} =\boldsymbol{M}_{wr}^{D}+\boldsymbol{\nabla}\alpha _{r}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{wr}\\ \boldsymbol{M}_{wi} =\boldsymbol{M}_{wi}^{D}+\boldsymbol{\nabla}\alpha _{i}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{wi}+\varGamma_{w}\boldsymbol{u}_{wi} \end{cases}. \end{gather}$$

3.2. Transfer conditions

With reference to the mass transfer condition, the existence of only two phases that exchange mass reduces (2.18) to the following relation:

(3.5)\begin{equation} \varGamma_{wi}={-}\varGamma_{iw} ,\end{equation}

and, therefore, thanks to (3.3b) and (3.3c), it becomes

(3.6)\begin{equation} \varGamma_{w}={-}\varGamma_{i}. \end{equation}

This condition states that the two mass transfers are equal to each other in absolute value and opposite in sign. The ice melting is responsible for a loss of ice mass and for an acquisition of water mass. Therefore, $\varGamma _{i}$ and $\varGamma _{w}$ are, respectively, negative and positive. To better highlight these features, in the following we set $\varGamma _{i}=-\vert \varGamma _{i}\vert$ and $\varGamma _{w}=+\vert \varGamma _{i}\vert$.

Concerning with the transfer conditions related to the momentum, the surface tension affects only the conditions involving water. Therefore, (2.21) and (2.23) become

(3.7a)$$\begin{gather} \boldsymbol{M}_{ri}+\boldsymbol{M}_{ir}=0, \end{gather}$$

(3.7b)$$\begin{gather}\boldsymbol{M}_{rw}+\boldsymbol{M}_{wr}=[\sigma\kappa\boldsymbol{n}]_{rw}, \end{gather}$$

(3.7c)$$\begin{gather}\boldsymbol{M}_{iw}+\boldsymbol{M}_{wi}=[\sigma\kappa\boldsymbol{n}]_{iw}. \end{gather}$$

In accordance with (3.4), we obtain

(3.8a)$$\begin{gather} \boldsymbol{M}_{ri}^{C}+\boldsymbol{M}_{ri}^{F}={-} (\boldsymbol{M}_{ir}^{C}+\boldsymbol{M}_{ir}^{F} ), \end{gather}$$

(3.8b)$$\begin{gather}( \boldsymbol{M}_{wr}^{D}+\boldsymbol{\nabla}\alpha_{r}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{wr} ) +(\boldsymbol{M}_{rw}^{D}-\boldsymbol{\nabla}\alpha_{r}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{rw} )=[\sigma\kappa\boldsymbol{n}]_{rw}, \end{gather}$$

(3.8c)$$\begin{gather}(\boldsymbol{M}_{wi}^{D}+\boldsymbol{\nabla}\alpha_{i}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{wi}+ \vert \varGamma_i \vert\boldsymbol{u}_{wi} ) + (\boldsymbol{M}_{iw}^{D}-\boldsymbol{\nabla}\alpha_{i} \boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{iw}- \vert \varGamma_i \vert\boldsymbol{u}_{iw} ) =[\sigma\kappa\boldsymbol{n}]_{iw} . \end{gather}$$

3.3. The final set

Introducing all the assumptions specified in the previous section, the final set of equations for the RIW model is composed of the macroscopic mass and momentum balance laws related to each phase (namely (2.11) and (2.13)), and of the basic relation (3.1). The system reads

(3.18a)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{{r}}\rho_{{r}} ) +\boldsymbol{\nabla\cdot} ( \alpha_{{r}}\rho_{{r}}\boldsymbol{u}_{{r} } ) = 0, \end{gather}

(3.18b)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{{i}}\rho_{{i}} ) +\boldsymbol{\nabla\cdot} ( \alpha_{{i}}\rho_{{i}}\boldsymbol{u}_{{i} } ) ={-} \vert \varGamma_{{i}} \vert, \end{gather}

(3.18c)\begin{gather}\frac{\partial}{\partial t} ( \alpha_{{w}}\rho_{{w}} ) +\boldsymbol{\nabla\cdot} ( \alpha_{{w}}\rho_{{w}}\boldsymbol{u}_{{w} } ) ={+} \vert \varGamma_{{i}} \vert, \end{gather}

(3.18d)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{{r}}\rho_{{r}}\boldsymbol{u}_{{r}} ) +\boldsymbol{\nabla\cdot} (\alpha_{{r}}\rho_{{r}}\boldsymbol{u}_{{r}}\otimes\boldsymbol{u}_{{r}} ) ={-}\boldsymbol{\nabla\cdot} ( \alpha_{{r}}{\boldsymbol{\mathsf{T}}}_{{r}} ) +\alpha_{{r}}\rho_{{r}}\boldsymbol{g}+\boldsymbol{\nabla}\alpha_{{r}} \boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{{w}}\nonumber\\ \hspace{8.3pc}\qquad\quad +\, \boldsymbol{M}_{{rw}}^{D}+\boldsymbol{M}_{{ri}}^{C}+\boldsymbol{M}_{{ri}}^{F}, \end{gather}

(3.18e)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{{i}}\rho_{{i}}\boldsymbol{u}_{{i}} ) +\boldsymbol{\nabla\cdot} ( \alpha_{{i}}\rho_{{i}}\boldsymbol{u}_{{i}}\otimes\boldsymbol{u}_{{i}} ) ={-}\boldsymbol{\nabla\cdot} ( \alpha_{{i}}{\boldsymbol{\mathsf{T}}}_{{i}} ) +\alpha_{{i}}\rho_{{i}}\boldsymbol{g}- \vert \varGamma_{{i}} \vert \boldsymbol{u}_{ss}+\boldsymbol{\nabla}\alpha_{{i}}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{{w}}\nonumber\\ \qquad\hspace{4.5pc}\,\, +\,\boldsymbol{M}_{{iw}}^{D}-\boldsymbol{M}_{{ri}}^{C}-\boldsymbol{M}_{{ri}}^{F}, \end{gather}

(3.18f)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{{w}}\rho_{{w}}\boldsymbol{u}_{{w}} ) +\boldsymbol{\nabla\cdot} ( \alpha_{{w}}\rho_{{w}}\boldsymbol{u}_{{w}}\otimes\boldsymbol{u}_{{w}} ) ={-}\boldsymbol{\nabla\cdot} (\alpha_{{w}}{\boldsymbol{\mathsf{T}}}_{{w}} )+\alpha_{{w}}\rho_{{w}}\boldsymbol{g}+ \vert \varGamma_{i} \vert \boldsymbol{u}_{{ss}}\nonumber\\ \hspace{14.5pc}+\boldsymbol{\nabla}\alpha_{{w}}\boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}_{{w}}-\boldsymbol{M}_{{iw}}^{D}-\boldsymbol{M}_{rw}^{D}, \end{gather}

(3.18g)\begin{gather} \alpha_{r}+\alpha_{i}+\alpha_{w} =1, \end{gather}

where the first three equations represent the mass balance equations related to rock, ice and water, respectively, while the next three equations identify the momentum balance equations for the three different phases (rock, ice and water, respectively). On the right-hand side of the mass and momentum balance equations $\pm \vert \varGamma _i\vert$ and $\pm \vert \varGamma _i\vert \boldsymbol {u}_{ss}$ correspond, respectively, to the mass and momentum transfers associated with the ice melting process; $\boldsymbol {\nabla \cdot }( \alpha _{{k}}{\boldsymbol{\mathsf{T}}}_{{k}})$ (with $k=r, i\,\text {and}\,w$) represents the internal stress term related to each phase; $\alpha _{{k}}\rho _{{k}}\boldsymbol {g}$ (with $k=r, i \text { and }w$) identifies the weight term; $\boldsymbol {\nabla }\alpha _{{k}}\boldsymbol {\cdot }{\boldsymbol{\mathsf{T}}}_{{w}}$ (with $k=r, i \text { and } w$) represents the static contribution of the solid–liquid interaction stress and corresponds to the buoyancy; $\boldsymbol {M}_{kw}^D$ (with $k=r\text { and }i$) takes account of the drag stress (drag, lift and virtual mass); $\boldsymbol {M}_{ri}^C$ and $\boldsymbol {M}_{ri}^F$ express the collisional and frictional stresses that arise from the contact between rock and ice.

It is worthy of note that the stress tensors of the three-phases in (3.18d)–(3.18f) have been multiplied by $-1$ to make positive the compressive stresses (Savage & Hutter Reference Savage and Hutter1989; Iverson & Denlinger Reference Iverson and Denlinger2001).

This system is composed of thirteen scalar equations and describes the flow in terms of thirteen unknowns. Twelve of them are present in the equations in an explicit way and correspond to the three phase concentrations $\alpha _{r},\alpha _{i},\alpha _{w}$ and the nine scalar velocity components of $\boldsymbol {u}_{{r}},\boldsymbol {u}_{{i}},\boldsymbol {u}_{{w}}$. The last unknown is the water pressure $p_{w}$, the scalar quantity that represents the isotropic part of the water stress tensor. All the other terms in the system must be referred to these thirteen unknowns through suitable closure relations (algebraic or differential).

It is worth noting that, by dividing (3.18a)–(3.18c) by the corresponding density, adding together the resulting equations and by using (3.18g), the following equation can be obtained:

(3.19)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} (\alpha_r\boldsymbol{u}_r+\alpha_i \boldsymbol{u}_i+\alpha_w\boldsymbol{u}_w )= \vert\varGamma_i \vert \left(\frac{1}{\rho_w}-\frac{1}{\rho_i} \right).\end{equation}

It expresses that the rate of the volumetric deformation of the mixture is connected with the mass transfer and with the change in density of the molten ice due to the phase transformation. Since $\rho _i<\rho _w$, the term on the right-hand side of (3.19) is negative. As a consequence, the mixture shows a decreasing volumetric deformation rate due to the ice melting.

3.4. The need for simplified models

The RIW model (3.18) is described by a complex system of equations, whose complexity arises from the high number of unknowns and from the interphase terms that couple the equations. Furthermore, the literature on the interaction stresses between rock and ice does not seem so well established. In fact, we think that there is some uncertainty surrounding both the closure relations and the way these interphase stresses are expressed in the balance principles. We suppose indeed that both melting and fragmentation of ice can occur during collisions and that the closure relations for the solid–solid interaction stresses should be affected by these processes. Moreover, we are not sure that the splitting into collisional ($\boldsymbol {M}_{ri}^{C}$) and frictional ($\boldsymbol {M}_{ri}^{F}$) stresses is the best way to represent the rock–ice interphase terms. Therefore, considering the current state of knowledge, the usage of simplified approaches seems to be the only reasonable option to tackle the RIW problem, both from a mathematical and a numerical point of view.

4.1. Partially isokinetic RIW model

The partially isokinetic RIW (PI-RIW) model can be obtained from the complete RIW approach by imposing the isokinetic condition between rock and ice,

(4.6)\begin{equation} \boldsymbol{u}_{r}=\boldsymbol{u}_{i}=\boldsymbol{u}_{s}, \end{equation}

where $\boldsymbol {u}_{s}$ is the velocity of the bulk ‘solid’ phase. From a physical point of view, this hypothesis means that some of the features that distinguish the two phases, such as the differences in grain sizes, are neglected. As pointed out in the introduction of this section, the bulk solid phase has a concentration and a density that are defined as

(4.7)$$\begin{gather} \alpha_{s}=\alpha_{r}+\alpha_{i}, \end{gather}$$

(4.8)$$\begin{gather}\rho_{s}=\frac{\alpha_{{r}}\rho_{{r}}+\alpha_{i}\rho_{i}}{\alpha_{s}}. \end{gather}$$

Differently from the densities of the RIW phases, $\rho _{s}$ is not constant in time and space.

For the derivation of the PI-RIW model, it is useful to define for the overall solid phase the drag stress vector $\boldsymbol {M} _{sw}^{D}$ and the internal stress tensor ${\boldsymbol{\mathsf{T}}}_{s}$ as follows:

(4.9)$$\begin{gather} \boldsymbol{M}_{sw}^{D}=\boldsymbol{M}_{rw}^{D}+\boldsymbol{M}_{iw}^{D}, \end{gather}$$

(4.10)$$\begin{gather}{\boldsymbol{\mathsf{T}}}_{s}=\frac{\alpha_{r}{\boldsymbol{\mathsf{T}}}_{r}+\alpha_{i} {\boldsymbol{\mathsf{T}}}_{i}}{\alpha_{s}}. \end{gather}$$

Assuming that both the phases have the same ‘solid’ mechanical behaviour, we do not distinguish rock from ice except for the phase transformation that affects only ice. Thus, the closure relations can be provided for these overall terms rather than for the stresses related to each phase.

In this class of models, the flow problem can be expressed as a function of $\boldsymbol {u}_{s}$, $\boldsymbol {u}_{w}$, $\alpha _{s}$, $\alpha _{w}$, $\rho _{s}$ and $p_{w}$. Therefore, the original unknowns $\alpha _{r}$ and $\alpha _{i}$ must be rewritten as functions of the new unknowns. By using (4.4a,b), these two concentrations can be specified as follows:

(4.11a,b)\begin{equation} \alpha_{i}=\alpha_{s}\frac{\rho_{s}-\rho_{r}}{\rho_{i}-\rho_{r}}; \quad \alpha_{r}=\alpha_{s}\frac{\rho_{s}-\rho_{i}}{\rho_{r} -\rho_{i}}.\end{equation}

Concerning the equations, the isokinetic condition allows us to replace the momentum balances related to rock and ice with the momentum equation for the bulk solid phase. This equation can be obtained by adding together the momentum equations for the rock and ice phases, namely (3.18d) and (3.18e). In addition, it is also useful to derive the mass balance for the bulk solid phase by summing the mass conservation equations related to rock and ice, namely (3.18a) and (3.18b), and to replace one of these original mass balances with the resulting equation. In this work, we have chosen to substitute the ice mass balance. Finally, in (3.18a) it is possible to replace $\alpha _{r}$ with the expression given in (4.11a,b), thus obtaining

(4.12)\begin{equation} \frac{\partial}{\partial t}\left( \rho_{{r}}\frac{(\rho_s\alpha_s-\rho_i\alpha_s)}{\rho_r-\rho_i}\right) +\boldsymbol{\nabla\cdot}\left( \rho_{{r}}\frac{(\rho_s\alpha_s- \rho_i\alpha_s)}{\rho_r-\rho_i}\boldsymbol{u}_{{s} }\right) = 0.\end{equation}

Since $\rho _r/(\rho _r-\rho _i)$ is constant in time and space due to the incompressibility assumption applied to rock and ice, we can neglect this term from (4.12). In this way, the resulting system reads as follows:

(4.13a)\begin{gather} \frac{\partial}{\partial t} ( \rho_s\alpha_s-\rho_i\alpha_s ) +\boldsymbol{\nabla\cdot} ( (\rho_s\alpha_s-\rho_i\alpha_s ) \boldsymbol{u}_{{s} } ) = 0, \end{gather}

(4.13b)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{s}\rho_{s} ) +\boldsymbol{\nabla\cdot} ( \alpha_{s}\rho_{s}\boldsymbol{u}_{s} ) ={-} \vert \varGamma_{i} \vert, \end{gather}

(4.13c)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{w}\rho_{w} ) +\boldsymbol{\nabla\cdot} ( \alpha_{w}\rho_{w}\boldsymbol{u}_{w} ) ={+} \vert \varGamma_{i} \vert, \end{gather}

(4.13d)\begin{align} \frac{\partial}{\partial t} ( \alpha_{s}\rho_{s}\boldsymbol{u}_{s} ) +\boldsymbol{\nabla\cdot} ( \alpha_{s}\rho_{s}\boldsymbol{u}_{s} \otimes\boldsymbol{u}_{s} ) &={-}\boldsymbol{\nabla\cdot} ( \alpha_{s}{\boldsymbol{\mathsf{T}}}_{s} ) +\alpha_{s}\rho _{s}\boldsymbol{g}- \vert \varGamma_{i} \vert \boldsymbol{u}_{ss}\nonumber\\ &\quad +\boldsymbol{\nabla}\alpha_{s} \boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{w}+\boldsymbol{M}_{sw}^{D}, \end{align}

(4.13e)\begin{align} \frac{\partial}{\partial t} ( \alpha_{w}\rho_{w}\boldsymbol{u} _{w} ) +\boldsymbol{\nabla\cdot} ( \alpha_{w}\rho_{w} \boldsymbol{u}_{w}\otimes\boldsymbol{u}_{w} ) &={-}\boldsymbol{\nabla\cdot} ( \alpha_{w}{\boldsymbol{\mathsf{T}}}_{w} )+\alpha_{w}\rho_{w}\boldsymbol{g} + \vert \varGamma_{i} \vert \boldsymbol{u}_{ss}\nonumber\\ &\quad +\boldsymbol{\nabla}\alpha_{w} \boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{w}-\boldsymbol{M}_{sw}^{D},\end{align}

(4.13f)\begin{align} \alpha_{s}+\alpha_{w}= 1 , \end{align}

where (4.13a) corresponds to the mass balance equation related to the rock phase and divided by a constant term, (4.13b) and (4.13c) represent the mass balance equations related to the solid bulk phase and to water, respectively, while (4.13d) and (4.13e) represent the momentum balance principles related to the solid bulk phase and water. With regard to (4.13f), this algebraic equation connects the solid bulk concentration with the water volume fraction in agreement with (3.18g). This class of simplified models maintains all the interphase terms except for the rock–ice interaction stresses, which have been deleted from the momentum equations as a consequence of the isokinetic condition. Furthermore, this condition does not produce any change in the melting process both in the way it is modelled (mass and momentum transfers) and in the volumetric deformation rate. With regard to this last aspect, by applying the isokinetic condition between rock and ice, (3.19) becomes

(4.14)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha_s\boldsymbol{u}_s+ \alpha_w\boldsymbol{u}_w)=\vert\varGamma_i\vert \left(\frac{1}{\rho_w}-\frac{1}{\rho_i}\right).\end{equation}

This relation implies that the rate of the mixture volumetric deformation reduces due to the ice melting, as it happens in the RIW model.

After providing proper closure relations, the equation system (4.13) can be solved in terms of the chosen unknowns, i.e. $\boldsymbol {u}_s$, $\boldsymbol {u}_w$, $\alpha _s$, $\alpha _w$, $\rho _s$ and $p_w$. It is worth noting that it is possible to retrieve the rock and ice concentrations, namely $\alpha _r$ and $\alpha _i$, from $\rho _s$ and $\alpha _s$ by applying (4.11a,b).

Finally, following the work of Pudasaini & Krautblatter (Reference Pudasaini and Krautblatter2014), an alternative expression for this system can be derived considering a new variable, here called ‘rock equivalent concentration’, that is defined as

(4.15)\begin{equation} \alpha_{{r}}^{{e}}=\alpha_{{r}}+\alpha_{{i}}\frac{\rho_{{i}}}{\rho_{{r}}},\end{equation}

where the ratio $\rho _{i}/\rho _{r}$ transforms the ice concentration into a rock equivalent term. Thanks to this new variable, the flow problem can be expressed as a function of $\boldsymbol {u}_{s}$, $\boldsymbol {u}_{w}$, $\alpha _{s}$, $\alpha _{r}^{e}$, $\alpha _{w}$ and $p_{w}$. To express the system (4.13) in terms of these new unknowns, it is necessary to rewrite the solid mass per unit volume, namely $\rho _s\alpha _s$, in terms of the rock equivalent concentration. This can be achieved by multiplying (4.8) by $\alpha _s$, by collecting $\rho _{r}$ and by using (4.15), as follows:

(4.16)\begin{equation} \rho_{s}\alpha_s=\rho_r\left( \alpha_{{r}}+\alpha_{i}\frac{\rho_{i}}{\rho_{{r}}}\right) =\rho_r\alpha_{{r}}^{{e}}.\end{equation}

The detailed expression of the PI-RIW system written in terms of the new unknowns can be found in Appendix B.

4.2. Fully isokinetic RIW model

The fully isokinetic RIW (FI-RIW) model can be obtained starting from the PI-RIW approach by imposing the isokinetic condition to the solid and water phases,

(4.17)\begin{equation} \boldsymbol{u}_{s}=\boldsymbol{u}_{w}=\boldsymbol{u}, \end{equation}

where $\boldsymbol {u}$ is the velocity of the overall bulk mixture. In this approach, the bulk concentration is equal to one,

(4.18)\begin{equation} \alpha=\alpha_{s}+\alpha_{w}=\alpha_{r}+\alpha_{i}+\alpha_{w}=1,\end{equation}

while the bulk density of the mixture becomes

(4.19)\begin{equation} \rho=\alpha_{s}\rho_{s}+\alpha_{w}\rho_{w}, \end{equation}

and it is variable in space and time.

For the derivation of the FI-RIW model, it is useful to define the overall stress tensor ${\boldsymbol{\mathsf{T}}}$ in a way consistent with the definition (4.10) provided for the solid stress tensor in the PI-RIW model. Thus, remembering that $\alpha _s+\alpha _w=1$, the overall stress tensor can be expressed as

(4.20)\begin{equation} {\boldsymbol{\mathsf{T}}}=\alpha_{s}{\boldsymbol{\mathsf{T}}}_{s}+\alpha_{w}{\boldsymbol{\mathsf{T}}}_{w}, \end{equation}

where the isotropic part of ${\boldsymbol{\mathsf{T}}}$ is $p{\boldsymbol{\mathsf{I}}}$ with $p$ representing the mixture pressure. By contrast with the PI-RIW model, the definition of the closure relation for ${\boldsymbol{\mathsf{T}}}$ translates into the description of the averaged mechanical behaviour of the mixture rather than into an equality of the mechanical behaviour of the solid and water phases.

In this class of models, the flow problem can be expressed as a function of $\boldsymbol {u}$, $\rho$, $p$ and two other unknown variables, for example $\alpha _{s}$ and $\alpha _{w}$. To do that, $\alpha _{s}\rho _{s}$ has to be expressed in terms of the chosen unknowns by using (4.19), namely

(4.21)\begin{equation} \alpha_{s}\rho_{s}=\rho-\alpha_{w}\rho_{w}. \end{equation}

Also for the FI-RIW model, the isokinetic condition allows us to replace the momentum balance principles related to water and the solid phase with the momentum conservation equation for the bulk mixture. This equation is obtained by summing (4.13d) with (4.13e) and by using the definition of $\rho$. The resulting expression corresponds to (4.23d), where all the interaction stresses between the phases disappear. It is also useful to sum the solid and water mass balances, namely equations (4.13b) and (4.13c). In this way, it is possible to derive the bulk mass conservation equation that replaces the solid mass balance in the PI-RIW model. Moreover, by subtracting (4.13a) to (4.13b), we get

(4.22)\begin{equation} \frac{\partial}{\partial t} ( \rho_{i}\alpha_{s} ) +\boldsymbol{\nabla\cdot} ( \rho_{i}\alpha_{s}\boldsymbol{u} ) ={-} \vert \varGamma_{i} \vert, \end{equation}

which can be further divided by the ice density $\rho _i$. As a result, the FI-RIW model reads as follows:

(4.23a)\begin{gather} \frac{\partial \alpha_{s}}{\partial t} +\boldsymbol{\nabla\cdot} ( \alpha_{s}\boldsymbol{u} ) ={-}\frac{ \vert\varGamma_{i} \vert}{\rho_{i}} , \end{gather}

(4.23b)\begin{gather}\frac{\partial\rho}{\partial t}+\boldsymbol{\nabla\cdot} ( \rho\boldsymbol{u} ) =0, \end{gather}

(4.23c)\begin{gather}\frac{\partial}{\partial t} ( \alpha_{w}\rho_{w} ) +\boldsymbol{\nabla\cdot} ( \alpha_{w}\rho_{w}\boldsymbol{u} ) ={+} \vert \varGamma_{i} \vert , \end{gather}

(4.23d)\begin{gather}\frac{\partial}{\partial t} ( \rho\boldsymbol{u} ) +\boldsymbol{\nabla\cdot} ( \rho\boldsymbol{u}\otimes\boldsymbol{u} ) ={-}\boldsymbol{\nabla\cdot}{\boldsymbol{\mathsf{T}}}+\rho\boldsymbol{g}, \end{gather}

(4.23e)\begin{gather}\alpha_{s}+\alpha_{w} =1, \end{gather}

where (4.23a) expresses the changes over time and space in the solid volume fraction. These changes are connected to the melting process of ice and this aspect is justified by the presence of $\rho _i$ in the denominator of the source term. It is worthy of note that the balance principle (4.23a) can be derived directly from the RIW model (3.18) by applying the isokinetic condition to rock and ice, by dividing (3.18a) and (3.18b) by $\rho _r$ and $\rho _i$, respectively, and by adding together the two resulting equations.

Finally, by dividing (4.23c) by $\rho _w$, adding the resulting equation to (4.23a) and by using (4.23e), it is possible to derive the expression for the velocity divergence,

(4.24)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=\vert\varGamma_i \vert \left(\frac{1}{\rho_w}-\frac{1}{\rho_i}\right),\end{equation}

which shows that the rate of the mixture volumetric deformation is negative as in the RIW model.

This class of models maintains only the mass transfers associated with the melting process, while it does not contain all the other interphase terms because of the isokinetic condition.

It is worth noting that the solution of (4.23) is expressed in terms of the unknown variables $\rho$, $\alpha _s$, $\alpha _w$, $\boldsymbol {u}$ and $p$. Knowing the mixture density $\rho$ and the solid and liquid concentrations ($\alpha _s$, $\alpha _w$), it is possible to retrieve the rock and ice concentrations by computing $\rho _s$ from (4.21) and by replacing this quantity into (4.11a,b).

4.3. Two-phase RIW model

In the two-phase RIW (TP-RIW) approach, the mixture is seen as composed of only two phases: one solid (rock plus ice) and one liquid. In this way, only four balance equations (two scalar and two vectorial) plus a relation for the concentrations are required to describe the flow. Therefore, starting from the PI-RIW approach, the number of unknown variables needs to be reduced by one scalar unit by imposing that a variable is constant. The choice of this variable is compulsory: since the concentrations and the velocities of the solid and liquid phases must be maintained, the variable to impose constant corresponds to $\rho _{s}$. This is consistent with what we have said in the introduction of this section: the incompressibility hypothesis must be applied to the bulk solid phase and, consequently the ratio $\alpha _{r}/\alpha _{i}$ is constant in time and space.

With regard to the equations, due to the incompressibility assumption, we need to neglect from the PI-RIW system one of the two mass balances that contain $\rho _{s}$. Nevertheless, the choice is not equivalent. Looking at these two equations, they do not boil down to the same expression under the incompressibility assumption. With regard to (4.13a), the term $\rho _{s}-\rho _{i}$ is constant and, since the source term is null, it can be simplified, thus giving rise to the following expression:

(4.25)\begin{equation} \frac{\partial\alpha_{s}}{\partial t}+\boldsymbol{\nabla\cdot}( \alpha_{s}\boldsymbol{u}_{s}) =0 .\end{equation}

On the other hand, in (4.13b) the source term is not null and thus, $\rho _{s}$ does not disappear from the equation when simplified. As a result, (4.13b) reads as follows:

(4.26)\begin{equation} \frac{\partial\alpha_{s}}{\partial t}+\boldsymbol{\nabla\cdot}( \alpha_{s}\boldsymbol{u}_{s}) ={-}\frac{\vert \varGamma_{i}\vert }{\rho_{s}}. \end{equation}

By combining (4.13c) with (4.26), it is easy to observe that the mass transfer condition existing between the solid and water phases is respected, since the source terms in these two balances are equal in absolute value but opposite in sign (see § 3.2 for more details). On the other hand, coupling (4.13c) and (4.25), the mass transfer condition is violated, since the source term in (4.13c) is not balanced by any counterpart. Therefore, the mass balance to maintain coincides with (4.26) multiplied by the constant $\rho _s$. This equation expresses the fact that the bulk solid phase is melted as a whole: both ice and rock are transformed into water at the same time, maintaining the ratio $\alpha _{r} /\alpha _{i}$ (and $\rho _{s}$) constant in time and space. It is quite easy to prove that the melting involves also the rock phase by deriving the mass balance principle related to rock. Firstly, by substituting (4.7) into (4.26), secondly, by expressing $\alpha _i$ in terms of $\alpha _r$ thanks to (4.5) and, lastly, by dividing the resulting equation by the constant term $(\rho _r-\rho _i)/(\rho _s-\rho _i)$ and by multiplying it by $\rho _r$, we get

(4.27)\begin{equation} \frac{\partial}{\partial t}(\alpha_{r}\rho_r)+\boldsymbol{\nabla\cdot}( \alpha_{r}\rho_r\boldsymbol{u}_{s}) ={-}\vert \varGamma_{i}\vert\frac{\rho_r }{\rho_{s}}\frac{(\rho_s-\rho_i)}{(\rho_r-\rho_i)}, \end{equation}

where the right-hand side term is a fraction of the mass transfer $\vert \varGamma _i\vert$.

The flow problem can be expressed as a function of $\boldsymbol {u}_{s}$, $\boldsymbol {u}_{w}$, $\alpha _{s}$, $\alpha _{w}$ and $p_{w}$. In this way, the TP-RIW model reads as follows:

(4.28a)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{s}\rho_{s} ) +\boldsymbol{\nabla\cdot} ( \alpha_{s}\rho_{s}\boldsymbol{u}_{s} ) ={-} \vert \varGamma_{s} \vert , \end{gather}

(4.28b)\begin{gather} \frac{\partial}{\partial t} ( \alpha_{w}\rho_{w} ) +\boldsymbol{\nabla\cdot} ( \alpha_{w}\rho_{w}\boldsymbol{u}_{w} ) ={+} \vert \varGamma_{s} \vert , \end{gather}

(4.28c)\begin{align} \frac{\partial}{\partial t} ( \alpha_{s}\rho_{s}\boldsymbol{u}_{s} ) +\boldsymbol{\nabla\cdot} ( \alpha_{s}\rho_{s}\boldsymbol{u}_{s}\otimes\boldsymbol{u}_{s} ) &={-}\boldsymbol{\nabla\cdot} ( \alpha_{s}{\boldsymbol{\mathsf{T}}}_{s} ) +\alpha_{s}\rho _{s}\boldsymbol{g}- \vert \varGamma_{s} \vert \boldsymbol{u}_{ss}\nonumber\\ &\quad +\boldsymbol{\nabla}\alpha_{s} \boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{w}+\boldsymbol{M}_{sw}^{D},\end{align}

(4.28d)\begin{align} \frac{\partial}{\partial t} ( \alpha_{w}\rho_{w}\boldsymbol{u} _{w} ) +\boldsymbol{\nabla\cdot} ( \alpha_{w}\rho_{w} \boldsymbol{u}_{w}\otimes\boldsymbol{u}_{w} ) &={-}\boldsymbol{\nabla\cdot} ( \alpha_{w}{\boldsymbol{\mathsf{T}}}_{w} )+\alpha_{w}\rho_{w}\boldsymbol{g} + \vert \varGamma_{s} \vert \boldsymbol{u}_{ss}\nonumber\\ &\quad +\boldsymbol{\nabla}\alpha_{w} \boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}_{w}-\boldsymbol{M}_{sw}^{D}, \end{align}

(4.28e)\begin{align} \alpha_{s}+\alpha_{w}= 1, \end{align}

where the subscript $s$ in the mass and momentum transfers is used to stress that the melting process involves the bulk solid phase and not the sole ice phase.

Finally, by dividing (4.28a) and (4.28b), respectively, by $\rho _s$ and $\rho _w$, adding together the resulting equations and by using (4.28e), it is possible to derive the equation for the velocity divergence,

(4.29)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha_s\boldsymbol{u}_s+\alpha_w \boldsymbol{u}_w)=\vert\varGamma_s \vert \left(\frac{1}{\rho_w}-\frac{1}{\rho_s}\right),\end{equation}

which states that the rate of the mixture volumetric deformation depends on the densities of water and the overall solid phase. Differently from the previous cases, the TP-RIW model allows both negative and positive volumetric deformation rates. If $\rho _s<\rho _w$, the behaviour is similar to that of the RIW model, while if $\rho _s>\rho _w$, the behaviour is opposed to that occurring in the RIW model. Since the solid density is defined as in (4.8), the condition $\rho _s<\rho _w$ occurs when the solid phase is mainly composed of ice. As a consequence, the physics of the volume deformations is in agreement with that of the RIW model only if the RIW mixture has a solid phase consisting mainly of ice. This aspect shows another severe limit of this model.

It is worth noting that, from the solution of (4.28), it is possible to retrieve the rock and ice concentrations by applying (4.11a,b) and by remembering that the solid density is a constant.

4.4. Isokinetic two-phase RIW model

The isokinetic two-phase RIW (ITP-RIW) model can be derived from the TP-RIW system (4.28) by applying the isokinetic condition to the water and solid phases. This condition is equivalent to that used in the FI-RIW approach and, therefore, the (4.17)–(4.21) hold. The flow problem can be expressed as a function of $\boldsymbol {u}$, $\rho$, $p$ and one unknown concentration, such as $\alpha _{w}$.

With regard to the equations, the momentum balance related to the bulk phase can be derived by adding together the momentum conservation equations for water and the solid phase, namely (4.28c) and (4.28d), and by using the definition of $\rho$. The resulting equation corresponds to (4.30c) and it is equal to (4.23d). In addition, the solid mass balance (4.28a) can be replaced by the bulk mass conservation equation, which can be obtained by adding (4.28a) to (4.28b). Finally, the solid concentration $\alpha _{s}$ can be expressed in terms of the chosen unknowns using (4.21). The final system for the ITP-RIW model reads as follows:

(4.30a)$$\begin{gather} \frac{\partial\rho}{\partial t}+\boldsymbol{\nabla\cdot} ( \rho\boldsymbol{u} ) =0, \end{gather}$$

(4.30b)$$\begin{gather}\frac{\partial}{\partial t} ( \alpha_{w}\rho_{w} ) +\boldsymbol{\nabla\cdot} ( \alpha_{w}\rho_{w}\boldsymbol{u} ) ={+} \vert \varGamma_{s} \vert, \end{gather}$$

(4.30c)$$\begin{gather}\frac{\partial}{\partial t} ( \rho\boldsymbol{u} ) +\boldsymbol{\nabla\cdot} ( \rho\boldsymbol{u}\otimes\boldsymbol{u} ) ={-}\boldsymbol{\nabla\cdot}{\boldsymbol{\mathsf{T}}}+\rho\boldsymbol{g}, \end{gather}$$

(4.30d)$$\begin{gather}\frac{\rho}{\rho_{s}}+\alpha_{w} \left( 1-\frac{\rho_{w}}{\rho_{s}} \right) =1 . \end{gather}$$

It is necessary to stress that the mass balances do not violate the mass transfer condition. Since $\rho$ is variable, the loss of mass expressed by the mass transfer $-\vert \varGamma _{s}\vert$ is indeed hidden in (4.30a). Finally, the equation for the velocity divergence can be derived by applying the isokinetic condition to the water and solid phases in (4.29), thus obtaining

(4.31)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}= \vert\varGamma_s\vert \left(\frac{1}{\rho_w}-\frac{1}{\rho_s}\right).\end{equation}

As in the TP-RIW model, this relation states that the rate of the mixture volumetric deformation is consistent with the RIW model only if $\rho _s<\rho _w$.

Since this class of models is derived from the TP-RIW approach, it is affected by the same approximation as the TP-RIW model and thus, the melting process is assigned to both rock and ice as a whole.

It is worth noting that, from the solution of (4.30), it is possible to retrieve the rock and ice concentrations by remembering that the solid density $\rho _s$ is a constant, computing the solid volume fraction from (4.21) and substituting $\alpha _s$ into (4.11a,b).

4.5. Mono-phase RIW model

In the mono-phase RIW (MP-RIW) model, the mixture is seen as composed of a single phase. In this way, only two balance equations are required to describe the flow and, therefore, the model can be obtained from the ITP-RIW approach by assuming that the bulk mixture $\rho$ is incompressible. Since the bulk mixture concentration is equal to one (see (4.18)), this assumption implies that even the solid and water concentrations, expressed by

(4.32a,b)\begin{equation} \alpha_{s}=\frac{\rho-\rho_{w}}{\rho_{s}-\rho_{w}},\quad \alpha_{w} =\frac{\rho-\rho_{s}}{\rho_{w}-\rho_{s}}, \end{equation}

are constant in time and space. As a consequence, one of the mass balances must be disregarded from the system (4.30), but, as in the TP-RIW approach, the choice is not arbitrary. Starting from the water mass balance (4.30b), the constancy of the water concentration allows us to simplify from the equation the term $\alpha _{w}\rho _{w}$, thus obtaining the following expression:

(4.33)\begin{equation} \boldsymbol{\nabla\cdot}\boldsymbol{u}={+}\frac{\vert \varGamma_{s}\vert }{\alpha_{w}\rho_{w}}. \end{equation}

Since the incompressibility condition, applied in this case to the mixture, reduces the required number of mass balances from two to one, the term appearing on the right-hand side of (4.33) is not balanced by any counterpart in another mass balance. As a consequence, this term appears to be linked to a mass production rather than to an internal mass transfer. As mass production is not considered in this work, this equation cannot be maintained in the system. Conversely, by applying the incompressibility assumption to the mixture mass balance (4.30a), we derive

(4.34)\begin{equation} \boldsymbol{\nabla\cdot}\boldsymbol{u}=0.\end{equation}

The monophase model needs to be derived choosing this equation, and in this way it is not able to describe the phase transformation associated with the melting process.

In this class of models, the flow problem can be expressed as a function of $\boldsymbol {u}$ and $p$. Therefore, the MP-RIW model reads as follows:

(4.35a)$$\begin{gather} \boldsymbol{\nabla\cdot}\boldsymbol{u} =0, \end{gather}$$

(4.35b)$$\begin{gather}\frac{\partial}{\partial t}( \rho\boldsymbol{u}) +\boldsymbol{\nabla\cdot}( \rho\boldsymbol{u}\otimes\boldsymbol{u} ) ={-}\boldsymbol{\nabla\cdot}{\boldsymbol{\mathsf{T}}}+\rho\boldsymbol{g}. \end{gather}$$

By contrast with the TP-RIW and ITP-RIW approaches, the MP-RIW model is based on a stronger approximation than the other two approaches. While the first two classes of models assign the melting to the overall solid phase, the MP-RIW approach overlooks it. Furthermore, due to (4.35a), no changes in the volumetric deformation rate occur within the mixture.

4.6. Trade-off between simplicity and completeness in the RIW framework

As described in the previous sections, the simplified RIW models are obtained reducing gradually the number of unknowns in the equation systems. Table 1 shows a synoptic panel concerning the unknowns used to describe each class of models. Furthermore, the simplified RIW models are derived by applying alternately or simultaneously the isokinetic and incompressibility conditions. As shown in table 2, these two assumptions have different effects on the original system of equations. While the isokinetic condition acts by reducing the complexity of the momentum balances, the incompressibility condition acts by simplifying the composition of the mixture and, in this way, deleting from the original system one or two mass balances.

Table 1. Classes of models and corresponding unknowns. The subscript $s$ denotes the overall solid phase (rock plus ice). No subscript indicates the overall bulk mixture. Details on the unknowns are given in the related model paragraphs.

Table 2. Effects of the simplifications used in the framework of simplified RIW models. The symbols ✓and X denote, respectively, whether the quantities appear in the class of models or not. Under ‘concentration variability’, the term ‘ratio’ stands for the concentration ratio and the term ‘single’ indicates the single concentration of the phases in each class of models.

With reference to the isokinetic condition, the classes based on this assumption simplify the dynamics of RIW mixtures deleting from the original system of equations the solid–solid (PI-RIW, TP-RIW) and solid–liquid (FI-RIW, ITP-RIW, MP-RIW) interactions. From a physical point of view, the isokinetic condition implies that the simplified models may be unable to simulate phenomena such as the segregation between two distinct solid phases and the phase separation between the solid and liquid phases (for specific references to these phenomena, see e.g. Drahun & Bridgwater (Reference Drahun and Bridgwater1983), Iverson (Reference Iverson1997), Gray & Thornton (Reference Gray and Thornton2005), Johnson et al. (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012) and Pudasaini & Fischer (Reference Pudasaini and Fischer2020)). More precisely, the mathematical approaches derived applying the isokinetic condition to rock and ice (PI-RIW, TP-RIW) can simulate the phase separation, but not the phenomenon of segregation. Conversely, the fully isokinetic models (FI-RIW, ITP-RIW, MP-RIW) are not able to simulate either the phase separation or the segregation processes, since all the phases move with the same velocity.

With regard to the incompressibility condition, this hypothesis imposes an important constraint on the phase concentrations, which translates into a strong simplification of the melting process. In the TP-RIW and ITP-RIW models, this constraint implies that the two solid phases (rock and ice) can change over time and space their concentrations, although maintaining their ratio constant. This translates into the fact that the concentration constraint assigns the melting process not only to the ice phase but to the overall bulk solid phase. Conversely, in the MP-RIW model, the constraint translates into a constancy of all the phase concentrations, thus overlooking the melting process. As a result, these three models (TP-RIW, ITP-RIW, MP-RIW) simplify the description of the melting process.

Since the hypothesis of incompressibility strongly approximates the melting process, it is advisable to model RIW mixtures with the PI-RIW and FI-RIW models. Unlike the FI-RIW approach, the PI-RIW description can detect differences between the solid and liquid velocities. For this reason, the PI-RIW model appears to be an appropriate trade-off between simplicity and completeness.

4.7. Classification of literature models

The mathematical models existing in the literature and briefly presented in the introduction can be classified on the basis of the RIW framework that we have developed. In this way, it is possible to catch the degree of simplification that these models show. In figure 2 we report a scheme in which the literature models are arranged within the corresponding class of simplified RIW approaches. It is worth noting that in this scheme the different simplified RIW models are grouped into two categories, called ‘reduced’ and ‘approximated’, respectively. The reason for this distinction lies in the fact that, while the ‘reduced’ models are derived imposing only a condition on the flow and maintaining the melting process as it appears in the RIW model, the ‘approximated’ approaches modify or even overlook the mechanism of melting.

Figure 2. Classes of simplified RIW models with their corresponding acronyms and classification of the models existing in the literature within the RIW framework. The differences in colours give, as in figure 1, the differences in the number of phases considered in each class of simplified RIW models.

The works of Evans et al. (Reference Evans, Tutubalina, Drobyshev, Chernomorets, McDougall, Petrakov and Hungr2009b) and Schneider et al. (Reference Schneider, Bartelt, Caplan-Auerbach, Christen, Huggel and McArdell2010), which are based on the mathematical models of McDougall (Reference McDougall2006), Hungr & McDougall (Reference Hungr and McDougall2009) and Christen et al. (Reference Christen, Kowalski and Bartelt2010), can be all classified as MP-RIW approaches since these models are described by the system (4.35). As a consequence, they overlook the melting process.

The model proposed by Pudasaini & Krautblatter (Reference Pudasaini and Krautblatter2014), known as the two-phase model of Pudasaini & Krautblatter (TPPK), is a two-phase approach, in which the mixture is composed of a solid phase (rock plus ice) and a liquid phase. The related system of equations consists of two mass and two momentum balance equations as in the TP-RIW model presented in § 4.3, but it shows some differences. The reason for these differences lies in the fact that the authors assumed a particular definition of the liquid volume fraction (Pudasaini & Krautblatter Reference Pudasaini and Krautblatter2014, p. 2279) which, using our notation, becomes expressed as

(4.36)\begin{equation} \alpha_w=1-\alpha_r^e,\end{equation}

rather than as $\alpha _w=1-\alpha _s$, as specified in our work. Considering the same liquid phase and equalizing the two definitions of the liquid concentration, it follows that $\alpha _s=\alpha _r^e$, and using (4.16), this condition implies that $\rho _s=\rho _r$. These equalities represent the extra hypotheses that we have to add to the isokinetic and incompressibility conditions in order to derive the TPPK model from the TP-RIW approach. Due to these extra hypotheses, the TPPK model is classified as ‘other simplified models’. Despite this aspect, the TPPK model shares the same physical limits as the TP-RIW approach, but it can be considered as a reference model for the description of the effects of basal lubrication and internal fluidization which, to the best of our knowledge, have never been considered before in multiphase flows.

Finally, the core layer of the model proposed by Bartelt et al. (Reference Bartelt, Christen, Bühler and Buser2018) can be classified as a FI-RIW model, since rock, ice and water propagate downslope with the same velocity. Thus, as shown in (4.23), this model is based on the isokinetic condition between the three different phases and the melting process is modelled through a mass transfer between ice and water.

From this analysis of the literature, it appears that the PI-RIW model represents a novel approach to tackle the rock–ice avalanche problem.

5.1. Solid and liquid rheologies

In this paragraph, we discuss the rheologies used for the solid and liquid phases. These relations depend on the properties of the material considered and affect strongly its dynamics.

For the water stress tensor, we assume a Newtonian closure

(5.1)\begin{equation} {\boldsymbol{\mathsf{T}}}_{w}=p_{w}{\boldsymbol{\mathsf{I}}}-\boldsymbol{\tau}_{w}, \end{equation}

where ${\boldsymbol{\mathsf{I}}}$ and $\boldsymbol {\tau }_{w}=2\rho _{w} \nu {\boldsymbol{\mathsf{D}}}_{w}$ represent, respectively, the identity tensor and the water deviatoric tensor, while ${\boldsymbol{\mathsf{D}}}_{w}$ and $\nu$ correspond to the shear velocity tensor and the kinematic viscosity. As pointed out in the introduction, this rheology can be used not only whether the liquid phase is water, but also whether the liquid phase consists of water with low concentrations of silty solid particles (Armanini Reference Armanini2013; Armanini et al. Reference Armanini, Larcher, Nucci and Dumbser2014). In cases of high concentrations of silty solid particles, for the liquid stress tensor, it is more appropriate to use a non-Newtonian closure relation, as proposed in Pudasaini & Krautblatter (Reference Pudasaini and Krautblatter2014) and Pudasaini & Mergili (Reference Pudasaini and Mergili2019).

With regard to the solid rheology, many constitutive relations exist in the literature (Savage & Hutter Reference Savage and Hutter1989; Iverson & Denlinger Reference Iverson and Denlinger2001; Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2006; Gray & Edwards Reference Gray and Edwards2014; Pudasaini & Mergili Reference Pudasaini and Mergili2019). Nevertheless, for the effective solid stress tensor $\mathcal {T}_{s}$, we consider a Mohr–Coulomb plasticity closure (Savage & Hutter Reference Savage and Hutter1989; Pudasaini & Hutter Reference Pudasaini and Hutter2007; Pudasaini & Mergili Reference Pudasaini and Mergili2019) in agreement with the literature models discussed in § 4.7:

(5.2)$$\begin{gather} \mathcal{T}_{s}^{xx}=K\mathcal{T}_{s}^{zz}, \end{gather}$$

(5.3)$$\begin{gather}\mathcal{T}_{s}^{xz}={-}{sgn}( u_{s}^{x}) \mu \mathcal{T}_{s}^{zz}, \end{gather}$$

where $K$ and $\mu$ are, respectively, the active-passive earth pressure coefficient and the friction coefficient, while $-{sgn}( u_{s}^{x})$ specifies that the shear stress opposes the motion. The stress $\mathcal {T}_{s}^{zz}$ is defined by

(5.4)\begin{equation} \mathcal{T}_{s}^{zz}=T_{s}^{zz}-p_{w}+\tau_{w}^{zz}. \end{equation}

5.2. Boundary conditions

The boundary conditions used to derive the depth-integrated versions of the models proposed in § 4 correspond to those used by Pitman & Le (Reference Pitman and Le2005) and Pudasaini (Reference Pudasaini2012) and can be summarized as follows.

(i) The free surface and the bottom are assumed to be material surfaces and thus, they need to satisfy the kinematic boundary conditions. (ii) The dynamic boundary condition is applied at the free surface. Hence, we assume that the action exerted on the mixture by the outside environment is negligible. Moreover, assuming small curvatures of the free surface, also the surface tension is negligible. (iii) The bottom is assumed to be fixed and impermeable. As a consequence, the $z$ components of the velocity vectors for the solid and water phases are null. (iv) The liquid phase needs to satisfy the no-slip condition which implies that (over a fixed bed) the whole liquid velocity vector at the bottom is null.

5.3. Integration along the flow depth

An important step used to derive the depth-integrated version of the PI-RIW model (omitted for brevity) consists in checking if rock–ice avalanches are characterized by the shallow flow condition. As pointed out in the work of Huggel et al. (Reference Huggel, Zgraggen-Oswald, Haeberli, Kääb, Polkvoj, Galushkin and Evans2005, p. 183) concerning the Kolka–Karmadon rock–ice avalanche, this flow had a value of the ratio flow height-to-longitudinal extent that was similar to those related to debris flows and lahars. Assuming this statement general, it is possible to affirm that also rock–ice avalanches are characterized by the shallow flow condition, as debris flows and lahars (Iverson Reference Iverson1997). As a consequence, it is possible to perform the scale analysis combined with the shallow flow condition. As a result of this analysis, the momentum balance principles along the $z$ direction simplify strongly, thus leading to the following equations:

(5.5a)$$\begin{gather} \frac{\partial p_{w}}{\partial z}={-}\rho_{w}g\cos\zeta , \end{gather}$$

(5.5b)$$\begin{gather}\frac{\partial( \alpha_{s}\mathcal{T}_{s}^{zz})}{\partial z} ={-}\rho_{r}\left(\alpha_{r}^{e}-\alpha_{s}\frac{\rho_{w}}{\rho_{r}}\right) g\cos\zeta . \end{gather}$$

With regard to (5.5a), it states that the water pressure is hydrostatic. Considering (5.5b), this equation might appear strange due to the presence of both $\alpha _r^e$ and $\alpha _s$. Nevertheless, by multiplying the term inside the round brackets by $\rho _r$, we get $\rho _r\alpha _r^e-\rho _w\alpha _s$, which becomes equal to $\rho _s\alpha _s-\rho _w\alpha _s$ thanks to (4.16). In this way, the momentum balance principle (5.5b) can be rewritten as follows:

(5.6)\begin{equation} \frac{\partial(\alpha_{s}\mathcal{T}_{s}^{zz})}{\partial z}={-} (\rho_s\alpha_{s}-\rho_w\alpha_{s}) g\cos\zeta, \end{equation}

and it states that the normal solid stress in the $z$ direction depends on the solid mass per unit volume reduced by the buoyancy. It is worthy of note that, in case of mixtures composed only of water and rock, (5.5b) coincides with that derived for two-phase debris flow models (Pitman & Le Reference Pitman and Le2005; Pudasaini Reference Pudasaini2012), since $\alpha _r^e=\alpha _s=\alpha _r$.

With reference to the momentum balance principles along the $x$ direction, the shallow flow condition allows the viscous stresses to be neglected (Iverson & Denlinger Reference Iverson and Denlinger2001; Pitman & Le Reference Pitman and Le2005), and this simplification can be considered a reasonable approximation, as demonstrated by the recent laboratory experiments of Armanini et al. (Reference Armanini, Larcher, Nucci and Dumbser2014).

Once the scale analysis has been performed, it is possible to integrate the equations between the free surface and the bottom and to express them as functions of the depth-averaged variables. On the basis of the flow depth $h$, the depth-averaged value of a generic function $b(x, z, t)$ can be defined as follows:

(5.7)\begin{equation} \bar{b}(x, t)=\frac{1}{h(x,t)}\int_{z_{b}(x)}^{\eta(x,t)} b(x, z, t)\, \text{d}z, \end{equation}

where the averaged variable is denoted here by an overline. In the following, to simplify the notation, we will disregard the explicit indication of the space–time variability of a quantity. In addition, it is worth noting that some terms appear in the equation system as a product of two or three quantities. In these cases, the depth-averaged terms read as follows:

(5.8)$$\begin{gather} \beta_{2}\bar{b}_{1}\bar{b}_{2} =\frac{1}{h}\int_{z_{b}}^{\eta }b_{1}\,b_{2} \,\text{d}z, \end{gather}$$

(5.9)$$\begin{gather}\beta_{3}\bar{b}_{1}\bar{b}_{2}\bar{b}_{3} =\frac{1}{h} \int_{z_{b}}^{\eta}b_{1}\,b_{2}\,b_{3} \,\text{d}z, \end{gather}$$

where $b_{1},b_{2}$ and $b_{3}$ identify three generic functions, while $\beta _{2}$ and $\beta _{3}$ represent suitable corrective factors. Although these parameters can affect the further analysis of the eigenvalues, all these corrective factors can be approximated to one in agreement with Pitman & Le (Reference Pitman and Le2005), Pelanti, Bouchut & Mangeney (Reference Pelanti, Bouchut and Mangeney2008) and Pudasaini (Reference Pudasaini2012). Thanks to these expressions, we define the following depth-integrated variables:

(5.10)$$\begin{gather} C_{k} =\frac{1}{h}\int_{z_{b}}^{\eta}\alpha_{k}\,{{\rm d}}z, \end{gather}$$

(5.11)$$\begin{gather}\boldsymbol{U}_{k} =\frac{1}{h}\int_{z_{b}}^{\eta}\boldsymbol{u} _{k} \,\text{d}z , \end{gather}$$

where $k\in \{s, w\}$,

(5.12)\begin{equation} C_{r}^{e}=\frac{1}{h}\int_{z_{b}}^{\eta}\alpha_{r}^{e} \,\text{d}z \end{equation}

and

(5.13)$$\begin{gather} \varUpsilon_{i} =\frac{1}{h}\int_{z_{b}}^{\eta}\vert \varGamma _{i}\vert \,\text{d}z, \end{gather}$$

(5.14)$$\begin{gather}\varUpsilon_{i}\boldsymbol{U}_{ss} =\frac{1}{h}\int_{z_{b}}^{\eta}\vert \varGamma_{i}\vert \boldsymbol{u}_{ss} \,\text{d}z, \end{gather}$$

(5.15)$$\begin{gather}\bar{M}_{sw}^{D} =\frac{1}{h}\int_{z_{b}}^{\eta}M_{sw} ^{D}|^{x} \,\text{d}z. \end{gather}$$

As a result of the depth-integration, the 1-D PI-RIW model reads as follows:

(5.16a)\begin{gather} \frac{\partial( C_{s}h) }{\partial t}+\frac{\partial( C_{s}U_{s}h) }{\partial x} ={-}\frac{\varUpsilon_i h}{\rho_{i} } , \end{gather}

(5.16b)\begin{gather} \frac{\partial( C_{r}^{e}h) }{\partial t}+\frac{\partial( C_{r}^{e}U_{s}h) }{\partial x} ={-}\frac{\varUpsilon_i h}{\rho_{r}} , \end{gather}

(5.16c)\begin{gather} \frac{\partial( C_{w}h) }{\partial t}+\frac{\partial( C_{w}U_{w}h) }{\partial x} ={+}\frac{\varUpsilon_i h}{\rho_{w}} , \end{gather}

(5.16d)\begin{align} \frac{\partial( C_{r}^{e}U_{s}h) }{\partial t}+\frac {\partial( C_{r}^{e}U_{s}^{2}h) }{\partial x}&={-}\frac{\partial}{\partial x}\left[ Kg\cos \zeta\left( C_{r}^{e}-C_{s}\frac{\rho_{w}}{\rho _{r}}\right) \frac{h^{2}}{2}\right]\nonumber\\ & \quad -gh\cos\zeta\left( K\frac{\partial z_{b}}{\partial x}+{sgn}(U_s)\mu\right) \left(C_{r}^{e}- C_{s}\frac{\rho_{w}}{\rho_{r}}\right)\nonumber\\ &\quad -g\frac{\rho_{w}}{\rho_{r}}C_{s}h\cos\zeta\frac{\partial h}{\partial x}-g\frac{\rho_{w}}{\rho_{r}}C_{s}h\cos\zeta\frac{\partial z_{b}}{\partial x}\nonumber\\ &\quad+\frac{1}{\rho_{r}}[\rho_rC_{r}^{e}gh\sin\zeta- \varUpsilon_{i}U_{ss}h+\bar{M}_{sw}^{D}h], \end{align}

(5.16e)\begin{align} \frac{\partial( C_{w}U_{w}h) }{\partial t}+\frac{\partial ( C_{w}U_{w}^{2}h) }{\partial x}&={-}g\cos\zeta C_{w} h\frac{\partial h}{\partial x}-g\cos\zeta C_{w}h\dfrac{\partial z_{b} }{\partial x}\nonumber\\ &\quad+\frac{1}{\rho_{w}}[\rho_wC_{w}gh\sin\zeta+ \varUpsilon_{i}U_{ss}h-\bar{M}_{sw}^{D}h], \end{align}

(5.16f)\begin{align} C_{s}+C_{w} =1, \end{align}

where (5.16a)–(5.16e) are divided by $\rho _i$, $\rho _r$, $\rho _w$, $\rho _r$ and $\rho _w$, respectively. The (5.16a)–(5.16c) derive from the mass balance principles related to the solid phase, the rock equivalent phase and water, while (5.16d) and (5.16e) correspond to the momentum balance principles related to the rock equivalent phase and water. Looking at the first three balance principles in (5.16), it might appear that (5.16a) and (5.16b) express the changes over time and space in the concentration of the same phase. Nevertheless, while $C_s$ describes the overall solid phase existing in the mixture, $C_r^e$ can be seen as a variable that describes the composition of the overall solid phase. In fact, by using (4.11a,b) and (4.16), it is possible to combine together $C_s$ and $C_r^e$ and to derive the depth-integrated concentrations of rock and ice,

(5.17a,b)\begin{equation} C_{i}=\rho_r\frac{C_{r}^e-C_{s}}{\rho_{i}-\rho_{r}}; \quad C_{r}=\frac{C_r^e\rho_{r}-C_s\rho_{i}}{\rho_{r}-\rho_{i}}. \end{equation}

It is worthy of note that it is possible to detect the changes over time and space in the rock and ice concentrations only thanks to both (5.16a) and (5.16b), and this detection can be computed without imposing any other assumption, such as the incompressibility of the solid phase as in the TP-RIW model.

Closure relations. The unknown variables of the 1-D PI-RIW model (5.16) correspond to $h$, $C_s$, $C_r^e$, $C_w$, $U_s$ and $U_w$, while the terms that remain unspecified consist in the mass and momentum transfers ($\varUpsilon _{i}$, $U_{ss}$) and the drag stress ($\bar {M} _{sw}^{D}$).

(i) The mass transfer $\varUpsilon _{i}$ can be estimated, as proposed in the model of Pudasaini & Krautblatter (Reference Pudasaini and Krautblatter2014), on the basis of the melting efficiency derived by Sosio et al. (Reference Sosio, Crosta, Chen and Hungr2012). More precisely, the molten ice is expressed in terms of a thickness that is computed by assuming that only a fraction of the heat produced by basal friction is used for the ice melting. In this way, the mass transfer $\varUpsilon _{i}$ can be defined as follows:

(5.18)\begin{equation} \varUpsilon_i=\frac{e\mu g\cos\zeta}{\lambda_f}\rho_rC_r^eU_s,\end{equation}

where $\lambda _f$ identifies the latent heat of ice, while $e$ is the parameter that estimates the fraction of heat produced by basal friction and used for the ice melting.

(ii) The velocity $U_{ss}$ appearing in the momentum transfer is defined in agreement with Kolev (Reference Kolev2007). Thus, $U_{ss}$ is assumed to be equal to the velocity of the ‘donor’ phase, namely the solid bulk phase,

(5.19)\begin{equation} U_{ss}=U_s.\end{equation}

(iii) With regard to the drag stress $\bar {M}_{sw}^D$, many authors provided closure relations for this term (Wang & Hutter Reference Wang and Hutter1999; Pitman & Le Reference Pitman and Le2005; Meruane, Tamburrino & Roche Reference Meruane, Tamburrino and Roche2010; Pudasaini Reference Pudasaini2012, Reference Pudasaini2020; Nucci, Armanini & Larcher Reference Nucci, Armanini and Larcher2019). A possible choice consists in using the expression defined in Pelanti et al. (Reference Pelanti, Bouchut and Mangeney2008):

(5.20)\begin{equation} \bar{M}_{sw}^D=\mathcal{D}(U_w-U_s),\end{equation}

where $\mathcal {D}$ is a function depending on the solid bulk volume fraction, on the difference in the solid and liquid velocities and on constant physical parameters (e.g. the characteristic grain size) and is based on the drag correlation proposed by Gidaspow (Reference Gidaspow1994). Thanks to this choice, quantitative comparisons between the debris-flow model of Pelanti et al. (Reference Pelanti, Bouchut and Mangeney2008) and the PI-RIW approach proposed in this work might be computed. However, since the drag stresses can play an important role in multiphase mass flows (Nucci et al. Reference Nucci, Armanini and Larcher2019; Pudasaini Reference Pudasaini2020), further analyses should be performed using or deriving different drag closure relations. In this way, it is possible to detect which of these closure relations may fit with the rock–ice avalanche case.

6.1. The PI-RIW eigenvalue set

The differential part of the system (5.16) can be written in compact form as

(6.1)\begin{equation} \frac{\partial\boldsymbol{U}(\boldsymbol{W})}{\partial t}+\frac{\partial \boldsymbol{F}(\boldsymbol{W})}{\partial x}+{\boldsymbol{\mathsf{H}}}(\boldsymbol{W} )\frac{\partial\boldsymbol{W}}{\partial x}=\boldsymbol{S}(\boldsymbol{W}), \end{equation}

where

(6.2)\begin{equation} \boldsymbol{W} =[ \begin{array} {ccccc} h & C_{r}^{e} & C_{w} & U_{s} & U_{w} \end{array} ]^\textrm{T} \end{equation}

is the vector of primitive variables, while $\boldsymbol {U}(\boldsymbol {W})$, $\boldsymbol {F}(\boldsymbol {W})$ and $\boldsymbol {S}(\boldsymbol {W})$ represent the vectors of the conserved variables, the conservative fluxes and of the source terms, respectively. Finally, ${\boldsymbol{\mathsf{H}}}(\boldsymbol {W})$ is the matrix of the non-conservative fluxes. All these terms read as follows:

(6.3) \begin{equation} \left.\begin{array}{c@{}} \boldsymbol{U}=\left[ \begin{array}{c} (1-C_{w})h\\ C_{r}^{e}h\\ C_{w}h\\ C_{r}^{e}U_{s}h\\ C_{w}U_{w}h\\ \end{array}\right] ;\quad\boldsymbol{F}=\left[ \begin{array}{c} (1-C_{w})U_{s}h\\ C_{r}^{e}U_{s}h\\ C_{w}U_{w}h\\ C_{r}^{e}U_{s}^{2}h+Kg\cos\zeta\left( C_{r}^{e}-\dfrac{\rho_{w}}{\rho_{r}}(1-C_{w})\right) \dfrac{h^{2}}{2}\\ C_{w}U_{w}^{2}h\\ \end{array}\right] \\ {\boldsymbol{\mathsf{H}}}=\left[ \begin{array}{ccccc} 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ g\dfrac{\rho_{w}}{\rho_{r}}\cos\zeta(1-C_{w})h & 0 & 0 & 0 & 0\\ g\cos\zeta C_{w}h & 0 & 0 & 0 & 0\\ \end{array}\right] ;\quad\boldsymbol{S}=\left[ \begin{array}{c} -\dfrac{\varUpsilon_i}{\rho_{i}}h \\ -\dfrac{\varUpsilon_i}{\rho_{r}}h \\ +\dfrac{\varUpsilon_i}{\rho_{w}}h \\ S_{s} \\ S_{w}\\ \end{array}\right] \end{array}\right\}, \end{equation}

where

(6.4) \begin{align} S_{s} &={-}g\cos\zeta\left( K\frac{\partial z_{b}}{\partial x}+\mu\right) \left( C_{r}^{e}- (1-C_{w})\frac{\rho_{w}}{\rho_{r}}\right) h-g\frac{\rho _{w}}{\rho_{r}}(1-C_{w})h\cos\zeta\frac{\partial z_{b}}{\partial x}\nonumber\\ &\quad+\frac{1}{\rho_{r}}[\rho_rC_{r}^{e}gh\sin\zeta- \varUpsilon_{i}U_{ss}h+\bar{M}_{sw}^{D}h], \end{align}

(6.5) $$\begin{gather}S_{w}={-}g\cos\zeta C_{w}h\frac{\partial z_{b}}{\partial x}+\frac{1} {\rho_{w}}[ \rho_{w}C_{w}gh\sin\zeta+\varUpsilon_{i}U_{ss}h- \bar{M}_{sw}^{D}h]. \end{gather}$$

It is worth noting that in these expressions we have used (5.16f) to substitute $C_{s}$ with $(1-C_{w})$.

The quasi-linear form of the homogeneous part of the system becomes

(6.6)\begin{equation} \frac{\partial\boldsymbol{U}}{\partial\boldsymbol{W}}\frac{\partial \boldsymbol{W}}{\partial t}+\frac{\partial\boldsymbol{F}}{\partial \boldsymbol{W}}\frac{\partial\boldsymbol{W}}{\partial x} +{\boldsymbol{\mathsf{H}}}(\boldsymbol{W})\frac{\partial\boldsymbol{W}}{\partial x}=0, \end{equation}

where $\partial \boldsymbol {U}/\partial \boldsymbol {W}$ and $\partial \boldsymbol {F}/ \partial \boldsymbol {W}$ are the Jacobians of the conserved vector and of the conservative flux vector with respect to the primitive variables. The characteristic polynomial is computed from the following expression:

(6.7)\begin{equation} P(\lambda)=\mbox{Det}\left[ \frac{\partial\boldsymbol{F}}{\partial \boldsymbol{W}}+{\boldsymbol{\mathsf{H}}}(\boldsymbol{W})-\lambda\frac{\partial \boldsymbol{U}}{\partial\boldsymbol{W}}\right], \end{equation}

thus becoming

(6.8)\begin{equation} P(\lambda)=( \lambda- U_{s}) \sum_{l=0}^{4}b_{l}\lambda^{l}.\end{equation}

The coefficients $b_{l}$ are written in an analogous way to those derived for two-phase flows by Pelanti et al. (Reference Pelanti, Bouchut and Mangeney2008):

(6.9a)$$\begin{gather} b_{4}=1, \end{gather}$$

(6.9b)$$\begin{gather}b_{3}={-}2 ( U_{w}+U_{s} ) , \end{gather}$$

(6.9c)$$\begin{gather}b_{2}= ( U_{s}+U_{w} )^{2}+2U_{s}U_{w}-a^{2} ( 1+\beta ^{2} ) , \end{gather}$$

(6.9d)$$\begin{gather}b_{1}={-}2U_{s}U_{w} ( U_{s}+U_{w} ) +2a^{2} ( U_{s} +\beta^{2}U_{w} ) -2a^{2} ( 1-C_{w} ) ( U_{s} -U_{w} ), \end{gather}$$

(6.9e)$$\begin{gather}b_{0}=U_{s}^{2}U_{w}^{2}-a^{2}(U_{s}^{2}+\beta^{2}U_{w}^{2})+a^{4} ( \beta^{2}-\delta^{2} ) +a^{2} ( 1-C_{w} ) ( U_{s}^{2}-U_{w}^{2} ), \end{gather}$$

where

(6.10a)$$\begin{gather} a =\sqrt{gh\cos\zeta}, \end{gather}$$

(6.10b)$$\begin{gather}\beta^2 ={\left[ C_{w}-1-\frac{K}{2}( C_{w}-2) \right] \left( 1-\frac{\rho_{w}}{\rho_{r}}\frac{1-C_{w}}{C^{e}_{r}}\right) }, \end{gather}$$

(6.10c)$$\begin{gather}\delta^2 ={( 1-C_{w}) ( K-1) \left( 1-\frac{\rho_{w}}{\rho_{r}}\frac{1-C_{w}}{C^{e}_{r}}\right) }. \end{gather}$$

To analyse the trend of the eigenvalues, it is useful to express the characteristic polynomial (6.8) in a non-dimensional form by dividing it by $a^{5}$. Therefore, its expression can be written in terms of the dimensionless eigenvalue $\hat {\lambda }=\lambda /a$ and of the following dimensionless quantities:

(6.11)\begin{equation} C_{r}^{e}, C_{w}, K, Fr_{w}=\frac{U_{w}}{a}, Fr_{s}=\frac{U_{s}}{a}, \end{equation}

where the last two terms represent the Froude numbers of the water and solid phases, respectively.

Following a classical way to represent the eigenvalues for free-surface flows (Lyn & Altinakar Reference Lyn and Altinakar2002; Garegnani et al. Reference Garegnani, Rosatti and Bonaventura2013), we plot the behaviour of the different $\hat {\lambda }$ as a function of the water Froude number $Fr_{w}$ by keeping fixed the values of the other parameters. A typical behaviour of these plots is presented in figure 4. It allows us to make three general considerations.

Figure 4. Trend of the dimensionless eigenvalues obtained considering $Fr_{s}/Fr_{w}=0.7$, $C_{w}=0.5$, $C_{r}^{e}=0.37$, $K=1$ and $\rho _{i} /\rho _{r}=0.353$.

(i) One eigenvalue, namely $\hat {\lambda }=Fr_{s}$, has a linear trend, whose slope depends on the values of the parameters used.

(ii) Two eigenvalues are nearly linear and, in this case, their slope depends on the values of the parameters used.

(iii) The eigenvalues are real and distinct for all values of $Fr_{w}$ except for a given range where two eigenvalues become complex conjugate, thus delimiting a region with a loss of hyperbolicity. The $Fr_{w}$ range in which this loss occurs changes as a function of the parameters used.

To understand how the changes in the parameters affect the features highlighted, a careful analysis is performed by analysing the changes in the plots due to the variation in a single parameter at a time.

6.1.1. Effect of the Froude ratio

Figure 5 shows three eigenvalue plots in which the ratio $Fr_{s}/Fr_{w}$, representing the ratio between the solid and water velocities, is increased from $0.7$ (solid slower than water) to $1.3$ (solid faster than water). It is possible to notice that: (i) the eigenvalues tend to increase more rapidly as the ratio of the Froude numbers increases; (ii) for $Fr_{s}/Fr_{w}=1$ (isokinetic condition), the eigenvalues remain distinct and real for all values of $Fr_{w}$; (iii) moving away from the isokinetic condition, the loss of hyperbolicity appears and is confined to a range whose boundaries tend to move towards higher values of the water Froude number as $Fr_{s}/Fr_{w}\rightarrow 1$.

Figure 5. Effect of the Froude ratio on the dimensionless eigenvalue trends. Results are obtained considering $C_{w}=0.5$, $C_{r}^{e}=0.37$, $K=1$ and $\rho _{i}/\rho _{r}=0.353$.

This aspect can be highlighted in figure 6, where the lower and upper limits of the non-hyperbolic region, namely $Fr_{w}^{low}$ and $Fr_{w}^{up}$, are plotted as functions of the Froude ratio $Fr_{s}/Fr_{w}$. This figure shows that: (i) the behaviour is symmetrical with respect to $Fr_{s}/Fr_{w}=1$, at least for $0< Fr_{s}/Fr_{w}<2$; (ii) $Fr_{w}^{low}\rightarrow \infty$ and $Fr_{w}^{up}\rightarrow \infty$ as $Fr_{s}/Fr_{w}\rightarrow 1$; (iii) the extent of the non-hyperbolic region increases as $Fr_{s} /Fr_{w}\rightarrow 1$.

Figure 6. Lower ($Fr_{w}^{low}$) and upper ($Fr_{w}^{up}$) limits of the non-hyperbolic region as functions of the Froude ratio $Fr_{s}/Fr_{w}$. Results are obtained, as in figure 4, considering $C_{w}=0.5$, $C_{r}^{e}=0.37$, $K=1$ and $\rho _{i}/\rho _{r}=0.353$.

To sum up, approaching the isokinetic condition, the region of the loss of hyperbolicity increases in size, but it starts at values of $Fr_{w}$ that become greater and greater, and that are far from what can be expected in real situations.

It is worthy of note that in the isokinetic condition the strict hyperbolicity of the system occurs independently from the value of all the other parameters. Therefore, in the following paragraphs, we focus the attention only on the effects of the other parameters on the non-hyperbolic range. Due to the symmetrical behaviour of this range, the analyses reported below are performed considering a single value for the Froude ratio.

6.1.2. Effect of $C_{w}$

As regards the water concentration, figure 7 shows the eigenvalue plots when $C_{w}$ changes from $0.4$ to $0.6$. The trends of these eigenvalues do not show an evident dependence on this quantity. The only feature that distinguishes the three conditions corresponds to the values of $Fr_{w}$ that delimit the region of the loss of hyperbolicity. By decreasing the water content in the mixture, both the limits of $Fr_{w}$ tend to move towards smaller values and the region of the loss of hyperbolicity tends to be enlarged (table 3). These results are obtained considering $Fr_{s}/Fr_{w}=0.7$, but similar results can be detected also with ratios greater than one.

Figure 7. Effect of the water concentration $C_{w}$ on the dimensionless eigenvalue trends. Results are obtained considering $Fr_{s}/Fr_{w} =0.7$, $C_{r}^{e}=0.34$, $K=1$ and $\rho _{i}/\rho _{r}=0.353$.

Table 3. Effects of the parameters $C_{w}$, $C_{r}^{e}$ and $K$ on the lower limit $Fr_{w}^{low}$, upper limit $Fr_{w}^{up}$ and range of the non-hyperbolic region. The values are referred to the different plots reported in the previous paragraphs and can be derived imposing both $Fr_{s} /Fr_{w}=0.7$ and $Fr_{s}/Fr_{w}=1.3$.

6.1.3. Effect of $C_{r}^{e}$

A change in the rock equivalent concentration keeping a fixed value for the water content (or equivalently, the solid concentration) translates into a change in the composition of the solid phase. Low values of $C_{r}^{e}$ imply that the overall solid phase is composed mainly of ice, while high values of $C_{r}^{e}$ implicate that the rock is predominant. As shown in figure 8, while three eigenvalues are not affected by a change in $C_{r}^{e}$, the lower limit of the region where the loss of hyperbolicity occurs moves significantly towards low values of $Fr_{w}$ by decreasing the rock equivalent concentration. Conversely, the upper limit of this region tends to move towards greater values of $Fr_{w}$ by reducing the rock equivalent concentration (table 3). Consequently, it seems that an increase in the ice content expands the region of the loss of hyperbolicity. Also in this case, a similar result can be detected for a Froude ratio greater than one.

Figure 8. Effect of the rock equivalent concentration $C_{r}^{e}$ on the dimensionless eigenvalue trends. Results are obtained considering $Fr_{s}/Fr_{w}=0.7$, $C_{w}=0.4$, $K=1$ and ${\rho }_{i}/\rho _{r}=0.353$.

6.1.4. Effect of $K$

As shown in figure 9, an increase in the earth-pressure coefficient $K$ does not produce evident changes in the trends of the dimensionless eigenvalues, except for the values of $Fr_{w}$ that delimit the region of the loss of hyperbolicity. By increasing the earth-pressure coefficient from an extensive ($K<1$) to a compressive ($K>1$) condition, the region of the loss of hyperbolicity tends to move towards higher values of $Fr_{w}$.

Figure 9. Effect of the earth-pressure coefficient $K$ on the dimensionless eigenvalue trends. Results are obtained considering $Fr_{s}/Fr_{w}=0.7$, $C_{w}=0.5$, $C^{e} _{r}=0.37$ and $\rho _{i}/\rho _{r}=0.353$.

6.2. Comparison with the other simplified 1-D models

In this section, we compare the eigenvalues related to the PI-RIW model with those referred to the FI-RIW and TP-RIW approaches with the aim to detect how the hypotheses introduced in § 4 affect the mathematical nature of the PDE systems. For the comparison of the eigenvalues of the PI-RIW model with those related to the ITP-RIW and MP-RIW approaches, we refer the reader to Appendix D.

6.2.1. FI-RIW versus PI-RIW

The characteristic polynomial of the FI-RIW model is

(6.12)\begin{equation} P(\lambda^{FI})= ( \lambda^{FI}-U )^{2} ( ( \lambda^{FI} )^{2}-2U\lambda^{FI}+U^{2}-Ka^{2} ),\end{equation}

where $U$ is the mixture velocity and $a$ is given by (6.10a). The corresponding eigenvalues are

(6.13a)$$\begin{gather} \lambda_{1}^{FI}=U-\sqrt{Ka^{2}}, \end{gather}$$

(6.13b)$$\begin{gather}\lambda_{2, 3}^{FI}=U, \end{gather}$$

(6.13c)$$\begin{gather}\lambda_{4}^{FI}=U+\sqrt{Ka^{2}}. \end{gather}$$

By defining the mixture Froude number as $Fr=U/a$, the dimensionless eigenvalues can be derived by dividing (6.13) by $a$. As pointed out in (6.13), the dimensionless eigenvalues depend only on the mixture Froude number $Fr$ and on the earth pressure coefficient $K$. In case of $K=1$, these eigenvalues coincide with those related to a free-surface water problem with two advected passive scalars.

In figure 10 the trends of the dimensionless eigenvalues as a function of $Fr$ are plotted and compared with the PI-RIW eigenvalues, where we have set $Fr_{w}=Fr$. It can be noticed that, unlike the PI-RIW model, the FI-RIW approach is always hyperbolic. The three dimensionless eigenvalues obviously coincide with those of the PI-RIW model for $Fr=Fr_{s}=Fr_{w}$, while they deviate significantly in the other cases. For $Fr_{s}/Fr_{w}<1$, the highest FI-RIW eigenvalue approximates the highest PI-RIW one, while for $Fr_{s}/Fr_{w}>1$ the lowest FI-RIW eigenvalue approximates the lowest PI-RIW one. The remaining intermediate eigenvalues are quite different from those related to the PI-RIW model for all values of $Fr_{w}$.

Figure 10. Comparison between the dimensionless eigenvalues of the FI-RIW and of the PI-RIW models. Parameters used in these plots are the same used in figure 5. The FI-RIW eigenvalues are marked by continuous lines with symbols, while those of the PI-RIW model are denoted by dashed lines without any type of distinction (see figure 5 for details in the PI-RIW eigenvalues).

6.2.2. TP-RIW versus PI-RIW

The characteristic polynomial of the TP-RIW model reads as follows:

(6.14)\begin{equation} P( \lambda^{TP}) =\sum_{l=0}^{4}b_{l}( \lambda^{TP})^{l}. \end{equation}

The coefficients $b_{l}$ are defined as in the PI-RIW model by (6.9) where $\beta ^2$ and $\delta ^2$ are substituted by the following expressions:

(6.15) \begin{equation} \left.\begin{gathered} \beta_{TP}^2={\left[ C_{w}-1-\dfrac{K}{2}( C_{w}-2) \right] \left( 1-\dfrac{\rho_{w}}{\rho_{r}}\right) }\\ \delta_{TP}^2={( 1-C_{w}) ( K-1) \left( 1-\dfrac{\rho_{w}}{\rho_{r}}\right) } \end{gathered}\right\}. \end{equation}

Imposing $K=1$, the coefficients of the characteristic polynomial related to the TP-RIW model coincide with those derived by Pelanti et al. (Reference Pelanti, Bouchut and Mangeney2008).

As performed for the PI-RIW model, the dimensionless eigenvalues can be plotted as a function of $Fr_{w}$. As shown in figure 11, the TP-RIW eigenvalues have the same trend as the dimensionless eigenvalues of the PI-RIW model. The only difference between the two approaches lies in the fact that the TP-RIW model does not show the eigenvalue $\hat {\lambda }=Fr_{s}$. This overlap attests that the TP-RIW approach is not only a consistent simplification of the PI-RIW model, but also that the PI-RIW approach converges to the two-phase model when the ice component disappears from the mixture (complete melting).

Figure 11. Comparison between the dimensionless eigenvalues of the TP-RIW and of the PI-RIW models. Parameters used in these plots are the same used in figure 5. The TP-RIW eigenvalues are marked by continuous lines with symbols, while those of the PI-RIW model are denoted by dashed lines without any type of distinction (see figure 5 for details in the PI-RIW eigenvalues).