2.1 Dimensional equations

Figure 1 shows a schematic diagram of the domain: a region of partially molten rock of height $H$ in the $z$ direction, composed of a solid phase ( $s$ , matrix, mantle rock) and a liquid phase ( $l$ , magma).

Figure 1.

Diagram of the problem. (a) A region of partially molten rock of depth $H$ with a flux of liquid phase from beneath and free-flux boundary condition above. (b) The gradient $\unicode[STIX]{x1D6FD}$ of the equilibrium composition of the liquid phase $c_{l}^{eq}$ . When a parcel of liquid (dashed blue circle) is raised (full blue circle), it has a concentration below the equilibrium, leading to reactive melting along the horizontal blue arrow. The composition of reactively produced melts $c_{\unicode[STIX]{x1D6E4}}$ is greater than equilibrium by a constant offset $\unicode[STIX]{x1D6FC}$ .

We account for conservation in both phases. Mass conservation in the solid and liquid is given by respectively

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(1-\unicode[STIX]{x1D719})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }((1-\unicode[STIX]{x1D719})\boldsymbol{v}_{s})=-\unicode[STIX]{x1D6E4}, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}\boldsymbol{v}_{l})=\unicode[STIX]{x1D6E4}, & \displaystyle\end{eqnarray}$$

$t$

$\unicode[STIX]{x1D719}$

$(1-\unicode[STIX]{x1D719})$

$\boldsymbol{v}_{l}$

$\boldsymbol{v}_{s}$

$\unicode[STIX]{x1D6E4}$

We use conservation of momentum to determine the solid and liquid velocities (McKenzie Reference McKenzie1984). In general, the solid phase (mantle) can deform viscously by both deviatoric shear and isotropic compaction. The latter is related to the pressure difference between the liquid and solid phases. We neglect deviatoric stresses on the solid phase and consider only the isotropic part of the stress and strain-rate tensors. The compaction rate $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{s}$ is related to the compaction pressure ${\mathcal{P}}$ according to the linear constitutive law

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{s}={\mathcal{P}}/\unicode[STIX]{x1D701},\end{eqnarray}$$

where $\unicode[STIX]{x1D701}$ is an effective compaction or bulk viscosity. The solid matrix behaves like a rigid porous matrix when the bulk viscosity is sufficiently large (an idea we will relate to a non-dimensional matrix stiffness later). The effective compaction or bulk viscosity $\unicode[STIX]{x1D701}$ can be estimated using micromechanical models of partially molten rocks, and may depend on the porosity (Sleep Reference Sleep1988). The most recent calculations show that the bulk viscosity depends only weakly on porosity (Rudge Reference Rudge2017). Therefore, we make the simplifying assumption that $\unicode[STIX]{x1D701}$ is constant. We discuss this issue further in appendix C.

Fluid flow is given by Darcy’s law,

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D719}(\boldsymbol{v}_{l}-\boldsymbol{v}_{s})=K[(1-\unicode[STIX]{x1D719})\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}g\hat{\boldsymbol{z}}-\unicode[STIX]{x1D735}{\mathcal{P}}].\end{eqnarray}$$

A Darcy flux $\unicode[STIX]{x1D719}(\boldsymbol{v}_{l}-\boldsymbol{v}_{s})$ is driven by gravity $g\hat{\boldsymbol{z}}$ associated with the density difference between the phases $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ and by compaction pressure gradients. Crucially, the mobility $K$ ( $\equiv$ permeability divided by liquid viscosity) of the liquid depends on the porosity,

(2.4) $$\begin{eqnarray}K=K_{0}(\unicode[STIX]{x1D719}/\unicode[STIX]{x1D719}_{0})^{n},\end{eqnarray}$$

where $K_{0}$ is a reference mobility at a reference porosity $\unicode[STIX]{x1D719}_{0}$ (equal to the porosity at the base of the column, $z=0$ ) and $n$ is a constant (we take $n=3$ in our numerical calculations). It is thought that $2\leqslant n\leqslant 3$ for the geological systems of interest (von Bargen & Waff Reference von Bargen and Waff1986; Miller et al. Reference Miller, Zhu, Montési and Gaetani2014; Rudge Reference Rudge2018).

Finally, we must determine the melting or reaction rate $\unicode[STIX]{x1D6E4}$ . The focus of this paper is the mechanics of the instability, so we adopt a fairly simple treatment of its chemistry, largely following Aharonov et al. (Reference Aharonov, Whitehead, Kelemen and Spiegelman1995). The reaction associated with the reaction-infiltration instability is one of chemical dissolution. At a simple level, this can be described as follows. As magma rises, its pressure decreases and it becomes undersaturated in silica. This, in turn, drives a reaction in which pyroxene is dissolved from the solid while olivine is precipitated (cf. figure 8 in Longhi Reference Longhi2002). Schematically, the dissolution reaction can be written as

(2.5) $$\begin{eqnarray}\text{magma}_{1}(l)+\text{pyroxene}(s)\rightarrow \text{magma}_{2}(l)+\text{olivine}(s),\end{eqnarray}$$

where ( $l$ ) denotes a component in the liquid phase and ( $s$ ) a component in the solid phase, and we use subscripts 1 and 2 to indicate magmas of slightly different compositions. Crucially, this reaction involves a net transfer of mass from solid to liquid (Kelemen Reference Kelemen1990), and hence it is typically called a melting reaction. Because the reaction replaces pyroxene with olivine, geological observations of tabular dunite bodies in exhumed mantle rock are interpreted as evidence for the reaction-infiltration instability (dunites are mantle rocks, residual after partial melting, that are nearly pure olivine) (Kelemen et al. Reference Kelemen, Dick and Quick1992).

We now formulate the reactive chemistry in terms of the simplest possible mathematics. We assume that $\unicode[STIX]{x1D6E4}$ is proportional to the undersaturation of a soluble component in the melt. The concentration of this component in the melt is denoted $c_{l}$ ; the equilibrium concentration is denoted $c_{l}^{eq}$ . Hence, the melting rate is written as

(2.6) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}=-R(c_{l}-c_{l}^{eq}),\end{eqnarray}$$

where $R$ is a kinetic coefficient with units 1/time. We assume that $R$ is a constant, independent of the concentration of the soluble component in the solid phase. This is valid for the purposes of studying the onset of instability if the soluble component is abundant and homogeneously distributed, both reasonable assumptions (Liang et al. Reference Liang, Schiemenz, Hesse, Parmentier and Hesthaven2010).

In this formulation, the chemical reaction rate depends on the composition of the liquid phase $c_{l}$ . Chemical species conservation in the liquid phase is given by

(2.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D719}c_{l})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}\boldsymbol{v}_{l}c_{l})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}D\unicode[STIX]{x1D735}c_{l})+\unicode[STIX]{x1D6E4}c_{\unicode[STIX]{x1D6E4}},\end{eqnarray}$$

where the effective diffusivity of chemical species is $\unicode[STIX]{x1D719}D$ (diffusivity in the liquid phase is written $D$ ; diffusion through the solid phase is negligible) and $c_{\unicode[STIX]{x1D6E4}}$ is the concentration of reactively produced melts. We then expand out the partial derivatives and simplify using (2.1b ) to obtain

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}c_{l}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D719}\boldsymbol{v}_{l}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c_{l}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}D\unicode[STIX]{x1D735}c_{l})+(c_{\unicode[STIX]{x1D6E4}}-c_{l})\unicode[STIX]{x1D6E4}.\end{eqnarray}$$

To close the system, we suppose that the equilibrium concentration has a constant gradient $\unicode[STIX]{x1D6FD}\hat{\mathbf{z}}$ , as shown in figure 1. If we define (without loss of generality) the equilibrium concentration at the base of the region ( $z=0$ ) to be zero, then

(2.9) $$\begin{eqnarray}c_{l}^{eq}=\unicode[STIX]{x1D6FD}z.\end{eqnarray}$$

We suppose further that the concentration $c_{\unicode[STIX]{x1D6E4}}$ of the reactively produced melts is offset from the equilibrium concentration by $\unicode[STIX]{x1D6FC}$ , a positive constant, so

(2.10) $$\begin{eqnarray}c_{\unicode[STIX]{x1D6E4}}=\unicode[STIX]{x1D6FD}z+\unicode[STIX]{x1D6FC}.\end{eqnarray}$$

A scaling argument clarifies the meaning of the compositional parameters: for a fast reaction ( $R\rightarrow \infty$ ), and hence for a liquid that is close to equilibrium, a vertical liquid flux $f_{0}$ would cause reactive melting at a characteristic rate $\unicode[STIX]{x1D6E4}_{0}\sim f_{0}\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC}$ , so $\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC}$ is the rate of reactive melting per unit of liquid flux. Our formulation of $c_{\unicode[STIX]{x1D6E4}}$ is slightly different from that of Aharonov et al. (Reference Aharonov, Whitehead, Kelemen and Spiegelman1995), who take $c_{\unicode[STIX]{x1D6E4}}=1$ . Their resulting simplified equations are equivalent to ours when $\unicode[STIX]{x1D6FC}=1$ (following the non-dimensionalization in our § 2.2).

At this point, we remark briefly on two simplifications inherent in the approach described above. First, we assume that the equilibrium chemistry of the liquid phase is a function of depth. A fuller treatment might consider the chemistry of the liquid as a function of pressure (Longhi Reference Longhi2002). However, to an excellent approximation, the liquid pressure is equal to the lithostatic pressure $\unicode[STIX]{x1D70C}_{s}g(H-z)$ , in which case pressure and depth are linearly related. Indeed, the dimensionless error in making this approximation is $O({\mathcal{S}}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{s})$ , where ${\mathcal{S}}$ is the matrix stiffness parameter introduced below. Thus, we neglect the difference relative to lithostatic pressure consistent with a Boussinesq approximation $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{s}\ll 1$ , where $\unicode[STIX]{x1D70C}_{s}$ is the density of the solid phase.

Second, we use a very simple treatment of melting that neglects, for example, latent heat and temperature variations. Hewitt (Reference Hewitt2010) developed a consistent thermodynamic model of melting and showed that latent heat may suppress instability because it reduces the melting rate. Such an effect can be represented within our simpler model by reducing the melting-rate factor $\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC}$ (see further discussion in appendix C).

2.2 Simplified non-dimensional equations

The governing equations (2.1a ), (2.1b ), (2.3) and (2.8) can be non-dimensionalized according to the characteristic scales

(2.11) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle [x,z]=H,\quad [\unicode[STIX]{x1D719}]=\unicode[STIX]{x1D719}_{0},\\ \displaystyle [\boldsymbol{v}_{l}]=w_{0}\equiv K_{0}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}g/\unicode[STIX]{x1D719}_{0},\quad [\boldsymbol{v}_{s}]=\unicode[STIX]{x1D719}_{0}w_{0},\quad [t]=\unicode[STIX]{x1D6FC}/(w_{0}\unicode[STIX]{x1D6FD}),\\ \displaystyle [{\mathcal{P}}]=\unicode[STIX]{x1D701}\unicode[STIX]{x1D719}_{0}w_{0}\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC},\quad [c_{l}]=\unicode[STIX]{x1D6FD}H,\quad [\unicode[STIX]{x1D6E4}]=\unicode[STIX]{x1D719}_{0}w_{0}\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC}.\end{array}\right\}\end{eqnarray}$$

The dimensionless parameters of the system are as follows. First, ${\mathcal{M}}=\unicode[STIX]{x1D6FD}H/\unicode[STIX]{x1D6FC}\ll 1$ , which is the change in solubility across the domain height and characterizes the reactivity of the system. Second, stiffness ${\mathcal{S}}={\mathcal{M}}\unicode[STIX]{x1D6FF}^{2}/H^{2}$ , which characterizes the rigidity of the medium, where $\unicode[STIX]{x1D6FF}=\sqrt{K_{0}\unicode[STIX]{x1D701}}$ is the dimensional compaction length, an emergent length scale (e.g. Spiegelman Reference Spiegelman1993). Third, $Da=\unicode[STIX]{x1D6FC}RH/(\unicode[STIX]{x1D719}_{0}w_{0})\gg 1$ is the Damköhler number, which characterizes the importance of reaction relative to advection. Fourth, $Pe=w_{0}H/D\gg 1$ is the Péclet number, which characterizes the importance of advection relative to diffusion.

Then, the equations can be simplified by taking the limit of small porosity, $\unicode[STIX]{x1D719}_{0}\ll {\mathcal{M}}\ll 1$ , and considering only horizontal diffusion (because we expect channelized features with a short horizontal wavelength compared with their vertical structure). We also assume that the reaction rate is fast, so we neglect terms of $O({\mathcal{M}}/Da)\ll 1$ . We also expand out the divergence term in (2.1a ) using (2.2). Thus, the four governing equations (2.1a ), (2.1b ), (2.3) and (2.8) become

(2.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}={\mathcal{P}}+\unicode[STIX]{x1D712}, & \displaystyle\end{eqnarray}$$

(2.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}\boldsymbol{v}_{l})={\mathcal{M}}\unicode[STIX]{x1D712}, & \displaystyle\end{eqnarray}$$

(2.12c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}\boldsymbol{v}_{l}=K[\hat{\boldsymbol{z}}-{\mathcal{S}}\unicode[STIX]{x1D735}{\mathcal{P}}], & \displaystyle\end{eqnarray}$$

(2.12d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}\boldsymbol{v}_{l}\boldsymbol{\cdot }\left[\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}}{Da}-\hat{\boldsymbol{z}}\right]=\frac{1}{DaPe}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}{\unicode[STIX]{x2202}x}\right)-\unicode[STIX]{x1D712}, & \displaystyle\end{eqnarray}$$

$K=\unicode[STIX]{x1D719}^{n}$

$\unicode[STIX]{x1D712}$

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D712}=Da(z-c_{l}).\end{eqnarray}$$

The dimensionless reactive melting rate is equal to the scaled undersaturation $\unicode[STIX]{x1D712}$ .

A set of appropriate boundary conditions is

(2.14a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D719}=1,\quad \unicode[STIX]{x1D712}=1,\quad \frac{\unicode[STIX]{x2202}{\mathcal{P}}}{\unicode[STIX]{x2202}z}=0,\quad (z=0),\end{eqnarray}$$

(2.14d ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}{\mathcal{P}}}{\unicode[STIX]{x2202}z}=0,\quad (z=1).\end{eqnarray}$$

The boundary conditions at $z=0$ combine with (2.12c ) to give an incoming vertical liquid velocity $w=1$ . At the upper boundary, there is no driving compaction pressure gradient (a ‘free-flux’ condition).

3.1 The base state

The leading-order flow is purely vertical. The conservation equations at this order are

(3.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle 0={\mathcal{P}}_{0}+\unicode[STIX]{x1D712}_{0}, & \displaystyle\end{eqnarray}$$

(3.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}z}(\unicode[STIX]{x1D719}_{0}w_{0})={\mathcal{M}}\unicode[STIX]{x1D712}_{0}, & \displaystyle\end{eqnarray}$$

(3.2c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{0}w_{0}=K_{0}\left[1-{\mathcal{S}}\frac{\text{d}{\mathcal{P}}_{0}}{\text{d}z}\right], & \displaystyle\end{eqnarray}$$

(3.2d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{0}w_{0}\left[\frac{1}{Da}\frac{\text{d}\unicode[STIX]{x1D712}_{0}}{\text{d}z}-1\right]=-\unicode[STIX]{x1D712}_{0}, & \displaystyle\end{eqnarray}$$

$K_{0}=\unicode[STIX]{x1D719}_{0}^{n}$

$Da$

${\mathcal{P}}_{0}=-\unicode[STIX]{x1D712}_{0}$

$\unicode[STIX]{x1D712}_{0}=\exp ({\mathcal{M}}z)$

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

$w_{0}$

${\mathcal{M}}\ll 1$

$\exp ({\mathcal{M}}z)\approx 1$

(3.3) $$\begin{eqnarray}-{\mathcal{P}}_{0}=\unicode[STIX]{x1D712}_{0}=\unicode[STIX]{x1D719}_{0}=w_{0}=1.\end{eqnarray}$$

The uniformity of the base state significantly simplifies the subsequent analysis.

3.2 Perturbation equations

The equations governing the perturbations can be written as

(3.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t}={\mathcal{P}}_{1}+\unicode[STIX]{x1D712}_{1}, & \displaystyle\end{eqnarray}$$

(3.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D719}_{0}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{1}+w_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}={\mathcal{M}}\unicode[STIX]{x1D712}_{1}, & \displaystyle\end{eqnarray}$$

(3.4c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{0}\boldsymbol{v}_{1}=-{\mathcal{S}}K_{0}\unicode[STIX]{x1D735}{\mathcal{P}}_{1}+(n-1)w_{0}\unicode[STIX]{x1D719}_{1}\hat{\boldsymbol{z}}, & \displaystyle\end{eqnarray}$$

(3.4d ) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D719}_{0}w_{1}+\unicode[STIX]{x1D719}_{1}w_{0})\left[\frac{1}{Da}\frac{\text{d}\unicode[STIX]{x1D712}_{0}}{\text{d}z}-1\right]+\frac{\unicode[STIX]{x1D719}_{0}w_{0}}{Da}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}_{1}}{\unicode[STIX]{x2202}z}=\frac{\unicode[STIX]{x1D719}_{0}}{DaPe}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D712}_{1}}{\unicode[STIX]{x2202}x^{2}}-\unicode[STIX]{x1D712}_{1}. & \displaystyle\end{eqnarray}$$

$K_{0}=\unicode[STIX]{x1D719}_{0}^{n}$

$K_{0}^{\prime }=nK_{0}/\unicode[STIX]{x1D719}_{0}$

We eliminate $\unicode[STIX]{x1D712}_{1}$ using (3.4a ) and $\boldsymbol{v}_{1}$ using (3.4c ). We also use (3.2d ) to simplify the expressions and obtain

(3.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle -{\mathcal{S}}K_{0}\unicode[STIX]{x1D6FB}^{2}{\mathcal{P}}_{1}+nw_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}=-{\mathcal{M}}{\mathcal{P}}_{1}, & \displaystyle\end{eqnarray}$$

(3.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left(-{\mathcal{S}}K_{0}\frac{\unicode[STIX]{x2202}{\mathcal{P}}_{1}}{\unicode[STIX]{x2202}z}+nw_{0}\unicode[STIX]{x1D719}_{1}\right)\left[\frac{-\unicode[STIX]{x1D712}_{0}}{\unicode[STIX]{x1D719}_{0}w_{0}}\right]=-\!\left[\frac{\unicode[STIX]{x1D719}_{0}w_{0}}{Da}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x1D719}_{0}}{DaPe}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}+1\right]\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t}-{\mathcal{P}}_{1}\right).\qquad & \displaystyle\end{eqnarray}$$

$O({\mathcal{M}})$

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

(3.6) $$\begin{eqnarray}\left[\frac{1}{Da}\unicode[STIX]{x2202}_{tz}-\frac{1}{DaPe}\unicode[STIX]{x2202}_{txx}+\unicode[STIX]{x2202}_{t}-n\right]\unicode[STIX]{x1D6FB}^{2}{\mathcal{P}}_{1}=\frac{n}{{\mathcal{S}}}\left[\left(\frac{1}{Da}-{\mathcal{S}}\right)\unicode[STIX]{x2202}_{z}-\frac{1}{DaPe}\unicode[STIX]{x2202}_{xx}+1\right]\unicode[STIX]{x2202}_{z}{\mathcal{P}}_{1}.\end{eqnarray}$$

For brevity in this equation, subscripts are used to denote partial derivatives.

We seek normal-mode solutions ${\mathcal{P}}_{1}\propto \exp (\unicode[STIX]{x1D70E}t+\text{i}kx+mz)$ of this linear equation, where $\unicode[STIX]{x1D70E}$ is the growth rate and $k$ is a horizontal wavenumber. Thus, we obtain the characteristic polynomial (dispersion relationship),

(3.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70E}}{Da}m^{3}+\left(\unicode[STIX]{x1D70E}{\mathcal{K}}-\frac{n}{Da{\mathcal{S}}}\right)m^{2}-\left(\frac{n{\mathcal{K}}}{{\mathcal{S}}}+\frac{\unicode[STIX]{x1D70E}}{Da}k^{2}\right)m+\left(n-\unicode[STIX]{x1D70E}{\mathcal{K}}\right)k^{2}=0,\end{eqnarray}$$

where ${\mathcal{K}}=1+k^{2}/DaPe$ . Equation (3.7) has three roots $m_{j}$ $(j=1,2,3)$ , and hence the compaction pressure perturbation will be given by

(3.8) $$\begin{eqnarray}{\mathcal{P}}_{1}=\mathop{\sum }_{j=1}^{3}A_{j}\exp (\unicode[STIX]{x1D70E}t+\text{i}kx+m_{j}z).\end{eqnarray}$$

The three unknown prefactors $A_{j}$ are determined by the boundary conditions.

3.3 Boundary conditions on the perturbation

We previously eliminated $\unicode[STIX]{x1D712}_{1}$ and $\unicode[STIX]{x1D719}_{1}$ in favour of the compaction pressure ${\mathcal{P}}_{1}$ . The corresponding boundary conditions on ${\mathcal{P}}_{1}$ , derived from (2.14a-c ) and (3.4a ), are

(3.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{P}}_{1}=0\quad \text{at }z=0, & \displaystyle\end{eqnarray}$$

(3.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{P}}_{1}}{\unicode[STIX]{x2202}z}=0\quad \text{at }z=0. & \displaystyle\end{eqnarray}$$

The upper boundary condition (2.14d ) is

(3.9c ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}{\mathcal{P}}_{1}}{\unicode[STIX]{x2202}z}=0\quad \text{at }z=1.\end{eqnarray}$$

The boundary conditions can be expressed in matrix form in terms of the coefficients of the normal-mode expansion (3.8) as

(3.10) $$\begin{eqnarray}\left(\begin{array}{@{}ccc@{}}1 & 1 & 1\\ m_{1} & m_{2} & m_{3}\\ m_{1}\text{e}^{m_{1}} & m_{2}\text{e}^{m_{2}} & m_{3}\text{e}^{m_{3}}\end{array}\right)\left(\begin{array}{@{}c@{}}A_{1}\\ A_{2}\\ A_{3}\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\\ 0\end{array}\right).\end{eqnarray}$$

A necessary (but not sufficient) condition for a non-trivial solution $A_{j}$ to exist is that the boundary-condition matrix $\unicode[STIX]{x1D648}$ has zero determinant.

3.4 Analysis of the dispersion relationship

We analyse the characteristic polynomial (3.7) for the case of real growth rate $\unicode[STIX]{x1D70E}$ (that is, we look for channel modes rather than compaction–dissolution waves, as discussed in § 1). The characteristic polynomial, a cubic, has three roots $m_{j}$ $(j=1,2,3)$ . The character of these roots is controlled by the cubic discriminant. If the discriminant is strictly positive, the roots are distinct and real. If the discriminant is zero, the roots are real but at least one root is repeated (degenerate). If the discriminant is strictly negative, then there is one real root ( $m_{1}$ , say) and a pair of complex conjugate roots ( $m_{2},m_{3}$ ).

For the case of real and distinct roots, the columns of $\unicode[STIX]{x1D648}$ are linearly independent, the determinant of $\unicode[STIX]{x1D648}$ is non-zero and the only solution has $A_{j}=0$ . When the roots are real but degenerate, $\det \unicode[STIX]{x1D648}=0$ , but there is no set of coefficients $A_{j}$ that can satisfy the boundary condition at $z=1$  (3.9c ). Hence, there are physically meaningful roots only when the cubic discriminant of (3.7) is strictly negative.

In this latter case, with one real root and two complex conjugate roots, $\det \unicode[STIX]{x1D648}$ is purely imaginary. A proof of this follows. Consider a $2\times 2$ matrix whose columns are complex conjugate, say $\unicode[STIX]{x1D654}=(\boldsymbol{X},\boldsymbol{X}^{\ast })$ , where $\boldsymbol{X}=[X_{1},X_{2}]^{\text{T}}$ . Then, $\det \unicode[STIX]{x1D654}=X_{1}X_{2}^{\ast }-X_{2}X_{1}^{\ast }$ , so $\det \unicode[STIX]{x1D654}+\det \unicode[STIX]{x1D654}^{\ast }=0$ , i.e.  $\det Y$ is pure imaginary. The boundary-condition matrix $\unicode[STIX]{x1D648}$ is $3\times 3$ , but $\det \unicode[STIX]{x1D648}$ can be written as the sum of purely imaginary determinants of $2\times 2$ submatrices, multiplied by purely real numbers; hence, $\det \unicode[STIX]{x1D648}$ is pure imaginary.

With $m_{1}$ real, and $m_{2}$ and $m_{3}=m_{2}^{\ast }$ complex, there are eigenvalues of $\unicode[STIX]{x1D70E}$ for which the imaginary part of $\det \unicode[STIX]{x1D648}$ vanishes. At these eigenvalues, $\det \unicode[STIX]{x1D648}=0$ , and there exists an eigenvector $A_{j}$ such that the boundary conditions are satisfied. We find these eigenvalues/vectors by numerically solving the coupled problem of the cubic polynomial (3.7) and $\det \unicode[STIX]{x1D648}=0$ .

3.5 Physical discussion of instability mechanism (part I: growth rate)

Figure 2(a) shows an example of the dispersion relationship $\unicode[STIX]{x1D70E}(k)$ . The curves are a series of valid solutions. The solutions on the uppermost dispersion curve have the largest growth rate $\unicode[STIX]{x1D70E}$ and are monotonic in $z$ . Curves below this fundamental mode are higher order, with increasing numbers of turning points in $z$ as $\unicode[STIX]{x1D70E}$ decreases at fixed $k$ . In this example, the instability is only present at $k\gtrsim 1$ , which roughly translates to channels that are narrower than the domain height. Hence, we expect that the lateral wavelength will always be smaller than the domain height.

Figure 2.

(a) Example dispersion relationship calculated at $Da=10^{2}$ , $Pe=10^{2}$ , ${\mathcal{S}}=1$ , $n=3$ . (b) Perturbation corresponding to the most unstable wavenumber (indicated by a diamond symbol in a). The background colour scale shows the porosity perturbation $\unicode[STIX]{x1D719}_{1}$ (normalized to have a maximum value of 1). The black curves are contours of liquid undersaturation $\unicode[STIX]{x1D712}_{1}$ , which is positively correlated with $\unicode[STIX]{x1D719}_{1}$ (solid $=$ positive, dashed $=$ negative). The magenta arrows show the perturbation liquid velocity $\boldsymbol{v}_{1}$ . One can note the flow into the proto-channels (regions of elevated porosity $\unicode[STIX]{x1D719}_{1}>0$ ). The compaction pressure ${\mathcal{P}}_{1}$ (not shown) is anti-correlated with $\unicode[STIX]{x1D719}_{1}$ , consistent with flow direction from high to low pressure.

We now explain the physical mechanism that gives rise to the instability. Figure 2(b) shows an example of the structure of the fastest-growing perturbation (most unstable mode). Regions of positive porosity perturbation $\unicode[STIX]{x1D719}_{1}$ (which we call proto-channels) create a positive perturbation of the vertical flux, according to (3.4c ). For didactic purposes, we consider the case of no compaction pressure (which is directly applicable to a rigid porous medium). Then,

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{0}w_{1}+\unicode[STIX]{x1D719}_{1}w_{0}=nw_{0}\unicode[STIX]{x1D719}_{1}.\end{eqnarray}$$

It should be noted that the positive vertical flux perturbation only occurs because the permeability increases with porosity ( $n>0$ ); this is a crucial aspect of the instability.

Positive vertical advection against the background equilibrium concentration gradient leads to positive liquid undersaturation $\unicode[STIX]{x1D712}_{1}$ , according to (3.4d ). In more physical terms, the enhanced vertical flux advects corrosive liquid from below. Thus, the equilibrium concentration gradient is the other crucial aspect of the instability, alongside the porosity-dependent permeability. For didactic purposes, we consider the case of very fast reaction ( $Da\gg 1$ ), in which the leading-order balance in (3.4d ) gives

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{1}=\unicode[STIX]{x1D719}_{0}w_{1}+\unicode[STIX]{x1D719}_{1}w_{0}=nw_{0}\unicode[STIX]{x1D719}_{1}.\end{eqnarray}$$

Positive liquid undersaturation in turn causes reactive melting and hence increasing porosity by (3.4a ), so the proto-channel emerges. Again, neglecting compaction pressure, replacing $\unicode[STIX]{x2202}_{t}\rightarrow \unicode[STIX]{x1D70E}$ and substituting (3.12), we find

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D70E}\unicode[STIX]{x1D719}_{1}=nw_{0}\unicode[STIX]{x1D719}_{1}\quad \Rightarrow \quad \unicode[STIX]{x1D70E}=n,\end{eqnarray}$$

where we used $w_{0}=1$ . It should be noted that the maximum growth rate in figure 2(a) is approximately $n=3$ . Recalling the non-dimensionalization of time in (2.11), we see that the time scale for channel growth is the time scale for reactive melting ( $\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD}w_{0}$ ) multiplied by the sensitivity of melt flux to porosity ( $n$ ).

Further consideration of (3.4d ) reveals two stabilizing mechanisms. The instability is weakened by diffusion, especially at high wavenumber, since diffusion acts to smooth out lateral gradients in the undersaturation. It is also weakened by advection of the liquid undersaturation, because the undersaturation in the proto-channel increases with height ( $\unicode[STIX]{x2202}\unicode[STIX]{x1D712}_{1}/\unicode[STIX]{x2202}z>0$ ). The subsequent analysis shows that this latter mechanism is also more important at large wavenumber, so both advection and diffusion of liquid undersaturation play a role in wavelength selection (see §§ 4.2 and 4.3 respectively). Indeed, figure 2(a) shows that the growth rate decreases at large $k$ .

Finally, we consider the effect of compaction, which is a further stabilizing mechanism at both large and small wavenumbers (Aharonov et al. Reference Aharonov, Whitehead, Kelemen and Spiegelman1995) (see §§ 4.1, 4.4 and appendix A). The instability only occurs if the matrix stiffness exceeds some critical value (see §§ 4.5 and 4.6). To leading order ( ${\mathcal{M}}\ll 1$ ), if we consider (3.4b ) governing liquid mass conservation, then

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{0}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{1}=-w_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}\quad \Rightarrow \quad K_{0}\unicode[STIX]{x1D6FB}^{2}{\mathcal{P}}_{1}=\frac{nw_{0}}{{\mathcal{S}}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z},\end{eqnarray}$$

where we substitute in (3.4c ) to achieve the last expression (cf. (3.5a )). Proto-channels are regions of increasing porosity perturbation ( $\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}/\unicode[STIX]{x2202}z>0$ ). Thus, by liquid mass conservation, they are regions of convergence of the perturbation velocity $\boldsymbol{v}_{1}$ . Therefore, proto-channels are regions of negative compaction pressure perturbation, which reduces the porosity perturbation, according to the equation of solid mass conservation (3.4a ). Again, this stabilizing mechanism is wavelength-dependent through the Laplacian in (3.14). It should be noted further that the perturbation to the compaction pressure decreases with increasing matrix stiffness ${\mathcal{S}}$ , so we recover the rigid porous medium case as ${\mathcal{S}}\gg 1$ . We return to the physical discussion of the instability in § 4.7 to explain the wavelength selection and the critical matrix stiffness.

4.1 Dependence on wavenumber $k$

Starting at small $k$ , $\unicode[STIX]{x1D716}$ initially decreases with $k$ , reaches some minimum value $\unicode[STIX]{x1D716}^{\ast }$ at $k=k^{\ast }$ (corresponding to the most unstable mode with maximum growth rate $\unicode[STIX]{x1D70E}^{\ast }=n(1-\unicode[STIX]{x1D716}^{\ast })$ ) and then increases as $k\rightarrow \infty$ .

Scaling arguments make these statements more precise. When $k\ll k^{\ast }$ ,

(4.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle a\sim \frac{1}{2{\mathcal{S}}}, & \displaystyle\end{eqnarray}$$

(4.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}\sim \frac{b^{2}+a^{2}}{k^{2}}. & \displaystyle\end{eqnarray}$$

${\mathcal{B}}({\mathcal{S}})=b^{2}+a^{2}$

$a=1/2{\mathcal{S}}$

$b$

$\unicode[STIX]{x1D716}=O(1)$

$\unicode[STIX]{x1D716}\ll 1$

$k\sim {\mathcal{B}}^{1/2}$

Conversely, at large $k$ , $k\gg k^{\ast }$ ,

(4.11a ) $$\begin{eqnarray}\displaystyle & \displaystyle a\sim \frac{1}{2{\mathcal{S}}}+\frac{k^{2}}{2Da}, & \displaystyle\end{eqnarray}$$

(4.11b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}\sim k^{2}\left(\frac{1}{(2Da)^{2}}+\frac{1}{DaPe}\right). & \displaystyle\end{eqnarray}$$

$Da\gg Pe$

$\unicode[STIX]{x1D716}\sim k^{2}/DaPe$

$k\sim (DaPe)^{1/2}$

$Pe\gg Da$

$\unicode[STIX]{x1D716}\sim k^{2}/(2Da)^{2}$

$k\sim 2Da$

It is also possible to determine strict minimum and maximum wavenumbers for instability, although this is more technical so we leave the details for appendix A. In summary, we find

(4.12a,b ) $$\begin{eqnarray}k_{min}\sim \frac{1.5171}{{\mathcal{S}}^{1/2}}\quad ({\mathcal{S}}\gg 1),\quad k_{min}\sim \frac{1}{{\mathcal{S}}}\quad ({\mathcal{S}}\ll 1),\end{eqnarray}$$

(4.13) $$\begin{eqnarray}k_{max}\sim {\mathcal{S}}DaPe.\end{eqnarray}$$

The dependence on the matrix stiffness ${\mathcal{S}}$ means that compaction stabilizes the system at both large and small wavenumbers (Aharonov et al. Reference Aharonov, Whitehead, Kelemen and Spiegelman1995). Indeed, in a rigid medium ( ${\mathcal{S}}\gg 1$ ), there is no minimum or maximum wavenumber.

4.2 Advection-controlled instability, $Pe\gg Da\gg 1$

We first consider the case of negligible diffusion. In this case, it is natural to introduce a change of variables: $\tilde{k}=kDa^{-1/2}$ , $\tilde{\unicode[STIX]{x1D716}}=\unicode[STIX]{x1D716}Da$ . Then, to leading order,

(4.14a ) $$\begin{eqnarray}\displaystyle & \displaystyle a\sim \frac{1}{2{\mathcal{S}}}+\frac{\tilde{k}^{2}}{2}, & \displaystyle\end{eqnarray}$$

(4.14b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D716}}\sim \frac{a^{2}+b^{2}}{\tilde{k}^{2}}. & \displaystyle\end{eqnarray}$$

$b$

$a$

$\tilde{\unicode[STIX]{x1D716}}$

$(\tilde{k},{\mathcal{S}})$

Figure 4.

Numerical calculations of the full cubic dispersion relation (3.7) (solid black curves) compared with power-law scalings (dashed red lines) for the maximum growth rate and the corresponding wavenumber as a function of $Da$ , at fixed $Pe=10^{12}$ (a,b) and $Pe=10^{2}$ (c,d). For all calculations, ${\mathcal{S}}=1$ .

We find the maximum growth rate by differentiating (4.14b ) and seeking the (unique) turning point, which satisfies

(4.15) $$\begin{eqnarray}b^{2}+(\tilde{k}^{2}-1)b\cos (b)\sin (b)-\tilde{k}^{2}\sin ^{2}(b)=0.\end{eqnarray}$$

We solve numerically to obtain the solution $\tilde{k}^{\ast }=\tilde{k}^{\ast }({\mathcal{S}})$ . The corresponding growth rate is $\tilde{\unicode[STIX]{x1D716}}^{\ast }({\mathcal{S}})$ . In summary, the most unstable wavelength $k^{\ast }\sim \tilde{k}^{\ast }Da^{1/2}$ (consistent with the numerical results of Aharonov et al. Reference Aharonov, Whitehead, Kelemen and Spiegelman1995) and the corresponding growth rate $\unicode[STIX]{x1D70E}^{\ast }\sim n[1-\tilde{\unicode[STIX]{x1D716}}^{\ast }Da^{-1}]$ . These scaling results are shown in figure 4(a,b). The dependence on compaction through the matrix stiffness $({\mathcal{S}})$ is shown in figure 5(a,b). The wavenumber is controlled by advection of liquid undersaturation (see § 4.7).

4.3 Diffusion-controlled instability, $Da\gg Pe$

The other limit occurs when diffusion is significant. For this case, it is natural to introduce a different change of variables: $\hat{k}=k(DaPe)^{-1/4}$ , $\hat{\unicode[STIX]{x1D716}}=\unicode[STIX]{x1D716}(DaPe)^{1/2}$ . Then,

(4.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle a\sim \frac{1}{2{\mathcal{S}}}, & \displaystyle\end{eqnarray}$$

(4.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D716}}\sim \frac{{\mathcal{B}}+\hat{k}^{4}}{\hat{k}^{2}}, & \displaystyle\end{eqnarray}$$

${\mathcal{B}}({\mathcal{S}})$

Figure 5.

Dependence on matrix compaction (stiffness ${\mathcal{S}}$ ) in the two regimes $Pe\gg Da$ (a,b) and $Da\gg Pe$ (c,d). Panels (a,c) show the maximum growth rate and (b,d) show the corresponding wavenumber. The solid black curves are numerical calculations of the full cubic dispersion relation (3.7). The dashed red curves (which are almost indistinguishable) are asymptotic results in the limit of large $Da$ . The blue dot-dashed lines are asymptotic results in the limit of small ${\mathcal{S}}$ .

This dispersion relation is simple enough to analyse by hand. The minimum of $\hat{\unicode[STIX]{x1D716}}^{\ast }=2{\mathcal{B}}^{1/2}$ and occurs when $\hat{k}^{\ast }={\mathcal{B}}^{1/4}$ . Thus, the maximum growth rate $\unicode[STIX]{x1D70E}^{\ast }$ that occurs at wavenumber $k^{\ast }$ satisfies

(4.17a ) $$\begin{eqnarray}\displaystyle & \displaystyle k^{\ast }\sim (PeDa{\mathcal{B}})^{1/4}, & \displaystyle\end{eqnarray}$$

(4.17b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}^{\ast }\sim n\left[1-\frac{2\sqrt{{\mathcal{B}}}}{\sqrt{DaPe}}\right]. & \displaystyle\end{eqnarray}$$

$k^{\ast }\sim Pe^{1/4}$

$Da$

${\mathcal{S}}$

${\mathcal{B}}({\mathcal{S}})$

${\mathcal{S}}$

4.4 Effect of compaction (dependence on ${\mathcal{S}}$ )

Asymptotic estimates of the dependence on ${\mathcal{S}}$ are obtained by analysing the roots of (4.5): $\tan b+b/a=0$ . The first non-trivial root $b$ of this equation occurs for $b\in (\unicode[STIX]{x03C0}/2,\unicode[STIX]{x03C0})$ . At small $a$ (large ${\mathcal{S}}$ ), the root $b\rightarrow \unicode[STIX]{x03C0}/2^{+}$ . At large $a$ (small ${\mathcal{S}}$ ), the root $b\rightarrow \unicode[STIX]{x03C0}^{-}$ .

Next, we determine the maximum growth rate for large and small ${\mathcal{S}}$ . First, we consider the case of advection-controlled growth ( $Pe\gg Da\gg 1$ ). At large ${\mathcal{S}}\gg 1$ , $a\sim \tilde{k}^{2}/2$ independent of ${\mathcal{S}}$ . Thus, $a,b,\tilde{k}^{\ast },\tilde{\unicode[STIX]{x1D716}}^{\ast }$ approach some limit that is independent of ${\mathcal{S}}$ . By solving (4.15) numerically, we find that

(4.18a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{k}^{\ast }\rightarrow 1.898, & \displaystyle\end{eqnarray}$$

(4.18b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D716}}^{\ast }\rightarrow 2.302. & \displaystyle\end{eqnarray}$$

${\mathcal{S}}\ll 1$

(4.19a ) $$\begin{eqnarray}\displaystyle & \displaystyle a^{\ast }\sim {\mathcal{S}}^{-1}, & \displaystyle\end{eqnarray}$$

(4.19b ) $$\begin{eqnarray}\displaystyle & \displaystyle b^{\ast }\sim \unicode[STIX]{x03C0}(1-{\mathcal{S}}), & \displaystyle\end{eqnarray}$$

(4.19c ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{k}^{\ast }\sim {\mathcal{S}}^{-1/2}, & \displaystyle\end{eqnarray}$$

(4.19d ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D716}}^{\ast }\sim {\mathcal{S}}^{-1}. & \displaystyle\end{eqnarray}$$

Second, we consider the case of diffusion-controlled growth ( $Da\gg Pe$ ). As before, the growth rate approaches a constant as ${\mathcal{S}}$ increases, namely

(4.20a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{k}^{\ast }\rightarrow (\unicode[STIX]{x03C0}/2)^{1/2}, & \displaystyle\end{eqnarray}$$

(4.20b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D716}}^{\ast }\rightarrow \unicode[STIX]{x03C0}. & \displaystyle\end{eqnarray}$$

${\mathcal{S}}\ll 1$

$a\sim {\mathcal{S}}^{-1}$

(4.21a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{k}^{\ast }\sim (2{\mathcal{S}})^{-1/2}, & \displaystyle\end{eqnarray}$$

(4.21b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D716}}^{\ast }\sim {\mathcal{S}}^{-1}. & \displaystyle\end{eqnarray}$$

${\mathcal{S}}\rightarrow \infty$

${\mathcal{S}}\gtrsim 1$

The scalings for wavenumber and growth rate are the same in terms of the power-law dependence on ${\mathcal{S}}$ . In either case, a compactible medium is less unstable than a rigid medium. That is, compaction stabilizes the system. We can interpret (4.19d ) and (4.21b ) in terms of a critical stiffness such that the instability occurs when ${\mathcal{S}}\geqslant {\mathcal{S}}_{crit}$ , where

(4.22a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}_{crit}\propto \frac{1}{Da}\quad (Pe\gg Da), & \displaystyle\end{eqnarray}$$

(4.22b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}_{crit}\propto \frac{1}{\sqrt{DaPe}}\quad (Da\gg Pe). & \displaystyle\end{eqnarray}$$

We can also estimate the aspect ratio ${\mathcal{A}}$ of the instability for the case ${\mathcal{S}}\ll 1$ by noting that $a\sim {\mathcal{S}}^{-1}$ . The ratio of horizontal to vertical length scale is approximately ${\mathcal{A}}\sim a/k^{\ast }$ . Substituting in the wavenumber scalings, we find

(4.23a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{A}}\propto ({\mathcal{S}}Da)^{-1/2}\quad (Pe\gg Da), & \displaystyle\end{eqnarray}$$

(4.23b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{A}}\propto ({\mathcal{S}}^{2}DaPe)^{-1/4}\quad (Da\gg Pe). & \displaystyle\end{eqnarray}$$

${\mathcal{S}}_{crit}$

${\mathcal{S}}_{crit}$

${\mathcal{S}}\gtrsim 1$

Figure 6.

Behaviour as ${\mathcal{S}}\rightarrow {\mathcal{S}}_{crit}^{+}$ (example with $Da=Pe=10$ ). (a) Series of closed loops in $(k,\unicode[STIX]{x1D70E})$ -space as ${\mathcal{S}}$ decreases towards the critical value (purple to light blue; black dot-dashed loop corresponds to the final iteration; method described in appendix B). (b) The porosity perturbation $\unicode[STIX]{x1D719}_{1}$ corresponding to the most unstable mode indicated by the diamond symbol in (a).

4.5 Numerical investigation of the critical stiffness

We next test these asymptotic predictions of a critical stiffness by numerically calculating the dispersion relationship at successive values of ${\mathcal{S}}\rightarrow {\mathcal{S}}_{crit}^{+}$ . Figure 6 shows that (a) the dispersion relationship forms closed loops whose size approaches zero, (b) the perturbation is localized in an $O({\mathcal{S}})$ boundary layer near the upper boundary. The latter observation is consistent with the asymptotic result that the vertical length scale $a^{-1}\sim {\mathcal{S}}$ . We estimate the critical value ${\mathcal{S}}_{crit}$ using the method described in appendix B, and map out the dependence on Damköhler number and Péclet number.

Figure 7.

(a) The critical stiffness ${\mathcal{S}}_{crit}$ as a function of $Da$ at $Pe=10,10^{2},10^{3},10^{4}$ (from dark purple to light blue). We also show power-law scalings $Da^{-1}$ (dash-dotted) and $Da^{-1/2}$ (dotted). (b) The corresponding wavenumber, which obeys the scaling (4.25). (c) The corresponding growth rate, which appears to approach a constant at large $Da$ .

Figure 7(a) shows the dependence of ${\mathcal{S}}_{crit}$ on $Da$ at $Pe=10,10^{2},10^{3},10^{4}$ . The calculations with high $Pe$ support the prediction of (4.22a ) that ${\mathcal{S}}_{crit}\propto Da^{-1}$ when $Pe\gg Da$ . The calculations with lower $Pe$ support the prediction of (4.22b ) that ${\mathcal{S}}_{crit}\propto Da^{-1/2}$ when $Da\gg Pe$ ; they are also consistent with the predicted $Pe^{-1/2}$ dependence. By estimating the prefactors numerically, we obtain the following scalings:

(4.24a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}_{crit}\sim \frac{1}{Da}\quad (Pe\gg Da), & \displaystyle\end{eqnarray}$$

(4.24b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}_{crit}\sim \frac{2}{\sqrt{DaPe}}\quad (Da\gg Pe). & \displaystyle\end{eqnarray}$$

${\mathcal{S}}=O(1)$

Figure 7(b) shows that, across the range of parameters considered, the wavenumber at ${\mathcal{S}}_{crit}$ obeys the scaling

(4.25) $$\begin{eqnarray}k({\mathcal{S}}_{crit})\sim (DaPe)^{1/2}.\end{eqnarray}$$

This is the same as the scaling for the large-wavenumber cutoff when $Da\gg Pe$ , identified in § 4.1. However, we see in § 4.6 that a different scaling eventually emerges at higher $Pe$ .

Finally, figure 7(c) shows that the growth rate at ${\mathcal{S}}_{crit}$ appears to approach a (non-zero) constant at large $Da$ , which might be independent of $Pe$ , a hypothesis we confirm in § 4.6.

4.6 Analysis of behaviour near the critical stiffness

We now analyse the structure of the bifurcation at ${\mathcal{S}}_{crit}$ . Our goal is to complement the numerical results obtained previously by mapping out the bifurcation structure and obtaining asymptotic results at very large $Da$ and $Pe$ , regimes that were hard to achieve numerically.

We proceed by rescaling equations (4.8a ) and (4.8b ). As we have seen previously, there are two distinguished limits depending on the relative magnitude of $Da$ to $Pe$ .

4.6.1 Case $Pe\gg Da\gg 1$

We first consider the case in which chemical diffusion is negligible. We use the rescaling

(4.26a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{x}=k^{2}Da^{-2}, & \displaystyle\end{eqnarray}$$

(4.26b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{a}=aDa^{-1}, & \displaystyle\end{eqnarray}$$

(4.26c ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{{\mathcal{S}}}={\mathcal{S}}Da, & \displaystyle\end{eqnarray}$$

(4.26d ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D70E}}=\unicode[STIX]{x1D70E}n^{-1} & \displaystyle\end{eqnarray}$$

(extending the scaling first introduced in § 4.2). Then, we note that ${\mathcal{K}}=1+k^{2}/DaPe=1+\tilde{x}(Da/Pe)$ , so ${\mathcal{K}}\sim 1$ . We also note that $\unicode[STIX]{x03C0}/2<b<\unicode[STIX]{x03C0}$ , so $b^{2}\ll a^{2}=O(Da^{2})$ . Thus, to leading order, (4.8a ) and (4.8b ) become respectively

(4.27a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{a}=\frac{(\tilde{{\mathcal{S}}}^{-1}+\tilde{\unicode[STIX]{x1D70E}}\tilde{x})}{2(\tilde{\unicode[STIX]{x1D70E}}-\tilde{{\mathcal{S}}}^{-1})}, & \displaystyle\end{eqnarray}$$

(4.27b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{a}^{2}=\frac{(1-\tilde{\unicode[STIX]{x1D70E}})\tilde{x}}{(\tilde{\unicode[STIX]{x1D70E}}-\tilde{{\mathcal{S}}}^{-1})}. & \displaystyle\end{eqnarray}$$

We can eliminate $\tilde{a}$ and rearrange into a quadratic for $\tilde{\unicode[STIX]{x1D70E}}$ ,

(4.28) $$\begin{eqnarray}\tilde{{\mathcal{S}}}^{2}\tilde{x}(4+\tilde{x})\tilde{\unicode[STIX]{x1D70E}}^{2}-2\tilde{{\mathcal{S}}}\tilde{x}(1+2\tilde{{\mathcal{S}}})\tilde{\unicode[STIX]{x1D70E}}+(1+4\tilde{x}\tilde{{\mathcal{S}}})=0.\end{eqnarray}$$

There are repeated roots when the discriminant of the quadratic is zero, corresponding to the left- and right-hand limits of the loops shown in figure 6(a). The discriminant is

(4.29) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D6E5}}=-16\tilde{x}\tilde{{\mathcal{S}}}^{2}[\tilde{x}^{2}\tilde{{\mathcal{S}}}+\tilde{x}(3\tilde{{\mathcal{S}}}-\tilde{{\mathcal{S}}}^{2})+1].\end{eqnarray}$$

Given that $\tilde{{\mathcal{S}}}>0$ and $\tilde{x}\propto k^{2}>0$ , the roots $\tilde{\unicode[STIX]{x1D6E5}}=0$ must satisfy the quadratic equation

(4.30) $$\begin{eqnarray}\tilde{x}^{2}\tilde{{\mathcal{S}}}+\tilde{x}(3\tilde{{\mathcal{S}}}-\tilde{{\mathcal{S}}}^{2})+1=0.\end{eqnarray}$$

The bifurcation (at the critical matrix stiffness) occurs when the discriminant of this quadratic is zero, when

(4.31) $$\begin{eqnarray}0=(3\tilde{{\mathcal{S}}}-\tilde{{\mathcal{S}}}^{2})^{2}-4\tilde{{\mathcal{S}}}=\tilde{{\mathcal{S}}}(\tilde{{\mathcal{S}}}-1)^{2}(\tilde{{\mathcal{S}}}-4).\end{eqnarray}$$

The root $\tilde{{\mathcal{S}}}=0$ is excluded because $\tilde{{\mathcal{S}}}>0$ . The roots $\tilde{{\mathcal{S}}}=1$ are excluded because they correspond to repeated roots $\tilde{x}=-1$ in (4.30). The only physically meaningful root is $\tilde{{\mathcal{S}}}=4$ , which corresponds to repeated roots $\tilde{x}=1/2$ in (4.30). We substitute back into (4.28) and find that the corresponding $\tilde{\unicode[STIX]{x1D70E}}=1/2$ .

This gives us the critical matrix stiffness and the corresponding properties of the solution (horizontal and vertical wavenumbers and growth rate). In summary, we find that

(4.32a-d ) $$\begin{eqnarray}{\mathcal{S}}_{crit}=4Da^{-1},\quad k_{crit}\sim Da/\sqrt{2},\quad a_{crit}\sim Da,\quad \unicode[STIX]{x1D70E}_{crit}\sim n/2\quad (Pe\gg Da\gg 1).\end{eqnarray}$$

In the numerical results (§ 4.5), we found the same ${\mathcal{S}}_{crit}\propto Da^{-1}$ scaling, albeit with a different prefactor. However, we did not observe the $k_{crit}\propto Da$ scaling (independent of $Pe$ ), which indicates that our numerical calculations were not performed at sufficiently high $Pe$ to observe the asymptotic regime. Our analysis in this section allows access to that regime. Furthermore, observations such as those in figure 6 of loops emerging at finite (non-zero) values of the growth rate $\unicode[STIX]{x1D70E}$ emerge as generic features of the bifurcation.

4.6.2 Case $Da\gg Pe$

We second consider the opposite case in which advection of the liquid undersaturation is negligible relative to diffusion. We apply the same type of methodology as before. We use the rescaling

(4.33a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{x}=k^{2}(DaPe)^{-1}, & \displaystyle\end{eqnarray}$$

(4.33b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{a}=a(DaPe)^{-1/2}, & \displaystyle\end{eqnarray}$$

(4.33c ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{{\mathcal{S}}}={\mathcal{S}}(DaPe)^{1/2}, & \displaystyle\end{eqnarray}$$

(4.33d ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D70E}}=\unicode[STIX]{x1D70E}n^{-1} & \displaystyle\end{eqnarray}$$

(the scaling extends that introduced in § 4.3). We note that ${\mathcal{K}}=1+k^{2}/DaPe=1+\hat{x}$ . Again, $b^{2}\ll a^{2}=O(DaPe)$ and $1/Da{\mathcal{S}}\sim (Pe/Da)^{1/2}\ll 1$ . Thus, to leading order, (4.8a ) and (4.8b ) become respectively

(4.34a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{a}=\frac{1}{2\hat{\unicode[STIX]{x1D70E}}\hat{{\mathcal{S}}}}, & \displaystyle\end{eqnarray}$$

(4.34b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{a}^{2}=\frac{(1-\hat{\unicode[STIX]{x1D70E}}(1+\hat{x}))\hat{x}}{\hat{\unicode[STIX]{x1D70E}}(1+\hat{x})}. & \displaystyle\end{eqnarray}$$

The corresponding quadratic for $\hat{\unicode[STIX]{x1D70E}}$ is,

(4.35) $$\begin{eqnarray}4\hat{{\mathcal{S}}}^{2}\hat{x}(1+\hat{x})\hat{\unicode[STIX]{x1D70E}}^{2}-4\hat{{\mathcal{S}}}^{2}\hat{x}\hat{\unicode[STIX]{x1D70E}}+(1+\hat{x})=0,\end{eqnarray}$$

whose repeated roots are zeros of the discriminant

(4.36) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D6E5}}=-16\hat{x}\hat{{\mathcal{S}}}^{2}[\hat{x}^{2}+\hat{x}(2-\hat{{\mathcal{S}}}^{2})+1].\end{eqnarray}$$

Again, we have $\hat{{\mathcal{S}}}>0$ and $\hat{x}\propto k^{2}>0$ , so roots $\hat{\unicode[STIX]{x1D6E5}}=0$ satisfy

(4.37) $$\begin{eqnarray}\hat{x}^{2}+\hat{x}(2-\hat{{\mathcal{S}}}^{2})+1=0.\end{eqnarray}$$

The bifurcation occurs when the discriminant of this quadratic is zero, when

(4.38) $$\begin{eqnarray}0=(2-\hat{{\mathcal{S}}}^{2})^{2}-4=\hat{{\mathcal{S}}}^{2}(\hat{{\mathcal{S}}}+2)(\hat{{\mathcal{S}}}-2).\end{eqnarray}$$

The only physically meaningful root is $\hat{{\mathcal{S}}}=2$ , which corresponds to repeated roots $\hat{x}=1$ in (4.37). We substitute back into (4.35) and find that the corresponding $\hat{\unicode[STIX]{x1D70E}}=1/4$ . In conclusion,

(4.39a-c ) $$\begin{eqnarray}{\mathcal{S}}_{crit}=2(DaPe)^{-1/2},\quad k_{crit},a_{crit}\sim (DaPe)^{1/2},\quad \unicode[STIX]{x1D70E}_{crit}\sim n/4\quad (Da\gg Pe).\end{eqnarray}$$

In the numerical results (§ 4.5), we found the same scaling relationships (with the same prefactors), so our numerics were able to access this regime adequately. The analysis in this section additionally obtained $\unicode[STIX]{x1D70E}_{crit}\sim n/4$ . Therefore, in both regimes, we find that the bifurcation at ${\mathcal{S}}_{crit}$ results in an instability with a finite growth rate.

4.7 Physical discussion of instability mechanism (part II: wavelength selection)

In § 3.5, we explained the basic structure of the physical instability mechanism. We found that there is an enhanced vertical flux in proto-channels (regions of positive porosity perturbation) caused by the porosity-dependent permeability. This vertical flux across a background equilibrium concentration gradient dissolves the solid matrix, increasing the porosity and establishing an instability that grows at a rate $\unicode[STIX]{x1D70E}\sim n$ . In this section, we use the insights gained from our asymptotic analysis to explain the physical controls on the vertical and horizontal length scales of the instability, and on the critical matrix stiffness. All of the following estimates are consistent with the results of our asymptotic analysis and numerical calculations.

We derive scalings focusing on the more interesting case of a compacting porous medium ( ${\mathcal{S}}\ll 1$ ). Results for a rigid porous medium (up to an unknown prefactor) can be obtained by substituting ${\mathcal{S}}=1$ into the subsequent scalings, consistent with our numerical results that the rigid-medium limit applies when ${\mathcal{S}}\gtrsim 1$ .

First, we consider the vertical length scale of the instability at fixed horizontal wavenumber. In a compacting porous medium, mass conservation implies that gradients in porosity are sources or sinks of compaction pressure, as expressed in (3.14), which we now rewrite by substituting expressions for the base state variables,

(4.40) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}{\mathcal{P}}_{1}=\frac{n}{{\mathcal{S}}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}.\end{eqnarray}$$

We next substitute in the balance between compaction and porosity change from (3.4a ), namely $\unicode[STIX]{x1D70E}\unicode[STIX]{x1D719}_{1}\sim {\mathcal{P}}_{1}$ , use $\unicode[STIX]{x1D70E}\sim n$ and scale $\unicode[STIX]{x2202}_{z}\sim m$ , neglecting horizontal derivatives at fixed $k$ , to obtain

(4.41) $$\begin{eqnarray}m\sim {\mathcal{S}}^{-1}.\end{eqnarray}$$

Therefore, the vertical structure is controlled by the matrix stiffness. In the rigid-medium case, $m\sim 1$ and the instability extends through the full depth of the melting region.

Second, we consider the horizontal length scale of the most unstable mode $k^{\ast }$ . We combine (3.4a ), (4.40) and (4.41) to obtain the estimate

(4.42) $$\begin{eqnarray}-k^{\ast 2}{\mathcal{P}}_{1}\sim \unicode[STIX]{x1D712}_{1}/{\mathcal{S}}^{2}.\end{eqnarray}$$

Physically, reactive dissolution in the channels requires a convergent flow of liquid into the proto-channels, which must be down a gradient in the compaction pressure. Then, by substituting (3.4a ) and (3.4c ) into the liquid concentration equation (3.4d ), we find that

(4.43) $$\begin{eqnarray}{\mathcal{S}}\frac{\unicode[STIX]{x2202}{\mathcal{P}}_{1}}{\unicode[STIX]{x2202}z}-{\mathcal{P}}_{1}\sim -\frac{1}{Da}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}_{1}}{\unicode[STIX]{x2202}z}+\frac{1}{DaPe}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D712}_{1}}{\unicode[STIX]{x2202}x^{2}}.\end{eqnarray}$$

Both terms on the left-hand side are $O({\mathcal{P}}_{1})$ . When $Pe\gg Da$ , diffusion on the right-hand side is negligible compared with advection of the liquid undersaturation. Then, by substituting (4.42), we find

(4.44) $$\begin{eqnarray}k^{\ast 2}\sim \frac{Da}{{\mathcal{S}}}\quad (Pe\gg Da).\end{eqnarray}$$

Conversely, if $Da\gg Pe$ , diffusion dominates the right-hand side of (4.43) and we find

(4.45) $$\begin{eqnarray}k^{\ast 2}\sim \frac{\sqrt{DaPe}}{{\mathcal{S}}}\quad (Da\gg Pe).\end{eqnarray}$$

Physically, the perturbed compaction-driven advection against the equilibrium concentration gradient is balanced by either advection or diffusion of liquid undersaturation. These results mean that the most unstable horizontal wavelength is proportional to the compaction length (it should be recalled that ${\mathcal{S}}\propto \unicode[STIX]{x1D6FF}^{2}$ ). However, it is much smaller than the compaction length since $Da\gg 1$ , as seen in the 2D numerical calculations of Spiegelman et al. (Reference Spiegelman, Kelemen and Aharonov2001).

Third, if the matrix stiffness ${\mathcal{S}}$ is reduced below some critical value, then the stabilizing influence of compaction is dominant over the destabilizing influence of reactive melting such that the instability is suppressed. We can obtain an estimate of this critical value as follows. At the critical value, (3.4a ) gives that compaction balances reaction, so $-{\mathcal{P}}_{1}\sim \unicode[STIX]{x1D712}_{1}$ . We next use (3.12) to obtain $-{\mathcal{P}}_{1}\sim nw_{0}\unicode[STIX]{x1D719}_{1}$ . We substitute into (3.14) to obtain

(4.46) $$\begin{eqnarray}{\mathcal{S}}_{crit}\sim m/k^{\ast 2}.\end{eqnarray}$$

We then substitute in our estimates of the vertical ( $m$ ) and horizontal ( $k^{\ast }$ ) wavenumbers to obtain

(4.47a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}_{crit}\sim \frac{1}{Da}\quad (Pe\gg Da), & \displaystyle\end{eqnarray}$$

(4.47b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{S}}_{crit}\sim \frac{1}{\sqrt{DaPe}}\quad (Da\gg Pe). & \displaystyle\end{eqnarray}$$

(4.48) $$\begin{eqnarray}k_{min}^{2}\sim m/{\mathcal{S}}\sim 1/{\mathcal{S}}^{2}\quad \Rightarrow \quad k_{min}\sim 1/{\mathcal{S}}.\end{eqnarray}$$

Thus, the maximum wavelength for the instability is proportional to $\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}^{2}/\unicode[STIX]{x1D6FC}$ , the product of the compaction length and the amount of reactive melting over a compaction length. In the rigid-medium limit ( ${\mathcal{S}}\gg 1$ ), the vertical wavelength is the full height of the domain ( $m\sim 1$ ). Thus,

(4.49) $$\begin{eqnarray}k_{min}^{2}\sim m/{\mathcal{S}}\sim 1/{\mathcal{S}}\quad \Rightarrow \quad k_{min}\sim 1/{\mathcal{S}}^{1/2}.\end{eqnarray}$$