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

In source regions of magmatic systems the temperature is above solidus and melt ascent is assumed to occur 5 predominantly by two-phase flow which includes a fluid phase (melt) and a porous deformable matrix. Since McKenzie (1984) introduced his 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 shearand bulk viscosity of the matrix, significantly overestimating the viscosity or completely neglecting the porositydependence of the bulk viscosity. Schmeling et al. (2012) suggested viscosity laws in which the viscosity decreases very 10 rapidly for small melt fractions. They are incorporated into a 2D 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 effect 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 15 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.

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 55 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 is the effective shear viscosity of the matrix, the bulk viscosity, 0 the intrinsic shear 60 viscosity of the matrix, a constant of order 1, the porosity, and = 0 for equal shear and bulk viscosities or = 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 65 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, (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 70 (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 equations (1) and (2) (Rudge, 2018).
The aim of this study is to model 2D-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. 75

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 80 governing equations for the melt-in-solid two-phase flow dynamics is based on McKenzie (1984), Spiegelman & McKenzie (1987) and Schmeling (2000) and is valid for infinite Prandtl number (i.e. The fluid pressure in equation (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 phase (i.e. fluid separation flow, or the segregation velocity) is driven by the buoyancy of the fluid with respect to the solid, and the 110 viscous stress in the matrix which includes the compaction pressure.
Following Šrámek et al. (2007) the matrix velocity, , can be written as the sum of the incompressible flow velocity, 1 , and the irrotational (compaction) flow velocity, 2 , as follows: with the incompressible velocity potential or stream function and the irrotational (compaction related) 115 velocity potential, . From eq. (11) it follows that the latter is given as the solution of the Poisson equation The divergence term ∇ ⃗ ⃗ • can be derived from summing up eq. (3) and eq. (4) to give 120 Eq. (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 2 , i.e. 2 , 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 sake of simplicity v 2n = 0 was chosen.
Taking the curl of the matrix momentum eq. (6) eliminates the pressure. Inserting the viscous stress 125 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 The governing equations are non-dimensionalized by the compaction length, 0 , (McKenzie, 1984) and a scaling separation velocity, 0 , 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 0 , 0 , and 0 , respectively. The compaction length is given by and is the length scale over which a variation in fluid flux gives a response on the compaction. The scaling separation velocity is given as and for the momentum equations we get with as unit vector in z-direction (positive upward).

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 155 melt porosity of 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 effective matrix shear and bulk viscosities for a melt network geometry consisting of 100% films for < 1 (23) 160 with 1 = 1 ( 2 + (1 − 2 )), 1 = 1 1+ 2 3 , 2 = 4 3 1 − 2 • ( 3 (1 − ) + ) where 1 = 0.97, 2 = 0.8, 1 = 2.2455, 2 = 3.45, 2 = 1.25, 3 = 1.29, 3 = 2, and is the aspect ratio of the ellipsoidal inclusions. At the disaggregation threshold found as = 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 165 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 Tab Figure 1 shows the effective shear and bulk viscosities for different aspect ratios together with the simplified previous laws (1) and (2). 175 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 of the order of a few 10% of the intrinsic matrix viscosity. In our approach we always have a finite melt porosity, thus we may identify the zero porosity viscosity 0 in our formulation with the initially weakened value of Takei and Holtzman (2009) or Rudge (2018). 180 (1) and (2)).

Methods and model setup
For the model we use a square box (1x1), which is initially partially molten to a certain degree, the background porosity. We place an initial porosity anomaly with a higher porosity centered at 0 = 0.5 and 0 = 0.2 from which a porosity wave will develop. As the shape and width of a solitary wave with a 190 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 of the wave. The influence of the boundaries on the ascending wave was investigated and found to be fairly small. In 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 of the center of the wave to the upper boundary is greater 200 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 corresponds to errors between 0.2 and 0.05%.  (3) and (4) and the 2 nd order error of the interpolation. These oscillations are smoothed by 220 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 amplitudephase 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 225 porosity and the rheology law have been varied. All models were carried out using n=2 and n=3 in the permeability-porosity law.

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 230 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 is 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 235 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 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 is 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 250 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 255 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 (c.f. equ. 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 equ. 8) which is 260 responsible for 1 (c.f. equ. 12) in the matrix momentum equation, and for an important contribution in the separation flow (equ. 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. Fig. 4 shows the time-dependent dispersion curves of models with 4 different initial amplitudes (4 to 10), 270 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 section 2, the 275 amplitude 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, m=1 dispersion curve determined semi-analytically by Simpson and Spiegelman (2011) (Fig. 4, dashed curve). However, given 280 the 10% numerical overestimation of phase velocities of our models (c.f. section 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 semianalytical 1-D solution of Simpson and Spiegelman (2011) with the full analytical 1-D solution of Yarushina et al. (2015) shows that for low porosities these solutions fit very well together. For higher 285 porosities the full solution becomes slower than the simplified one. Tentatively transferring this result to 2D 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.

Fig. 4. Dispersion curves for 44 models with 4 different initial amplitudes (4 to 10) and 11 different 290
initial widths each. All models were carried out for a melt network geometry consisting of 100% films and an aspect ratio of 0.1. The background porosity is 0.005 and n = 3.

Effect of different viscosity laws for n=2 and 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/ellipsoidal melt pockets and 50% tubes, the second of 100% films/ellipsoidal 295 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 (c.f. Fig.   1).
Waves with these different viscosity laws give only minor differences in the dispersion curves (Fig. 5a, b).
Especially with the films & tubes case the curves for different aspect ratios (Fig. 5a) are not 300 distinguishable, both during the transient and final stage. In contrast, the analytic viscosity case (equ. 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 in 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 4 (at 4 % melt 305 corresponding to a porosity amplitude 8). It is also to be noted that the steady state part of our dispersion curve calculated with the analytical viscosity (eq. 1 and 2) excellently agrees with the semianalytical 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 (c.f. section 2.2). Thus, their neglect of shear stresses and other simplifications have only a very minor effect compared to the effect of 310 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) 315 solutions even if the numerical 10 % overestimation is considered. As for the n=3 case, porosity waves with the stronger analytical viscosity case (equ. 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 viscositiy laws, the width of the wave changes more strongly. In Figure 6 the half-widths of the solitary waves of 320 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 325 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 compact melt pockets with large aspect ratios significantly broaden the waves. For the stiff case of analytical viscosity (equ. 1 and 2) the half width of the wave is comparable to the weaker 0.2 films, but the velocities are larger (Fig. 6a,b, light brown symbols). 330 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 an aspect ratio between 0.2 and 0.3. For n = 2 shear viscosity was increased from 1 to 9 for n = 3. We see this switch around  = 0.25 corresponding to 340 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 1 -related, incompressible, and upward in the center of the wave, with respect to the compaction flow, which is 2 -related, irrotational, and downward in the center of the wave (c.f. equ. 12). Weakening of the bulk viscosity within the porosity wave relative to the shear viscosity allows stronger decompaction and compaction rates which 345 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, of order 10% for e.g. the aspect ratio 0.1 models and of order 5% for the stiffer analytical viscosity laws (1), (2). This explains the rather mild rheology effect when comparing the effect of the different viscosity 360 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 365 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 370 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 same amplitude, i.e. less buoyancy, which 375 slows down the rising phase velocity.

Discussion
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 385 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 equation (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 equation (11), which is equivalent to assuming zero shear viscosity. Surprisingly it turned out that separation velocity changed only insignificantly while the 390 amplitudes of matrix divergence and convergence increases by about 25%, and the compaction related term driving the separation velocity in equation (11) increases 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 395 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 he obtains a somewhat stronger weakening of the shear viscosity at smaller melt 400 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 (=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 Rudge's (2018) analyses. 405 It should be noted that in our study the viscosity law has been varied by assuming various melt geometries of melt films and films/melt pockets superimposed with tubes, while the permeabilityporosity has been varied independently between = 2 corresponding to the ideal case of only interconnected tubes and = 3 corresponding to the ideal case of interconnected thin films. Threedimensional melt distributions of partially molten mantle rocks have been studied e.g. by serial 410 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 3D-fluid flow modelling and found an exponent of 2.6. Thus, our 415 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 420 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. 1b, Suppl. Mat.). 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. 1a, Suppl. Mat.). As Richard et al. (2012) used a 1-D model, we suspect that 2-D effects such as including the 425 incompressible flow velocity, 1 , are responsible for the different shapes of the wave at different background porosities.

Conclusion
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, 430 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 435 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 is of the order of the intrinsic shear viscosity. Then, our models predict significant narrowing of the porosity waves and 440 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 of the compaction contribution. 445