Magma ascent mechanisms in the transition regime from solitary porosity waves to diapirism

In partially molten regions inside the Earth, melt buoyancy may trigger upwelling of both solid and fluid phases, i.e., diapirism. If the melt is allowed to move separately with respect to the matrix, melt perturbations may evolve into solitary porosity waves. While diapirs may form on a wide range of scales, porosity waves are restricted to sizes of a few times the compaction length. Thus, the size of a partially molten perturbation in terms of compaction length controls whether material is dominantly transported by porosity waves or by diapirism. We study the transition from diapiric rise to solitary porosity waves by solving the two-phase flow equations of conservation of mass and momentum in 2D with porosity-dependent matrix viscosity. We systematically vary the initial size of a porosity perturbation from 1.8 to 120 times the compaction length.

If the perturbation is of the order of a few compaction lengths, a single solitary wave will emerge, either with a positive or negative vertical matrix flux. If melt is not allowed to move separately to the matrix a diapir will emerge. In between these end members we observe a regime where the partially molten perturbation will split up into numerous solitary waves, whose phase velocity is so low compared to the Stokes velocity that the whole swarm of waves will ascend jointly as a diapir, just slowly elongating due to a higher amplitude main solitary wave.

Only if the melt is not allowed to move separately to the matrix will no solitary waves build up, but as soon as two-phase flow is enabled solitary waves will eventually emerge. The required time to build them up increases nonlinearly with the perturbation radius in terms of compaction length and might be too long to allow for them in nature in many cases.

1 Introduction

In geodynamic settings such as mid-ocean ridges, hotspots, subduction zones, or orogenic belts partial melts are generated within the asthenosphere or lower continental crust and ascend by fluid migration within deforming rocks (e.g., Sparks and Parmentier, 1991; Katz, 2008; Keller et al., 2017; Schmeling et al., 2019). Inherent tectonic or rock heterogeneities in such systems may result in spatially varying melt fractions on length scales varying over several orders of magnitudes. These length scales play an important role in determining whether melt anomalies may rise as porous waves (Jordan et al., 2018) or by other mechanisms such as diapirs (Rabinowicz et al., 1987), focused channel networks (Spiegelman et al., 2001), or dikes (Rivalta et al., 2015). Here we focus on the effect of the length scale on the formation and evolution of buoyancy-driven porous waves or diapirs.

The physics of fluid moving relatively to a viscously deformable porous matrix were first described by McKenzie (1984), and it was later shown by several authors that these equations allow for the emergence of solitary porosity waves (Scott and Stevenson, 1984; Barcilon and Lovera 1989; Wiggins and Spiegelman, 1995). Porosity waves are regions of localized excess fluid that ascend with permanent shape and constant velocity, controlled by compaction and decompaction of the surrounding matrix. They have been extensively studied as mechanisms transporting geochemical signatures or magma through the asthenosphere, lower crust, and middle crust (e.g., Watson and Spiegelman, 1994; McKenzie, 1984; Connolly, 1997; Connolly and Podladchikov, 2013, Jordan et al., 2018, Richard et al., 2012). It has been shown that the dynamics of porous waves strongly depend on the porosity dependence of the matrix rheology (e.g., Connolly and Podladchikov, 1998, 2015; Yarushina et al., 2015; Omlin et al., 2017; Dohmen et al., 2019). However, one open question is how the length scale of solitary porosity waves relates to an arbitrary length scale of a possible porosity anomaly in given geodynamic settings.

The size of a solitary porosity wave is usually of the order of a few compaction lengths (McKenzie, 1984; Scott and Stevenson, 1984; Simpson and Spiegelman, 2011), but this length scale varies over a few orders of magnitude, depending on the shear and volume viscosity of the matrix, fluid viscosity, and permeability (see Eq. 19) with typical values of 100–10 000 m (McKenzie, 1984; Spiegelman, 1993a, b). However, partially molten regions in the lower crust or upper mantle are prone to gravitational instabilities such as Rayleigh–Taylor instabilities or diapirism (e.g., Griffith, 1986; Bittner and Schmeling, 1995; Schmeling et al., 2019). Originating from the Greek “diapeirein”, i.e., “to pierce through”, diapirism describes the “buoyant upwelling of relatively light rock” (Turcotte and Schubert, 1982) through and into a denser overburden. In the general definition, the rheology of the diapir and ambient material is not specified and both can be ductile, as in our case. Buoyancy may be of compositional or phase-related origin, e.g., due to the presence of non-segregating partial melt (Wilson, 1989). In this model we describe a diapir as a partially molten perturbation, whose rising velocity, characterizable by the Stokes velocity, is lower than the corresponding solitary-wave phase velocity.

As characteristic wavelengths of Rayleigh–Taylor instabilities may be similar but also of significantly different order to those of porosity waves, and the Stokes velocity is strongly affected by the spatial expansion, the question arises as to how these two mechanisms interact and what the transition from a porosity wave to a rising partially molten diapir looks like. Scott (1988) already investigated a similar scenario. He calculated porosity waves changing the compaction length by altering the constant shear to volume viscosity ratio. In contrast, we vary the radius of a partially molten perturbation in terms of compaction lengths while keeping the porosity-dependent viscosity law the same. While Scott (1988) was not able to reach the single-phase flow endmember due to his setup, we can reach this endmember with our description and can explore the transition.

In this work we will address the question of length scale of a partially molten region with respect to the length scale of a solitary porosity wave by varying the sizes of initial porosity perturbations. We further focus on the numerical implications on modeling magma transport.

2 Methods

2.1 Governing equations

The formulation of the governing equations for the melt-in-solid two-phase flow dynamics is based on McKenzie (1984), Spiegelman and McKenzie (1987), and Schmeling (2000), assuming an infinite Prandtl number, a low fluid viscosity with respect to the effective matrix viscosity, zero surface tension, and the Boussinesq approximation. In the present formulation, the Boussinesq approximation assumes the same constant density for the solid and fluid except for the buoyancy terms of the momentum equations for the solid and fluid. In the following all variables associated with the pore fluid (melt) have the subscript f and those associated with the solid matrix have the subscript s. The equation for the conservation of the mass of the melt is

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\phi }}{\partial t}}+\mathbf{\nabla }\cdot \left(\mathit{\phi }{\mathit{v}}_{\text{f}}\right)=\mathrm{0},\end{array}$$

and the mass conservation of the solid is

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial \left(\mathrm{1}-\mathit{\phi }\right)}{\partial t}}+\mathbf{\nabla }\cdot \left((\mathrm{1}-\mathit{\phi }){\mathit{v}}_{\text{s}}\right)=\mathrm{0}.\end{array}$$

φ is the volumetric rock porosity (often called melt fraction), v_f and v_s are the fluid and solid velocities, respectively. The momentum equations are given as a generalized Darcy equation for the fluid separation flow

$$\begin{array}{}\text{(3)}& {\mathit{v}}_{\text{f}}-{\mathit{v}}_{\text{s}}=-{\displaystyle \frac{k_{\mathit{\phi }}}{\mathit{\mu }\mathit{\phi }}}\left(\mathbf{\nabla }{P}_{\text{f}}-{\mathit{\rho }}_{\text{f}}\mathit{g}\right),\end{array}$$

where ρ_f is the fluid density and P_f is the fluid pressure (including the lithostatic pressure), whose gradient is driving the fluid segregation by porous flow, μ is the melt dynamic viscosity, and g is the gravitational acceleration. k_φ is the permeability that depends on the rock porosity

$$\begin{array}{}\text{(4)}& k_{\mathit{\phi }}=k_{\mathrm{0}}{\mathit{\phi }}^n,\end{array}$$

with n being the power law exponent constant, usually equal to 2 or 3. This relation is known as the Kozeny–Carman relation (e.g., Costa, 2006). The Stokes equation for the mixture is given as

$$\begin{array}{}\text{(5)}& \stackrel{\mathrm{‾}}{\mathit{\rho }}\mathit{g}-\mathbf{\nabla }{P}_{\text{f}}+\mathbf{\nabla }\cdot \mathit{\tau }=\mathrm{0}.\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ is the density of the melt–solid mixture and τ is the effective viscous stress tensor of the matrix, including both shear and compaction components

$$\begin{array}{}\text{(6)}& \mathit{\tau }=\mathit{\eta }\left({\displaystyle \frac{\partial v_{si}}{\partial x_j}}+{\displaystyle \frac{\partial v_{sj}}{\partial x_i}}\right)+\left(\mathit{\zeta }-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\mathit{\eta }\right){\mathit{\delta }}_{ij}\mathbf{\nabla }\cdot {\mathit{v}}_{\text{s}},\end{array}$$

where ζ is the volume viscosity, and v_si and v_sj are the ith and jth component of v_s. The linearized equation of state for the mixture density is given as

$$\begin{array}{}\text{(7)}& \stackrel{\mathrm{‾}}{\mathit{\rho }}={\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{1}-c_{\text{f}}\mathit{\phi }\right),\end{array}$$

with ρ₀ as the solid density and $c_{\text{f}}=\frac{{\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\text{f}}}{{\mathit{\rho }}_{\mathrm{0}}}$ . The shear and volume viscosity are given by the equations

$$\begin{array}{}\text{(8)}& \mathit{\eta }={\mathit{\eta }}_{\mathrm{0}}\left(\mathrm{1}-\mathit{\phi }\right),\end{array}$$

and

$$\begin{array}{}\text{(9)}& \mathit{\zeta }={\mathit{\eta }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{1}-\mathit{\phi }}{\mathit{\phi }}},\end{array}$$

where η₀ is the constant intrinsic shear viscosity of the matrix.

As in both Eqs. (3) and (5), P_f is the fluid pressure (see McKenzie, 1984, Appendix A), these equations can be merged to eliminate the pressure resulting in

$$\begin{array}{}\text{(10)}& {\mathit{v}}_{\text{f}}-{\mathit{v}}_{\text{s}}=-{\displaystyle \frac{k_{\mathrm{0}}{\mathit{\phi }}^{n-\mathrm{1}}}{\mathit{\mu }}}\left({\mathit{\rho }}_{\mathrm{0}}{c}_{\text{f}}\mathit{g}\left(\mathrm{1}-\mathit{\phi }\right)+\mathbf{\nabla }\cdot \mathit{\tau }\right).\end{array}$$

This equation states that the fluid separation flow (i.e., melt segregation velocity) is driven by the buoyancy of the fluid with respect to the solid and the viscous stress in the matrix including compaction and decompaction.

Following Šrámek et al. (2010), the Stokes equation (Eq. 3) can be rewritten by expressing the matrix velocity, v_s, as the sum of the incompressible flow velocity, v₁, and the irrotational (compaction) flow velocity, v₂, as follows:

$$\begin{array}{}\text{(11)}& {\mathit{v}}_{\text{s}}={\mathit{v}}_{\mathrm{1}}+{\mathit{v}}_{\mathrm{2}}=\left(\begin{array}{c}\frac{\partial \mathit{\psi }}{\partial z}\\ -\frac{\partial \mathit{\psi }}{\partial x}\end{array}\right)+\left(\begin{array}{c}\frac{\partial \mathit{\chi }}{\partial x}\\ \frac{\partial \mathit{\chi }}{\partial z}\end{array}\right),\end{array}$$

with ψ as stream function and χ as the irrotational velocity potential, given as the solution of the Poisson equation

$$\begin{array}{}\text{(12)}& {\mathbf{\nabla }}^{\mathrm{2}}\mathit{\chi }=\mathbf{\nabla }\cdot {\mathit{v}}_{\text{s}}.\end{array}$$

The divergence term  ∇⋅v_s can be derived from Eqs. (1) and (2) to give

$$\begin{array}{}\text{(13)}& \mathbf{\nabla }\cdot {\mathit{v}}_{\text{s}}=-\mathbf{\nabla }\cdot \left[\mathit{\phi }\left({\mathit{v}}_{\text{f}}-{\mathit{v}}_{\text{s}}\right)\right].\end{array}$$

In the small fluid viscosity limit the viscous stresses within the fluid phase are neglected, resulting in a viscous stress tensor in the Stokes equation of the mixture (Eq. 5), in which only the stresses in the solid phase are relevant. This is evident from the definition of the viscous stress tensor, which only contains matrix and not fluid viscosities. Melt viscosities of carbonatitic, basaltic, or silicic wet or dry melts span a range from < 1 Pa s to extreme values up to 10¹⁴ Pa s (see the discussion in Schmeling et al., 2019), while effective viscosities of mafic or silicic partially molten rocks may range between 10¹⁶ and 10²⁰ Pa s, depending on melt fraction, stress, and composition. Thus, in most circumstances the small fluid viscosity limit is justified.

In the limit of this small viscosity assumption, inserting the above solid velocity Eq. (11) into the viscous stress Eq. (6), this in turn into the Stokes equation (Eq. 5), and then taking the curl of the x and z equations, the pressure is eliminated and one gets

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& \left({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right)\left[{\mathit{\eta }}_{\text{s}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial z^{\mathrm{2}}}}\right)\right]\\ & +\mathrm{4}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x\partial z}}\left[{\mathit{\eta }}_{\text{s}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial x\partial z}}\right]=-g{\displaystyle \frac{\partial \mathit{\rho }}{\partial x}}+A\left(\mathit{\chi }\right),\end{array}\end{array}$$

with

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& A\left(\mathit{\chi }\right)=-\mathrm{2}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x\partial z}}\left[{\mathit{\eta }}_{\text{s}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\chi }}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\chi }}{\partial z^{\mathrm{2}}}}\right)\right]\\ & +\mathrm{2}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right)\left[{\mathit{\eta }}_{\text{s}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\chi }}{\partial x\partial z}}\right].\end{array}\end{array}$$

To describe the transition from solitary waves to diapirs it is useful to non-dimensionalize the equations. As scaling quantities, we use the radius r of the anomaly, the reference viscosity η₀, and the scaling Stokes sphere velocity (e.g., Turcotte and Schubert, 1982) based on the maximum porosity of the anomaly φ_max

$$\begin{array}{}\text{(16)}& v_{\text{St}}=C_{\text{St}}{\displaystyle \frac{{\mathit{\phi }}_{\text{max}}\mathrm{\Delta }\mathit{\rho }g{r}^{\mathrm{2}}}{{\mathit{\eta }}_{\mathrm{0}}}},\end{array}$$

resulting in the following non-dimensionalization where non-dimensional quantities are primed:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& \left(x,z\right)=\left(x{}^{\prime },z^{\prime }\right)\cdot r,{\mathit{v}}_{\text{s,f}}={\mathit{v}}_{\text{s,f}}^{\prime }\cdot v_{\text{St}},t=t^{\prime }\cdot {\displaystyle \frac{r}{v_{\text{St}}}},\\ & \left({\mathit{\tau }}_{ij},P\right)=\left({\mathit{\tau }}^{\prime },P^{\prime }\right)\cdot {\displaystyle \frac{{\mathit{\eta }}_{\mathrm{0}}{v}_{\text{St}}}{r}},\\ & \left(\mathit{\eta },\mathit{\zeta }\right)=\left(\mathit{\eta }{}^{\prime },{\mathit{\zeta }}^{\prime }\right)\cdot {\mathit{\eta }}_{\mathrm{0}},\left(\mathit{\psi },\mathit{\chi }\right)=\left({\mathit{\psi }}^{\prime },{\mathit{\chi }}^{\prime }\right)\cdot r{v}_{\text{St}}.\end{array}\end{array}$$

For r the half width of the prescribed initial perturbation, consisting of a 2D Gaussian bell, is chosen. This is reasonable as the rising velocity in our code is best described by the Stokes velocity using this radius. The exact shape of the perturbation is given later in the model setup.

C_St is calculated by using the analytic solution of an infinite Stokes cylinder within another cylinder (Popov and Sobolev, 2008, based on the drag force derived by Slezkin, 1955) because the cylinder gets effectively slowed due to boundary effects. C_St is calculated using $C_{\text{St}}=\ln \left(k\right)-\frac{k^{\mathrm{2}}-\mathrm{1}}{k^{\mathrm{2}}+\mathrm{1}}$, where k is the ratio of outer cylinder radius to inner cylinder radius. For our model setup, C_St is equal to 0.17.

With these rules, the Darcy equation (Eq. 10) is given in non-dimensional form

$$\begin{array}{}\text{(18)}& {\mathit{v}}_{\text{f}}^{\prime }-{\mathit{v}}_{\text{s}}^{\prime }=-{\displaystyle \frac{{\mathit{\delta }}_{\text{c}}^{\mathrm{2}}}{r^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\eta }}{}^{\prime }\mathit{\phi }}}\left({\mathit{e}}_z{\displaystyle \frac{\left(\mathrm{1}-\mathit{\phi }\right)}{{\mathit{\phi }}_{\text{max}}}}+{\mathbf{\nabla }}^{\mathbf{\prime }}\cdot {\mathit{\tau }}^{\prime }\right),\end{array}$$

where e_z is the unit vector in the z direction and $\stackrel{\mathrm{̃}}{\mathit{\eta }}{}^{\prime }$ is equal to ${\mathit{\zeta }}^{\prime }+\frac{\mathrm{4}}{\mathrm{3}}{\mathit{\eta }}^{\prime }$. The momentum equation of the mixture Eq. (12) is given by

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& \left({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial {x^{\prime }}^{\mathrm{2}}}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial {z^{\prime }}^{\mathrm{2}}}}\right)\left[{\mathit{\eta }}_{\text{s}}^{\prime }\left({\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}^{\prime }}{\partial {x^{\prime }}^{\mathrm{2}}}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}^{\prime }}{\partial {z^{\prime }}^{\mathrm{2}}}}\right)\right]\\ & +\mathrm{4}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\prime }\partial z^{\prime }}}\left[{\mathit{\eta }}_{\text{s}}^{\prime }{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}^{\prime }}{\partial x^{\prime }\partial z^{\prime }}}\right]={\displaystyle \frac{\mathrm{1}}{{\mathit{\phi }}_{\text{max}}}}{\displaystyle \frac{\partial \mathit{\phi }}{\partial x^{\prime }}}+A^{\prime }\left({\mathit{\chi }}^{\prime }\right).\end{array}\end{array}$$

${\mathit{\delta }}_{\text{c}}^{\mathrm{2}}/r^{\mathrm{2}}$ in Eq. (18) is the squared ratio of compaction length δ_c to the system length scale r, which is the main parameter describing our system. The compaction length is a natural length scale emerging from the problem and of particular importance in our context because 2D porosity waves have half-width radii of the order of 3⋅δ_c to 5⋅δ_c (Simpson and Spiegelman, 2011). It is defined as follows:

$$\begin{array}{}\text{(20)}& {\mathit{\delta }}_{\text{c}}=\sqrt{{\displaystyle \frac{\mathit{\zeta }+\frac{\mathrm{4}}{\mathrm{3}}\mathit{\eta }}{\mathit{\mu }}}{k}_{\mathit{\phi }}}.\end{array}$$

All quantities in the other equations are simply replaced by their non-dimensional primed equivalents (Eqs. 1, 2, 6, 11, 12, 13, and 15).

We now compare the two limits, where segregation or two-phase flow dominates (solitary-wave regime), and where fluid and solid rise together with the same velocity as partially molten bodies, which we identify with the diapir regime. We compare the characteristic segregation velocity within solitary waves, which scales as

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& v_{\text{sgr}}\approx {\displaystyle \frac{k_{\mathrm{0}}{\mathit{\phi }}_{\text{max}}^{n-\mathrm{1}}}{\mathit{\mu }}}\left(\mathrm{\Delta }\mathit{\rho }g\left(\mathrm{1}-{\mathit{\phi }}_{\text{max}}\right)-{\mathbf{\nabla }}^{\mathbf{\prime }}\cdot \mathit{\tau }\right)\\ & =C_{\text{sgr}}{\displaystyle \frac{k_{\mathrm{0}}{\mathit{\phi }}_{\text{max}}^{n-\mathrm{1}}\mathrm{\Delta }\mathit{\rho }g\left(\mathrm{1}-{\mathit{\phi }}_{\text{max}}\right)}{\mathit{\mu }}},\end{array}\end{array}$$

where C_sgr is of the order 0.5 for 2D solitary waves (Schmeling, 2000), with the characteristic Stokes sphere rising velocity given by Eq. (15). The ratio of these is given by

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{v_{\text{sgr}}}{v_{\text{St}}}}={\displaystyle \frac{C_{\text{sgr}}}{C_{\text{St}}}}{\displaystyle \frac{{\mathit{\delta }}_{c\mathrm{0}}^{\mathrm{2}}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\phi }}_{\text{max}}^{n-\mathrm{2}}\left(\mathrm{1}-{\mathit{\phi }}_{\text{max}}\right)}{\stackrel{\mathrm{̃}}{{\mathit{\eta }}_{\mathrm{0}}}{}^{\prime }{\mathit{\phi }}_{\mathrm{0}}^n}}.\end{array}$$

Here $\stackrel{\mathrm{̃}}{{\mathit{\eta }}_{\mathrm{0}}}{}^{\prime }$ refers to $\stackrel{\mathrm{̃}}{\mathit{\eta }}{}^{\prime }$ for the background porosity φ₀ and δ_c0 to the compaction length of the background porosity. In contrast to Scott (1988), who varies the volume viscosity in his model series, we vary the ratio of initial Stokes radius to compaction length.

Thus, in the solitary-wave limit

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{C_{\text{sgr}}}{C_{\text{St}}}}{\displaystyle \frac{{\mathit{\delta }}_{c\mathrm{0}}^{\mathrm{2}}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\phi }}_{\text{max}}^{n-\mathrm{2}}\left(\mathrm{1}-{\mathit{\phi }}_{\text{max}}\right)}{\stackrel{\mathrm{̃}}{{\mathit{\eta }}_{\mathrm{0}}}{}^{\prime }{\mathit{\phi }}_{\mathrm{0}}^n}}\gg \mathrm{1}.\end{array}$$

Darcy's law (18) results in large segregation velocity, which scales as

$$\begin{array}{}\text{(24)}& v_{\text{sgr}}^{\prime }={\displaystyle \frac{C_{\text{sgr}}}{C_{\text{St}}}}{\displaystyle \frac{{\mathit{\delta }}_{c\mathrm{0}}^{\mathrm{2}}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\phi }}_{\text{max}}^{n-\mathrm{2}}\left(\mathrm{1}-{\mathit{\phi }}_{\text{max}}\right)}{\stackrel{\mathrm{̃}}{{\mathit{\eta }}_{\mathrm{0}}}{}^{\prime }{\mathit{\phi }}_{\mathrm{0}}^n}}.\end{array}$$

From Eq. (13), it follows that the irrotational part of the matrix velocity scales with

$$\begin{array}{}\text{(25)}& v_{\mathrm{1}}\approx -{\mathit{\phi }}_{\text{max}}{v}_{\text{sgr}},\end{array}$$

while the rotational part is given by Eq. (19). In that equation, A^′ scales with χ^′, which, via Eqs. (12) and (13), scale with v_sgr, i.e., with ${\mathit{\delta }}_{c\mathrm{0}}^{\mathrm{2}}/r^{\mathrm{2}}$. In other words, the second term on the right-hand side of Eq. (19) dominates for small $r^{\mathrm{2}}/{\mathit{\delta }}_{c\mathrm{0}}^{\mathrm{2}}$ as the first term is of the order of 1. Thus, the rotational matrix velocity has the same order as the irrotational compaction velocity and serves to accommodate the compaction flow. In this limit, the buoyancy term in Eq. (19), $\frac{\mathrm{1}}{{\mathit{\phi }}_{\text{max}}}\frac{\partial \mathit{\phi }}{\partial x^{\prime }}$, is of vanishing importance for the matrix velocity, and the matrix velocity, v₁+v₂, is of the order of φ_maxv_sgr. In the small porosity limit, matrix velocities are negligible with respect to fluid velocities.

In the diapir limit,

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{C_{\text{sgr}}}{C_{\text{St}}}}{\displaystyle \frac{r^{\mathrm{2}}}{{\mathit{\delta }}_{\text{c}}^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\phi }}_{\text{max}}^{n-\mathrm{2}}\left(\mathrm{1}-{\mathit{\phi }}_{\text{max}}\right)}{\stackrel{\mathrm{̃}}{{\mathit{\eta }}_{\mathrm{0}}}{}^{\prime }{\mathit{\phi }}_{\mathrm{0}}^n}}\ll \mathrm{1},\end{array}$$

and Eq. (18) predicts vanishing segregation velocities. As A^′ and χ^′ scale with $r^{\mathrm{2}}/{\mathit{\delta }}_{c\mathrm{0}}^{\mathrm{2}}$ , both vanish in the diapir limit, no irrotational matrix velocity occurs, and Eq. (19) reduces to the classical biharmonic equation (i.e., Stokes equation) driven by melt buoyancy and describing classical diapiric ascent. Segregation velocities are negligible with respect to matrix velocities.

Figure 1 The segregation to Stokes velocity ratio, following Eq. (22), is given as a function of initial perturbation radius r in terms of compaction length δ_c. Each colored line refers to different values of perturbation amplitude φ_max, which are given in the legend.

In Fig. 1, the results of this simple analysis are shown, where we calculated the velocity ratios as a function of initial perturbation radius for several perturbation radii. In our models we use a φ_max of 2 %, for which we get a switch from solitary wave to diapir dominant behavior at $r=\mathrm{48}\cdot {\mathit{\delta }}_{\text{c}}$. Smaller amplitudes lead to a switch at a smaller radius and larger amplitudes to a switch at a larger radius.

2.2 Model setup

The model consists of a $L^{\prime }\times L^{\prime }$ box with a background porosity, φ₀, of 0.5%. L^′ is the non-dimensional side length of the box and equal to 6 times the initial radius of the perturbation. As initial condition, a non-dimensional Gaussian bell-shaped porosity anomaly is placed in the middle of the model at $x_{\mathrm{0}}{}^{\prime }=\mathrm{3}$ and $z_{\mathrm{0}}{}^{\prime }=\mathrm{3}$. The Gaussian wave is given by

$$\begin{array}{}\text{(27)}& \mathit{\phi }={\mathit{\phi }}_{\text{max}}\cdot \exp \left(-{\left({\displaystyle \frac{x^{\prime }-x_{\mathrm{0}}^{\prime }}{w^{\prime }}}\right)}^{\mathrm{2}}-{\left({\displaystyle \frac{z^{\prime }-z_{\mathrm{0}}^{\prime }}{w^{\prime }}}\right)}^{\mathrm{2}}\right),\end{array}$$

where φ_max is the amplitude equal to 0.02 in our models and w^′ corresponds to the width where φ has reached ${\mathit{\phi }}_{\text{max}}/e$. In our case w^′ is equal to 1.2.

In our model series we vary the ratio of Stokes radius to compaction length from 1.8 to 48 to explore the transition from solitary wave towards diapiric regime. The resolution of the models is chosen to be at least 201×201 grid points and was increased for higher ratios of Stokes radius to compaction length so that the compaction length is resolved by at least 3–4 grid points.

At the top and the bottom domain boundaries, we prescribe an outflow and inflow for both melt and solid, respectively, to prevent melt accumulations at the top. The segregation velocity of the background porosity φ₀ is calculated using Eq. (18) without the viscous stress term. The corresponding matrix velocity is calculated using the conservation of mass.

At the sides we enforce no horizontal flux boundary conditions. The permeability–porosity relation exponent in our models is always n=3.

To run models for a longer, practically infinite, amount of time we let the models coordinate system follow the maximum melt fraction.

2.3 Numerical approach

We discretize the set of equations using finite differences on a staggered grid and solve the system using the code FDCON (Schmeling et al., 2019). Starting from the prescribed initial condition for φ and assuming $A^{\prime }\left({\mathit{\chi }}^{\prime }\right)=\mathrm{0}$ at time 0, the time loop is entered and the biharmonic Eq. (19) is solved for ψ^′ by Cholesky decomposition, from which ${\mathit{v}}_{\mathrm{1}}^{\prime }$ is derived. Together with ${\mathit{v}}_{\mathrm{2}}^{\prime }$ the resulting solid velocity is used to determine the viscous stress term in the segregation velocity Eq. (18). This equation and the melt mass Eq. (1) are solved iteratively with strong under-relaxation for φ and ${\mathit{v}}_{\text{f}}^{\prime }-{\mathit{v}}_{\text{s}}^{\prime }$ for the new time step using the upwind scheme and an implicit formulation of Eq. (1). During this internal iteration these quantities are used, via Eq. (13), to give ∇⋅v_s, the divergence of the matrix velocity, which is needed in the viscous stress term (Eq. 6). After convergence, ∇⋅v_s is used via Eq. (12) to determine χ by LU decomposition and then to get ${\mathit{v}}_{\mathrm{2}}^{\prime }$. Now $A^{\prime }\left({\mathit{\chi }}^{\prime }\right)$ can be determined to be used on the right-hand side of Eq. (19). The procedure is then repeated upon entering the next time step.

Time steps are dynamically adjusted by the Courant criterion times 0.2 based on the fastest velocity (either melt or solid).

The model resolution is a critical parameter in this kind of numerical calculation and should always be kept in mind. With increasing length scale ratio, the compaction length in the model gets smaller, and the resolution needs to be increased to keep it equally resolved.

According to several authors (e.g., Räss et al., 2019; Keller et al., 2013), the compaction length should be resolved by at least 4–8 grid points to solve for waves with sufficient accuracy. For small length scale ratios this is no problem, where up to nearly 30 grid points per compaction length can be achieved with a model resolution of 201×201. The highest resolution our code can run is 601×601, which is enough to resolve the compaction length by three grid points for the model with a length scale ratio of 48. Everything above that cannot be sufficiently resolved with respect to studying solitary waves.

Figure 2 The six panels depict a model with an initial perturbation radius of 12 times the compaction length but with different numerical grid resolutions: (a) 13×13 (b) 26×26 (c) 51×51 (d) 101×101, (e) 201×201, and (f) 401×401. The size of the compaction length in terms of grid length is given in the lower-left corner in each panel.

Figure 2 shows the resulting models for a length scale ratio of 12 for six different resolutions. The model states where φ_max has risen to approximately 0.25 times the initial Stokes radius ($t{}^{\prime }=\mathrm{0.25})$ are shown. With increasing resolution, the maximum melt fraction increases strongly from 101×101 to 401×401 by approximately 20 %, but the velocity of φ_max decreases by 7 % (not shown in the figure). Both values converge for resolutions higher than 51×51, corresponding to ${\mathit{\delta }}_{\text{c}}/\mathrm{d}x=\mathrm{1}$. Even though the compaction length is not sufficiently resolved in Fig. 2d, one can still observe the main features of the model: a main solitary wave has emerged from the original Gaussian perturbation, and secondary porosity waves are beginning to emerge within its wake. Even with ${\mathit{\delta }}_{\text{c}}/\mathrm{d}x=\mathrm{1}$, these features can be observed but are clearly underresolved. With even lower resolutions, accumulations at the top of the perturbation can be seen, which can be broadly interpreted as the attempt of a solitary wave to build up. With ${\mathit{\delta }}_{\text{c}}/\mathrm{d}x=\mathrm{0.24}$, the model is too coarse and the results cannot be trusted anymore.

The solitary waves modeled with our code have been compared to the semi-analytical solution of Simpson and Spiegelman (2011), and more benchmarking was carried out in Dohmen et al. (2019).

In a single-phase flow case, where the melt is not allowed to move relative to the solid, the initial perturbation ascends, shortly after beginning, with a velocity of 0.95 times the calculated Stokes velocity and then slowly decreases as the original Gauss-shaped wave deforms and loses its amplitude.

3 Results

3.1 The transition from porosity wave to diapirism: Varying the initial wave radius

In this model series we vary the initial wave radius to cover the transition from porosity waves towards diapirism. As a reminder, due to our scaling the initial wave has always the same size with respect to the model box, and “increasing the initial wave radius” is equivalent to decreasing the compaction length or the size of the emerging solitary waves with respect to the model box. In Fig. 3, the models are shown at $t{}^{\prime }=\mathrm{0.2}$. For small radii ($r\le \mathrm{12}\cdot {\mathit{\delta }}_{\text{c}})$, a single porosity wave emerges from the original perturbation. The melt that is not situated within the emerging wave is left behind and has (for the most part) already left the model region. For $r=\mathrm{2.4}\cdot {\mathit{\delta }}_{\text{c}}$, the emerged solitary wave is about the size of the initial perturbation, and even smaller radii would lead to too big waves that would not fit into the model. With increasing radius, the emerging solitary wave gets smaller. With $r=\mathrm{12}\cdot {\mathit{\delta }}_{\text{c}}$, the resulting wave only has a size that is ∼ 20 % of the initial perturbation size.

Figure 3 Melt ascent morphology as function of initial perturbation radius in terms of compaction length. (a) Initial conditions of the model are valid for all cases apart of the change in compaction length. (b–j) Melt fraction distribution after $t{}^{\prime }=\mathrm{02}$ for length scale ratios varying between 2.4 and 48. (k) Diapiric rise resulting from a compaction length of zero at $t{}^{\prime }=\mathrm{9}$. (l) Model transition time as a function of length scale ratios varying between 1.8 and 120. The transition time gives the time after which the main wave has reached solitary-wave status.

We compare the observed solitary-wave velocities of Fig. 3b–e to equivalent Stokes velocities for a diapir based on Eq. (16). While the dimensional Stokes velocity of a porosity anomaly is proportional to the amplitude of porosity and the square of the radius, the non-dimensional Stokes velocity is always equal to 1. In Fig. 4, this non-dimensional Stokes velocity is indicated by the dashed line with the value 1. The colored lines give 2D solitary-wave velocities with their appropriate radii, given by Simpson and Spiegelman (2011), normalized by the Stokes velocity corresponding to different initial perturbation radii. These semi-analytical solutions are in good agreement with our solitary-wave models and differ only by 3 %–5 % percent in velocity, as already shown in Dohmen et al. (2019). The velocities in this figure correspond to ratios of solitary-wave velocity to initial perturbation Stokes velocity. Inspection of Fig. 4 reveals that for the first four cases of Fig. 3b–e with radii smaller than or equal to 12⋅δ_c, the phase velocities are always larger than the Stokes velocity. For example, for $r=\mathrm{12}\cdot {\mathit{\delta }}_{\text{c}}$, an emerging solitary wave with a typical radius of 4.5⋅δ_c has a higher phase velocity than a $r=\mathrm{12}\cdot {\mathit{\delta }}_{\text{c}}$ melt anomaly rising by Stokes flow. Thus, the cases are always in the solitary-wave regime.

For greater radii (e.g., $r=\mathrm{18}\cdot {\mathit{\delta }}_{\text{c}}-\mathrm{30}\cdot {\mathit{\delta }}_{\text{c}}$, Fig. 3e–g) the phase velocities of solitary waves are of the order of the Stokes velocity (see Fig. 4), and they therefore need more time to separate from the remaining melt of the initial perturbation, still rising with the order of the Stokes velocity. The amount of melt accommodated within the main solitary wave is just a small percentage of the original perturbation, and secondary waves evolve in its remains. With further ascending, more and more solitary waves build up, and the former perturbation will sooner or later consist of solitary waves in an ordered cluster or a formation. This formation elongates during ascent as the main wave has a larger amplitude than all of the following waves, whose amplitudes are also decreasing with depth, as a higher proportion of melt accumulated at the top of the perturbation. Similar formations of strongly elongated fingers can be also observed in 3D, as shown by Räss et al. (2019), who used decompaction weakening. In the models with smaller radii, the main solitary wave consisted of the majority of melt originally situated within the perturbation, and just few or none secondary waves are able to emerge. With increasing radius, more melt is left behind, which allows for second and higher generations of solitary waves.

Figure 4 The dashed line marks the velocity of the Stokes sphere ($v{}^{\prime }=\mathrm{1})$. The colored lines refer to the velocity of a 2D solitary wave, calculated semi-analytically by Simpson and Spiegelman (2011), in our non-dimensionalization, based on the radii shown in the legend.

For greater radii (e.g., $r=\mathrm{24}\cdot {\mathit{\delta }}_{\text{c}}-\mathrm{48}\cdot {\mathit{\delta }}_{\text{c}}$, Fig. 3f–j), the phase velocities of solitary waves are almost equal to the Stokes velocity (See Fig. 4). This leads to almost no separation after $t{}^{\prime }=\mathrm{0.2}$. While for $r=\mathrm{36}\cdot {\mathit{\delta }}_{\text{c}}$ a solitary wave has already built up and is rising just ahead of the perturbation, for $r=\mathrm{42}\cdot {\mathit{\delta }}_{\text{c}}$ and $r=\mathrm{48}\cdot {\mathit{\delta }}_{\text{c}}$ just the accumulation of melt at the top of the perturbation can be observed, which will eventually lead to a solitary wave. Secondary waves also build up with higher run times, as can already be seen for $r=\mathrm{36}\cdot {\mathit{\delta }}_{\text{c}}$.

For even greater radii, the compaction length cannot be sufficiently resolved with our approach, but tests using models that are not sufficiently resolved have shown that solitary waves can be observed for $r\ge \mathrm{48}\cdot {\mathit{\delta }}_{\text{c}}$. At some point they no longer appear, probably due to lack of sufficient resolution, but our tests show that solitary waves should always emerge, even if its phase velocity is way below the Stokes velocity. As long as the ascending time is long enough and melt is able, independent of segregation velocity, to move separately to the matrix, a diapir will evolve into a swarm of a certain number of solitary waves based on the compaction length. Because the phase velocities of each small solitary wave are small compared to the Stokes velocity of the full swarm, we consider such a rising formation of melt a large-scale diapir.

Figure 3l shows the required time for the initial perturbation to build up a solitary wave. This status is achieved after the dispersion relation of the main wave reaches a point from where it follows the solitary wave dispersion relation. This time increases nearly linearly for small radii ($r\le \mathrm{48}\cdot {\mathit{\delta }}_{\text{c}})$ but increases nonlinearly for greater radii. This might be due to lack of proper resolution, but a nonlinear trend can be already observed for small radii. The transition time for radii smaller than 30⋅δ_c is smaller than 0.2, the time at which the models in Fig. 3b–j are shown. The other models already show solitary-wave-like blobs but did not yet reach their final form.

A classical diapir will evolve only in cases with zero compaction length ($r=\mathrm{\infty }\cdot {\mathit{\delta }}_{\text{c}})$, i.e., melt is not able to move with respect to the matrix (Fig. 3k). Here, no focusing into solitary waves can be observed and transition time is infinity.

Figure 5 The upper row depicts the solid and fluid mass fluxes of a horizontal line cutting through the maximum melt fraction at time steps where the main wave has just reached the status of a solitary wave. These time steps are $t{}^{\prime }=\mathrm{002}$, 0068, 0155, and 0416 from left to right, respectively. The bottom row depicts the corresponding melt porosity fields. All quantities shown are non-dimensional.

Summarizing Fig. 4, the comparison of Stokes and porosity wave velocities explains our observations shown in Fig. 3 well. For small initial radii the solitary wave velocity is clearly higher and will therefore build up and separate from the melt left behind quickly. For cases with approximately equal perturbation to solitary-wave radius only one solitary wave will build up, which includes most of the melt of the initial perturbation. With increasing perturbation radius, the velocity ratio decreases and multiple solitary waves, requiring more time, will emerge, each including only a fraction of the melt originally situated in the initial perturbation. However, even with velocity ratios smaller than 1, solitary waves emerge and, not being able to separate, rise just ahead of the remains, slowly elongating the initial perturbation.

3.2 Effects on the mass flux

It is important to study the partitioning between rising melt and solid mass fluxes in partially molten magmatic systems because melts and solids are carriers of different chemical components. Within our Boussinesq approximation, we may neglect the density differences between solid and melt. Thus, our models allow the evaluation of vertical mass fluxes of solid or fluid by quantifying the vertical velocity components multiplied with the melt or solid fractions, respectively:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}& q_{sz}{}^{\prime }=\left(\mathrm{1}-\mathit{\phi }\right)\cdot v{{}^{\prime }}_{sz}\\ & \text{and}\\ & q_{fz}{}^{\prime }=\mathit{\phi }\cdot {v^{\prime }}_{fz}.\end{array}\end{array}$$

Horizontal profiles of the mass fluxes through rising melt bodies are calculated at the vertical positions of maximum melt fraction at time steps where the main wave has just reached the status of a solitary wave (Fig. 5).

The mass fluxes of solid and fluid are strongly affected by the change of the initial radius from the solitary-wave regime towards the diapiric regime. For $r=\mathrm{2.4}\cdot {\mathit{\delta }}_{\text{c}}$, where we observe a solitary wave, the fluid has its peak mass flux in the middle of the wave and the solid is going downwards, i.e., against the phase velocity. In the center, the fluid flux is about 10 times higher than the solid net flux. The upward flow in the center is balanced by the matrix-dominated downward flow inside and outside the wave. For $r=\mathrm{12}\cdot {\mathit{\delta }}_{\text{c}}$, the wave area is much smaller and the ratio between solid and fluid flux is still around the order of 10. At the boundary of the wave, the solid is nearly not moving at all, but a minimum can be observed within the center of it. For $r{}^{\prime }=\mathrm{24}\cdot {\mathit{\delta }}_{\text{c}}$, the solid flux is just above zero in the center and increases to a maximum towards the flanks of the wave that is still about 10 times smaller than the maximum fluid flux.

With $r{}^{\prime }=\mathrm{48}\cdot {\mathit{\delta }}_{\text{c}}$, the solid flux is about 3 times smaller than the fluid flux, but most of the material ascent is accomplished by the solid. This suggests that diapiric rise begins to dominate.

The transition from solitary waves towards diapirism on qualitative model observations has so far only been based on observations. We now invoke a more quantitative criterion. In a horizontal line passing through the anomaly's porosity maximum, we define the total vertical mass flux of the rising magma body by ${\int }_{\mathit{\phi }>{\mathit{\phi }}_{\mathrm{0}}}\left(q_{\text{f}}+q_{\text{s}}\right)\mathrm{d}x$, where the integration is carried out only in the region of increased porosity φ>φ₀. This mass flux is partitioned between the fluid mass flux, ${\int }_{\mathit{\phi }>{\mathit{\phi }}_{\mathrm{0}}}{q}_{\text{f}}\mathrm{d}x$, and the solid mass flux, ${\int }_{\mathit{\phi }>{\mathit{\phi }}_{\mathrm{0}}}{q}_{\text{s}}\mathrm{d}x$. With these we define the partition coefficients

$$\begin{array}{}\text{(29)}& C_{\mathrm{soli}}={\displaystyle \frac{{\int }_{\mathit{\phi }>{\mathit{\phi }}_{\mathrm{0}}}{q}_{\text{f}}\mathrm{d}x}{{\int }_{\mathit{\phi }>{\mathit{\phi }}_{\mathrm{0}}}\left(q_{\text{f}}+q_{\text{s}}\right)\mathrm{d}x}},\end{array}$$

and

$$\begin{array}{}\text{(30)}& C_{\mathrm{dia}}={\displaystyle \frac{{\int }_{\mathit{\phi }>{\mathit{\phi }}_{\mathrm{0}}}{q}_{\text{s}}\mathrm{d}x}{{\int }_{\mathit{\phi }>{\mathit{\phi }}_{\mathrm{0}}}\left(q_{\text{f}}+q_{\text{s}}\right)\mathrm{d}x}}.\end{array}$$

The sum C_soli+C_dia is always 1, and if C_soli>C_dia, then the solitary wave proportion is dominant, while for C_soli<C_dia diapirism is dominant. In Fig. 6a these partition coefficients for several initial radii are shown. In red the diapir is shown, and in blue the solitary wave partition coefficients are shown.

For $r=\mathrm{1.8}\cdot {\mathit{\delta }}_{\text{c}}$, C_soli is equal to 5 and C_dia is equal to −4, i.e., we have a downward solid flux. With increasing radius, C_dia increases until it changes its sign and the matrix flows upward, at $r\approx \mathrm{20}\cdot {\mathit{\delta }}_{\text{c}}$. It eventually becomes bigger than C_soli at $r=\mathrm{36}\cdot {\mathit{\delta }}_{\text{c}}$ and then approaches 1 for bigger radii. C_soli changes so that the sum of both is always equal to 1. Even though diapirism is dominant for $r>\mathrm{36}\cdot {\mathit{\delta }}_{\text{c}}$, we still observe solitary waves, yet their phase velocities are much smaller than the large-scale rising velocities of the full melt formation.

The ratio of maximum fluid velocity (i.e., v_f) to absolute matrix velocity (Fig. 6b) shows that for small radii, where C_soli≫C_dia, this ratio is approximately constant with a high value of about 100. The absolute velocity maxima themselves are not constant but decrease with the same rate until the switch of negative to positive matrix mass flux, where the absolute matrix velocity starts to increase, while the fluid velocity keeps decreasing. At this zero crossing we would expect a ratio of infinity, but while the zero crossing takes place within the center of the solitary wave, other regions near the wave still have finite vertical velocities. This switch from negative to positive mass flux was already observed by Scott (1988), but while they changed the viscosity ratio as an independent constant model parameter, we change the radius and keep the viscosity law the same while still evolving with φ. Both describe the transition from a two-phase limit towards the Stokes limit, but in our formulation we are able to reach the Stokes limit, whereas Scott's formulation (1988) is restricted to two-phase flow. With even greater radii the velocity ratio will eventually converge towards 1, where melt is no longer able to move relatively to the matrix (i.e., v_f=v_s), and material will be transported collectively as in single-phase flow. These last models are not sufficiently resolved to obtain leading and secondary solitary waves but still show the expected behavior in terms of macroscopically rising partially molten diapir.

Figure 6 Quantitative parameters as function of initial perturbation radius in terms of compaction length. (a) Solitary wave (blue) and diapir (red) partition coefficients for several initial perturbation radii. (b) Ratio of maximum fluid velocity to maximum absolute solid velocity in the entire model.

Based on these observations, the evolution of these models can be divided into three regimes. (i) In the solitary wave regime ($r\le \mathrm{36}\cdot {\mathit{\delta }}_{\text{c}})C_{\mathrm{soli}}$ is larger than C_dia, and the initial perturbation emerges into waves that have the properties of solitary waves and ascend with constant velocity and staying in shape. This regime can be further divided into (ia) ($r<\mathrm{20}\cdot {\mathit{\delta }}_{\text{c}})$, where the solid mass flux is negative, and (ib) ($\mathrm{20}\cdot {\mathit{\delta }}_{\text{c}}\le r<\mathrm{36}\cdot {\mathit{\delta }}_{\text{c}})$, where the solid moves upwards with the melt. Waves in these regimes are very similar, but the further we are in regime (ia), the fewer solitary waves will emerge out of the initial perturbation. For radii smaller than about 4.8⋅δ_c only one wave will merge. In regime (ib) the perturbation will always emerge into multiple solitary waves.

In the diapirism-dominated regime of (ii) ($r\ge \mathrm{36}\cdot {\mathit{\delta }}_{\text{c}})$, C_dia is larger than C_soli, but as the fluid melt is still able to move relatively to the solid matrix, solitary waves build up and the whole partially molten region will evolve into a swarm of them. The phase velocities of these waves are very small compared to the Stokes velocity of the perturbation and the whole swarm will rise as a diapir, whose buoyancy is still comparable to the buoyancy of the initial perturbation.

The endmember of the second regime can be reached by prohibiting the relative movement of fluid ($r=\mathrm{\infty }\cdot {\mathit{\delta }}_{\text{c}})$ for which the compaction length has not to be sufficiently resolved. In this regime the initial perturbation will not disintegrate into solitary waves but instead will rise as a well-formed partially molten diapir. In every other case in the present model where fluid is able to move with respect to the solid, at some point all diapirs will evolve into a swarm of solitary waves that can be infinitely small compared to the initial perturbation. However, this is expected to happen only after a long distance of diapiric rise. In cases where the size of solitary waves is comparable to the perturbation (e.g., regime i) this will occur sooner, and in cases where solitary waves are much smaller this will occur later. Their observation is mostly limited by resolution. For models that allow for the diapir to grow (e.g., Keller et al., 2013), they may not dissolve into solitary waves as they approach the single-phase limit.

4 Discussion

4.1 Application to nature

While in our models the perturbation size in terms of compaction lengths was systematically varied but kept constant within in each model, our results might also be applicable to natural cases in which the compaction length varies vertically. In the case of compaction length decreasing with ascent, a porosity anomaly might start rising as a solitary wave, but it may at some point enter the second regime where diapiric rise is dominant. If this boundary is sharp, the solitary wave might disintegrate into several smaller solitary waves that rise as a diapiric swarm. If the boundary is a continuous transition, the wave should slowly shrink and become slower. The melt left behind might also evolve into secondary solitary waves.

A decreasing compaction length could be accomplished by decreasing the matrix viscosity or the permeability or by increasing the fluid viscosity. Decreasing matrix viscosity might be explainable by, for example, local heterogeneities, temperature anomalies due to secondary convective overturns in the asthenosphere, or a vertical gradient of water content, which may be the result of melt-segregation-aided volatile enrichment at shallow depths in magmatic systems. This could lead to the propagation of magma-filled cracks (Rubin, 1995), as already pointed out in Connolly and Podladchikov (1998). The latter authors have looked at the effects of rheology on compaction-driven fluid flow and had similar results for an upward weakening scenario. A decrease in permeability due to a decrease in background porosity might be an alternative explanation. In the hypothetical case of a porosity wave reaching the top of partially molten region within the Earth's upper mantle or lower crust, the background porosity might decrease, which would most certainly lead to focusing because the compaction length will decrease, and eventually, when reaching melt-free rocks, the solitary waves might be small enough and have a high enough amplitude to trigger the initiation of dikes.

Even though most diapirs should, according to our models, disintegrate into numerous solitary waves, not all will inevitably do so. Within regime (i), solitary waves are possible and most probably expected, but the deeper we are in regime (ii) the less expected the disintegration is because a long time is needed for it to build up. In nature, in contrast with our models, they cannot rise for an infinite amount of time. The time needed to build up a solitary wave increases nonlinearly with r (see Fig. 3l). For example, while for $r=\mathrm{4.8}\cdot {\mathit{\delta }}_{\text{c}}$ a solitary wave is completely evolved after $t{}^{\prime }=\mathrm{0.02}$, for $r=\mathrm{48}\cdot {\mathit{\delta }}_{\text{c}}$ it needs until $t{}^{\prime }=\mathrm{0.4}$, i.e., equivalent to the diapiric rise time necessary to ascend the distance of approximately half the initial radii. Additionally, as already pointed out, if a model setup allows for the diapir to grow, it could approach the single-phase flow, prohibiting the emergence of solitary waves (see Keller et al., 2013).

4.2 Model limitations

The introduced partition coefficients help to distinguish whether solitary-wave or diapiric rise is dominant but cannot be solely consulted as to whether a solitary wave or a diapir can be expected. As the fluid velocity and flux are still very high in the wave center for diapiric-dominant cases, small solitary waves will build up. However, the net mass flux is dominated by the large-scale rising solid, and the formation time of small solitary waves might be long. Additionally, the internal circulation of diapirs can be faster than the phase velocity, which would smear out the emergence of solitary waves and would not allow for them to emerge. Due to the limitations of our model, we are not able to reach regions where solitary waves are small enough and their phase velocity slow enough to observe this.

While the minimum size of solitary waves in nature might be in some way limited by the grain size, in numerical models the minimum size is limited by the model's resolution. We restrict our models in this study to cases where the compaction length is at least resolved by three grid lengths dx (i.e., ${\mathit{\delta }}_{\text{c}}\ge \mathrm{3}\cdot \mathrm{d}x)$ to get fairly resolved solitary waves, but they can be also observed for compaction lengths that are resolved much more poorly. The resolution test (Fig. 2) shows that, even though they are not solved decently, probable solitary waves can be observed for cases with δ_c=dx. Smaller resolutions can show indications of solitary waves but should not be trusted, as other tests (not shown here) with similar resolutions result in spurious channeling. For very poorly resolved compaction lengths (${\mathit{\delta }}_{\text{c}}<\mathrm{0.25}\cdot \mathrm{d}x$ for our models), no indications of solitary waves can be observed, and the partially molten perturbation ascends as a diapir. The deeper we are in regime 2, the more dominant the dynamics of diapirism are on a length scale of r compared to Darcy flow or solitary waves on the unresolved length scale of δ_c. Thus, two-phase flow, either Darcy flow or solitary waves, becomes negligible for r≫δ_c, and partially molten diapirs can be regarded as well resolved.

5 Conclusions

This work shows that, depending on the extent of a partially molten region within the Earth, the resulting ascent of melt may not only occur by solitary waves or by diapirs but instead by a composed mechanism where a diapir splits up into numerous solitary waves. Their phase velocities might become so slow that the whole swarm will ascend as a diapir, only slowly elongating due to the main solitary wave having a higher amplitude and therefore higher phase velocity than the following ones. Depending on the ratio of the melt anomaly size to the compaction length (or rather the model length scale to compaction length ratio), we can classify the ascent behavior into two different regimes using mass flux and velocity of matrix and melt: (ia + b) solitary wave a and b regimes, and the (ii) diapirism-dominated regime. In regime (ia), the matrix sinks with respect to the rising melt, and in regime (ib) the matrix also rises but only very slowly. The further we are in this regime, the fewer solitary waves will emerge out of the initial perturbation until eventually only one solitary wave will emerge. On the first order, these regimes can be explained by comparing Stokes velocity of the rising perturbation with the solitary-wave phase velocity. If the solitary-wave velocity is higher than the Stokes velocity, a solitary wave will evolve, and if it is lower it means that diapirism is dominant but that solitary waves will still build up if the ascending time is long enough. The deeper we are in regime (ii), the more time is needed to build up solitary waves and the less likely it is that they will appear in nature. The endmember regime (ii), i.e., pure diapirism, can be reached if fluid is not allowed to move separately from the matrix.

Numerical resolution plays an important role, especially around the transition of the regimes, as the compaction length may be under-resolved to allow for the emergence of solitary waves. Hence, it should be generally important for two-phase flow models to inspect whether solitary waves are expected and, if so, whether they have a major influence on the conclusions made.