**Research article**
03 Dec 2019

**Research article** | 03 Dec 2019

# The effect of effective rock viscosity on 2-D magmatic porosity waves

Janik Dohmen, Harro Schmeling, and Jan Philipp Kruse

**Janik Dohmen et al.**Janik Dohmen, Harro Schmeling, and Jan Philipp Kruse

- Institute for Geoscience, Goethe University Frankfurt, Frankfurt, Germany

- Institute for Geoscience, Goethe University Frankfurt, Frankfurt, Germany

**Correspondence**: Janik Dohmen (dohmen@geophysik.uni-frankfurt.de)

**Correspondence**: Janik Dohmen (dohmen@geophysik.uni-frankfurt.de)

Received: 26 Jun 2019 – Discussion started: 08 Jul 2019 – Revised: 24 Oct 2019 – Accepted: 31 Oct 2019 – Published: 03 Dec 2019

In source regions of magmatic systems the temperature is above solidus, and melt ascent is assumed to occur predominantly by two-phase flow, which includes a fluid phase (melt) and a porous deformable matrix. Since McKenzie (1984) introduced equations for two-phase flow, numerous solutions have been studied, one of which predicts the emergence of solitary porosity waves. By now most analytical and numerical solutions for these waves used strongly simplified models for the shear- and bulk viscosity of the matrix, significantly overestimating the viscosity or completely neglecting the porosity dependence of the bulk viscosity. Schmeling et al. (2012) suggested viscosity laws in which the viscosity decreases very rapidly for small melt fractions. They are incorporated into a 2-D finite difference mantle convection code with two-phase flow (FDCON) to study the ascent of solitary porosity waves. The models show that, starting with a Gaussian-shaped wave, they rapidly evolve into a solitary wave with similar shape and a certain amplitude. Despite the strongly weaker rheologies compared to previous viscosity laws, the effects on dispersion curves and wave shape are only moderate as long as the background porosity is fairly small. The models are still in good agreement with semi-analytic solutions which neglect the shear stress term in the melt segregation equation. However, for higher background porosities and wave amplitudes associated with a viscosity decrease of 50 % or more, the phase velocity and the width of the waves are significantly decreased. Our models show that melt ascent by solitary waves is still a viable mechanism even for more realistic matrix viscosities.

- Article
(1802 KB) -
Supplement
(213 KB) - BibTeX
- EndNote

Magmatic phenomena such as volcanic eruptions on the earth's surface show, among others, that melt is able to ascend from partially molten regions in the earth's mantle. The melt initially segregates through the partially molten source region and then ascends through the unmolten lithosphere until it eventually reaches the surface. Within supersolidus source regions at low melt fractions, melt is assumed to slowly percolate by two-phase porous flow within a deforming matrix (McKenzie, 1984; Schmeling, 2000; Bercovici et al., 2001), followed by melt accumulation within rising high-porosity waves (Scott and Stevenson, 1984; Spiegelman, 1993; Wiggins and Spiegelman, 1995; Richard et al., 2012) or focusing into channels which can possibly penetrate into subsolidus regions. Stevenson (1989) carried out a linear stability analysis and found conditions at which flow instabilities may arise, which may result in different 3-D shapes like elongated pockets, channels or porosity waves (Richardson, 1998; Wiggins and Spiegelman, 1995). Formation of 3-D channels within a deforming matrix has been demonstrated in Omlin et al. (2018) or Räss et al. (2014). Here we focus on the supersolidus source region and in particular on the dynamics of porosity waves. An essential parameter controlling the width and phase velocity of porosity waves is the effective shear and bulk matrix viscosity (Simpson and Spiegelman, 2011; Richard et al., 2012). Most of the porosity wave model approaches used either equal bulk and shear viscosities or simple laws in the form of

where *η*_{s} is the effective shear viscosity of the matrix, *η*_{b} the bulk viscosity, *η*_{s0} the intrinsic shear viscosity of the
matrix, *C* a constant of order 1, *φ* the porosity, and *m*=0 for
equal shear and bulk viscosities or *m*=1 otherwise. There are also recent
models that use more complex pressure-dependent weakening viscosities, but
they still use the simplified equations mentioned above for the porosity
dependence of the viscosity (Omlin et al., 2018; Yarushina et al., 2015).
Schmeling et al. (2012) developed an effective viscosity model depending on
a simplified geometry of the fluid phase within a viscous matrix. Possible
melt geometries include flat, ellipsoid-shaped melt inclusions with an
aspect ratio *α* and melt tubes with circular or triangular cross
sections with tapered edges. Comparison of the previous viscosity laws, Eqs. (1)
and (2), with the ones by Schmeling et al. (2012) clearly shows that for
aspect ratio 1, and particularly for smaller *α*, the effective matrix
viscosities are significantly weaker, and disaggregation of the solid occurs
at melt fractions significantly smaller than 100 % as predicted by laws Eqs. (1) and (2). Recent viscosity models based on microscopic diffusion through
grains, grain faces and the melt phase confirm the significance of weakening
with respect to Eqs. (1) and (2) (Rudge, 2018). The aim of this study
is to model 2-D porosity waves with the viscosity laws by Schmeling et al. (2012) and test the influence of the weaker rheology on their shape and
ascent velocity in the absence of melting or freezing.

## 2.1 Governing equations

The mathematical formulation of differential movement between solid matrix
and melt basically builds on that described in Schmeling (2000) and
Schmeling et al. (2019) and is applied here to a porosity wave. We solve the
equations for mass and momentum conservation for melt and solid. 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), and it is valid for infinite Prandtl number (i.e., neglecting
inertia terms in the momentum equations), and small fluid-to-matrix-viscosity ratios. In the following all variables associated with the fluid
(melt) have the subscript *f* and those associated with the solid have the
subscript *s*. Without melting and freezing, the equation for the
conservation of the mass of the melt is

and the mass conservation of the solid is

*v*_{f} and *v*_{s} are the fluid and solid velocities,
respectively. The velocities are derived from the momentum equations, which
is a generalized Darcy equation for the fluid separation flow,

and the Stokes equation for the solid–fluid mixture in the limit of zero fluid viscosity is

*k*_{φ} is the permeability that depends on the rock porosity (i.e.,
melt fraction) with the power *n*,

where *η*_{f} is the dynamic melt viscosity; ** g** is the gravitational
acceleration;

*ρ*is the density of the melt–solid mixture,

*ρ*

_{f}is the density of the melt;

*P*is the fluid pressure, whose gradient is driving the motion; and

*τ*

_{ij}is the effective viscous stress tensor of the matrix.

with the effective shear viscosity *η*_{s} and the effective volumetric
or bulk viscosity *η*_{b} of the porous matrix. The term $\left({\mathit{\eta}}_{\mathrm{b}}-\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\eta}}_{\mathrm{s}}\right)\mathbf{\nabla}\cdot {\mathit{v}}_{\mathrm{s}}$
is often referred to as compaction pressure. The linearized equation of
state for the mixture density is given as

with the density of the matrix *ρ*_{s}.

The fluid pressure in Eqs. (5) and (6) is the same and can be eliminated by merging the two equations. Inserting the density of the mixture, and using Eq. (7), Eq. (5) is recast into

This equation states that the velocity difference between fluid and solid phases (i.e., fluid separation flow or the segregation velocity) is driven by the buoyancy of the fluid with respect to the solid and the viscous stress in the matrix which includes the compaction pressure.

Following Šrámek et al. (2007) the matrix velocity, *v*_{s}, can be written as the sum of the incompressible flow
velocity, *v*_{1}, and the irrotational (compaction) flow velocity,
*v*_{2}, as follows:

with the incompressible velocity potential or stream function *ψ* and
the irrotational (compaction related) velocity potential, *χ*. From
Eq. (11) it follows that the latter is given as the solution of the Poisson
equation

The divergence term **∇**⋅*v*_{s}
can be derived from summing up Eqs. (3) and (4) to give

Equation (12) represents a Poisson equation which can be solved for *χ* once
the melt porosity and segregation velocity are known. As boundary condition
the normal velocity of *v*_{2}, i.e., *v*_{2n}, can be prescribed,
which is equivalent to a normal derivative of *χ*, i.e., a Neuman
boundary condition. If the normal velocity is constant along the boundary,
it automatically fulfills free slip. For the sake of simplicity,
v_{2n}=0 was chosen.

Taking the curl of the matrix momentum (Eq. 6) eliminates the pressure.
Inserting the viscous stress tensor (Eq. 8), the density (Eq. 9) and the
matrix velocity (Eq. 11) into the resulting equation gives the momentum
equation in terms of the stream function *ψ* and the irrotational
velocity potential *χ*

with

The governing equations are non-dimensionalized by the compaction length,
*δ*_{c0} (McKenzie, 1984), and a scaling separation velocity,
*v*_{sc0}, both of which are taken at a reference state which assumes a
constant background porosity *φ*_{0}. The corresponding scaling
viscosities and the scaling permeability are denoted by *η*_{b0}, *η*_{s0} and *k*_{φ0}. The compaction length is
given by

and is the length scale over which a variation in fluid flux gives a response to the compaction. The scaling separation velocity is given as

This defines the scaling law, where the primes denote non-dimensional values and the subscript 0 refers to the background porosity:

The resulting governing equations for the mass are

and for the momentum equations we get

with *e*_{z} as unit vector in the *z* direction (positive upward).

## 2.2 The effective viscosity of a porous matrix

The effective viscosity laws proposed by Schmeling et al. (2012) assume ellipsoidal melt inclusions, or melt films if the inclusions are flat, or melt tubules embedded within an effective viscous medium. This self-consistent assumption is able to predict viscous weakening of a solid matrix with a disaggregation melt porosity on the order of 50 % or less depending on the assumed melt geometry. From their numerical models, Schmeling et al. (2012) derive approximate formulas for the porosity dependence of the dimensional effective matrix shear and bulk viscosities for a melt network geometry consisting of 100 % films,

with ${k}_{\mathrm{1}}={a}_{\mathrm{1}}\left({a}_{\mathrm{2}}+\mathit{\alpha}\left(\mathrm{1}-{a}_{\mathrm{2}}\right)\right)$,
${c}_{\mathrm{1}}=\frac{{b}_{\mathrm{1}}\mathit{\alpha}}{\mathrm{1}+{b}_{\mathrm{2}}{\mathit{\alpha}}^{{k}_{\mathrm{3}}}}$ and
${c}_{\mathrm{2}}=\frac{\mathrm{4}}{\mathrm{3}}\mathit{\alpha}{c}_{\mathrm{1}}^{-{k}_{\mathrm{2}}}\cdot \left({c}_{\mathrm{3}}\left(\mathrm{1}-\mathit{\alpha}\right)+\mathit{\alpha}\right)$, where *a*_{1}=0.97, *a*_{2}=0.8, *b*_{1}=2.2455, *b*_{2}=3.45, *k*_{2}=1.25, *k*_{3}=1.29, *c*_{3}=2
and *α* is the aspect ratio of the ellipsoidal inclusions. At the
disaggregation threshold found as *φ*=*c*_{1}, the partially molten
material loses its cohesiveness and both viscosities approach zero.

For a melt network consisting of 50 % tubes and 50 % films, the following approximate equations have been derived from the model of Schmeling et al. (2012):

The parameters needed to calculate these viscosities for different aspect
ratios between 0.2 and 0.5 are given in Table 1. *k* is given by $k={a}_{\mathrm{1}}\mathit{\phi}+{b}_{\mathrm{1}}$.

Figure 1 shows the effective shear and bulk viscosities for different aspect ratios together with the simplified previous laws (1) and (2).

Takei and Holtzman (2009) and Rudge (2018) suggest that in the presence of
an infinitesimal amount of connected melt the effective viscosity undergoes
a finite drop on the order of a few tens of percents of the intrinsic matrix
viscosity. In our approach we always have a finite melt porosity, and thus we
may identify the zero porosity viscosity *η*_{s0} in our formulation
with the initially weakened value of Takei and Holtzman (2009) or Rudge (2018).

## 2.3 Methods and model setup

For the model we use a square box (1×1), which is initially partially molten
to a certain degree, the background porosity. We place an initial porosity
anomaly with a higher porosity centered at *x*_{0}=0.5 and *z*_{0}=0.2 from
which a porosity wave will develop. As the shape and width of a solitary
wave with a certain rheology law and amplitude is not known a priori, we use
a Gaussian wave of the form

for the perturbation and vary the initial width *w* of the wave.

At the sides of the symmetric box boundaries and at the top and the bottom, free slip boundaries are used. The in- and outflow velocities of matrix and melt at the top and bottom are prescribed in terms of the analytical solution for the background porosity.

The influence of the boundaries on the ascending wave was investigated and found to be fairly small. In Fig. 3 one can see the effect of the upper boundary on the phase velocity. At the end, as the waves approach the upper boundary, the dispersion curves slightly deviate from the supposed line. This error is smaller than 0.5 % as long as the distance from the center of the wave to the upper boundary is greater than 1.5 times its 10 % radius. This radius is defined as the radius at which the porosity has decreased to 10 % of the amplitude of the wave. For the side boundaries this distance has to be larger. For distances greater than 3 times the 10 % radius, this error is smaller than 1 %. In our models the waves have distances of 7–10 times the 10 % radius, which correspond to errors between 0.2 % and 0.05 %.

The equations are solved on a 201×201 grid by finite differences (FDs) using the code FDCON (e.g., Schmeling et al., 2019). Resolution tests have been made with grids varying from 101×101 to 401×401. They show that after a short transient time the phase velocity and amplitude of the evolved porosity wave approach constant values for very high resolutions for all viscosity laws used. The subsequently observed slow variations in the phase velocity and amplitude of the wave along a quasi-steady-state dispersion curve can be attributed to numerical diffusion at finite grid resolution. The resolution test shows that (1) the quasi-steady-state phase velocity and amplitude of the wave are of error order 1, and (2) the dispersion curves obtained on a 201×201 grid overestimate the extrapolated phase velocity values by about 10 %. Time step resolution tests show that the long-term temporal behavior of the porosity waves is significantly improved if the time steps are chosen smaller than approximately 0.2 times the Courant criterion.

The amplitude and phase velocity of the evolving porosity wave are obtained at every time step by quadratic interpolation of the porosity values on the FD grid and determining the value and velocity of the position of the maximum of the quadratic function. The resulting phase velocity shows small oscillations in time, which are probably due to the interaction of the 1st-order error in time when solving Eqs. (3) and (4) and the 2nd-order error of the interpolation. These oscillations are smoothed by applying a moving average including 50 neighboring points. The resulting time series of porosity amplitude and phase velocity can be plotted as a curve with time as curve parameter in an amplitude–phase velocity plot. This curve can be understood as a dispersion curve because the phase velocity depends on amplitude and thus implicitly on the width or wavelength of the porosity wave.

For the model series presented below the width and the amplitude of the
initial wave, the background porosity, and the rheology law have been varied.
All models were carried out using *n*=2 and *n*=3 in the
permeability–porosity law.

## 3.1 Dispersion curves for varied widths and amplitudes

As the shape of a two-dimensional porosity wave for a certain wave amplitude
is not known, the initial width is varied. In Fig. 2a we show a porosity
wave of amplitude 8 initially positioned at *x*=0.5 and *z*=0.2 (left)
as it rises through the model box. In Fig. 2b a horizontal cross section
through the maximum of an initial wave and the resulting solitary wave at a
late stage are shown. During the early stage the wave gains some amplitude as
the volume of an equivalent solitary wave with the same amplitude would be
smaller for this example. Then the amplitude of the ascending wave slowly
decreases again due to numerical diffusion and the evolving phase-velocity–amplitude curve describes the quasi-steady-state dispersion relation. At
this point the wave is expected to be a solitary wave. The shape of this
wave resembles a Gaussian bell-shaped curve quite well but does not fit exactly.
The upper part of the wave in this example fits very well, while the lower
part is slightly wider.

To analyze the evolution of the ascending solitary wave the phase velocity
and the amplitude are tracked over the full rising time and plotted into a
dispersion diagram. In Fig. 3 the dispersion curves of a model with a
starting wave width which is initially larger than the resulting solitary
wave, a model with a similar width and a model with a smaller initial width
are shown. The curves start with high velocities for the Gaussian bell-shaped wave and then rapidly slow down until they approach a specific point
visible as a sharp kink from which they slowly follow a line. For the bigger
and optimal width models, after this kink, the wave is expected to have
reached the solitary wave stage. For the bigger initial width this stage is
reached at a higher amplitude than initially assumed. It is important to
note that, independent of the initial wave width, after reaching a solitary
wave stage the velocities and shapes of waves of a certain amplitude are
always equal, i.e., the three curves merge on one dispersion curve. For
comparison with semi-analytic 2-D solitary porosity wave solutions, the
dashed curves in Fig. 3 and later figures show dispersion curves with
different power law *n* of the permeability–porosity relation and different
bulk viscosity laws, with *m*=0 assuming a constant bulk viscosity and *m*=1
for a 1∕*φ* proportionality (see Eq. 2) (Simpson and Spiegelman,
2011). In contrast to our models, these solutions (a) use a stiff rheology
(“Analytic viscosity” in Fig. 1); (b) neglect solid shear (first term of the
right-hand side of Eq. 8), which is responsible for *v*_{1} (see Eq. 12) in the matrix momentum equation and for an important contribution in
the separation flow (Eq. 11); and (c) apply the small porosity limit.

Based on this result one can carry out many models with different initial wave widths and different initial amplitudes and get one empirical steady-state solitary wave dispersion curve for one viscosity law for a wide range of amplitudes.

Figure 4 shows the time-dependent dispersion curves of models with four different initial amplitudes (4 to 10) and 11 different initial widths each. Depending on the initial widths, they either gain amplitudes as they approach the solitary wave stage or they monotonously loose amplitude. Depending on the initial amplitude and width, each case is characterized by a certain total melt volume, corresponding to a specific steady-state solitary wave with a specific amplitude. Therefore the 44 models finally reach one steady-state solitary wave dispersion curve at different amplitudes. As discussed in Sect. 2, the amplitudes of the waves slowly continue to decrease due to some small amount of numerical diffusion. Yet, they continue following the steady-state solitary wave dispersion curve.

Although we use a different rheology law and do not apply the
simplifications mentioned above, the steady-state dispersion curve of our
model is in general agreement with the *n*=3 and *m*=1 dispersion curve
determined semi-analytically by Simpson and Spiegelman (2011) (Fig. 4,
dashed curve). However, given the 10 % numerical overestimation of phase
velocities of our models (see Sect. 2.2), for high amplitudes our
dispersion curve shows a significantly smaller slope and correspondingly
smaller phase velocities than the semi-analytical curve by Simpson and
Spiegelman (2011). Comparison of the simplified semi-analytical 1-D solution
by Simpson and Spiegelman (2011) with the full analytical 1-D solution by
Yarushina et al. (2015) shows that for low porosities these solutions fit
very well together. For higher porosities the full solution becomes slower
than the simplified one. Tentatively transferring this result to 2-D, our
decrease in the slope can probably be explained by the low-porosity
limitation of the Simpson and Spiegelman (2011) solution, which overestimates
the velocity at high porosities.

## 3.2 Effect of different viscosity laws for *n*=2 and *n*=3 on dispersion curves

To investigate the effect of different viscosity laws, two melt network
geometries are chosen. The first one consists of 50 % films, or ellipsoidal
melt pockets, and 50 % tubes; the second of 100 % films or ellipsoidal melt
pockets. Furthermore, the aspect ratio *α* is varied, whereby a higher
aspect ratio corresponds to compact melt pockets and leads to stronger
viscosities and to a higher disaggregation threshold (see Fig. 1).

Waves with these different viscosity laws give only minor differences in the dispersion curves (Fig. 5a, b). Especially with the films and tubes case, the curves for different aspect ratios (Fig. 5a) are not distinguishable, both during the transient and final stages. In contrast, the analytic viscosity case (Eqs. 1 and 2) propagates along a different path and converges to a 4 %–6 % higher final phase velocity curve. With 100 % films the differences among curves with the different viscosity laws in the final velocity are higher and lie on the order of 6 %. These differences are surprisingly small if compared to the actual differences in effective shear viscosities of about 13 % and bulk viscosity of about a factor of 4 (at 4 % melt corresponding to a porosity amplitude of 8). It is also to be noted that the steady-state part of our dispersion curve calculated with the analytical viscosity (Eqs. 1 and 2) excellently agrees with the semi-analytical solution (dashed) by Simpson and Spiegelman (2011) for the same viscosity law, if we account for the 10 % numerical overestimation of our model phase velocity (see Sect. 2.2). Thus, their neglect of shear stresses and other simplifications have only a very minor effect compared to the effect of different viscosity laws. The overall effect of weakening of matrix viscosity due to decreasing aspect ratio is to slow down the phase velocity slightly.

Changing *n* of the permeability–porosity relation to 2 decreases the wave
velocities significantly (Fig. 5c, d). This drop is consistent with the
simplified semi-analytical solitary wave solutions (*n*=2, *m*=1, dashed
curves). In contrast to the *n*=3 cases, the n=2 velocities are above the
Simpson and Spiegelman (2011) solutions even if the numerical 10 %
overestimation is considered. As for the *n*=3 case, porosity waves with the
stronger analytical viscosity case (Eqs. 1 and 2) are slightly faster than
the new weaker viscosity cases.

While the ascending phase velocity of the wave is only slightly affected by
the different viscosity laws, the width of the wave changes more strongly.
In Fig. 6 the half widths of the solitary waves of amplitude of 8 are
plotted against the corresponding wave velocities for the different
viscosity laws. For *n*=2 (Fig. 6a) and 100 % films the wave gets wider
for higher aspect ratios, while for the mixed geometry the widths stay more
or less constant. The velocity increases only slightly with the aspect
ratio. For *n*=3 (Fig. 6b) and 100 % films the width increases with aspect
ratio, but in contrast to *n*=2 the phase velocity decreases with increasing
aspect ratio. For the mixed geometry the velocity and half width variations
are minor again. These results show that as long as melt tubes represent a
significant portion of the total melt volume (here 50 %), they control the
porosity wave dynamics and keep the porosity wave properties rather fixed.
Only in the absence of tubes do compact melt pockets with large aspect ratios
significantly broaden the waves. For the stiff case of analytical viscosity
(Eqs. 1 and 2), the half width of the wave is comparable to the weaker films with aspect 0.2, but the velocities are larger (Fig. 6a, b, light-brown symbols).

Another interesting phenomenon to observe is the matrix velocity in the
center of the wave, which increases for all geometries with aspect ratio
(Fig. 7). While for 100 % films this increase is stronger, for both
geometries the velocities are approximately equal at aspect ratios between
0.2 and 0.3. For *n*=2 (Fig. 7a), the matrix velocities are always
positive, meaning that despite a slow negative background velocity of the
matrix it rises in the center of the wave (together with the melt).
Interestingly, for *n*=3 (Fig. 7b) and small aspect ratios (0.1 and 0.2,
i.e., weaker effective matrix viscosities), the direction of flow of the
matrix is changed and the matrix in the center flows downwards, i.e., against the
direction of melt flow. Assuming constant matrix shear and bulk viscosities,
Scott (1988) observed a similar switch from negative to positive matrix
velocities in the center of a 2-D solitary wave when the ratio of the bulk-to-shear viscosity was increased from 1 to 9 for *n*=3. We see this switch
around *α*=0.25, corresponding to a bulk-to-shear viscosity at the
center of the porosity wave of about 16, and higher elsewhere. Such a switch
can be explained by an increasing role of diapiric flow, which is
*v*_{1} related, incompressible and moving upward in the center of the wave,
with respect to the compaction flow, which is *v*_{2} related,
irrotational and moving downward in the center of the wave (see Eq. 12).
Weakening of the bulk viscosity within the porosity wave relative to the
shear viscosity allows for stronger decompaction and compaction rates, which
amplify the downward compaction flow with respect to the upward diapiric
flow.

In the previous models the scaling background porosity of 0.005 and maximum
wave amplitudes of 10 to 12 imply maximum melt fractions of 5 % to 6 %.
Thus, the matrix shear viscosity decrease was only small, on the order of 10 %
for the aspect ratio 0.1 models and on the order of 5 % for the stiffer
analytical viscosity laws (1) and (2). This explains the rather mild rheology
effect when comparing the effect of the different viscosity laws. With the
aim to reach higher maximum melt fractions associated with stronger
rheological effects, we carried out a model series with increased background
porosities, both applying the analytical viscosity law (*m*=1) and our
weaker matrix viscosities with 100 % films with an aspect ratio 0.1 (Fig. 8). The increase in the background porosity from 0.5 % to 1.5 % has
only a minor influence on the behavior of the solitary wave for models which
use the analytical viscosity law (*m*=1): the half width of the wave is
almost completely unaffected (by ∼1 %), while the phase
velocity is increased by only approximately 2.5 %. Using a viscosity law
based on a melt geometry consisting of 100 % films and an aspect ratio of
0.1, the differences become significant. The half width decreases to
∼70 % of its initial value, and the phase velocity decreases
by up to 20 % with increasing background porosity, i.e., with an increased
maximum porosity within the wave. Thus, the half widths and phase velocities
show a significant difference to the analytical viscosity law (Fig. 8). In
fact, the phase velocities show the opposite behavior to the analytical
viscosity law (see Fig. 8b). These models suggest that the high melt
fractions within the waves, which are associated with a significant local
matrix weakening, both for shear and bulk viscosity, lead to effectively
shortened compaction lengths within the wave, i.e., to a narrowing and
focusing of the wave. Such narrower waves contain less melt than broader
waves of the same amplitude, i.e., less buoyancy, which slows down the rising
phase velocity.

It is interesting to note that although the semi-analytic solutions of
Simpson and Spiegelman (2011) neglect the shear term in the matrix momentum
equation and in the separation flow equation, they are in good agreement with
the low *φ*_{0} models which include this term. To understand this we
made a test with a model with 100 % films and aspect ratio 0.1 and found
that in the separation flow Eq. (11) the shear term has a significant
amplitude of about 50 % compared to the compaction term. We then switched
off this term in the separation flow Eq. (11), which is equivalent to
assuming zero shear viscosity. Surprisingly, it turned out that separation
velocity changed only insignificantly while the amplitudes of matrix
divergence and convergence increased by about 25 %, and the compaction-related term driving the separation velocity in Eq. (11) increased by
about 50 %, i.e., by the same amount the shear term had before. Obviously,
the buoyancy forces of the solitary wave are partitioned between the
decompaction pressure controlled by the bulk viscosity and the shear
stresses, namely the vertical normal shear stresses. If these stresses are
neglected by assuming a zero shear viscosity, the buoyancy forces are
balanced by the compaction pressure alone, and the shear contribution of the
downward segregation flow is taken over by the increased compaction
contribution.

Recently, Rudge (2018) developed a diffusion creep model based on
microscopic diffusion calculations in the presence of melt in textural
equilibrium with truncated octahedrons. Assuming infinite diffusivity in the
melt phase, Rudge (2018) obtains a somewhat stronger weakening of the shear viscosity
at smaller melt fractions than in our model but comparable disaggregation
porosities as in Fig 1. However, due to the infinite diffusivity assumption,
the bulk viscosity remains finite (equal to 5∕3 of the effective shear viscosity)
even at very small melt porosities, while in our model it increases
infinitely in the limit of zero porosity. We expect that our results with
increased weakening effect (*φ*_{0} increased to 1.5 %) might be
applicable also to the rheology based on the analyses of Rudge (2018).

It should be noted that in our study the viscosity law has been varied by
assuming various melt geometries of melt films and films or melt pockets
superimposed with tubes, while the permeability–porosity relation has been varied
independently between *n*=2 corresponding to the ideal case of only
interconnected tubes and *n*=3 corresponding to the ideal case of
interconnected thin films. Three-dimensional melt distributions of partially
molten mantle rocks have been studied by serial sectioning (Garapić
et al., 2013), identifying a network of melt tubes and films, and by
microtomography (Zhu et al., 2011), suggesting the predominance of melt tubes
along grain edges. Yet, at higher melt fractions the latter distributions
are characterized by tapered edges of the melt tubes partly or completely
wetting grain faces between adjacent grains. From the latter experiments
Miller et al. (2014) determined the permeability by 3-D fluid flow modeling
and found an exponent of 2.6. Thus, our simplified melt viscosities and
permeabilities cover quite well observed partially molten olivine-basalt
systems in textural equilibrium.

In Richard et al. (2012) it was observed that with increasing background
porosities the waves will widen and the phase velocities will slow down. In
our models we observe faster velocities with increasing background porosity
if the analytical viscosity is used. This can be explained by the different
scaling which was used by Richard et al. (2012). They used just the shear
viscosity to calculate the compaction length and not the sum of shear and
bulk viscosity. If the same scaling is used, we get the same behavior for
the phase velocity (Fig. S1b, the Supplement). In contrast to Richard et al. (2012), we observe a narrowing effect of the waves for larger background
porosities, which cannot be explained by scaling (Fig. S1a). As
Richard et al. (2012) used a 1-D model, we suspect that 2-D effects such as
including the incompressible flow velocity, *v*_{1}, are responsible
for the different shapes of the wave at different background porosities.

As the shape of a solitary wave in our models cannot be described analytically, we start with a Gaussian wave which develops quite rapidly into a solitary wave with a similar shape and a certain amplitude, depending on the initial width of the wave.

Even though the rheologies used are much weaker than the simplified analytical ones, the effect on dispersion curves and wave shape are only moderate as long as the shear viscosity does not drop below about 80 % of the intrinsic shear viscosity. This value corresponds to a melt fraction of 5 %, equivalent to 20 % of the disaggregation value. At this porosity the bulk viscosity is approximately 5–7 times the intrinsic shear viscosity. In this case the phase velocity changes just slightly for all cases, while the waves broaden in the absence of tubes with increasing aspect ratio.

In contrast, for higher melt fractions of about 12 %, equivalent to 50 % of the disaggregation values, the shear viscosity decreases to 50 % of the intrinsic viscosity, and the bulk viscosities are on the order of the intrinsic shear viscosity. Then, our models predict significant narrowing of the porosity waves and slowing down of the phase velocities. For such conditions a strong discrepancy in solitary wave behavior between our viscosity law and the analytical ones is found.

For low melt fractions our models are in good agreement with semi-analytic solutions which neglect the shear stress term, because the matrix shear contribution of the downward segregation flow is taken over by the increase in the compaction contribution.

The used code and all model results, namely the melt porosity, segregation velocity and matrix velocity fields, are available on request.

The supplement related to this article is available online at: https://doi.org/10.5194/se-10-2103-2019-supplement.

The idea for this project came from HS. The models were carried out by JD. JPK carried out the calculations for the viscosity laws. JD and HS prepared the article. All authors have contributed to the discussion and article review.

The authors declare that they have no conflict of interest.

We would like to thank the reviewers Guillaume Richard and Viktoriya Yarushina for their detailed and thoughtful reviews, which helped to improve the article.

This paper was edited by Susanne Buiter and reviewed by Guillaume Richard and Viktoriya Yarushina.

Bercovici, D., Ricard, Y., and Schubert G., A two phase model for compaction and damage, 1: General theory, J. Geophys. Res., 106 , 8887–8906, 2001.

Garapić, G., Faul, U. H., and Brisson, E.: High resolution imaging of the melt distribution in 1 partially molten upper mantle rocks: evidence for 2 wetted two-grain boundaries, G-Cubed, 14, 556–566 https://doi.org/10.1029/2012GC004547, 2013.

McKenzie, D.: The generation and compaction of partially molten rock, J. Petr., 25, 713–765, 1984.

Miller, K. J., Zhu, W. I., Montési, L. G. J., and Geatani, G. A.: Experimental quantification of permeability of partially molten mantle rock Earth Planet. Sc. Lett., 388, 273–282, 2014.

Omlin, S., Räss, L., and Podladchikov, Y. Y.: Simulation of three-dimensional viscoelastic deformation coupled to porous fluid flow, Tectonophysics, 746, 695–701, 2018.

Räss, L., Yarushina, V. M., Simon, N. S., and Podladchikov, Y. Y.:
Chimneys, channels, pathway flow or water conducting features-an explanation
from numerical modelling and implications for CO_{2} storage, Energy Proced., 63, 3761–3774,
2014.

Richard, G. C., Kanjilal, S., and Schmeling, H.: Solitary-waves in geophysical two-phase viscous media: a semi-analytical solution, Phys. Earth Planet. Int., 198/199, 61–66, 2012.

Richardson, C. N.: Melt flow in a variable viscosity matrix, Geophys. Res. Lett., 25, 1099–1102, 1998.

Rudge, J. F.: The viscosities of partially molten materials undergoing diffusion creep, J. Geophy. Res., 123, 10534–10562, https://doi.org/10.1029/2018JB016530, 2018.

Schmeling, H.: Partial melting and melt segregation in a convecting mantle, in: Physics and Chemistry of Partially Molten Rocks, edited by: Bagdassarov, N., Laporte, D., and Thompson, A. B., Kluwer Academic Publ., Dordrecht, 141–178, 2000.

Schmeling, H., Marquart, G., Weinberg, R., and Wallner, H.: Modelling melting and melt segregation by two-phase flow: new insights into the dynamics of magmatic systems in the continental crust, Geophys. J. Int., 217, 422–450, 2019.

Schmeling, H., Kruse, J. P., and Richard, G.: Effective shear and bulk viscosity of partially molten rock based on elastic moduli theory of a fluid filled poroelastic medium, Geophys. J. Int., 190, 1571–1578, https://doi.org/10.1111/j.1365-246X.2012.05596.x, 2012.

Scott, D. R.: The competition between percolation and circulation in a deformable porous medium, J. Geophys. Res., 93, 6451–6462, 1988.

Scott, D. R. and Stevenson, D. J.: Magma solitons, Geophys. Res. Lett., 11, 1161–1164, 1984.

Simpson, G. and Spiegelman, M.: Solitary wave benchmarks in magma dynamics, J. Sci. Comput., 49, 268–290, https://doi.org/10.1007/s10915-011-9461-y, 2011.

Spiegelman, M.: Physics of melt extrraction: theory, implications and applications, Phil. Trans. R. Soc. Lond. A, 342, 23–41, 1993.

Spiegelman, M. and McKenzie, D. : Simple 2-D models for melt extraction at mid-oceanic ridges and island arcs, Earth Planet. Sc. Lett., 83, 137–152, https://doi.org/10.1016/0012-821X(87)90057-4, 1987.

Šrámek, O., Ricard, Y., and Bercovici, D.: Simultaneous melting and compaction in deformable two-phase media, Geophys. J. Int., 168, 964–982, https://doi.org/10.1111/j.1365-246X.2006.03269.x, 2007.

Stevenson, D. J.: Spontaneous small-scale segregation in partial melts undergoing deformation, Geophys. Res. Lett., 9, 1064–1070, 1989.

Takei, Y. and Holtzman, B. K.: Viscous constitutive relations of solid-liquid composites in terms of grain boundary contiguity: 2. Compositional model for small melt fractions, J. Geophys. Res., 114, B06206, https://doi.org/10.1029/2008JB005851, 2009.

Wiggins, C. and Spiegelman, M.: Magma migration and magmatic solitary waves in 3-D, Geophys. Res. Lett., 22, 1289–1292, 1995.

Yarushina, V. M., Podladchikov, Y. Y., and Connolly, J. A.: (De) compaction of porous viscoelastoplastic media: Solitary porosity waves, J. Geophys. Res.-Sol. Ea., 120, 4843–4862, 2015.

Zhu, W., Gaetani, G. A., Fusseis, F., Montesi, L. G. J., and De Carlo, F.: Microtomography of partially molten rocks: Three-dimensional melt distribution in mantle peridotite, Science, 332, 88–91, 2011.