2.1 Combined equations

We combine all of the above equations into a set of two governing equations:

(2.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}\frac{\text{D}p}{\text{D}t}=\frac{k_{0}}{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }[(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{c})_{+}^{b}p\unicode[STIX]{x1D735}p]-\frac{1}{\unicode[STIX]{x1D701}}p(p-p^{lava}) & \displaystyle\end{eqnarray}$$

(2.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\unicode[STIX]{x1D719}}{\text{D}t}=\frac{1}{\unicode[STIX]{x1D701}}(1-\unicode[STIX]{x1D719})(p-p^{lava}), & \displaystyle\end{eqnarray}$$

$\text{D}/\text{D}t=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$

To non-dimensionalize this system, we use the transformation

(2.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}p^{\ast }=p/p_{0};\quad \boldsymbol{V}^{\ast }=\boldsymbol{V}/V_{c};\quad \boldsymbol{x}^{\ast }=\boldsymbol{x}/h_{0};\quad t^{\ast }=t/t_{c};\\ h^{\ast }=h/h_{0};\quad f=p^{lava}/p_{0}=f_{a}+h^{\ast }-z^{\ast },\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $h_{0}$ is the maximum dome height, $p_{0}=\unicode[STIX]{x1D70C}^{lava}gh_{0}$ is the maximum lava pressure at the base of the dome, $f_{a}=p_{a}/p_{0}$ is the non-dimensional ambient air pressure, $f$ represents the non-dimensional lava pressure and $V_{c}$ and $t_{c}$ are scale factors for the lava velocity and time respectively. This yields a non-dimensional system (with stars dropped) where $p$ , $\unicode[STIX]{x1D719}$ , $\boldsymbol{x}$ , $t$ , $\boldsymbol{V}$ and $f$ are now $O(1)$ dimensionless variables:

(2.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}(\unicode[STIX]{x2202}_{t}p+\unicode[STIX]{x1D712}\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p)=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D735}\boldsymbol{\cdot }[(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{c})_{+}^{b}p\unicode[STIX]{x1D735}p]-\unicode[STIX]{x1D6FD}p(p-f) & \displaystyle\end{eqnarray}$$

(2.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}+\unicode[STIX]{x1D712}\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D6FD}(1-\unicode[STIX]{x1D719})(p-f) & \displaystyle\end{eqnarray}$$

(2.6a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D712}=\frac{V_{c}t_{c}}{h_{0}};\quad \unicode[STIX]{x1D6FC}=\frac{k_{0}p_{0}t_{c}}{\unicode[STIX]{x1D707}h_{0}^{2}};\quad \unicode[STIX]{x1D6FD}=\frac{p_{0}t_{c}}{\unicode[STIX]{x1D701}}. & & \displaystyle\end{eqnarray}$$

Typical lava dome volumetric flow rates are ${\sim}1{-}10~\text{m}^{3}~\text{s}^{-1}$ , implying that for typical lava dome volumes ( $O(10^{6}){-}O(10^{8})~\text{m}^{3}$ ), the mean velocity throughout a dome is only $O(10^{-5}){-}O(10^{-4})~\text{m}~\text{s}^{-1}$ (Sheldrake et al. Reference Sheldrake, Sparks, Cashman, Wadge and Aspinall2016). For typical values of the parameters in $\unicode[STIX]{x1D712}$ , $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ the diffusion and compaction terms are much stronger than the advection term (table 1). Therefore we neglect the advection and take the limit $\unicode[STIX]{x1D712}\rightarrow 0$ throughout the present study. This limit simplifies (2.5a ) and (2.5b ) to a system of two equations for two unknown fields: the pressure $p$ and porosity $\unicode[STIX]{x1D719}$ .

Table 1.

Values of parameters used including ranges where uncertain.

2.2 One-dimensional analysis

We will consider a simplified version of these combined equations: purely vertical gas flow, casting the equations over only the vertical coordinate and time as

(2.7a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}\,\unicode[STIX]{x2202}_{t}p=\unicode[STIX]{x1D6FC}\,\unicode[STIX]{x2202}_{z}[K(\unicode[STIX]{x1D719})p\,\unicode[STIX]{x2202}_{z}p]-\unicode[STIX]{x1D6FD}p(p-f) & \displaystyle\end{eqnarray}$$

(2.7b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D6FD}(1-\unicode[STIX]{x1D719})(p-f), & \displaystyle\end{eqnarray}$$

$K(\unicode[STIX]{x1D719})=(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{c})_{+}^{b}$

(2.8a ) $$\begin{eqnarray}\displaystyle & \displaystyle p=f_{a}\quad \text{on }\{z=1\}\times [0,\infty ] & \displaystyle\end{eqnarray}$$

(2.8b ) $$\begin{eqnarray}\displaystyle & \displaystyle K(\unicode[STIX]{x1D719})p\,\unicode[STIX]{x2202}_{z}p=-\unicode[STIX]{x1D6EC}_{s}Q(t)\quad \text{on }\{z=0\}\times [0,\infty ], & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}_{s}=\frac{{\dot{m}}_{s}\unicode[STIX]{x1D707}RT^{dome}}{\unicode[STIX]{x1D70C}^{2}g^{2}h_{0}k_{0}A}, & & \displaystyle\end{eqnarray}$$

which represent a constant ambient pressure at the dome surface and some flux of gas from the lower boundary (volcanic conduit) which is $\unicode[STIX]{x1D6EC}_{s}$ during steady-state conditions and varies from steady state according to $Q(t)$ . The lower boundary condition is given here in a general form to be detailed later. These simplified equations represent a coupled system in which the diffusion of gas pressure is modified by changes in the porosity and permeability due to expansion and compression of the porous lava medium. This compaction is itself dependent on the gas–lava pressure difference. Equation (2.7a ) is a nonlinear parabolic PDE which becomes singular in both of the limits $p\rightarrow 0^{+}$ and $\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{c}^{+}$ . The first limit is a common degeneracy and is related to the similarity between (2.7a ) and the equation governing porous flow of gases in a rigid network, the so-called porous medium equation (PME). Generally, the PME yields nonlinear diffusion due to pressure-proportional diffusivity resulting in steeper pressure gradients than linear diffusion and finite propagation speed (Vazquez Reference Vazquez2007). In the present context, the $p\rightarrow 0$ limit is unphysical; however, porosities below the percolation threshold are found widely in lava dome samples (Boudon et al. Reference Boudon, Balcone-Boissard, Villemant and Morgan2015) and so the second limit ( $\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{c}$ ) is physically relevant. If $\unicode[STIX]{x1D719}\leqslant \unicode[STIX]{x1D719}_{c}$ , equations (2.7a ) and (2.7b ) can be combined as

(2.10) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}_{t}p}{p}=\frac{\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}(1-\unicode[STIX]{x1D719})} & & \displaystyle\end{eqnarray}$$

and integrated to find a phase-plane solution:

(2.11) $$\begin{eqnarray}\displaystyle p(z,t)=p_{c}(z)\frac{\unicode[STIX]{x1D719}_{c}}{1-\unicode[STIX]{x1D719}_{c}}\left[\frac{1-\unicode[STIX]{x1D719}(z,t)}{\unicode[STIX]{x1D719}(z,t)}\right], & & \displaystyle\end{eqnarray}$$

where $p_{c}(z)$ is the local pore pressure when the porosity drops below the percolation threshold or $p_{c}(z)=p(z;\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{c})$ . Although there is no closed-form solution for $p$ and $\unicode[STIX]{x1D719}$ in terms of $z$ and $t$ only, this relation indicates that in a hydraulically unconnected region, the pressure and porosity are inversely related, that is, when the porosity drops below the percolation threshold due to a negative pressure anomaly ( $p-f<0$ ), the porosity will continue to contract, causing the pressure to rise back to $p=f$ .

Additionally, equations (2.7a ) and (2.7b ) have three steady-state solutions ( $p_{s}(z)$ , $\unicode[STIX]{x1D719}_{s}(z)$ ), although only one is physical and allows porous gas flow:

(2.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle p_{s}(z)=f(z)=f_{a}+1-z & \displaystyle\end{eqnarray}$$

(2.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{s}(z)=\unicode[STIX]{x1D719}_{c}+\left[\frac{\unicode[STIX]{x1D6EC}_{s}}{f(z)}\right]^{1/b}, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6EC}_{s}$

${\dot{m}}_{s}$

$A$

$\unicode[STIX]{x1D6EC}_{s}$

$0\leqslant \unicode[STIX]{x1D719}\leqslant 1$

$\unicode[STIX]{x1D719}_{a}$

(2.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{a}=\unicode[STIX]{x1D719}_{c}+\left(\frac{\unicode[STIX]{x1D6EC}_{s}}{f_{a}}\right)^{1/b}=\unicode[STIX]{x1D719}_{c}+\left(\frac{{\dot{m}}_{s}\unicode[STIX]{x1D707}RT^{dome}}{p_{a}\unicode[STIX]{x1D70C}^{lava}gk_{0}A}\right)^{1/b}. & & \displaystyle\end{eqnarray}$$

To find physically reasonable values of the flux $\unicode[STIX]{x1D6EC}_{s}$ , the parameters must be constrained such that $\unicode[STIX]{x1D719}_{c}<\unicode[STIX]{x1D719}_{a}(\unicode[STIX]{x1D6EC}_{s})<1$ .

The exact form of the porosity steady-state solution is governed by the $k-\unicode[STIX]{x1D719}$ relationship. Although $K(\unicode[STIX]{x1D719}_{s})\propto 1/f$ in general, the exact form of $\unicode[STIX]{x1D719}_{s}(z)$ will depend upon inverting the particular $K(\unicode[STIX]{x1D719})$ for the problem. Consideration of other porosity–permeability relationships is beyond the scope of this work.

These solutions indicate a magmastatic steady-state pore gas pressure and that the steady-state porosity profile accommodates this linear pressure increase by increased compaction with depth (figure 2 a). Furthermore, the dimensionless pressure gradient is always normalized as $\text{d}f/\text{d}z=-1$ regardless of the choice of the parameters $f_{a}$ and $\unicode[STIX]{x1D719}_{a}$ ; however, the porosity profile changes significantly for different choices of these parameters (figure 2 b). Indeed, for small $f_{a}$ , the porosity profiles are most heavily curved, which can be seen as phase plane curves with the lowest values of $\text{d}p/\text{d}\unicode[STIX]{x1D719}$ (figure 2). The steady-state pressure plays a particularly important role for understanding how the transient equations evolve since the expansion and compaction terms of (2.7a ) and (2.7b ) are functions of the deviation from steady-state pressure, that is, perfect gas–lava pressure coupling.

Figure 2.

(a) One-dimensional steady-state (equilibrium) profiles for gas pressure (non-dimensional) (grey) and porosity (black). (b) Variation of these profiles for different parameter values represented in the $p-\unicode[STIX]{x1D719}$ phase plane. $\unicode[STIX]{x1D719}_{a}=0.4$ (solid), $\unicode[STIX]{x1D719}_{a}=0.6$ (dashed), $\unicode[STIX]{x1D719}_{a}=0.8$ (dotted).

3.1 Uniform background state

To examine the effects of the background state on the behaviour of perturbations, we first consider the assumption of a uniform background state that has typically been used in linear stability analyses (LSA) of similar gas–lava mixtures (Manga Reference Manga1996; Jellinek & Bercovici Reference Jellinek and Bercovici2011; Yarushina et al. Reference Yarushina, Bercovici and Michaut2015). In this formulation, $f$ , $\unicode[STIX]{x1D719}_{s}$ and $k_{s}$ are taken as constants. This implies either of the unphysical conditions: depth-invariant lava pressure or zero lava dome thickness. However, we make this assumption to simplify the analysis and give insight into the contrast between the behaviour of the system under uniform and non-uniform background states which we examine in greater detail below. Following the linearization procedure yields

(3.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{s}\,\unicode[STIX]{x2202}_{t}p^{(1)}=\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}\unicode[STIX]{x1D719}^{b}f_{a}\,\unicode[STIX]{x2202}_{zz}p^{(1)}-\unicode[STIX]{x1D6FD}fp^{(1)} & \displaystyle\end{eqnarray}$$

(3.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}^{(1)}=\unicode[STIX]{x1D6FD}(1-\unicode[STIX]{x1D719}_{s})p^{(1)}, & \displaystyle\end{eqnarray}$$

$p^{(1)},\unicode[STIX]{x1D719}^{(1)}$

$(p^{(1)},\unicode[STIX]{x1D719}^{(1)})=\exp [\text{i}\unicode[STIX]{x1D705}z+\unicode[STIX]{x1D706}(\unicode[STIX]{x1D705})t]$

$\unicode[STIX]{x1D705}$

$\unicode[STIX]{x1D706}(\unicode[STIX]{x1D705})$

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}(\unicode[STIX]{x1D705})=-\frac{f}{\unicode[STIX]{x1D719}_{s}}(\unicode[STIX]{x1D6FC}k_{s}\unicode[STIX]{x1D705}^{2}+\unicode[STIX]{x1D6FD})\leqslant 0, & & \displaystyle\end{eqnarray}$$

which implies that for all physical values of the parameters, the background state is stable and perturbations will decay. Additionally, because $\unicode[STIX]{x1D706}$ is purely real, perturbations decay in place, eliminating any wave motion. However, the simplicity of a constant background state does not capture the reality of background gradients present in the steady-state solutions.

3.2 Non-uniform background state

Using the properties of the steady states (2.12a ), (2.12b ), the linearization ansatz yields (after much simplification)

(3.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{\unicode[STIX]{x1D6FC}}\,\unicode[STIX]{x2202}_{t}p^{(1)}=w_{p}\,\unicode[STIX]{x2202}_{z}p^{(1)}+w_{\unicode[STIX]{x1D719}}\,\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}^{(1)}+D_{p}\,\unicode[STIX]{x2202}_{zz}p^{(1)}+R_{pp}p^{(1)}+R_{p\unicode[STIX]{x1D719}}\unicode[STIX]{x1D719}^{(1)} & \displaystyle\end{eqnarray}$$

(3.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{\unicode[STIX]{x1D6FC}}\,\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}^{(1)}=R_{\unicode[STIX]{x1D719}p}p^{(1)}, & \displaystyle\end{eqnarray}$$

(3.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{z}p^{(1)}=0\quad \text{on }\{z=0\}\times [0,\infty ] & \displaystyle\end{eqnarray}$$

(3.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle p^{(1)}=0\quad \text{on }\{z=1\}\times [0,\infty ] & \displaystyle\end{eqnarray}$$

(3.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle w_{p}(z)=-k_{s}(z){\unicode[STIX]{x1D719}_{s}(z)}^{-1} & \displaystyle\end{eqnarray}$$

(3.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle w_{\unicode[STIX]{x1D719}}(z)=-b\unicode[STIX]{x0394}\unicode[STIX]{x1D719}^{b-1}f_{a}^{(b-1)/b}f(z)^{1/b}{\unicode[STIX]{x1D719}_{s}(z)}^{-1} & \displaystyle\end{eqnarray}$$

(3.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle D_{p}(z)=\unicode[STIX]{x0394}\unicode[STIX]{x1D719}^{b}f_{a}{\unicode[STIX]{x1D719}_{s}(z)}^{-1} & \displaystyle\end{eqnarray}$$

(3.6d ) $$\begin{eqnarray}\displaystyle & \displaystyle R_{pp}(z)=-[k_{s}(z)f(z)^{-1}+\unicode[STIX]{x1D6FE}f(z)]{\unicode[STIX]{x1D719}_{s}(z)}^{-1} & \displaystyle\end{eqnarray}$$

(3.6e ) $$\begin{eqnarray}\displaystyle & \displaystyle R_{p\unicode[STIX]{x1D719}}(z)=-k_{s}(z)^{(b-1)/b}{\unicode[STIX]{x1D719}_{s}(z)}^{-1} & \displaystyle\end{eqnarray}$$

(3.6f ) $$\begin{eqnarray}\displaystyle & \displaystyle R_{\unicode[STIX]{x1D719}p}(z)=\unicode[STIX]{x1D6FE}[1-\unicode[STIX]{x1D719}_{s}(z)], & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}=\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6FC}}=\frac{\unicode[STIX]{x1D707}h_{0}^{2}}{k_{0}\unicode[STIX]{x1D701}}=\frac{h_{0}^{2}}{\unicode[STIX]{x1D6FF}^{2}}. & & \displaystyle\end{eqnarray}$$

For a particular dome height $h_{0}$ , the parameter $\unicode[STIX]{x1D6FE}$ balances the permeability $k_{0}$ and volume viscosity $\unicode[STIX]{x1D701}$ of the lava medium or alternately, it represents a comparison of the dome thickness and compaction length. Additionally, we define a new scaled time $t_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D6FC}t$ and write the system as a matrix equation:

(3.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t_{\unicode[STIX]{x1D6FC}}}\boldsymbol{u}=\unicode[STIX]{x1D652}\,\unicode[STIX]{x2202}_{z}\boldsymbol{u}+\unicode[STIX]{x1D63F}\,\unicode[STIX]{x2202}_{zz}\boldsymbol{u}+\unicode[STIX]{x1D64D}\boldsymbol{u}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}=(p^{(1)},\unicode[STIX]{x1D719}^{(1)})^{\intercal }$ , $\unicode[STIX]{x1D652}=\binom{w_{p}\;\;w_{\unicode[STIX]{x1D719}}}{0\;\;\;\;0}$ , $\unicode[STIX]{x1D63F}=\binom{D_{p}\;\;\;0}{0\;\;\;\;0}$ and $\unicode[STIX]{x1D64D}=\binom{R_{pp}\;\;R_{p\unicode[STIX]{x1D719}}}{R_{\unicode[STIX]{x1D719}p}\;\;\;\;0}$ . In this case, the plane-wave ansatz is unjustified because all of the coefficients in (3.8) depend on position as the waves are propagating through a non-uniform background. Here, we are concerned mainly with the growth or decay of waves and so we convert this problem to an eigenvalue problem for the wave growth rates.

3.2.1 LSA as an eigenproblem

We seek a separated solution of the form: $\boldsymbol{u}(z,t)=\text{e}^{\unicode[STIX]{x1D706}t_{\unicode[STIX]{x1D6FC}}}\hat{\boldsymbol{u}}(z)$ , where the growth constant for the unscaled non-dimensional time can be obtained from the solution as $\boldsymbol{u}(z,t)=\text{e}^{\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FC}t}\hat{\boldsymbol{u}}(z)$ . This leads to the generalized eigenvalue problem

(3.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}\hat{\boldsymbol{u}}(z)={\mathcal{L}}\hat{\boldsymbol{u}}(z), & & \displaystyle\end{eqnarray}$$

where the linear operator ${\mathcal{L}}=\unicode[STIX]{x1D652}(z)\,\unicode[STIX]{x2202}_{z}+\unicode[STIX]{x1D63F}(z)\,\unicode[STIX]{x2202}_{zz}+\unicode[STIX]{x1D64D}(z)$ . We solve this eigenproblem numerically because of the variable coefficients on all derivatives, meaning that the eigenfunctions ( $\hat{\boldsymbol{u}}(z)$ ) will not be plane waves. Instability of the steady-state flow (wave growth) is given if any eigenvalues exist with $\text{Re}\{\unicode[STIX]{x1D706}\}>0$ , whereas linear stability is given if no such eigenvalues exist.

3.2.2 Numerical solution to eigenvalue problem

The above eigenvalue problem is solved numerically by discretizing the linear operator with central finite differences (FD) and the coefficients are evaluated at $z=j\unicode[STIX]{x0394}z$ :

(3.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}\hat{\boldsymbol{u}}_{j}=\frac{1}{2\unicode[STIX]{x0394}z}\unicode[STIX]{x1D652}_{j}(\hat{\boldsymbol{u}}_{j+1}-\hat{\boldsymbol{u}}_{j-1})+\frac{1}{(\unicode[STIX]{x0394}z)^{2}}\unicode[STIX]{x1D63F}_{j}(\hat{\boldsymbol{u}}_{j+1}-2\hat{\boldsymbol{u}}_{j}+\hat{\boldsymbol{u}}_{j-1})+\unicode[STIX]{x1D64D}_{j}\hat{\boldsymbol{u}}_{j}, & & \displaystyle\end{eqnarray}$$

with boundary conditions at $z=0$ ( $j=1$ ) and $z=1$ ( $j=J$ ) as specified by (3.5a ) and (3.5b ). Because the discrete representation of ${\mathcal{L}}$ is a $2(J-2)\times 2(J-2)$ block-tridiagonal matrix, there will be $J-2$ pairs of eigenvalues corresponding to the growth rates of the discretized eigenfunctions.

Figure 3.

(a) Distribution of a subset of eigenvalues (growth parameters) in the complex plane for $f_{a}=0.05,\unicode[STIX]{x1D719}_{a}=0.5$ . The colours represent different values of $\unicode[STIX]{x1D6FE}$ . The eigenvalue with largest real part at each value of $\unicode[STIX]{x1D6FE}$ is represented by a larger, brighter circle, corresponding to the dominant wave that will grow for that particular $\unicode[STIX]{x1D6FE}$ . Inset figure shows the distribution of all of the eigenvalues from which those in the main figure are subset. The approximate range of the inset figure is $-2500<\text{Re}\{\unicode[STIX]{x1D706}\}<1,-2<\text{Im}\{\unicode[STIX]{x1D706}\}<2$ . (b) Eigenfunctions associated with the dominant eigenvalues over six values of $\unicode[STIX]{x1D6FE}$ showing the real part (blue, solid), imaginary part (red, dotted) and the region bounded by the modulus (green).

3.2.3 Distribution and variation of eigenvalues and eigenfunctions

For a given set of free parameters in the eigenvalue problem $(\unicode[STIX]{x1D6FE},f_{a},\unicode[STIX]{x1D719}_{a})$ , the eigenvalues ( $\unicode[STIX]{x1D706}_{j}\in \mathbb{C}$ ) are mainly distributed along the negative half – $\text{Re}\{\unicode[STIX]{x1D706}\}$ line; however, a small fraction have non-zero imaginary part ( $\text{Im}\{\unicode[STIX]{x1D706}_{j}\}\neq 0$ ) and are near the stability–instability threshold (figure 3). For some ranges of the parameters, there exists eigenvalues with $\text{Re}\{\unicode[STIX]{x1D706}_{j}\}>0$ . Focusing on these positive eigenvalues, for a given value of $\unicode[STIX]{x1D6FE}$ , they are distributed in two symmetric branches (in the upper and lower half-planes) each of which generally begin near the origin and increase in imaginary and real parts (figure 3 a). In many cases, there is one eigenvalue on each branch with $\text{Re}\{\unicode[STIX]{x1D706}_{j}\}$ significantly larger than nearby eigenvalues. This larger positive eigenvalue corresponds to a dominating growth rate of its associated eigenfunction wave profiles. Because of this, these scenarios represent unstable flow regimes, causing particular excited oscillation to grow with time.

Varying the ratio of expansion rate to diffusion rate ( $\unicode[STIX]{x1D6FE}$ ) and examining the eigenvalues with largest real parts ( $\unicode[STIX]{x1D706}_{max}$ ), it is clear that $\unicode[STIX]{x1D706}_{max}(\unicode[STIX]{x1D6FE})$ increases with $\unicode[STIX]{x1D6FE}$ , meaning that for greater expansion and contraction and lesser diffusion, waves grow more quickly and the instability is stronger (figure 3). Stated in terms of physical parameters, less viscous or less permeable lavas allow greater wave growth than more viscous or permeable lavas. Where permeability is relatively low, pressure anomalies cannot relax sufficiently quickly by porous flow and as a consequence, the porosity will expand or contract as the pressure anomaly grows. Inversely, very permeable lavas allow efficient gas escape, smoothing out steep pressure gradients, limiting the growth of high and low pressure anomalies. Similarly, if the matrix viscosity is relatively low, even a small pressure anomaly will cause a significant change in porosity since over-pressurized regions are allowed to expand and under-pressurized regions allowed to contract against a weaker viscous resisting force in the matrix. These effects become self-reinforcing since the waves travel ( $\text{Im}\{\unicode[STIX]{x1D706}\}\neq 0$ ). In particular, the porosity governing equation predicts that a travelling pressure wave will induce a porosity wave that lags behind as it travels. As a result, the negative phase of a pressure wave lowers the porosity and permeability of the lava and the next arriving positive pressure wave phase will not diffuse through the low-permeability region easily, causing the pressure to rise behind the constriction, resulting in overall wave growth.

Figure 4.

(a, left) The variation of $\unicode[STIX]{x1D706}_{max}(\unicode[STIX]{x1D6FE})$ for different values of $f_{a}$ showing the dependence of the stability–instability threshold on $\unicode[STIX]{x1D6FE}$ and (a, right) the associated steady-state porosity profile corresponding to each value of $f_{a}$ . (b) Contours of $\unicode[STIX]{x1D706}_{max}(\unicode[STIX]{x1D6FE},f_{a})$ for two values of $\unicode[STIX]{x1D719}_{a}$ showing the modification of the instability threshold (bold lines) with surface porosity.

The dominant eigenfunctions at each value of $\unicode[STIX]{x1D6FE}$ show that the maxima of the eigenfunction moduli or the location of the most principally excited waves gets closer to the surface for increasing $\unicode[STIX]{x1D6FE}$ (figure 3 b). Additionally, as $\unicode[STIX]{x1D6FE}$ increases, the pressure and porosity eigenfunctions have the same general shape; however, for the pressure, the maximum of the modulus decreases whereas for the porosity, the maximum increases, meaning that for larger $\unicode[STIX]{x1D6FE}$ , excited pressure waves are proportionally smaller in magnitude than their associated porosity waves. The structure of these eigenfunctions shows that the number of waves in the domain also increases, that is, the mean wavelength of the excited waves decreases with increasing $\unicode[STIX]{x1D6FE}$ . This trend corresponds with the notion that for small $\unicode[STIX]{x1D6FE}$ , the system is diffusion dominated, resulting in larger wavelengths and smoother gradients; however, for larger $\unicode[STIX]{x1D6FE}$ , the system experiences very little diffusion and so high frequency waves are allowed to grow. This instability is not universal for all values of $f_{a}$ and $\unicode[STIX]{x1D719}_{a}$ . For increasing values of $f_{a}$ and smaller $\unicode[STIX]{x1D6FE}$ , $\text{Re}\{\unicode[STIX]{x1D706}_{max}\}$ decreases until all even the largest eigenvalue becomes negative, implying linear stability of the steady state for those parameters (figure 4 a). For each set of stable $f_{a}$ and $\unicode[STIX]{x1D719}_{a}$ , there are stable and unstable values of $\unicode[STIX]{x1D6FE}$ . The critical value at which this threshold occurs ( $\unicode[STIX]{x1D6FE}_{c}$ ) increases with both $f_{a}$ and $\unicode[STIX]{x1D719}_{a}$ (figure 4 b). The complexity of this bifurcation contrasts with the unconditional stability predicted from the assumption of a uniform background state. An important consideration is the fact that as $f_{a}$ increases, the pressure and porosity background states become more approximately uniform and linear stability is predicted over a wider range of $\unicode[STIX]{x1D6FE}$ (figure 4 b). This mirrors the conditions of the LSA for a uniform background and agrees with that analysis. Critically, when $f_{a}$ is small, significant gradients form in the steady-state porosity field resulting from gravitational compaction of the column (figure 4 b). When $f_{a}$ is large, the dome is short and there is less overburden on the base of the dome which experiences less gravitational compaction. From these considerations, it is clear that increased dome height causes the porous lava to compact more, creating larger gradients in the dome’s porosity, which in turn initiate instability onset. Although this generalized eigenvalue problem yields insight into the generation of fluid instabilities in this system, significant simplifications have been made in reducing the highly nonlinear governing equations to a linear form. Consequently, the LSA above does not provide insight into higher-order processes which may limit the initial wave growth nor does it yield any insight into the dynamics near the $\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{c}^{+}$ limit. Building from these simple analyses, we will next solve the one-dimensional (1-D) system ((2.7a ) and (2.7b )) numerically to study the nonlinear effects that follow from the linear wave growth analysed here and gain insight into general characteristics of solutions.

4.1 Finite-difference scheme

In solving this system numerically, we employ a modified backward Euler finite-difference (BEFD) scheme based on the well-established methods of Huyakorn & Pinder (Reference Huyakorn and Pinder1983) for modelling multiphase subsurface flow problems. To solve the coupled equations for pressure and porosity, equations (4.1a ) and (4.1b ) are discretized on a staggered, uniform grid where the pressures are located at the nodes ( $z=j\unicode[STIX]{x0394}z$ ) and the porosity and fluxes evaluated between the nodes ( $z=(j\pm 1/2)\unicode[STIX]{x0394}z$ ) on the interval $z\in [0,1]=\unicode[STIX]{x0394}z\{1,\ldots ,J\}$ . This approach respects the fact that variables like permeability and porosity control the flux of gas through the dome and thus should be located between nodes. Consequently, all pressure terms in the porosity equation are necessarily averaged between the two nearest nodes to the porosity half-node sites. For the pressure equation, since the flux is evaluated at the half-nodes, the pressure in the flux is similarly averaged. This local averaging in the numerical scheme can be thought of as a local volume element averaging as in Whitaker (Reference Whitaker1999). Using the staggered grid and local volume element averaging ensures that in simulating of the steady-state conditions, the diffusivity in the pressure equation ( $K(\unicode[STIX]{x1D719})p$ ) remains constant and reproduces the steady state upon updating.

The main modifications to the BEFD scheme is that the compaction terms are evolved according to forward differencing and the time steps are dynamic and are performed with a mid-time predictor–corrector step. The dynamic time stepping is based on the Jacobian of the compaction term in (4.1a ) at the last time step, that is

(4.3) $$\begin{eqnarray}\displaystyle t^{n+1}-t^{n}=\unicode[STIX]{x0394}t^{n+1}=0.1\left[\max _{j\in \{2,\ldots ,J-1\}}|(J_{p})_{j}^{n}|\right]^{-1},\quad \text{where }(J_{p})_{j}^{n}=\unicode[STIX]{x1D6FD}\frac{f_{j}-2p_{j}^{n}}{\unicode[STIX]{x1D719}_{j}^{n}}. & & \displaystyle\end{eqnarray}$$

Due to the nonlinearity in the diffusion term, for numerical solution to be tractable, this scheme is partially linearized by evaluating the non-constant diffusivity ( $K(\unicode[STIX]{x1D719})p$ ) at the previous time step. In the discretization of (4.1a ), the derivatives of the fluxes are evaluated as the difference across the edges of the cells:

(4.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \unicode[STIX]{x1D719}_{j+1/2}^{\ast }\rangle _{-}\frac{p_{j}^{n+1}-p_{j}^{n}}{\unicode[STIX]{x0394}t^{n+1}}=\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x0394}z}(F_{j+1/2}^{n+1}-F_{j-1/2}^{n+1})-\unicode[STIX]{x1D6FD}p_{j}^{\ast }(p_{j}^{\ast }-f_{j}) & \displaystyle\end{eqnarray}$$

(4.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D719}_{j\pm 1/2}^{n+1}-\unicode[STIX]{x1D719}_{j\pm 1/2}^{n}}{\unicode[STIX]{x0394}t^{n+1}}=\unicode[STIX]{x1D6FD}(1-\unicode[STIX]{x1D719}_{j\pm 1/2}^{\ast })\langle p_{j}^{\ast }-f_{j}\rangle _{\pm }, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle F_{j\pm 1/2}^{n+1}=K_{j\pm 1/2}^{\ast }\langle p_{j}^{\ast }\rangle _{\pm }\left(\pm \frac{p_{j\pm 1}^{n+1}-p_{j}^{n+1}}{\unicode[STIX]{x0394}z}\right) & & \displaystyle\end{eqnarray}$$

and the local volume averages are $\langle (\cdot )_{j}\rangle _{\pm }=1/2[(\cdot )_{j}+(\cdot )_{j\pm 1}]$ , centred at $z=(j\pm 1/2)\unicode[STIX]{x0394}z$ . Boundary values for (4.4a ) are discretized as

(4.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{p_{2}^{n+1}-p_{1}^{n+1}}{\unicode[STIX]{x0394}z}=-\frac{\unicode[STIX]{x1D6EC}_{s}}{K_{1+1/2}^{\ast }\langle p_{1}^{\ast }\rangle _{+}}[1+q\sin (2\unicode[STIX]{x03C0}t^{n+1})] & \displaystyle\end{eqnarray}$$

(4.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle p_{J}^{n+1}=f_{a}. & \displaystyle\end{eqnarray}$$

$\ast =n$

$n+1/2$

$p_{j}^{n+1/2}$

$\unicode[STIX]{x1D719}_{j\pm 1/2}^{n+1/2}$

4.2 Wave-diffusion behaviour

Solutions of (4.4a ) and (4.4b ) oscillate around the steady-state solutions $p_{s}(z)$ and $\unicode[STIX]{x1D719}_{s}(z)$ in a series of diffusive waves propagating toward the surface. As pressure waves pass, there are concordant porosity waves that lag just behind the pressure pulses (figure 5). Early in simulations, waves grow in a manner similar to the eigenfunctions of the LSA above; however, if the solution is continued, these waves grow and steepen as they propagate toward the surface (figure 5). This is an indication that nonlinear processes in the diffusion term of the pressure equation are dominating the simple linear wave growth predicted by the LSA. In general, the propagating pressure waves are steepest in the deepest troughs of the porosity waves. The amplitude of these waves continues to grow until $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{c}$ , where the next arriving pressure wave steepens into a discontinuity (shock) (figure 5). In the regions of $\unicode[STIX]{x1D719}\leqslant \unicode[STIX]{x1D719}_{c}$ , the pressure and porosity evolve as predicted by (2.11), with the pressure increasing asymptotically back to the local magmastatic pressure ( $f(z)$ ) from $p=p_{c}(z)$ . Once multiple shocks develop, the excess pressure ( $p-f$ ) takes the form of a series of ‘N-waves’ separated by regions of magmastatic pressure in which $\unicode[STIX]{x1D719}<\unicode[STIX]{x1D719}_{c}$ (figure 5). The shocks generally slow and decay slightly as they propagate to shallower levels in the domain. Because the N-waves are separated by regions which are impermeable, there is no flow from one shock-bounded region to the next. Although no systematic study of the Rankine–Hugoniot relations for shock propagation speed has been made here owing to the complexity of the flux term, it is noted that in the simulations, successive N-waves rarely interact. These rare interactions always result in the up-gradient N-wave quickly diffusing into the down-gradient N-wave as the impermeable region separating the converging shocks becomes narrower until it vanishes (figure 5 a,b (right panels)). This behaviour is merely a consequence of the sudden hydraulic connection of two isolated permeable domains with different mean pressures. Although shock development is common across a wide set of free parameters in the problem, this wave behaviour is not universal. Still, in cases without shock development, waves steepen as they propagate toward the atmospheric boundary.

Figure 5.

Numerical solution $(f_{a}=0.05,\unicode[STIX]{x1D719}_{a}=0.5,\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FD}=10,q=0.1)$ showing wave propagation, steepening and shock formation for porosity (a) and pressure (b) profiles. (a, left) Dashed lines on first and last porosity profiles represent steady-state porosity, the dotted line on first profile represents the percolation threshold $\unicode[STIX]{x1D719}_{c}$ , and the black regions represent the areas with $\unicode[STIX]{x1D719}<\unicode[STIX]{x1D719}_{c}$ . (b, left) Light and dark grey shaded regions represent positive and negative pressure anomalies respectively. (a, right) and (b, right) show the behaviour of the porosity and pressure respectively late in the simulation during an interaction of two shocks between $0\leqslant z/h_{0}\leqslant 0.25$ .

5.1 Porosity-dependent fragmentation threshold

From experimental studies of lava fragmentation, Spieler et al. (Reference Spieler, Kennedy, Kueppers, Dingwell, Scheu and Taddeucci2004) determined an empirical fragmentation threshold relationship that relates the pressure difference required to fragment a sample ( $\unicode[STIX]{x0394}p$ ) to the sample porosity as

(5.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}p=\frac{\unicode[STIX]{x1D70E}_{m}}{\unicode[STIX]{x1D719}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{m}$ is taken as 1 MPa, which they noted should not be expected to hold for highly porous rocks. This relationship is an alternative to earlier studies which had concluded that porosity was the main control on fragmentation, considering fragmentation thresholds of porosities between 65 and 85 % (Sparks Reference Sparks1978; Sparks et al. Reference Sparks, Barclay, Jaupart, Mader and Phillips1994). As a result we modify the threshold of Spieler et al. (Reference Spieler, Kennedy, Kueppers, Dingwell, Scheu and Taddeucci2004) to respect the failure of this relationship for large porosity:

(5.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}p(\unicode[STIX]{x1D719})=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{1~\text{MPa}}{\unicode[STIX]{x1D719}}\quad & \unicode[STIX]{x1D719}<0.8\\ 0\quad & \unicode[STIX]{x1D719}\geqslant 0.8.\end{array}\right. & & \displaystyle\end{eqnarray}$$

We incorporate this threshold into the numerical model presented above by predicting failure at time $t_{f}=t^{n+1}$ if

(5.3) $$\begin{eqnarray}\displaystyle p_{j}^{n+1}-f_{j}\geqslant \unicode[STIX]{x0394}p(\langle \unicode[STIX]{x1D719}_{j+1/2}^{n+1}\rangle _{-}). & & \displaystyle\end{eqnarray}$$

In doing so, we can examine the predicted failure times across a range of representative parameters for intermediate-composition lava domes (table 1).

In the development of the LSA, we showed that for initial wave growth, the diffusion rate parameter $\unicode[STIX]{x1D6FC}$ mainly scales time during wave growth whereas the ratio of expansion rate to diffusion rate ( $\unicode[STIX]{x1D6FE}$ ) controls the behaviour of these waves. In the nonlinear case, $\unicode[STIX]{x1D6FC}$ is expected to not only scale time, but also control some nonlinear effects due to the severe nonlinearity of the diffusion in (2.7a ) and so $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are varied independently in the ranges $0\leqslant \log (\unicode[STIX]{x1D6FC})\leqslant 4$ and $-1\leqslant \log (\unicode[STIX]{x1D6FD})\leqslant 2$ . Additionally, the boundary flux oscillation strength is varied in the range $0.1\leqslant q\leqslant 1$ . As a result, we consider the failure time as a function of the three-dimensional sample space: $t_{f}:S\subset \mathbb{R}^{3}\rightarrow \mathbb{R}$ .

Figure 6.

(a–c) The variation of $t_{f}$ with $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ and $q$ for several values of the other two parameters. (d) Logarithmic contour plot of $t_{f}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ generated by averaging values of $t_{f}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},q)$ through $q$ and a single convolution with a $3\times 3$ running average filter. Red dotted contours represent values of $\unicode[STIX]{x1D6FE}$ .

5.2 Parametric dependence of dome failure times

As anticipated by the derivation of the LSA, the diffusion rate parameter $\unicode[STIX]{x1D6FC}$ has the strongest control on the failure time followed by the compaction rate parameter $\unicode[STIX]{x1D6FD}$ and relative boundary oscillation amplitude $q$ . Although $q\neq 0$ is the only condition generating the non-equilibrium conditions, variation of $q$ explains almost none of the variability in computed survival times ( $t_{f}$ ) compared with $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ (figure 6 a–c). For example, the full domain of $q$ is only able to change $t_{f}(q)$ by less than one order of magnitude (figure 6 a) compared with ${\sim}2$ orders of magnitude for $\unicode[STIX]{x1D6FD}$ (figure 6 b) and $1-2$ orders of magnitude for $\unicode[STIX]{x1D6FC}$ (figure 6 c). Because of this, we principally analyse the variation of $t_{f}$ through the parameter space of $\{\unicode[STIX]{x1D6FC}\}\times \{\unicode[STIX]{x1D6FD}\}$ only, averaging $t_{f}$ through the full space along $q$ (figure 6 d). The computed function $t_{f}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ is only defined on a subset of the parameter space approximately where $\unicode[STIX]{x1D6FE}\geqslant O(10^{-2})$ (figure 6 d). The excluded region includes parameter values for which initial instability causes wave growth as predicted by the LSA; however, higher order terms limit the growth and the pressure and porosity waves are bounded below the fragmentation threshold. Analysis of these solutions is beyond the scope of the present work. Within the domain of finite fragmentation times, increased $\unicode[STIX]{x1D6FD}$ decreases the survival time with an approximately exponential decay. By contrast, increasing $\unicode[STIX]{x1D6FC}$ decreases $t_{f}$ on approximately $\unicode[STIX]{x1D6FC}<10\unicode[STIX]{x1D6FD}$ ; however, with higher $\unicode[STIX]{x1D6FC}$ , the survival time grows quickly until it meets the undefined region in which the higher-order effects prevent fragmentation. Although the value of $\unicode[STIX]{x1D6FE}$ plays a large role in the initial wave growth and in the behaviour of the excited waves, it does not seem to exert a strong control over the timing of failure by itself.