2.1. Pitman and Le’s model

The assumption of shallow flow, used in deriving (2.3a ) and (2.3b ), may also be used to infer from the slope-normal component of (2.1b ) that at leading order the normal stresses are in equilibrium with the interphase forces and gravity:

(2.11) \begin{equation} \frac {\partial \sigma _i^{zz}}{\partial z} = -f_i^z + \rho _i \varphi _i g^z. \end{equation}

In deriving their debris-flow model, Pitman & Le (Reference Pitman and Le2005) use this to obtain expressions for the stresses. The fluid tensor is assumed to be isotropic and the slope-normal interphase forces are considered to be dominated by buoyancy, so $d_s^z = 0$ and from (2.5), $f_s^z = -\varphi _s \partial p / \partial z = -f_f^z$ . Substituting this into (2.11), depth-integrating twice and using (2.6), gives

(2.12a,b) \begin{align}\sigma _f^{zz} = -\rho _f g^z(h-z) \quad \mathrm{and}\quad \overline {\sigma _f^{xx}}=\overline {\sigma _f^{zz}} = -\scriptsize{\frac{1}{2}}\rho _f g^z h. \end{align}

Note that the direction of the buoyancy force and gravity coincide to make the effective stress for the fluid phase equal to the intrinsic pressure $p$ of the fluid. Conversely, for the solids phase, buoyancy acts against gravity to reduce the effective normal stress to

(2.13) \begin{equation} \sigma _s^{zz}= -\varphi _s (\rho _s-\rho _f) g^z (h-z). \end{equation}

Since the flow is anticipated to be densely packed with grains, principles of soil mechanics are invoked to infer a proportional relationship between lateral and normal stresses, via an Earth pressure coefficient $K$ :

(2.14) \begin{equation} \sigma _s^{xx} = K \sigma _s^{zz}. \end{equation}

On depth-averaging and using (2.14), one may therefore deduce that

(2.15) \begin{equation} -\frac {\partial \,}{\partial x}\left ( h\overline {\sigma _s^{xx}} \right ) = \frac {\partial \,}{\partial x}\left [ \frac {1}{2} K (1-\gamma ) \rho _s g^z \overline {\varphi _s} h^2 \right ]\!. \end{equation}

The model may be expressed in full by substituting (2.7), (2.12b ) and (2.15) into (2.3a ) and (2.3b ) and algebraically simplifying. It is convenient at this stage to define variables that express the proportion of the mixture depth occupied by each phase:

(2.16) \begin{equation} H_i = \overline {\varphi _i} h. \end{equation}

Using these variables and noting in particular that $h = H_s + H_f$ , the following set of equations are obtained:

(2.17a) \begin{align} \frac {\partial H_s}{\partial t} + \frac {\partial \,}{\partial x}(H_s \overline {u_s}) & = 0, \end{align}

(2.17b) \begin{align} \frac {\partial \overline {u_s}}{\partial t} + \overline {u_s}\frac {\partial \overline {u_s}}{\partial x} + g^z \left [ \gamma + K(1-\gamma )\left (\! 1 + \frac {H_f}{2H_s} \right ) \right ] \frac {\partial H_s}{\partial x} + g^z \left [ \gamma + \frac {K}{2} (1 - \gamma ) \right ] \frac {\partial H_f}{\partial x} = S_s, \end{align}

(2.17c) \begin{align} \frac {\partial H_f}{\partial t} + \frac {\partial \,}{\partial x}(H_f \overline {u_f}) & = 0, \end{align}

(2.17d) \begin{align} \frac {\partial \overline {u_f}}{\partial t} + \overline {u_f}\frac {\partial \overline {u_f}}{\partial x} + g^z \frac {\partial H_s}{\partial x} + g^z \frac {\partial H_f}{\partial x} & = S_f. \end{align}

The particular case of $K = 1$ was studied in detail by Pelanti et al. (Reference Pelanti, Bouchut and Mangeney2008).

2.2. Meng et al.’s model

The model of Meng et al. (Reference Meng, Johnson and Gray2022) is derived using a conceptually different description of the flow that posits separate free surfaces for the depth of solid particles $h_s$ and depth of fluid $h_f$ . When $h_f \gt h_s$ , the particles are ‘oversaturated’ with fluid and assumed to have settled into a layer at the bottom of the flow, within which they occupy a constant volume fraction $\varphi _c$ . We consider this case only, since the analysis of Meng et al. (Reference Meng, Taylor-Noonan, Johnson, Take, Bowman and Gray2024) (in their Appendix A) establishes that their model equations in the ‘undersaturated’ regime $h_f\,{\lt}\,h_s$ are hyperbolic, with a differential operator whose structure decouples into separate shallow-layer terms for each phase, thereby leading to well-posed initial-value problems.

The solids stresses take the same form as in the Pitman & Le (Reference Pitman and Le2005) model’s (2.13), except they are only present up to the height $h_s$ of the solids layer, implying that the term inside the pressure derivative of (2.15) differs by a factor of $h_s/h_f$ . Moreover, $K = 1$ is assumed. Therefore,

(2.18) \begin{equation} -\frac {\partial \,}{\partial x}\left ( h\overline {\sigma _s^{xx}} \right ) = \frac {\partial \,}{\partial x}\left [ \frac {1}{2} (1-\gamma ) \rho _s g^z \overline {\varphi _s} h_s h_f \right ]\!. \end{equation}

Additionally, the viscous component of the fluid stress tensor is retained. Therefore, rather than appealing to (2.11), the constitutive relation

(2.19) \begin{equation} {\boldsymbol{\sigma }}_f = -p {\unicode{x1D644}} + \varphi _f \eta _f \left [\nabla {\boldsymbol{u}}_f + (\nabla {\boldsymbol{u}}_f)^T\right ] \end{equation}

is proposed, where $\eta _f$ is the dynamic viscosity of the fluid. The intrinsic pore fluid pressure is hydrostatic as before, so (2.6) applies and consequently,

(2.20) \begin{equation} -\frac {\partial \,}{\partial x}\left ( h \overline {\sigma _f^{xx}} \right ) = \rho _f g^z h_f \frac {\partial h_f}{\partial x} - \frac {\partial \,}{\partial x}\left ( 2 \eta _f h_f \overline {\varphi _f} \frac {\partial \overline {u_f}}{\partial x} \right )\!. \end{equation}

To obtain the final term on the right, $\overline {\partial u_f/\partial x} \approx \partial \overline {u_f}/\partial x$ is used, which follows from an assumption of low shear in the velocity profile $\overline {u_f} \approx u_f(h)$ , and is consistent with the approximation made in (2.4).

Averaging the solids volume fraction over the full depth gives $\overline {\varphi _s} = \varphi _c h_s / h_f$ . This implies that the equivalent partial depths (2.16) in this model are

(2.21) \begin{equation} H_s = \varphi _c h_s, \quad H_f = h_f - \varphi _c h_s. \end{equation}

On making these transformations, the derivative terms in the Meng et al. (Reference Meng, Johnson and Gray2022) model equations are the same as the Pitman & Le (Reference Pitman and Le2005) model’s (2.17a – d ), save for the components related to the different formulations for internal stresses. Therefore, we report only the solid and fluid momentum equations, which may be obtained by substituting (2.18) and (2.20), along with the buoyancy forces (2.7), into (2.3b ), using (2.21) and simplifying, leading to

(2.22a) \begin{align} \frac {\partial \overline {u_s}}{\partial t} + \overline {u_s}\frac {\partial \overline {u_s}}{\partial x} + g^z \! \left [ \gamma + \frac {1-\gamma }{\varphi _c} \right ] \! \frac {\partial H_s}{\partial x} + g^z\gamma \frac {\partial H_f}{\partial x} = S_s, \end{align}

(2.22b) \begin{align} \frac {\partial \overline {u_f}}{\partial t} + \overline {u_f}\frac {\partial \overline {u_f}}{\partial x} + g^z \frac {\partial H_s}{\partial x} + g^z \frac {\partial H_f}{\partial x} = \frac {2\eta _f}{\rho _f H_f} \frac {\partial \,}{\partial x}\left ( H_f \frac {\partial \overline {u_f}}{\partial x} \right ) + S_f. \end{align}

A typical choice for the solids fraction constant in the regimes relevant to this model might be expected to lie somewhere in the range $0.5 \lesssim \varphi _c \lesssim 0.75$ (Pierson Reference Pierson1995). Nevertheless, it should be noted that in the limit $\varphi _c \to 1$ (where there are no saturated gaps between particles) and assuming also that $\eta _f = 0$ , (2.22a ) and (2.22b ) together with (2.17a ) and (2.17c ) reduce to a system of depth-averaged equations for the motion of two immiscible fluids of different densities, whose properties have been widely studied (see e.g. Ovsyannikov Reference Ovsyannikov1979; Vreugdenhil Reference Vreugdenhil1979; Castro et al. Reference Castro, Macías and Parés2001; Abgrall & Karni Reference Abgrall and Karni2009; Kurganov & Petrova Reference Kurganov and Petrova2009; Chiapolino & Saurel Reference Chiapolino and Saurel2018). A model of this latter type has also been proposed by Meyrat et al. (Reference Meyrat, McArdell, Ivanova, Müller and Bartelt2022), for use in debris-flow modelling.

2.3. Pudasaini’s model

Pudasaini (Reference Pudasaini2012) uses an approach that is consistent with Pitman & Le (Reference Pitman and Le2005), but extends their framework in various ways. Of relevance to our analysis are the inclusion of the added mass term given previously in (2.9) and a fluid stress tensor that incorporates a non-Newtonian component.

The inclusion of added mass augments the inertial terms in the momentum equations. The coefficient $\overline {C}$ in (2.9) is assumed to be a constant. Furthermore, in order to simplify the conservative form of the equations Pudasaini (Reference Pudasaini2012) makes the assumption that $\overline {C'} \equiv \overline {C}\,\overline {\varphi _s}/\overline {\varphi _f}$ may be absorbed into the time and space derivatives of (2.9) without explicitly holding it constant. This does not appear to be justified in our view. Nevertheless, summing the added mass force terms for each phase with the corresponding inertial terms from (2.3b ) and converting to quasilinear form (i.e. by dividing through by $\rho _i H_i$ and simplifying) leads to

(2.23a) \begin{align} & (1+\gamma \overline {C})\left ( \frac {\partial \overline {u_s}}{\partial t} + \overline {u_s} \frac {\partial \overline {u_s}}{\partial x} \right ) -\gamma \overline {C}\left ( \frac {\partial \overline {u_f}}{\partial t} + \overline {u_f} \frac {\partial \overline {u_f}}{\partial x} \right ) -\underbrace {\frac {\gamma \overline {C}\,\overline {u_f}}{H_s} \left [ \frac {\partial H_{s}}{\partial t} + \frac {\partial \,}{\partial x}\left (H_s \overline {u_f}\right ) \right ]},\nonumber\\& \ \ \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad {}_{\,\,\mathrm{extra\,terms}}\\[-5pt]\nonumber \end{align}

(2.23b) \begin{align}& \left (1+ \frac {\overline {C}H_s}{H_f}\right ) \left ( \frac {\partial \overline {u_f}}{\partial t} + \overline {u_f}\frac {\partial \overline {u_f}}{\partial x} \right ) -\frac {\overline {C}H_s}{H_f}\left ( \frac {\partial \overline {u_s}}{\partial t} + \overline {u_s}\frac {\partial \overline {u_s}}{\partial x} \right ) +\overbrace {\frac {\overline {C}\,\overline {u_f}}{H_f} \left [ \frac {\partial H_s}{\partial t} + \frac {\partial \,}{\partial x}\left (H_s\overline {u_f}\right ) \right ]} \end{align}

for the inertia of the solids and fluid phases, respectively. The extra terms, highlighted by the braces, do not appear if (2.9) is depth-averaged directly and could arguably be omitted, since they correspond to a force between the phases whose physical origin is unclear. However, in order to analyse the model as it has appeared in prior publications, we retain them.

The assumed form of the fluid stress tensor is equal to the expression used by Meng et al. (Reference Meng, Johnson and Gray2022), given in (2.19), plus an additional phenomenological component

(2.24) \begin{equation} -\eta _f {\mathcal{A}}\left [ \nabla \varphi _s \otimes ({\boldsymbol{u}}_f - {\boldsymbol{u}}_s) + ({\boldsymbol{u}}_f - {\boldsymbol{u}}_s) \otimes \nabla \varphi _s \right ], \end{equation}

where ${\mathcal{A}}$ is a parameter that depends on the solids fraction. After adding on the Newtonian component, depth-averaging $\sigma _f^{xx}$ gives

(2.25) \begin{equation} -\frac {\partial \,}{\partial x}\left ( h \overline {\sigma _f^{xx}} \right ) = \rho _f g^z h \frac {\partial h}{\partial x} - \frac {\partial \,}{\partial x}\left ( 2 \eta _f h \overline {\varphi _f} \frac {\partial \overline {u_f}}{\partial x} - 2\eta _f \overline {{\mathcal{A}}} h (\overline {u_f}-\overline {u_s})\frac {\partial \overline {\varphi _s}}{\partial x} \right ) \end{equation}

for this model, where we have used $\overline {\partial \varphi _s/\partial x} \approx \partial \overline {\varphi _s}/\partial x$ , which is consistent with the assumption of negligible variation in volume fraction over the depth, $\overline {\varphi _s} \approx \varphi _s(h)$ . In the original derivation, Pudasaini (Reference Pudasaini2012) goes further, following an approach of Iverson & Denlinger (Reference Iverson and Denlinger2001) for averaging diffusive stresses by bringing $h$ outside the spatial derivatives of (2.25). This introduces extra terms, which, under the stress-free boundary condition, reduce to expressions that do not contain derivatives and may be modelled separately as source terms (Pudasaini Reference Pudasaini2012). These extra steps do not affect the forthcoming analysis of the model structure (since the linearised diffusion operator remains the same). Therefore, we leave (2.25) as it is.

3.1. Two-phase models

Given some putative model solution with fields $\boldsymbol{q} = (H_s, \overline {u_s}, H_f, \overline {u_f})^T$ , we would like to understand the local behaviour of the governing equations at an arbitrary space–time location $(x_0,t_0)$ . Denote a state vector there by

(3.1) \begin{equation} \boldsymbol{q}_0 = \boldsymbol{q}(x_0, t_0) = \left ( H_{s}^{(0)}\!\!,\,\overline {u_s}^{(0)}\!,\,H_{f}^{(0)}\!\!,\,\overline {u_f}^{(0)} \right )^T. \end{equation}

We assume non-vanishing fluid depth $H_f^{(0)} \gt 0$ and velocity $\overline {u_f}^{(0)} \neq 0$ , so that the governing equations for each model may be non-dimensionalised with respect to these scales. States may then be fully characterised by three dimensionless quantities:

(3.2a–c) \begin{equation} R_H = H_s^{(0)} / H_f^{(0)}, \quad R_u = \overline {u_s}^{(0)} / \overline {u_f}^{(0)}, \quad {\textit{Fr}} = \frac {\overline {u_f}^{(0)}}{\sqrt {g^z H_f^{(0)}}}, \end{equation}

where $\textit{Fr}$ is the local Froude number for the fluid phase. Therefore, hereafter the transformations

(3.3a–e) \begin{eqnarray} && x\mapsto x / H_f^{(0)}, \quad t \mapsto t \overline {u_f}^{(0)} / H_f^{(0)}, \quad H_i \mapsto H_i / H_f^{(0)}, \quad \overline {u_i} \mapsto \overline {u_i} / \overline {u_f}^{(0)}, \nonumber\\[4pt]&& \qquad\qquad\qquad\qquad\qquad S_i \mapsto S_i H_f^{(0)} / (\overline {u_f}^{(0)})^2 \end{eqnarray}

are made to the two-phase models analysed. Furthermore, note that the systems detailed in §§ 2.1–2.3 may be collectively cast in the general form

(3.4) \begin{equation} {\unicode{x1D63C}}(\boldsymbol{q}) \frac {\partial \boldsymbol{q}}{\partial t} + {\unicode{x1D63D}}(\boldsymbol{q}) \frac {\partial \boldsymbol{q}}{\partial x} = \boldsymbol{S}(\boldsymbol{q}) + {\unicode{x1D63F}}_1(\boldsymbol{q})\frac {\partial \,}{\partial x}\left ( {\unicode{x1D63F}}_2(\boldsymbol{q}) \frac {\partial \boldsymbol{q}}{\partial x} \right ), \end{equation}

where $\boldsymbol{S}(\boldsymbol{q}) = (0, S_s, 0, S_f)^T$ and $\unicode{x1D63C}$ , $\unicode{x1D63D}$ , ${\unicode{x1D63F}}_1$ , ${\unicode{x1D63F}}_2$ are matrices of (dimensionless) variable coefficients that that may be readily specified for each model.

We now ‘freeze’ the solution, by assuming that $\boldsymbol{q}(x,t) = \boldsymbol{q}_0$ within a local neighbourhood of $(x_0, t_0)$ and consider the evolution of a normal-mode perturbation $\boldsymbol{r} \exp ({\mathrm{i}} k x + \sigma t)$ to this state, where $k$ is a real-valued wavenumber, $\sigma$ a complex growth rate and $\boldsymbol{r}$ a vector constant with $|\boldsymbol{r}| \ll |\boldsymbol{q}_0|$ . Linearising (3.4) around the frozen base state $\boldsymbol{q}_0$ leads to the following eigenproblem for the pair $(\sigma , \boldsymbol{r}$ ):

(3.5) \begin{gather} \sigma {\unicode{x1D63C}}(\boldsymbol{q}_0)\boldsymbol{r} + {\mathrm{i}} k {\unicode{x1D63D}}(\boldsymbol{q}_0) \boldsymbol{r} = {\unicode{x1D63E}}(\boldsymbol{q}_0)\boldsymbol{r} - k^2 {\unicode{x1D63F}}(\boldsymbol{q}_0) \boldsymbol{r}, \end{gather}

where ${\unicode{x1D63E}} \equiv \partial \boldsymbol{S} / \partial \boldsymbol{q}$ and ${\unicode{x1D63F}} = {\unicode{x1D63F}}_1 {\unicode{x1D63F}}_2$ . If $\boldsymbol{q}_0$ happens to represent a state for which the model equations admit steady uniform flow, i.e. $\boldsymbol{S}(\boldsymbol{q}_0)=\boldsymbol{0}$ , then the solutions to (3.5) dictate the linear stability of such a flow, for which $\textrm {Re}(\sigma ) \gt 0$ indicates an unstable mode and the case of Hadamard instability occurs if $\textrm {Re}(\sigma ) \to \infty$ as $k\to \infty$ , indicating that the governing equations (3.4) are ill posed at $\boldsymbol{q}_0$ . Otherwise, the procedure of freezing the base state is justified insofar as it may be used to identify this latter pathology on the grounds that any candidate solution $\boldsymbol{q}$ must be effectively constant near $(x_0,t_0)$ , when measured with respect to the infinitesimal length and time scales over which the Hadamard instability develops (Joseph & Saut Reference Joseph and Saut1990; Joseph Reference Joseph1990).

Since diffusion can be small, relative to other terms in the equations, it is sometimes desirable to neglect its effects. Therefore, we first consider the case where $\unicode{x1D63F}$ is the zero matrix. It may then straightforwardly be determined from (3.5) that in the asymptotic limit of high $k$ , the leading-order components of the growth rates for the four eigenmodes are $\sigma = -{\mathrm{i}} k \lambda _j$ , where $\lambda _j$ ( $j=1,\ldots ,4$ ) denote the characteristic wave speeds of the problem, given by the solutions to the generalised eigenproblem ${\unicode{x1D63D}}\boldsymbol{r} = \lambda _j {\unicode{x1D63C}}\boldsymbol{r}$ . If all four are real and distinct then the system is said to be ‘strictly hyperbolic’ at $\boldsymbol{q}_0$ and is well posed as an initial-value problem. On the other hand, any complex characteristics must arise in conjugate pairs. Since one of the pair must have $\textrm {Im}(\lambda _j)\gt 0$ , the corresponding real part of $\sigma$ is positive and scales as $O(k)$ for $k \gg 1$ , giving rise to a Hadamard instability. Repeated real characteristics can also lead to growth rate blow-up, but the reasons for this are more subtle. This case is covered later in § 4.1.

3.1.1. Emergence of ill posedness

The inclusion of added mass leads to complications, which we address shortly, in §§ 3.1.2 and 4.3.2. If it is neglected, then $\unicode{x1D63C}$ simplifies to the identity matrix $\unicode{x1D644}$ and the problem reduces to computing the eigenvalues of the Jacobian $\unicode{x1D63D}$ , which has the same essential form for each of the models. Bearing in mind our transformation to dimensionless variables in (3.3a – e ), this matrix is

(3.6) \begin{gather} {\unicode{x1D63D}}(\boldsymbol{q}) = \begin{pmatrix} \overline {u_s} & H_s & 0 & 0 \\[4pt] (\gamma + \beta _1){\textit{Fr}}^{-2} & \overline {u_s} & (\gamma + \beta _2){\textit{Fr}}^{-2} & 0 \\[4pt] 0 & 0 & \overline {u_f} & H_f \\[4pt] {\textit{Fr}}^{-2} & 0 & {\textit{Fr}}^{-2} & \overline {u_f} \end{pmatrix} , \end{gather}

where $\beta _1 = K(1-\gamma )[1 + H_f/(2H_s)]$ , $\beta _2 = K(1-\gamma )/2$ for Pitman & Le (Reference Pitman and Le2005) and Pudasaini (Reference Pudasaini2012); $\beta _1 = (1-\gamma )/\varphi _c$ , $\beta _2 = 0$ in Meng et al. (Reference Meng, Johnson and Gray2022); and $\beta _1 = 1-\gamma$ , $\beta _2 = 0$ for two-fluid models (e.g. Ovsyannikov Reference Ovsyannikov1979), as well as the debris-flow model of Meyrat et al. (Reference Meyrat, McArdell, Ivanova, Müller and Bartelt2022). Note that at $\boldsymbol{q}=\boldsymbol{q}_0$ , $H_s = R_H$ , $\overline {u_s} = R_u$ and $H_f = \overline {u_f} = 1$ .

The possibility for ${\unicode{x1D63D}}(\boldsymbol{q}_0)$ to have complex characteristics arises due to the coupling between the momentum equations provided by the entries ${B}_{23}$ and ${B}_{41}$ . Physically, these terms arise because the buoyancy and solids stresses depend on the total depth $H_s + H_f$ . For systems without this coupling, i.e. ${B}_{23} = {B}_{41}=0$ , the eigenvalues of ${\unicode{x1D63D}}(\boldsymbol{q}_0)$ are

(3.7a,b) \begin{equation} \lambda _s^\pm \equiv R_u \pm \frac {\sqrt {R_H(\gamma + \beta _1)}}{{\textit{Fr}}} \quad \mathrm{and}\quad \lambda _f^\pm \equiv 1 \pm \frac {1}{{\textit{Fr}}}. \end{equation}

These are real provided $\beta _1 + \gamma \gt 0$ . For the model closures described above, this is certainly the case, since both $\beta _1$ and $\gamma$ are strictly positive. While the corresponding expressions for the eigenvalues of ${\unicode{x1D63D}}(\boldsymbol{q}_0)$ in the general case, ${B}_{23},{B}_{41}\neq 0$ , can be computed via the quartic formula, these are are too complicated to be especially useful (Pitman & Le Reference Pitman and Le2005; Pelanti et al. Reference Pelanti, Bouchut and Mangeney2008; Pudasaini Reference Pudasaini2012). Nevertheless, since ${\unicode{x1D63D}}(\boldsymbol{q}_0)$ is almost block diagonal, its characteristic polynomial is amenable to further analysis.

In particular, one can generalise an approach followed by Ovsyannikov (Reference Ovsyannikov1979) for the simpler two-fluid case ( $\beta _1 = 1-\gamma$ , $\beta _2 = 0$ ) and notice that the eigenvalues $\lambda _j$ are determined by an equation of the form

(3.8) \begin{equation} f(P_1, P_2) \equiv (P_1^2-1)(P_2^2-1) = c, \end{equation}

where

(3.9a–c) \begin{align} P_1^2 = \frac {(\lambda _j - R_u)^2{\textit{Fr}}^2}{R_H(\gamma + \beta _1)},\quad P_2^2 = (\lambda _j - 1)^2{\textit{Fr}}^2\quad \mathrm{and}\quad c = \frac {\gamma + \beta _2}{\gamma + \beta _1}.\end{align}

For a particular point in parameter space, characterised by the triple $(R_H, R_u, {\textit{Fr}}\,)$ , we can eliminate $\lambda _j$ from (3.9a , b ) to determine that the characteristics lie on the intersection of the line

(3.10) \begin{equation} P_2 = P_1\sqrt {R_H(\gamma + \beta _1)} + {\textit{Fr}}(R_u - 1) \end{equation}

with the level set given by the contour of the surface $f(P_1,P_2)$ (3.8) at the value $c$ . This is depicted graphically in figure 2(a). Coloured contours in the figure show the surface $f(P_1,P_2)$ , with an example level at $c = 1/3$ given by the solid black curves. Three dash–dot black lines illustrate possibilities for the characteristics. The line labelled I represents a strictly hyperbolic case, since it possesses four distinct intersections with the solid black contour. On shifting the line upwards to II (by increasing ${\textit{Fr}}(R_u-1)$ ) the characteristics associated with the central contour merge to form a complex conjugate pair and only two real solutions to (3.8) and (3.10) remain. Shifting the line further up recovers strict hyperbolicity, since at position III, it makes two additional intersections with the portion of the level set that is confined to $\{(P_1, P_2) : P_1 \lt -1, P_2 \gt 1\}$ . Provided that $c\gt 0$ and that $\beta _1$ , $\beta _2$ are either constants or a functions of $R_H$ only, as is the case for the models considered herein, we can see that there will always be an ill-posed region associated with the loss of strict hyperbolicity (i.e. regions without four distinct real eigenvalues). This is because a given $R_H$ fixes the level set determined by $c$ . Then, varying ${\textit{Fr}}(R_u - 1)$ shifts the dash–dotted lines in the $P_2$ direction, guaranteeing that they pass through a region with only two intersections. Indeed, by symmetry, there must be two such regions.

Figure 2. Geometric analysis of the characteristics for two-phase models. Filled contours of the surface $f(P_1,P_2)$ are plotted, spaced at intervals $\pm 10^m$ for $m = 0, \ldots , 4$ . The zero contour is marked separately (white dashed), as is the level set at $c = 1/3$ (black solid). Dash–dotted lines are $P_2 = P_1 - 0.75 + 2n$ , for $n = 0,1,2$ .

This framework encapsulates the analysis by Pitman & Le (Reference Pitman and Le2005) who showed for their model that cases close to $R_u = 1$ are always strictly hyperbolic. This is a consequence of the fact that the (3.10) lines pass through the origin at this point. Moreover, Pelanti et al. (Reference Pelanti, Bouchut and Mangeney2008) later gave bounds on $|R_u - 1|$ that guarantee well posedness for sufficiently small and sufficiently large values.

The white dashed lines in figure 2 show the level set contours at $c=0$ , given by $P_1 = \pm 1$ , $P_2 = \pm 1$ . In this special case, the characteristics are everywhere real and it is straightforward to see that they must be the same as the values for uncoupled systems, given in (3.7a , b ). Moreover, each point $(\pm 1, \pm 1)$ may be linked to one of the four possible intersections between the solid ( $\lambda _s^\pm$ ) and fluid ( $\lambda _f^\pm$ ) characteristics. For example, let $R_u$ , $R_H$ be fixed and suppose that $R_u \gt 1$ , implying that the dash–dotted lines of figure 2 (3.10) intercept the $P_2$ axis at positive values. From examination of the expressions for the characteristics in (3.7a , b ), it may determined that there always exists a $\textit{Fr}$ such that $\lambda _s^- = \lambda _f^+$ in this case. Moreover, depending on whether the gradients of the dash–dot lines $\sqrt {R_H(\gamma +\beta _1)}$ are greater or less than unity, the respective intersections $\lambda _s^-=\lambda _f^-$ and $\lambda _s^+ = \lambda _f^+$ are possible. By considering (geometrically) the corresponding options for the lines to pass through $(\pm 1, \pm 1)$ in this case, we infer that $\lambda _s^- = \lambda _f^+$ corresponds to the point $(-1, 1)$ and likewise that $\lambda _s^- = \lambda _f^-$ corresponds to $(-1, -1)$ and $\lambda _s^+=\lambda _f^+$ to $(1,1)$ . Symmetric reasoning for the case $R_u \lt 1$ determines the final intersection, $\lambda _s^+ = \lambda _f^-$ at $(1, -1)$ . The important points are $(-1, 1)$ and $(1, -1)$ , when the positive and negative branches coincide. This occurs when

(3.11) \begin{equation} R_u = 1 \pm \frac {1}{{\textit{Fr}}}\left [ 1 + \sqrt {R_H(\gamma + \beta _1)} \right ]\!. \end{equation}

It is from these intersections that the complex eigenvalues of ${\unicode{x1D63D}}(\boldsymbol{q}_0)$ emerge when the system is fully coupled. Therefore, the consequent blow-up in growth rate in these regions can be thought of as stemming from a resonant interaction between the characteristic wave speeds of the solid and fluid phases.

In figures 3(a) and 3(b), we plot the regions where ill posedness occurs for the Pitman & Le (Reference Pitman and Le2005) and Meng et al. (Reference Meng, Johnson and Gray2022) models, respectively (without diffusion), in terms of $R_H$ and $|{\textit{Fr}}(R_u-1)|$ . Note that these two parameters fully determine whether the characteristics are real-valued or not. As already inferred from the geometric analysis, the models are unconditionally well posed when $R_u = 1$ and at sufficiently high values of $|{\textit{Fr}}(R_u-1)|$ . Furthermore, the bands of ill posedness are organised around the condition in (3.11) (black dashed lines). The width of the bands is contingent on the model parameters, which select the level set(s) in figure 2, and the qualitative differences in shape between the bands for the two models are explained by the different dependence. Specifically, while the level set value $c$ for the Meng et al. (Reference Meng, Johnson and Gray2022) model is fixed, for Pitman & Le (Reference Pitman and Le2005), $c\equiv c(R_H)$ with $c\to 0$ as $R_H \to 0$ . This implies that the width of the figure 3(a) band approaches zero in this limit. Conversely, when $c$ is constant, the band has a finite width as $R_H\to 0$ , determined by the minimum distance between the central piece of the figure 2 level set and any of the lines $P_1, P_2 = \pm 1$ , which are asymptotically approached by the other sections of the level set. A brief calculation shows that this is $1-\sqrt {1-c}$ for $c \in [0, 1]$ , or $1$ for $c \gt 1$ (where in this latter case there is no central piece of the level set). Additionally, the uncoupled characteristic intersections must lie at the upper limit of the band as $R_H \to 0$ . Combining these observations with (3.11) determines that the interval $(\sqrt {1-c}, 1)$ remains ill posed in this case, as $R_H \to 0$ . This property is demonstrated for the Meng et al. (Reference Meng, Johnson and Gray2022) model with $c = 1/3$ , by examining figure 3(b).

Figure 3. Regions of parameter space which contain complex characteristics, indicated by the red shaded regions, for the models of (a) Pitman & Le (Reference Pitman and Le2005) with $\gamma = 0.5$ , $K=1$ and (b) Meng et al. (Reference Meng, Johnson and Gray2022) with $\gamma = 0.5$ , $\varphi _c = 0.5$ . In the case of (b), the parameter choices correspond to the solid black level set in figure 2. Outside the shaded regions, the characteristics are real and distinct. The black dashed curves are where positive and negative branches of characteristics from the corresponding uncoupled problems intersect, as given in (3.11).

3.1.2. Added mass effect

When the added mass effect is included in the two-phase model of Pudasaini (Reference Pudasaini2012), many additional terms are introduced that cause the equations to be more strongly coupled. Though this model also contains diffusive terms, it is informative to investigate first how the incorporation of this additional physics affects the model’s eigenstructure in the absence of diffusion. One reason for this is that, at least in some cases, ill posedness in some non-depth-averaged two-phase flow systems without diffusion can be regularised by including added mass terms (Drew Reference Drew1983).

On introducing added mass effects by generalising the inertial terms in the solids and fluid momentum equations to the expressions given previously in (2.23a ) and (2.23b ), respectively, the matrices $\unicode{x1D63C}$ and $\unicode{x1D63D}$ become

(3.12) \begin{align} {\unicode{x1D63C}} &= \begin{pmatrix} 1 & 0 & 0 & 0 \\[4pt] -\gamma \overline {C}\,\overline {u_f}H_s^{-1} & 1+\gamma \overline {C} & 0 & -\gamma \overline {C} \\[4pt] 0 & 0 & 1 & 0 \\[4pt] \overline {C}\,\overline {u_f}H_f^{-1} & -\overline {C}H_sH_f^{-1} & 0 & 1+\overline {C}H_sH_f^{-1} \end{pmatrix},\end{align}

(3.13) \begin{align} {\unicode{x1D63D}} &= \begin{pmatrix} \overline {u_s} & H_s & 0 & 0 \\[4pt] (\gamma + \beta _1){\textit{Fr}}^{-2}-\frac {\gamma \overline {C}\,\overline {u_f}^2}{H_s} & (1+\gamma \overline {C})\overline {u_s} & (\gamma + \beta _2){\textit{Fr}}^{-2} & -2\gamma \overline {C}\,\overline {u_f} \\[4pt] 0 & 0 & \overline {u_f} & H_f \\[4pt] {\textit{Fr}}^{-2}+\overline {C}\,\overline {u_f}^2H_f^{-1} & -\overline {C}\,\overline {u_s}H_sH_f^{-1} & {\textit{Fr}}^{-2} & \overline {u_f} + \frac {2 \overline {C}\,\overline {u_f}H_s}{H_f} \end{pmatrix} .\end{align}

When the added mass coefficient $\overline {C}$ is non-zero, the corresponding characteristic polynomial $p_c(\lambda ) = \det ({\unicode{x1D63D}}-\lambda {\unicode{x1D63C}})$ for this system lacks the advantageous structure that was leveraged in the previous section to analyse the eigenvalues geometrically. Nevertheless, they are straightforward to obtain numerically at any point in parameter space. On doing so, it was found (as might well be expected) that the boundaries of the well-posed regions do not collapse neatly onto curves in terms of the parameters $R_H$ and ${\textit{Fr}}(R_u-1)$ , as before. However, it is possible to observe the qualitative effect of increasing $\overline {C}$ from zero.

The plots in figure 4 show an illustrative example, in which $\gamma = 0.5$ , $R_H = 1$ and $\overline {C}$ is incremented up to the value of $0.5$ suggested by Pudasaini (Reference Pudasaini2012). When $\overline {C}=0$ , the system reduces to the structure of the Pitman & Le (Reference Pitman and Le2005) model. The ill-posed regions lie either side of $R_u = 1$ and take the form of bands around the curves given previously in (3.11). Increasing $\overline {C}$ to $0.1$ results in a slight narrowing of the ‘upper’ band with $R_u \gt 1$ and a slight thickening of the lower $R_u \lt 1$ band. Additionally, a new region of complex characteristics emerges beneath the lower band at higher Froude numbers. This region extends further towards lower $\textit{Fr}$ when $\overline {C}$ is increased to $0.5$ (figure 4 c), leaving most of the $R_u \lt 1$ half-plane ill posed. Furthermore, the upper band separates into two pieces, leaving a well-posed region in between them. Moreover, a small region of ill posedness appears at an $R_u$ closer to unity. It is approximately centred around the point $({\textit{Fr}}, R_u) = (3.04,1.05)$ . We inspected equivalent plots for other choices of $\gamma$ and $R_H$ in the ranges $0.3\lt \gamma \lt 0.8$ , $0.2\lt R_H\lt 1.5$ and found them to be qualitatively similar.

These results indicate that the added mass force in this case does little to ameliorate the problem of ill posedness on its own and arguably seems to make matters worse, especially when the fluids velocity greatly exceeds the solids velocity ( $R_u\lt 1$ ) – a situation which could be encountered when a less concentrated debris flow entrains a static pile of grains, for example. Some analytical insight into the emergence of the large ill-posed region for $R_u \lt 1$ is gained later in § 4.3.2.

Figure 4. Effect of added mass terms in the Pudasaini (Reference Pudasaini2012) model without diffusion. Regions of parameter space that possess complex characteristics are shaded red, for $R_H = 1$ , $\gamma = 0.5$ and $\overline {C} =$ (a) $0$ , (b) $0.1$ and (c) $0.5$ .

3.2. Three-phase models

For three-phase models that share the same essential structure as the two-phase models we have analysed, it is possible to generalise the geometric reasoning of § 3.1.1 to identify regions of parameter space that must contain complex characteristics, provided the added mass effect is negligible.

Therefore, we return to the case $\overline {C}=0$ , ${\unicode{x1D63C}}={\unicode{x1D644}}$ and consider models that possess a Jacobian of the form

(3.14) \begin{equation} {\unicode{x1D63D}} = \begin{pmatrix} \overline {u_1} & H_1 & 0 & 0 & 0 & 0 \\[4pt] \beta _{11} & \overline {u_1} & \beta _{12} & 0 & \beta _{13} & 0 \\[4pt] 0 & 0 & \overline {u_2} & H_2 & 0 & 0 \\[4pt] \beta _{21} & 0 & \beta _{22} & \overline {u_2} & \beta _{23} & 0 \\[4pt] 0 & 0 & 0 & 0 & \overline {u_3} & H_3 \\[4pt] \beta _{31} & 0 & \beta _{32} & 0 & \beta _{33} & \overline {u_3} \end{pmatrix}, \end{equation}

where $H_i$ denote partial heights for each phase, $\overline {u_i}$ the corresponding downstream velocities and $\beta _{ij}$ represent arbitrary functions of these flow variables. This generalises the essential structure of the two-phase Jacobian in (3.6) to three phases. Denote the characteristic polynomial of this matrix by $p_c$ . By direct computation, it may be shown that $p_c(\lambda ) = 0$ simplifies to

(3.15) \begin{align} f(P_1,P_2,P_3) &\equiv \prod _{i=1}^3 (P_i^2-1) - (P_1^2-1)\frac {\beta _{23}\beta _{32}}{\beta _{22}\beta _{33}} - (P_2^2-1)\frac {\beta _{13}\beta _{31}}{\beta _{11}\beta _{33}} \nonumber\\& \quad - (P_3^2-1)\frac {\beta _{12}\beta _{21}}{\beta _{11}\beta _{22}} = c, \end{align}

where

(3.16) \begin{equation} P^2_i = \frac {(\overline {u_i} - \lambda )^2}{\beta _{ii} H_i} \quad \mathrm{and} \quad c= \frac {\beta _{13}\beta _{21}\beta _{32} + \beta _{12}\beta _{23}\beta _{31}}{\beta _{11}\beta _{22}\beta _{33}}\gt 0, \end{equation}

for $i = 1,2,3$ . We retain the convention adopted previously, by non-dimensionalising with respect to the depth and velocity of the third phase, represented by the final two governing equations and hereafter assumed to represent the carrier fluid. There are now two pairs of relevant dimensionless quantities associated with the relative heights and velocities of the phases:

(3.17a–d) \begin{equation} R_{H_1} = H_1^{(0)} / H_3^{(0)}, \quad R_{H_2} = H_2^{(0)} / H_3^{(0)}, \quad R_{u_1} = \overline {u_1}^{(0)} / \overline {u_3}^{(0)}, \quad R_{u_2} = \overline {u_2}^{(0)} / \overline {u_3}^{(0)}, \end{equation}

alongside the Froude number for the fluid phase ${\textit{Fr}} = \overline {u_3}^{(0)}/\sqrt {g^z H_3^{(0)}}$ . On making the appropriate non-dimensionalising transformations and eliminating $\lambda$ from among the defining relations for the $P_i$ coordinates in (3.16), it may be concluded that the number of real roots of $p_c$ at a given point in parameter space is determined by the intersections of the level surface defined in (3.15) and (3.16), with the line given by the map

(3.18) \begin{equation} P_3 \mapsto \left ( \frac {R_{u_1}-1}{\sqrt {\beta _{11} R_{H_1}}} + P_3\sqrt {\frac {\beta _{33} }{\beta _{11} R_{H_1}}}, \frac {R_{u_2}-1}{\sqrt {\beta _{22} R_{H_2}}} + P_3\sqrt {\frac {\beta _{33} }{\beta _{22} R_{H_2}}}, P_3 \right )\!. \end{equation}

To illustrate the resulting geometric picture, we use the model of Pudasaini & Mergili (Reference Pudasaini and Mergili2019), which extends the two-phase system of Pudasaini (Reference Pudasaini2012) to incorporate an intermediate fraction of fine solid particles. When added mass effects are neglected, the Jacobian for this model matches the structure given in (3.14). If the equations are organised such that the first two rows denote the solid phase, the second two the fine-solid phase and the final two the fluid phase, then the (non-dimensionalised) $\beta _{ij}$ closure terms are

(3.19a) \begin{align} \beta _{11} &= \frac {1}{{\textit{Fr}}^2}\left [ 1 + \frac {1}{2}(1-\gamma _1)\left ( \frac {R_{H_1} + R_{H_2} + 1}{R_{H_1}} \right ) \right ],\\[-10pt]\nonumber\end{align}

(3.19b) \begin{align} \beta _{12} &= \beta _{13} = \frac {1}{2{\textit{Fr}}^2}(1+\gamma _1)\end{align}

and $\beta _{2i} = \gamma _2/{\textit{Fr}}^2$ , $\beta _{3i} = 1/{\textit{Fr}}^2$ , for $i = 1,2,3$ , where $\gamma _1$ is the ratio of fluid to solid densities and $\gamma _2$ is the ratio of fluid to fine-solid densities. These latter two parameters are fixed material constants. On substituting the expressions for $\beta _{ij}$ into (3.15), (3.16) and (3.18), it may be deduced that both the level surface $f(P_1,P_2,P_3)=c$ and the gradient of the line in (3.18) depend only on the flow via the relative heights $R_{H_1}$ and $R_{H_2}$ . Consequently, for a given $(R_{H_1}, R_{H_2})$ pair, the number of intersections between the line and the level set is determined by the remaining degrees of freedom for the line, namely the terms

(3.20a,b) \begin{equation} K_1 \equiv \frac {R_{u_1} - 1}{\sqrt {\beta _{11}R_{H_1}}} \quad \mathrm{and}\quad K_2 \equiv \frac {R_{u_2} - 1}{\sqrt {\beta _{22}R_{H_2}}}. \end{equation}

By substituting in the appropriate values for $\beta _{11}$ , $\beta _{22}$ , it may be seen that $K_i \propto {\textit{Fr}}(R_{u_i} - 1)$ .

In figure 5, we plot the surface corresponding to the case $R_{H_1} = R_{H_2} = 1$ and $\gamma _1 = \gamma _2 = 0.5$ . It consists of nine disjoint pieces, comprising eight surfaces in each corner octant, which we label $1$ – $8$ for later reference, and a central ‘cross-shaped’ surface. Far from the origin, the corner surfaces asymptote to the planes $P_i = \pm 1$ . This is a consequence of the more general property that in the limit $|P_i|\to \infty$ , (3.15) reduces to the two-dimensional level set corresponding to the equivalent two-phase problem with phase $i$ removed. This also explains the extended stems of the central cross, since slices of the surface in the far-field limits $|P_2|\to \infty$ and $|P_3|\to \infty$ (removing either of the fluid phases) may be compared with the two-phase level set in figure 2(a), with the stems of the cross giving rise to the closed curve around the origin.

Figure 5. The surface $f(P_1,P_2,P_3) = c$ , for the three-phase model of Pudasaini & Mergili (Reference Pudasaini and Mergili2019), with $\gamma _1 = \gamma _2 = 0.5$ and $R_{H_1} = R_{H_2} = 1$ . For visual clarity, the disjoint pieces of the surface are rendered with a triangular mesh and coloured from blue to red according to the value of the $P_3$ coordinate. Also plotted is the line defined by (3.18) for $R_{u_1}=R_{u_2}=1$ . This intersects with the surface at the four points marked with circles and at the origin (marked with a cross), which is an additional isolated solution of $f(P_1,P_2,P_3)=c$ , in this case. Supplementary movie 2 shows an animated view of the surface.

Also plotted in figure 5 is the corresponding line with $K_1 = K_2 = 0$ , which passes through $(0,0,0)$ . Using the fact that $\beta _{21}=\beta _{22}=\beta _{23}$ and $\beta _{31}=\beta _{32}=\beta _{33}$ for this model, it may be verified that the origin is an additional isolated point on the level set. Since the line also necessarily passes through the central cross surface and the corner surfaces $1$ and $8$ , it intersects with the level set five times in total. Therefore, this case corresponds to a repeated real root of $p_c$ . More generically, we should expect an even number of purely real eigenvalues, determined by the number of intersections of the plotted line through the origin after undertaking an appropriate translation in the $(P_1,P_2)$ plane by $(K_1, K_2)$ , depending on the values of ${\textit{Fr}}(R_{u_i}-1)$ at a given point in parameter space. The different possibilities are summarised in figure 6. In particular, figure 6(a) plots, as a function of $(K_1, K_2)$ , whether there are six real roots (white), four real roots (pink) or only two (red). Cases where there are repeated real roots are associated with either tangential intersections between the line and level set (the borders of each shaded region in figure 6 a) or isolated points such as the origin. Regions with complex eigenvalues cover a substantial part of the plane. Notably, the model is ill posed as an initial-value problem in this case for all $(K_1,K_2)\in [-1,1]^2$ , i.e. when

(3.21a) \begin{equation} |R_{u_1} - 1| \leqslant \sqrt {\beta _{11}R_{H_1}} \quad \mathrm{and} \quad |R_{u_2} - 1| \leqslant \sqrt {\beta _{22}R_{H_2}}. \end{equation}

From the geometric picture in figure 5, we see that this is because there are only four available intersections in these cases (excepting the special case $R_{u_1}=R_{u_2}=1$ , already discussed). This observation contrasts with the two-phase models analysed above, which are always well posed when $R_u$ is sufficiently close to unity.

Figure 6. (a) Regions of the $(K_1,K_2)$ plane, for which the ( $R_{H_1} = R_{H_2}=1$ ) surface geometry in figure 5 gives rise to six (white), four (pink) or two (red) real eigenvalues. (b–d) Intersections between the (3.18) line and the (3.15) level surface (solid lines), for $(K_1, K_2)$ values along the corresponding dashed lines plotted in (a). These are: (b) $K_2 = -K_1$ , (c) $K_2 = -K_1-6$ and (d) $K_2 = -K_1+4$ . The shaded bands indicate the number of intersections, in accordance with the colouring in (a). Labels denote regions enclosed by the numbered corner surfaces (see figure 5).

There are various other possibilities when $R_{u_1}$ or $R_{u_2}$ are large enough to lie outside the intervals in (3.21a ). The diagrams in figure 6(b–d) are useful for visualising them. These plots show how the intersections of the line and level set change as the line is translated along different trajectories in the $(K_1,K_2)$ plane, represented by dashed lines in figure 6(a). The first of these, in figure 6(b), considers translations with $K_2 = -K_1$ . When $|K_1| \lt 1.06$ (3 s.f.), four intersections are identified, as already discussed. However, when $1.06 \lt |K_1| \lt 2.30$ (3 s.f.), two intersections are lost, since the line no longer passes through the central cross-shaped surface. On increasing $|K_1|$ further, two pairs of intersections are created with corner surfaces (2 and 6 for $K_1 \gt 2.30$ , 3 and 7 for $K_1 \lt -2.30$ ), ultimately leading to well-posed regions when $|K_1| \gt 3.43$ (3 s.f.). Figure 6(c) shows translations with $K_2 = -K_1 -6$ . In this case, when $K_1 = 0$ , the line passes out of region 1, through the arm of the cross that extends along the $P_1 = 0$ plane and clips the sixth corner surface before passing into region 8, leading to six intersections. Larger $K_1$ values lead to a band of complex characteristics ( $1.04 \lt K_1 \lt 2.26$ , 3 s.f.), where the line misses the cross arm. When $K_1$ is lowered from zero, it misses region 6 and the cross arm in turn, leaving only two intersections for $-1.48 \lt K_1 \lt -1.04$ (3 s.f.). In the interval $-2.69 \lt K_1 \lt -1.48$ (3 s.f.), the line again intersects with the cross surface, this time through the arm extending along the $P_3$ axis. Lowering $K_1$ further leads to intersections with the seventh and third corner surfaces for $K_1 \lt -3.42$ (3 s.f.) and $K_1 \lt -9.40$ (3 s.f.) respectively. Finally, the intersections depicted in figure 6(d), which cover translations along $K_2 = -K_1 + 4$ , are similar, but highlight an additional case: for $0.786 \lt K_1 \lt 1.05$ (3 s.f ), the line clips through both arms of the cross, leading to a small well-posed band. Translations farther from the origin can also lead to intersections with regions $4$ and $5$ .

Varying $R_{H_1}$ and $R_{H_2}$ alters both the line and the level set. However, our earlier observation that the limits $|P_i|\to \infty$ reduce to two-phase models implies that the resulting parameter space must always contain ill-posed regions. This is likewise true for any model of the form given by (3.14). Therefore, while there may exist other three-phase models that possess more favourable properties near the origin (a fully general analysis would require us to classify all surfaces of the form given in (3.15)), none of these systems can be unconditionally well posed. Returning to the example case of the Pudasaini & Mergili (Reference Pudasaini and Mergili2019) model, an investigation of different values in the ranges $R_{H_1},R_{H_2}\in [0.1,10]$ led to regions of complex characteristics that qualitatively match the plot in figure 6(a), suggesting that the observations made thus far are robust across parameter space. It should be noted that in its full generality, this model also includes the option to include added mass forces and diffusive stresses. Though the presence of former terms alters the system’s characteristic structure, the § 3.1.2 analysis of the corresponding two-phase case does little to suggest that they will substantially improve matters. The effect of diffusion is dealt with in the next section.

4.1. A general framework for finding Hadamard instabilities

We return to the linear stability problem given in (3.5). A general procedure for detecting the presence or absence of Hadamard instabilities is developed. Since it is cast as an arbitrary matrix equation, there is no restriction on the dimensionality $N$ of the system, so our analysis in this subsection is applicable to models with any number of phases $n=N/2$ . Readers that would rather skip the linear algebra may proceed to the final paragraph of this subsection, where the method for determining posedness is recapitulated.

First, we bring (3.5) into a simpler form for analysis. The matrix $\unicode{x1D63C}$ must be invertible, in order for there to be $N$ independent time-evolving fields. Furthermore, we assume that the matrix ${\unicode{x1D63C}}^{-1}{\unicode{x1D63F}}$ is diagonalisable, since this covers all the specific cases in this paper. Then, the problem may be reformulated in terms of a basis $\{\hat {\boldsymbol{e}}_1,\ldots ,\hat {\boldsymbol{e}}_N\}$ with respect to which ${\unicode{x1D63C}}^{-1}{\unicode{x1D63F}}$ is diagonal. Therefore, for each matrix ${\unicode{x1D648}} \in \{{\unicode{x1D63D}}, {\unicode{x1D63E}}, {\unicode{x1D63F}}\}$ , we define

(4.1a,b) \begin{equation} \skew5\hat {{\unicode{x1D648}}} = {\unicode{x1D64B}}^{-1}{\unicode{x1D63C}}^{-1}{\unicode{x1D648}}{\unicode{x1D64B}}\quad \mathrm{and}\quad \hat {\boldsymbol{v}} = {\unicode{x1D64B}}^{-1}\boldsymbol{v}, \end{equation}

for any vector $\boldsymbol{v}$ , where $\unicode{x1D64B}$ is a basis change matrix that diagonalises ${\unicode{x1D63C}}^{-1}{\unicode{x1D63F}}$ . With respect to this transformation, (3.5) becomes

(4.2) \begin{equation} \sigma \hat {\boldsymbol{r}} + {\mathrm{i}} k \skew5\hat {{\unicode{x1D63D}}}\hat {\boldsymbol{r}} = \skew5\hat {{\unicode{x1D63E}}}\hat {\boldsymbol {r}} - k^2 \skew5\hat {{\unicode{x1D63F}}}\hat {\boldsymbol{r}}. \end{equation}

At high wavenumber $k\gg 1$ , we make the following asymptotic expansions:

(4.3a,b) \begin{equation} \sigma = -\sigma _2 k^2 - {\mathrm{i}}\sigma _1 k + \sigma _0 + \ldots , \quad \hat {\boldsymbol{r}} = \hat {\boldsymbol{r}}_0 + k^{-1}\hat {\boldsymbol{r}}_{-1} + \ldots , \end{equation}

substitute them into (4.2) and look for the leading-order terms. Therefore, at $O(k^2)$ , the problem reduces to

(4.4) \begin{equation} \skew5\hat {{\unicode{x1D63F}}}\hat {\boldsymbol{r}}_0 = \sigma _2 \hat {\boldsymbol{r}}_0. \end{equation}

Noting the sign convention in (4.3a , b ), the eigenvalues $\sigma _2$ , which represent diffusion coefficients for the linear problem, must each have non-negative real part in order to avoid blow-up of $\textrm {Re}(\sigma )$ . The growth of modes with $\sigma _2 = 0$ is determined beyond this leading-order balance. If $\skew5\hat {{\unicode{x1D63F}}}$ is not full rank, it has $i \in \{1,\ldots , N\}$ zero eigenvalues. Without loss of generality, we locate these in the first $i$ diagonal values of $\skew5\hat {{\unicode{x1D63F}}}$ . The corresponding eigenvectors are determined only up to an $i$ -dimensional subspace ( $\hat {\boldsymbol{r}}_0 \in \text{span}\{\hat {\boldsymbol{e}}_1,\ldots ,\hat {\boldsymbol{e}}_i\}$ ) by (4.4).

Therefore, we proceed to the $O(k)$ part of the asymptotic expansion of (4.2). When $\sigma _2 = 0$ , this is

(4.5) \begin{equation} (\skew5\hat {{\unicode{x1D63D}}} - \sigma _1 {\unicode{x1D644}}\,)\hat {\boldsymbol{r}}_0 = {\mathrm{i}}\skew5\hat {{\unicode{x1D63F}}}\hat {\boldsymbol{r}}_{-1}. \end{equation}

Since $\hat {\boldsymbol{r}}_0 \in \text{span}\{\hat {\boldsymbol{e}}_1,\ldots ,\hat {\boldsymbol{e}}_i\}$ , only the first $i$ columns of $\skew5\hat {{\unicode{x1D63D}}}-\sigma _1{\unicode{x1D644}}$ enter into this system of equations on the left-hand side. Furthermore, only the first $i$ rows of (4.5) are needed to determine $\hat {\boldsymbol{r}}_0$ and these are rows for which the right-hand side is zero. Consequently, the $\sigma _1$ values are the eigenvalues of the matrix $\skew5\hat {{\unicode{x1D63D}}}$ with the last $N - i$ rows and columns removed. We write ${\unicode{x1D648}}_{\textit{red}}$ to denote any matrix $\unicode{x1D648}$ reduced in this way by deleting rows and columns associated with the nullspace of the diagonal matrix $\skew5\hat {{\unicode{x1D63F}}}$ . Referring back to (4.3a , b ), we obtain a second criterion that must be met to avoid Hadamard instability: the eigenvalues $\sigma _1$ of $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ must be real. If these values are also distinct, then the growth rates stay bounded as $k \to \infty$ .

However, $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ may have repeated eigenvalues, which can also lead to blow-up of $\textrm {Re}(\sigma )$ . To see why, we proceed to the $O(1)$ equation with $\sigma _2 = 0$ , which reads

(4.6) \begin{equation} (\sigma _0 {\unicode{x1D644}} - \skew5\hat {{\unicode{x1D63E}}}\,)\hat {\boldsymbol{r}}_0 + {\mathrm{i}}(\skew5\hat {{\unicode{x1D63D}}} - \sigma _1 {\unicode{x1D644}}\,)\hat {\boldsymbol{r}}_{-1} = -\skew5\hat {{\unicode{x1D63F}}}\hat {\boldsymbol{r}}_{-2}. \end{equation}

To eliminate dependence of the left-hand side on the unknown vectors $\hat {\boldsymbol{r}}_{-1}$ and $\hat {\boldsymbol{r}}_{-2}$ , the left eigenvectors, corresponding to the eigenproblem adjoint to (4.2), may be used. By repeating the arguments used to determine $\hat {\boldsymbol{r}}_0$ , these may be expanded as $\hat {\boldsymbol{l}} = \hat {\boldsymbol{l}}_0 + k^{-1}\hat {\boldsymbol{l}}_{-1}+\cdots$ and inferred to satisfy $\hat {\boldsymbol{l}}_0^T\skew5\hat {{\unicode{x1D63F}}} = \boldsymbol{0}$ and $\hat {\boldsymbol{l}}_0^T(\skew5\hat {{\unicode{x1D63D}}} - \sigma _1 {\unicode{x1D644}}\,) = {\mathrm{i}} \hat {\boldsymbol{l}}_{-1}^T{\unicode{x1D63F}}$ (when $\sigma _2 = 0$ ). For any of the $i$ modes, the dot product of the leading-order left eigenvector $\hat {\boldsymbol{l}}_0$ may be taken with (4.6) and, on rearranging the result, the formula

(4.7) \begin{equation} \sigma _0 = \frac {\hat {\boldsymbol{l}}_0 \cdot \skew5\hat {{\unicode{x1D63E}}}\hat {\boldsymbol{r}}_0}{\hat {\boldsymbol{l}}_0 \cdot \hat {\boldsymbol{r}}_0} \end{equation}

is obtained. The corresponding left and right eigenvectors for $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ are the vectors $\hat {\boldsymbol{l}}_0$ , $\hat {\boldsymbol{r}}_0$ with the last $N - i$ entries (which are all zeros) deleted. When $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ is diagonalisable, these vectors form a biorthonormal set, with the left and right eigenvectors for each mode satisfying $\hat {\boldsymbol{l}}_0\cdot \hat {\boldsymbol{r}}_0 = 1$ , so the $O(1)$ growth rate in (4.7) is well defined. However, if $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ is not diagonalisable, at least one of its eigenvalues is defective. Therefore, $\sigma _1$ is a repeated eigenvalue associated with one or more Jordan chains of length at least two. Then for the full matrix $\skew5\hat {{\unicode{x1D63D}}}$ there are two pairs of corresponding generalised left and right eigenvectors $\hat {\boldsymbol{l}}_{0,1}$ , $\hat {\boldsymbol{l}}_{0,2}$ and $\hat {\boldsymbol{r}}_{0,1}$ , $\hat {\boldsymbol{r}}_{0,2}$ , respectively, satisfying

(4.8a,b) \begin{align} \begin{cases} \hat {\boldsymbol{l}}_{0,1}^T(\skew5\hat {{\unicode{x1D63D}}}-\sigma _1 {\unicode{x1D644}}\,) = \boldsymbol{0},\\[6pt] \hat {\boldsymbol{l}}_{0,2}^T(\skew5\hat {{\unicode{x1D63D}}}-\sigma _1 {\unicode{x1D644}}\,) = \hat {\boldsymbol{l}}_{0,1}^T \end{cases} \quad \mathrm{and} \quad \begin{cases} (\skew5\hat {{\unicode{x1D63D}}}-\sigma _1 {\unicode{x1D644}}\,)\hat {\boldsymbol{r}}_{0,1} = \boldsymbol{0},\\[4pt] (\skew5\hat {{\unicode{x1D63D}}}-\sigma _1 {\unicode{x1D644}}\,)\hat {\boldsymbol{r}}_{0,2} = \hat {\boldsymbol{r}}_{0,1}, \end{cases}\end{align}

where $\hat {\boldsymbol{r}}_{0,1} \equiv \hat {\boldsymbol{r}}_0$ and $\hat {\boldsymbol{l}}_{0,1}\equiv \hat {\boldsymbol{l}}_0$ . In this case, the formula in (4.7) is always singular, since projecting any left eigenvector onto (4.8b ) shows that $\hat {\boldsymbol{l}}_0 \cdot \hat {\boldsymbol{r}}_0 = 0$ . Physically, this singularity can be thought to emerge from a resonance between two or more modes that collapse onto one another when $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ becomes defective. Examples of this are given in § 4.3.

The failure of (4.7) in these cases suggests the need for an alternative asymptotic expansion. Anticipating growth of some intermediate order between $O(k)$ and $O(1)$ , we replace the expansions in (4.3a , b ) with

(4.9a,b) \begin{equation} \sigma = -{\mathrm{i}}\sigma _1 k + \sigma _{1/2} k^{1/2} + \sigma _0 + \cdots , \quad \hat {\boldsymbol{r}} = \hat {\boldsymbol{r}}_{0,1} + k^{-1/2}\hat {\boldsymbol{r}}_{-1/2} + k^{-1}\hat {\boldsymbol{r}}_{-1}+\cdots .\end{equation}

This leaves the analysis at $O(k)$ unchanged and introduces the following equation at $O(k^{1/2})$ :

(4.10) \begin{equation} \sigma _{1/2} \hat {\boldsymbol{r}}_{0,1} + {\mathrm{i}} (\skew5\hat {{\unicode{x1D63D}}} - \sigma _1 {\unicode{x1D644}}\,) \hat {\boldsymbol{r}}_{-1/2} + \skew5\hat {{\unicode{x1D63F}}} \hat {\boldsymbol{r}}_{-3/2} = \boldsymbol{0}. \end{equation}

We project this onto $\hat {\boldsymbol{l}}_{0,2}$ and use (4.8), along with the fact that $\hat {\boldsymbol{l}}_{0,2}$ is orthogonal to the range of $\skew5\hat {{\unicode{x1D63F}}}$ , to conclude that

(4.11) \begin{equation} \sigma _{1/2} \hat {\boldsymbol{l}}_{0,2}\cdot \hat {\boldsymbol{r}}_{0} + \hat {\boldsymbol{l}}_{0,1}\cdot \hat {\boldsymbol{r}}_{-1/2} = 0. \end{equation}

The unknown vector $\hat {\boldsymbol{r}}_{-1/2}$ is eliminated by proceeding to the $O(1)$ equation. With the new expansion, this is

(4.12) \begin{equation} \sigma _{1/2}\hat {\boldsymbol{r}}_{-1/2} + \sigma _0 \hat {\boldsymbol{r}}_{0,1} +{\mathrm{i}}(\skew5\hat {{\unicode{x1D63D}}} - \sigma _1 {\unicode{x1D644}}\,)\hat {\boldsymbol{r}}_{-1} - \skew5\hat {{\unicode{x1D63E}}}\hat {\boldsymbol{r}}_{0,1} + \skew5\hat {{\unicode{x1D63F}}}\hat {\boldsymbol{r}}_{-2} = \boldsymbol{0}. \end{equation}

Then, we project this onto $\hat {\boldsymbol{l}}_{0,1}$ . Since $\hat {\boldsymbol{l}}_{0,1}\cdot \hat {\boldsymbol{r}}_{0,1} = 0$ , the term containing $\sigma _0$ vanishes, along with the diffusive term which lies in an orthogonal subspace. After rearranging and using (4.11), we obtain a formula for the $O(k^{1/2})$ part of the growth rate:

(4.13) \begin{equation} \sigma _{1/2} = \pm \frac {1 + {\mathrm{i}}}{2} \left ( \frac {2 \hat {\boldsymbol{l}}_{0,1}\cdot \skew5\hat {{\unicode{x1D63E}}}\hat {\boldsymbol{r}}_{0,1}}{\hat {\boldsymbol{l}}_{0,2} \cdot \hat {\boldsymbol{r}}_{0,1}} \right )^{\!1/2}. \end{equation}

For Jordan chains of length two, $\hat {\boldsymbol{l}}_{0,2}\cdot \hat {\boldsymbol{r}}_{0,1}=|\hat {\boldsymbol{l}}_{0,2}||\hat {\boldsymbol{r}}_{0,1}|\neq 0$ , provided both the left and right vectors correspond to the same Jordan block. Consequently, (4.13) implies that there is a mode such that $\textrm {Re}(\sigma ) \sim k^{1/2}$ , provided the source terms in $\skew5\hat {{\unicode{x1D63E}}}$ do not interact with the eigenvectors of $\skew5\hat {{\unicode{x1D63D}}}$ in a way that causes the numerator to vanish.

Conversely, for longer Jordan chains, the denominator in the (4.13) formula is also guaranteed to be singular. Different asymptotic expansions are needed, depending on the length of the the chain. However, to avoid these further complications, we terminate our analysis here, since cases where three or more modes intersect at high wavenumber are far less commonly encountered.

To summarise the analysis above, models up to second order that may be cast in the general form of (3.4) are ill posed as initial-value problems if any of the following conditions are met:

1. (i) Any eigenvalue of $\skew5\hat {{\unicode{x1D63F}}}$ is negative, where $\skew5\hat {{\unicode{x1D63F}}}$ denotes a diagonalisation of ${\unicode{x1D63C}}^{-1}{\unicode{x1D63F}}$ .

2. (ii) Any eigenvalue of $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ is complex, where $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ denotes the matrix formed by representing ${\unicode{x1D63C}}^{-1}{\unicode{x1D63D}}$ in the basis used to diagonalise ${\unicode{x1D63C}}^{-1}{\unicode{x1D63F}}$ in (i) and deleting each row and column $j$ such that the $j$ th diagonal entry of $\skew5\hat {{\unicode{x1D63F}}}$ is zero. We refer to $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ as a ‘reduced Jacobian’ in later analysis.

3. (iii) Repeated real eigenvalues of $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ of algebraic multiplicity $2$ share the same left and right eigenvectors $\hat {\boldsymbol{l}}_{0,1}$ and $\hat {\boldsymbol{r}}_{0,1}$ (up to normalisation), and $\hat {\boldsymbol{l}}_{0,1}\cdot \skew5\hat {{\unicode{x1D63E}}}_{\textit{red}}\hat {\boldsymbol{r}}_{0,1}\neq 0$ , where $\skew5\hat {{\unicode{x1D63E}}}_{\textit{red}}$ is defined via ${\unicode{x1D63C}}^{-1}{\unicode{x1D63E}}$ in the same way as $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ . (More generally, the expectation following from (4.7) is that repeated real eigenvalues of any algebraic multiplicity $m\geqslant 2$ imply ill posedness if the dimension of their associated eigenspace is strictly less than $m$ , but this is not explicitly proven above.)

For the remainder of this section, we apply these steps to different example systems.

4.2. Velocity diffusion in every momentum equation

Before analysing individual models, we highlight a generic case, which is guaranteed to be well posed. Consider an $n$ -phase model, with each phase $i$ characterised by height $H_i$ and velocity $\overline {u_i}$ , organised into pairs of mass and momentum equations of the form

(4.14a,b) \begin{equation} \frac {\partial H_i}{\partial t} + \frac {\partial \,}{\partial x} \left ( H_i \overline {u_i} \right ) = 0, \quad \frac {\partial \overline {u_i}}{\partial t} + F_i(H_1,\overline {u_1},\ldots ,H_n,\overline {u_n}) = 0, \end{equation}

where the functions $F_i$ contain no dependence on time derivatives or spatial derivatives of first order or higher. To each of the the $j = 1,\ldots ,n$ momentum equations, add a term of the form $({\partial \,}/{\partial x})(\nu _j(\boldsymbol{q}) ({\partial \overline {u_j}}/{\partial x}))$ , where $\nu _j(\boldsymbol{q})$ denotes a general diffusivity coefficient function that stays strictly positive. When casting the linearised problem in matrix form, the equations are ordered so that the mass and momentum equations respectively lie on odd and even rows, as before. The corresponding diffusion matrix is already diagonal, so at any $\boldsymbol{q}=\boldsymbol{q}_0$ ,

(4.15) \begin{equation} \skew5\hat {{\unicode{x1D63F}}} = \begin{pmatrix} 0 & 0 & \ldots & 0 & 0 \\[3pt] 0 & \nu _1(\boldsymbol{q}_0) & \ldots & 0 & 0 \\[3pt] \vdots & \vdots & \ddots & \vdots & \vdots \\[3pt] 0 & 0 & \ldots & 0 & 0 \\[3pt] 0 & 0 & \ldots & 0 & \nu _n(\boldsymbol{q}_0) \end{pmatrix}, \end{equation}

which is positive semidefinite and damps out growth at high wavenumber for $n$ of the $2n$ stability modes. The reduced Jacobian matrix is simply

(4.16) \begin{equation} \skew5\hat {{\unicode{x1D63D}}}_{\textit{red}} = \begin{pmatrix} \overline {u_1} & 0 & \ldots & 0 \\[3pt] 0 & \overline {u_2} & \ldots & 0 \\[3pt] \vdots & \vdots & \ddots & \vdots \\[3pt] 0 & 0 & \ldots & \overline {u_n} \\[3pt] \end{pmatrix}, \end{equation}

which clearly possesses $n$ real eigenvalues and $n$ linearly independent eigenvectors. Therefore, debris-flow models with ${\unicode{x1D63C}}={\unicode{x1D644}}$ can always be regularised by adding positive diffusion to every momentum equation.

4.3. Existing models

4.3.1. Meng et al. (Reference Meng, Johnson and Gray2022)

As detailed in § 2.2, in the model of Meng et al. (Reference Meng, Johnson and Gray2022), ${\unicode{x1D63C}}={\unicode{x1D644}}$ and only diffusion in the fluid phase is included. The equations are (2.17a ), (2.22a ), (2.17c ) and (2.22b ), rendered dimensionless as described in § 3.1. The diffusion matrix is already diagonal and is given by

(4.17) \begin{equation} \skew5\hat {{\unicode{x1D63F}}} = \begin{pmatrix} 0 & 0 & 0 & 0 \\[4pt] 0 & 0 & 0 & 0 \\[4pt] 0 & 0 & 0 & 0 \\[4pt] 0 & 0 & 0 & 2\nu _f \end{pmatrix} , \end{equation}

where $\nu _f \equiv \eta _f / (\rho _f H_f^{(0)}\overline {u_f}^{(0)})\gt 0$ is a dimensionless kinematic viscosity coefficient (though it could equally be viewed as an eddy diffusivity if the flow is turbulent). Therefore, the corresponding reduced Jacobian is formed by removing the fourth row and column from the full matrix $\skew5\hat {{\unicode{x1D63D}}} = {\unicode{x1D63D}}$ , given in (3.6). At $\boldsymbol{q}=\boldsymbol{q}_0$ ,

(4.18) \begin{equation} \skew5\hat {{\unicode{x1D63D}}}_{\textit{red}} = \begin{pmatrix} R_u & R_H & 0 \\[4pt] (\gamma +\beta _1){\textit{Fr}}^{-2} & R_u & (\gamma + \beta _2){\textit{Fr}}^{-2} \\[4pt] 0 & 0 & 1 \end{pmatrix}. \end{equation}

Its eigenvalues are

(4.19) \begin{equation} \sigma _1 = 1, R_u \pm \frac {\sqrt {R_H(\gamma + \beta _1)}}{{\textit{Fr}}} \end{equation}

with corresponding eigenvectors

(4.20) \begin{equation} \begin{pmatrix} R_H \\[4pt] 1 - R_u \\[4pt] \displaystyle\frac {(R_u-1)^2{\textit{Fr}}^2-R_H(\gamma +\beta _1)}{\gamma +\beta _2} \end{pmatrix}, \,\, \begin{pmatrix} R_H \\[4pt] \pm \frac {1}{\textit{Fr}}\sqrt {R_H(\gamma +\beta _1)} \\[4pt] 0 \end{pmatrix}. \end{equation}

Firstly, note that the latter pair of $\sigma _1$ values equal the characteristics for the solids phase of the ‘decoupled’ problem, given in (3.7a , b ). Hence, all values are expected to be real. However, there is the opportunity for repeated eigenvalues, which occurs when $\sigma _1=1$ matches either of the other two growth rates, i.e. when

(4.21) \begin{equation} R_u = 1 \pm \frac {1}{{\textit{Fr}}}\sqrt {R_H(\gamma + \beta _1)}. \end{equation}

This is similar, but not equivalent, to the condition for intersecting decoupled characteristics, given previously in (3.11). By substituting (4.21) into (4.20), it may be verified that the corresponding eigenvectors are equal when this condition is satisfied. Consequently, the equations feature an instability with order $k^{1/2}$ growth rate in the high-wavenumber asymptotic limit and are ill posed wherever (4.21) is satisfied.

Diffusion of momentum is often also included in shallow models of dry granular flows (e.g. Gray & Edwards Reference Gray and Edwards2014; Baker et al. Reference Baker, Johnson and Gray2016). The general argument in § 4.2 implies that adding a diffusive term of the form $ ({\partial \,}/{\partial x})(\nu _s ({\partial \overline {u_s}}/{\partial x}))$ , where $\nu _s \equiv \nu _s(\boldsymbol{q}_0) \gt 0$ , to the solids momentum equation is sufficient to regularise this model. Moreover, using analogous arguments to those above, it may be verified that diffusive terms in both momentum terms are required, in order to guarantee that the model stays unconditionally well posed.

Specifically, if $\nu _f = 0$ , but diffusion in the solids momentum equation is included, then the reduced Jacobian is formed by eliminating the second row and column of $\skew5\hat {{\unicode{x1D63D}}}$ , to leave

(4.22) \begin{equation} \skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}= \begin{pmatrix} R_u & 0 & 0 \\[4pt] 0 & 1 & 1 \\[4pt] {\textit{Fr}}^{-2} & {\textit{Fr}}^{-2} & 1 \end{pmatrix}. \end{equation}

This matrix is defective when

(4.23) \begin{equation} R_u = 1 \pm 1/{\textit{Fr}}, \end{equation}

giving rise to a family of $O(k^{1/2})$ instabilities at these points in parameter space, similar to the case where only fluid diffusion is included.

Note that since we used the general form of $\unicode{x1D63D}$ from (3.6) to construct the reduced Jacobians, these assessments apply also to the case of adding simple diffusive terms to regularise the models of Pitman & Le (Reference Pitman and Le2005), Pelanti et al. (Reference Pelanti, Bouchut and Mangeney2008) and Meyrat et al. (Reference Meyrat, McArdell, Ivanova, Müller and Bartelt2022).

4.3.2. Pudasaini (Reference Pudasaini2012)

This model incorporates two diffusive stresses for the fluids phase: a Newtonian component, equivalent to the term used by Meng et al. (Reference Meng, Johnson and Gray2022), and a non-Newtonian closure defined in § 2.3. The relevant contributions to the depth-averaged downslope momentum equation are the second-order terms of (2.25). When the model is converted to the quasilinear form that was used for the local analysis, the fluid momentum equation is non-dimensionalised (as per § 3.1) and divided through by $\rho _f H_f$ , and the diffusive terms become

(4.24) \begin{equation} \frac {2\nu _f}{H_f}\frac {\partial \,}{\partial x}\left ( H_f\frac {\partial \overline {u_f}}{\partial x} \right ) + \frac {2 \nu _f {\mathcal{N}}}{H_f} \frac {\partial \,}{\partial x} \left [ H_f(\overline {u_s} - \overline {u_f}) \frac {\partial \,}{\partial x}\left ( \frac {H_s}{H_s + H_f} \right ) \right ]\!. \end{equation}

The parameter ${\mathcal{N}} = \overline {{\mathcal{A}}}/\overline {\varphi _f}$ is a ratio of the effective diffusion coefficients for the Newtonian and non-Newtonian parts (see (2.25)) and is assumed to be constant by Pudasaini (Reference Pudasaini2012).

After linearising around $\boldsymbol{q}=\boldsymbol{q}_0 = (R_H, R_u, 1, 1)^T$ , the diffusion matrix becomes

(4.25) \begin{equation} {\unicode{x1D63F}}= 2\nu _f \begin{pmatrix} 0 & 0 & 0 & 0 \\[4pt] 0 & 0 & 0 & 0 \\[4pt] 0 & 0 & 0 & 0 \\[4pt] \displaystyle\frac {\mathcal{N}(R_u-1)}{(1+R_H)^2} & 0 & \displaystyle\frac {{\mathcal{N}}R_H(1-R_u)}{(1+R_H)^2} & 1 \end{pmatrix} . \end{equation}

The matrices $\unicode{x1D63C}$ and $\unicode{x1D63D}$ were given previously, in (3.12) and (3.13), respectively. The basis change matrix

(4.26) \begin{equation} {\unicode{x1D64B}}= \begin{pmatrix} 1 & 0 & 0 & 0 \\[4pt] 0 & 1 & 0 & \gamma \overline {C} \\[4pt] 0 & 0 & 1 & 0 \\[4pt] \displaystyle\frac {{\mathcal{N}}(1-R_u)}{(1+R_H)^2} & 0 & \displaystyle\frac {{\mathcal{N}}R_H(R_u-1)}{(1+R_H)^2} & \gamma \overline {C}+1 \end{pmatrix} \end{equation}

diagonalises ${\unicode{x1D63C}}^{-1}{\unicode{x1D63F}}$ , so that $\skew5\hat {{\unicode{x1D63F}}}={\unicode{x1D64B}}^{-1}{\unicode{x1D63C}}^{-1}{\unicode{x1D63F}}{\unicode{x1D64B}}$ is the matrix of all zeros, save for its only eigenvalue $\sigma _2$ , located on the bottom right entry:

(4.27) \begin{equation} {{D}}_{44} = \sigma _2 = \frac {2\nu _f(\gamma \overline {C} + 1)}{1 + \overline {C}(\gamma +R_H)}. \end{equation}

Since $\sigma _2 \gt 0$ , there is no blow-up at $O(k^2)$ and it remains to check the properties of $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ , which is formed by removing the fourth row and column of ${\unicode{x1D64B}}^{-1}{\unicode{x1D63C}}^{-1}{\unicode{x1D63D}}{\unicode{x1D64B}}$ .

In the general case, analytical expressions for this matrix are cumbersome and it is better to compute its eigenvalues numerically. However, two limiting cases are tractable. Firstly, when the non-Newtonian viscosity is not included, ${\mathcal{N}}=0$ and the reduced Jacobian is

(4.28) \begin{equation} \skew5\hat {{\unicode{x1D63D}}}_{\textit{red}} = \begin{pmatrix} R_u & R_H & 0 \\[4pt] \displaystyle\frac {1}{\gamma \overline {C}+1}\left [\displaystyle\dfrac {\gamma + \beta _1}{{{{\textit{Fr}}}}^{2}} + \displaystyle\frac {\overline {C}(R_u-1)}{R_H}\right ]& R_u + \displaystyle\frac {\gamma \overline {C}}{\gamma \overline {C}+1} & \displaystyle\frac {\gamma + \beta _2}{(\gamma \overline {C}+1){{{\textit{Fr}}}}^2} \\[4pt] 0 & 0 & 1 \end{pmatrix}. \end{equation}

The eigenvalues $\sigma _1$ of this matrix are

(4.29) \begin{equation} \sigma _1 = 1, R_u + \frac {1}{\gamma \overline {C}+1}\left ( \frac {\gamma \overline {C}}{2} \pm \frac {\sqrt {\varDelta }}{{\textit{Fr}}} \right ), \end{equation}

where

(4.30) \begin{equation} \varDelta = (\gamma \overline {C}+1) \left [ {\textit{Fr}}^2\gamma \overline {C}(R_u-1) +(\gamma +\beta _1)R_H \right ] +\left (\frac {{\textit{Fr}}\hspace {0.1em}\gamma \overline {C}}{2}\right )^2\!\!. \end{equation}

The latter pair are complex conjugate if and only if $\varDelta \lt 0$ . Rearranging this inequality leads to

(4.31) \begin{equation} R_u - 1 \lt -\frac {\gamma \overline {C}}{4(\gamma \overline {C}+1)} - \frac {R_H(\gamma +\beta _1)}{\gamma \overline {C}{\textit{Fr}}^2}. \end{equation}

This describes a region of complex eigenvalues that is constrained to lie within $R_u \lt 1$ . Note that in the $\overline {C}\to 0$ limit, this region entirely recedes and inequality (4.31) is never satisfied. In addition to these complex eigenvalues, there is the opportunity for $O(k^{1/2})$ blow-up if $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ is defective, which can happen if $\sigma _1=1$ intersects with either of the other two eigenvalues in (4.29). The condition for this simplifies to

(4.32) \begin{equation} R_u = 1 \pm \frac {1}{{\textit{Fr}}}\sqrt {\frac {R_H(\gamma +\beta _1)}{\gamma \overline {C}+1}}, \end{equation}

which generalises (4.21) for cases where $\overline {C}\geqslant 0$ . It may be separately verified that only one eigenvector corresponding to $\sigma _1=1$ exists when $R_u$ satisfies (4.32), implying that $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ is defective here.

Figure 7. Regions where the reduced Jacobian for the model of Pudasaini (Reference Pudasaini2012) with diffusive terms possesses complex eigenvalues (red shading), for $R_H = 1$ , $\gamma = 0.5$ , ${\mathcal{N}}=0$ and $\overline {C} =$ (a) $0.02$ , (b) $0.1$ and (c) $0.5$ . The boundaries of these regions are given analytically by inequality (4.31) (dotted black). Along the black dashed lines, given by (4.32), the reduced Jacobian is defective. The model is ill posed as an initial-value problem for flow states that pass through either the dashed line or the red region.

In figure 7, we show the regions where the model is ill posed for ${\mathcal{N}}=0$ , $\overline {C}=$ $0.02$ , $0.1$ and $0.5$ and the same illustrative parameters used in figure 4(a). Dashed curves show the lines given by (4.32). The ill-posed region that emerges at low $R_u$ values via inequality (4.31), whose border is given by the dotted line, may be compared with similar regions present in the problem without diffusion, plotted in figure 4.

If instead ${\mathcal{N}}\gt 0$ and the limit of vanishing added mass $\overline {C}\to 0$ is taken, then

(4.33) \begin{equation} \skew5\hat {{\unicode{x1D63D}}}_{\textit{red}} = \begin{pmatrix} R_u & R_H & 0 \\[4pt] (\gamma + \beta _1){\textit{Fr}}^{-2} & R_u & (\gamma + \beta _2){\textit{Fr}}^{-2} \\[4pt] \displaystyle\frac {{\mathcal{N}}(1-R_u)}{(1+R_H)^2} & 0 & 1 + \displaystyle\frac {{\mathcal{N}}R_H(R_u-1)}{(1+R_H)^2} \end{pmatrix}. \end{equation}

Note that this is a generalisation of the reduced Jacobian in (4.18). When the non-Newtonian terms are included, ${\mathcal{N}}$ is expected to be a large number compared with the other parameters (Pudasaini (Reference Pudasaini2012) uses ${\mathcal{N}} = 5000$ ). By solving for roots of the characteristic polynomial via series expansion, when ${\mathcal{N}}\gg 1$ , the eigenvalues of $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ may be obtained:

(4.34) \begin{equation} \sigma _1 = \frac {(R_u-1)R_H}{(1+R_H)^2}{\mathcal{N}} + 1 + O({\mathcal{N}}^{-1}), \,\, \pm \frac {1}{{\textit{Fr}}}\sqrt {(\gamma +\beta _1)R_H + \gamma + \beta _2} + O({\mathcal{N}}^{-1}). \end{equation}

These expressions are real-valued and remain so in the limit. However, either of the second and third branches merges with the first when

(4.35) \begin{equation} R_u = 1 \pm \frac {(1+R_H)^2}{{\textit{Fr}} R_H{\mathcal{N}}}\sqrt {(\gamma +\beta _1)R_H + \gamma + \beta _2} + O({\mathcal{N}}^{-2}). \end{equation}

As we have seen previously, the merging of branches can give rise to complex eigenvalues. In this case, the merger originates in the limit ${\mathcal{N}}\to \infty$ . For large but finite ${\mathcal{N}}$ , we compute the eigenvalues of $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ numerically and summarise their type in figure 8. The two parametric lines given in (4.35) are flanked by bands where $\sigma _1$ takes complex values. Furthermore, the shape of these bands is self-similar in the asymptotic high- ${\mathcal{N}}$ regime.

Figure 8. Regions where the model of Pudasaini (Reference Pudasaini2012) with diffusive terms is ill posed as an initial-value problem (red shading), for $R_H = 1$ , $\gamma = 0.5$ , $\overline {C}=0$ and high values of the ratio ${\mathcal{N}}$ between non-Newtonian and Newtonian diffusion coefficients: ${\mathcal{N}} =$ (a) $50$ , (b) $500$ and (c) $5000$ . Note that each vertical axis is scaled with respect to $50/{\mathcal{N}}$ and that the shaded regions are near identical under this rescaling. Asymptotic expansions for the eigenvalues at high ${\mathcal{N}}$ intersect along the black dashed lines, whose formulae are given in (4.35).

4.3.3. Pudasaini & Mergili (Reference Pudasaini and Mergili2019)

To conclude this section, we touch upon the three-phase model of Pudasaini & Mergili (Reference Pudasaini and Mergili2019), which was introduced in § 3.2. It was shown previously that omitting the diffusive terms in this model can lead to ill-posed initial-value problems. However, it remains to be seen whether including the terms can eliminate this issue. As before, we neglect the complications of the added mass effect, though as we have just seen, this can be analysed using the same methods.

Diffusion of momentum in the Pudasaini & Mergili (Reference Pudasaini and Mergili2019) model appears in the equations for both fluid phases (which are labelled $2$ and $3$ in § 3.2). Each contains a Newtonian and non-Newtonian component similar to the terms in (4.24) for the Pudasaini (Reference Pudasaini2012) system. The diffusion matrix $\unicode{x1D63F}$ for the non-dimensionalised and linearised model equations is given explicitly in Appendix C. Due to the non-Newtonian terms, it has off-diagonal entries. It possesses two non-zero eigenvalues, which are $2\nu _2 \gt 0$ and $2\nu _3\gt 0$ , where $\nu _2$ and $\nu _3$ are the Newtonian diffusion coefficients associated with the second and third phases. A suitable basis change matrix $\unicode{x1D64B}$ that diagonalises $\unicode{x1D63F}$ was determined using computational algebra and is also specified in Appendix C. This allows us to form the reduced Jacobian (a $4\times 4$ matrix in this case) numerically and compute its eigenvalues.

Since there are five independent dimensionless variables $(R_{H_1}, R_{H_2}, R_{u_1}, R_{u_2}, {\textit{Fr}}\,)$ that specify a particular state (in addition to several fixed model parameters), we do not attempt an exhaustive study. Instead, we fix $R_{H_1}=1$ and investigate the effect of introducing the intermediate fluid phase by increasing $R_{H_2}$ from zero. Guided by our analysis in § 3.2, we shift $R_{u_2}$ slightly away from unity, setting $R_{u_2}=1.01$ , to allow for richer interactions between the phases. Figure 9 shows the results of these computations, using illustrative model parameters, given in Appendix C. These parameters were selected to match our choices for computations relating to the Pudasaini (Reference Pudasaini2012) model with ${\mathcal{N}}=5000$ (§ 4.3.2), so that when $R_{H_2}\to 0$ , the system collapses to the this two-phase case. When $R_{H_2} = 10^{-3}$ (figure 9 a), there are two bands where the reduced Jacobian features a pair of complex eigenvalues, either side of $R_{u_1}=1$ . As expected, these closely match the corresponding regions plotted in figure 8(c). However, the reflection symmetry of these bands about $R_{u_1} = 1$ is broken for any $R_{H_{2}}\gt 0$ . This becomes apparent at higher $R_{H_2}$ values. Figures 9(b) and 9(c) show the cases $R_{H_2} = 0.01$ and $ 0.1$ , respectively. The upper band drops below $R_{u_1} \lt 1$ and overlaps the lower band, with the Froude numbers at which this occurs decreasing as $R_{H_2}$ increases. Where the bands overlap, there are two pairs of complex eigenvalues. Additionally, a second upper band appears at higher $\textit{Fr}$ and draws towards lower $\textit{Fr}$ as $R_{H_2}$ increases to $0.5$ and $1$ , in figures 9(d) and 9(e), respectively. In the final plot (figure 9 f), at $R_{H_2} = 4$ , the two upper bands have merged, though in this case the merger does not double the number of complex eigenvalues present.

Figure 9. Illustrative computations of the eigenvalues of $\skew5\hat {{\unicode{x1D63D}}}_{\textit{red}}$ for the model of Pudasaini & Mergili (Reference Pudasaini and Mergili2019), with $R_{H_1} = 1$ and $R_{u_2} = 1.01$ . In regions shaded pink, the model possesses a single pair of complex eigenvalues, while red shading covers areas where two complex pairs were found. Elsewhere, all eigenvalues are real. The parameters for these computations are given in Appendix C.

Other choices for the flow variables lead to plots that are similar to figure 9, at least in the sense that they are constructed from complicated tangles of complex eigenvalue regions. While it may be possible to make sense of these diagrams in detail, this is perhaps beside the point. It is clear, even from this cursory investigation, that this three-phase model suffers from the same issues as the two-phase models.