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

5 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 a diapir or a porosity wave is 10 dominant. 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 15 with a positive or negative vertical matrix flux and 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.


Introduction
In many scenarios inside the earth the process of a fluid moving relatively to a viscously deformable porous matrix is an important transport mechanism.The physics of these scenarios were firstly described by McKenzie (1984) and it was later shown by several authors that these equations allow for the emergence of solitary porosity waves (Scott & Stevenson, 1984;Barcilon & Lovera 1989;Wiggins & 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.
Porosity waves were of vast interest for many authors over the last decades and the possible consequences on geochemistry and fluid flow in lower and middle crust in general (e.g. Watson & The formulation of the governing equations for the melt-in-solid two-phase flow dynamics is based on McKenzie (1984), Spiegelman & McKenzie (1987) and Schmeling (2000) assuming an infinite Prandtl number, a low fluid viscosity w.r.t. 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 fluid (melt) have the subscript  and those associated with the solid have the subscript .Applying the Boussinesq approximation the equation for the conservation of the mass of the melt is and the mass conservation of the solid is (2) is the volumetric rock porosity (often called melt fraction),   and   are the fluid and solid velocities, respectively.The momentum equations are given as a generalized Darcy equation for the fluid separation flow where   is the fluid density and  is the fluid pressure (including the lithostatic pressure), whose gradient is driving the motion.The Stokes equation for the mixture is given as is the permeability that depends on the rock porosity is the melt dynamic viscosity,  is the gravitational acceleration,  is the density of the meltsolid mixture and   is the viscous stress tensor is the bulk viscosity.The linearized equation of state for the mixture density is given as with  0 as the solid density and   =  0 −   0 . The shear and bulk viscosity are given by the simple equations and where  0 is the constant intrinsic shear viscosity of the matrix.
As in both equations ( 3) and ( 4)  is the fluid pressure (see McKenzie, 1984, Appendix A), these equations can be merged to eliminate the pressure resulting in ). (10) 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 ( 3) can be rewritten by expressing the matrix velocity,   , as the sum of the incompressible flow velocity,  1 , and the irrotational (compaction) flow velocity,  2 , as: with  as stream function and  as the irrotational velocity potential, given as the solution of the Poisson equation The divergence term ∇ ⃗ ⃗ ⋅   ⃗⃗⃗ can be derived from eqs. 1 and 2 to give 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 (equ.4), 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 14 Pa s (see the discussion in Schmeling et al., 2019), while effective viscosities of mafic or silicic partially molten rocks may range between 10 20 Pa s and 10 16 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 (11) into the viscous stress (6), this into the Stokes equation ( 4), and taking the curl of the x-and z equations the pressure is eliminated and one gets with To describe the transition from solitary waves to diapirs it is useful to non-dimensionalize the equations.As scaling quantities we use the radius  of the anomaly, the reference viscosity  0 , and the scaling Stokes sphere velocity (e.g.Turcotte & Schubert, 1982) based on the maximum porosity of the resulting to the following non-dimensionalization where non-dimensional quantities are primed: For  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.
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, due to boundary effects, the cylinder gets effectively slowed.  is calculated using , where k is the ratio of outer cylinder's to inner cylinder's radius.For our model setup   is equal to 0.17.
With these rules the Darcy equation ( 10) is given in non-dimensional form where   is the unit vector in z-direction and  ̃′ is equal to   ′ + 4 3   ′ , which is occasionally referred to as compaction viscosity.The momentum equation of the mixture ( 12) is given by 2 / 2 in equation ( 17) is the squared ratio of compaction length   to the system length scale , which is the main parameter describing our system.The compaction length is a typical length scale used in two-phase flow problems and of particular importance in our context, because 2D porosity waves have half width radii of the order of 3 ⋅   to 5 ⋅   (Simpson and Spiegelman, 2011).It is defined as: In the other equations ( 1), ( 2), ( 6), ( 11), ( 12), (13), and (14a) all quantities are simply replaced by their non-dimensional primed equivalents.
We now can 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 (batch melting), which we identify with the diapir regime.This can be done by comparing the characteristic segregation velocity within solitary waves, which scales as where   is of the order ½ for 2D solitary waves (Schmeling, 2000), with the characteristic Stokes sphere rising velocity given by (15).The ratio of these is given by Here  0 ̃′ refers to  ̃′ for the background porosity  0 .In contrast to Scott (1988), who varies the compaction viscosity in his model series, we vary the ratio of initial Stokes radius to compaction length.
Thus, in the solitary wave limit Darcy's law (17) results in large segregation velocity, which scales as From equation (13) it follows that the irrotational part of the matrix velocity scales with while the rotational part is given by ( 18  2 as the first term is of the order 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 equation ( 18), 1    ′ , is of vanishing importance for the matrix velocity and the matrix velocity,  1 +  2 , is of the order of     .In the small porosity limit, matrix velocities are negligible with respect to fluid velocities.
In the diapir limit, and equation ( 17) predict vanishing segregation velocities.As ′ and ′ scale with 2 , both vanish in the diapir limit, no irrotational matrix velocity occurs and equ.( 18) 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.
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   of 2%, for which we get a switch from solitary wave to diapir dominant behavior at  = 48 ⋅   .For bigger amplitudes this switch happens at a bigger radius and for smaller amplitudes the other way around.

Model setup
The model consists of a ′ × ′ box with a background porosity,  0 , of 0.5%.′ 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 wave porosity anomaly is placed in the middle of the model at  0 ′ = 3 and  0 ′ = 3.The Gaussian wave is given by Where   is the amplitude equal to 0.02 in our models and  ′ corresponds to the width where  has reached   /.In our case ′ 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 parameter range towards the diapiric regime.The resolution of the models is chosen to be at least 201 × 201 and was increased for higher length scale ratios so that the compaction length is resolved by at least 3-4 grid points.
At the top and the bottom, we prescribe an out-and inflow for both melt and solid, respectively, which is calculated analytically for the background porosity.This is necessary because we have a background melt fraction  0 that would lead to melt accumulation at the top of the model.We therefore calculate the segregation velocity of the background porosity  0 using equation ( 17) without the viscous stress term.The corresponding matrix velocity is calculated using the conservation of mass.
At the sides we use mirroring boundary conditions, which corresponds to a symmetry axis, where no horizontal flow is allowed.The permeability-porosity relation exponent in our models is always  = 3.
To run models for a longer, practically infinite, amount of time we let the models coordinate system follow the maximum melt fraction.This is done by shifting the whole model down one grid length after   has risen above + .This procedure allows us to zoom into the perturbation and follow it, not knowing its velocity and without carrying out any interpolations, which would strongly influence the model.

Numerical approach
The above equations in non-dimensional form are solved by the finite differences code FDCON developed essentially by one of the authors (Schmeling et al., 2019).Starting from the prescribed initial condition for , and assuming ′(′) = 0 at time 0, the time loop is entered and the biharmonic equation ( 19) is solved for ′ by Cholesky decomposition, from which  1 ′ is derived.Together with  2 ′ the resulting solid velocity is used to determine the viscous stress term in the segregation velocity equation ( 17).This equation and the melt mass equation ( 1) are solved iteratively with strong underrelaxation for  and   ′ −   ′ for the new time step using upwind and an implicit formulation of equ.(1).During this internal iteration these quantities are used, via equ.( 13), to give ∇ ⃗ ⃗ •   , the divergence of the matrix velocity, which is needed in the viscous stress term (equ.6).After convergence ∇ ⃗ ⃗ •   is used via equ.
(12) to determine  by LU-decomposition and then to get  2 ′.Now ′(′) can be determined to be used on the RHS of equ ( 18).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 calculations 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 at least resolved by 4-8 grid points to solve for waves sufficiently accurately.For small length scale ratios this is no problem, where, with a model resolution of 201 × 201, up to nearly 30 grid points per compaction length can be achieved.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.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   / = 1 these features can be observed but are much worse emphasized.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   / = 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 & 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 relatively 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 in amplitude.

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 w.r.For greater radii (e.g. = 18 ⋅   − 30 ⋅   , 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 order of Stokes velocity.The amount of melt accommodated within the main solitary wave is just a small percentage of the original perturbation and secondary waves can 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 ascend as the main wave has a larger amplitude than all 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 the emergence of secondary waves turns out zero or small, but with greater radii enough melt is left behind to observe the emergence of second and higher generations of solitary waves.
For greater radii (e.g. = 24 ⋅   − 48 ⋅   , Fig. 3 fj) the phase velocities of solitary waves are almost equal to the Stokes velocity (See Fig. 4).This leads to almost no separation after  ′ = 0.2.While for  = 36 ⋅   a solitary wave has already build up and is rising just ahead of the perturbation, for  = 42 ⋅   and  = 48 ⋅   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 runtimes, as can be already seen for  = 36 ⋅   .
For even greater radii the compaction length cannot be sufficiently resolved with our approach, but tests with not sufficiently resolved models have shown that solitary waves can be observed for  ≥ 48 ⋅   .
At some point they do 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 to move separately to the matrix, independently of segregation velocity, 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 is small compared to the Stokes velocity of the full swarm we consider such a rising formation of melt as a large scale diapir.
Fig. 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 ( ≤ 48 ⋅   ) but increases non-linearly for greater radii.This might be due to lack of proper resolution, but a non-linear trend can be already observed for small radii.The transition time for radii smaller than 30 ⋅   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.
Just in the case where the compaction length is zero ( = ∞ ⋅   ), i.e. melt is not able to move w.r.t. the matrix, a classical diapir, as shown in Fig. 3k, will evolve.Here, no focusing into solitary waves can be observed and the transition time in infinity.
Summarizing Fig. 4, the comparison of Stokes and porosity wave velocities correlates nicely with our observations shown in Fig. 3: For small initial radii the solitary wave velocity is clearly higher and will therefore build up and separate from the melt left behind quickly.With increasing radius, the velocity ratio decreases, and the waves need more time to build up and separate.But even with velocity ratios smaller than 1, solitary waves emerge and, not able to separate, rise just ahead of the remains, slowly elongating the initial perturbation.

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.Then our models allow to evaluate vertical mass fluxes of solid or fluid by depicting the vertical velocity components multiplied with the melt or solid fractions, respectively: Fig. 5 shows horizontal profiles through rising melt bodies at the vertical positions of maximum melt fraction at timesteps where the main wave has just reached the status of a solitary wave.
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  = 2.4 ⋅   , 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, 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  = 12 ⋅   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  ′ = 24 ⋅   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 ten times smaller than the maximum fluid flux.
With  ′ = 48 ⋅   the solid flux is just about three times smaller than the fluid flux, but most of the material ascend is accomplished by the solid.This suggests that diapiric rise begins to dominate.
So far, we have based our discussion of the transition from solitary waves towards diapirism on qualitative model observations.We now try to invoke a more quantitative criterion.In a horizontal line .With these we define the partition coefficients and The sum   +   is always 1 and if   >   then the solitary wave proportion is dominant, while for   <   diapirism is dominant.In Fig. 6a these partition coefficients for several initial radii are shown.In red are the diapir and in blue the solitary wave partition coefficients.
For  = 1.8 ⋅   ,   is equal to 5 and   is equal to -4, i.e. we have a downward solid flux.With increasing radius   increases until it changes its direction at  ≈ 20 ⋅   , eventually becomes bigger than   at  = 36 ⋅   and then approaches 1 for bigger radii.  changes so that the sum of both is always equal to 1.Even though diapirism is dominant for  > 36 ⋅   , we still observe solitary waves, yet their phase velocities are much smaller than the large scale rising velocities of the full melt formation.
In Fig. 6b the ratio of maximum fluid velocity (i.e.  ⃗⃗⃗⃗ ) to absolute matrix velocity is shown.For small radii, where   ≫   , this ratio is approximately constant with a high value of about 100.The absolute velocity maxima itself 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 he changed the viscosity ratio as an independent constant model parameter, we change the radius and keep the viscosity law the same, 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 while Scott (1988) is still in the two-phase flow regime.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.  ⃗⃗⃗⃗ =   ⃗⃗⃗ ) and material will be transported by a single phase.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.
Based on these observations, the evolution of these models can be divided into two regimes: (1) In the solitary wave regime ( ≤ 36 ⋅   )   is larger than   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 1a ( < 20 ⋅   ), where the solid mass flux is negative, and 1b (20 ⋅   ≤  < 36 ⋅   ), where the solid moves upwards with the melt.Waves in these regimes are very similar and differ only in the matrix flux.
In the solitary wave composed diapiric uprise regime (2) ( ≥ 36 ⋅   ),   is larger than   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's.
A third regime can be reached by prohibiting the relative movement of fluid ( = ∞ ⋅   ), for which the compaction length has not to be sufficiently resolved.In this regime the initial perturbation will not disintegrate into solitary waves but rise as a well-formed partially molten diapir.In every other case where fluid is able to move w.r.t. the solid, at some point all diapirs will evolve into a swarm of solitary waves which 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 (1)) this will occur sooner and in cases, where solitary waves are much smaller, later.Their observation is mostly limited by resolution.

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 then at some point might 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 diapiric as a 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 for example explainable by local heterogeneities, temperature anomalies for example due to secondary convective overturns in the asthenosphere or by 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 & Podladchikov (1998).The latter authors have looked at the effects of rheology on compaction-driven fluid flow and came to similar results for an upward weakening scenario.The decrease of permeability due to decrease in background porosity might be an alternative explanation.In the hypothetic case of a porosity wave reaching the top of a magma chamber, 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 its amplitude might be high enough to trigger the initiation of dykes.
While we propose that all partially molten diapirs which allow for two-phase flow will inevitably disintegrate into numerous solitary waves when two-phase flow is allowed, it doesn't mean that there

Other issues
The introduced partition coefficients help to distinguish whether solitary wave or diapiric rise is dominant but cannot be solely consulted whether a solitary wave or a diapir can be expected.As the fluid velocity and flux is still very high in the waves 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.
Nevertheless, we propose that all partially molten perturbations allowing for two-phase flow emerge into solitary waves, while those with batch melting rise as diapirs.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 3 grid lengths  (i.e.  ≥ 3 ⋅ ) to get fairly resolved solitary waves, but they can be also observed for much worse resolved compaction lengths.The resolution test (Fig. 2) shows that, even though they are not solved decently, probable solitary waves can be observed for cases with   = .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 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 are the dynamics of diapirism on a length scale of r compared to Darcy flow or solitary waves on the unresolved length scale of   .Thus, two-phase flow, either Darcy flow or solitary waves, becomes negligible for  ≫   and partially molten diapirs can be regarded as well resolved.

Conclusion
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 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, just 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 anomalies size to the compaction length, or rather the models length scale to compaction length ratio, we can classify the ascent behavior into three different regimes using mass flux and velocity of matrix and melt: (1) Solitary wave a and b, (2) solitary wave composite diapirism and (3) diapirism.In regime 1a the matrix sinks with respect to the rising melt, in 1b also the matrix rises, but very slowly.
On first order these regimes can be explained by comparing Stokes velocity of the rising perturbation with the solitary waves phase velocity.If the solitary wave velocity is higher than the Stokes velocity a solitary wave will evolve and, if lower, diapirism is dominant, but still solitary waves will build up if the ascending time is long enough.The deeper we are in regime 2, the more time is needed to build up solitary waves and the less likely it is that they will appear in nature.The third regime can be only reached if fluid is not allowed to move separately to the matrix.
Especially in the second regime numerical resolution plays an important role as the compaction length might be no longer resolved properly 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, do they have a major influence on the conclusions made.

Fig. 2
Fig. 2 shows the resulting models for a length scale ratio of 12 for six different resolutions.The pictures were taken after   has risen approximately 0.25 times the initial Stokes radius ( ′ = 0.25).With increasing resolution, the maximum melt fraction increases strongly from 101 × 101 to 401 × 401 by approximately 20% but the velocity of   decreases by 7% (not shown in the figure).Both values converge.Even though the compaction length is not sufficiently resolved in Fig. 2d, one can still observe t. the model box, and "increasing the initial wave radius" is equivalent to decreasing the compaction length or the size of the emerging solitary waves w.r.t. the model box.In Fig.3the models are shown at  ′ = 0.2.For small radii ( ≤ 12 ⋅   ) 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  = 2.4 ⋅   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  = 12 ⋅   , the resulting wave has just a size of ~20% the initial perturbation size.We can compare the observed rising velocities of these solitary waves of Fig.3b-e with hypothetical Stokes velocities of an equivalent diapir based on equation (15).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.4this 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 bySimpson & Spiegelman (2011), normalized by the Stokes velocity corresponding to different initial perturbation radii.These semi analytical solutions fit quite nicely to our solitary wave models, as already shown inDohmen et al. (2019).The velocities in this figure can be understood as ratios of solitary wave velocity to initial perturbation Stokes velocity.Inspection of Fig.4reveals that for the first four cases of Fig.3b-e with radii smaller or equal 12 ⋅   the phase velocities are always larger than the Stokes velocity.For example, for  = 12 ⋅   , an emerging solitary wave with a typical radius of 4.5 ⋅   has a higher phase velocity than a  = 12 ⋅   melt anomaly rising by Stokes flow.Thus, the cases are always in the solitary wave regime.
passing through the anomalies porosity maximum we define the total vertical mass flux of the rising magma body by ∫ (  +   ) > 0where the integration is carried out only in the region of increased porosity  >  0 .This mass flux is partitioned between the fluid mass flux, are no classical diapirs in nature.Within regime (1) solitary waves are possible and most probably expected but the deeper we are in regime (2) the less expected is the disintegration because they need a long time to build up.In nature compared to the time scale of diapiric rise, different from our models, they cannot rise for an infinite amount of time.The time needed to build up a solitary wave increases non-linearly with  (c.f.Fig.3 l).For example, while for  = 4.8 ⋅   a solitary wave is completely evolved after  ′ = 0.02, for  = 48 ⋅   it needs until  ′ = 0.4, i.e. equivalent to the diapiric rise time necessary to ascend the distance approximately half the initial radii.Our results show that large partially molten bodies with sizes of about  > 20 ⋅   are expected to rise as diapirs but have the potential to split up into a number of solitary waves.Such sizes translate to > 2 km -200 km for typical compaction length within the earth.If such swarms of solitary waves reach the base of the lithosphere each solitary wave may trigger a melt extraction and volcanic event.Given the size of the original diapir, its rising velocity and the number of solitary waves, we might speculate that the episodicity of melt extraction may be related to the time-dependent arrival of the solitary waves.For example, assuming a compaction length of 100 m to1 km, a 3D diapir with  = 40 ⋅   = 40  would possibly split up into several hundred to several 100 thousand solitary waves each having a radius of order 4 ⋅   .With a typical Stokes velocity of a 10 cm/a for the whole body it would release magmas from each solitary wave with a rate of order one per 1 ka to 1 per year.Such extraction rates are in good agreement with observed eruption rates e.g. in Hawaii (e.g.Schmeling, 2006, and references therein).

Fig. 1 :
Fig. 1: The segregation to Stokes velocity ratio, following equation (21), is given as a function of 565 initial perturbation radius  in terms of compaction length   .Each colored line refers to different values of perturbation amplitude   , given in the legend.

Fig. 2 :Fig. 4 :
Fig. 2: All six figures show a model with an initial perturbation radius of 12 times the compaction length but with different resolutions: a) 13x13 b) 26x26 c) 51x51 d) 101x101, e) 201x201, f) 570 401x401.In the lower left corner in each figure the size of the compaction length in terms of grid length is given. 585

Fig. 5 :
Fig. 5: The upper row gives the solid and fluid mass fluxes of a horizontal line cutting through the maximum melt fraction at timesteps where the main wave has just reached the status of a solitary wave.These timesteps are  ′ = .; .; .; ., respectively from left to right.The bottom row gives the corresponding melt porosity fields.All quantities shown are nondimensional.590 ): In that equation ′ scales with ′, which, via equ (12) and (13), scale with   , i.e. with