- Email: [email protected]
A mechanism for the lateral transport of gas bubbles in silicic lava rising in a vertical conduit
EPSL ELSEVIER
Earth and Planetary
Science Letters
I56 ( 1998)
13-18
A mechanism for the lateral transport of gas bubbles in silicic lava rising in a vertical conduit S.D.R. Wilson Department
of Mathematics.
University
Manchester
Received
13 June 1997: revised version received
Ml3
of Manchester:
0,rford
Road.
9PL. UK
10 December
1997: accepted
26 December
1997
Abstract It is well-known that the shear viscosity of silicate-rich melts is very sensitive to water content. Thus, as such, a material rises towards the earth’s surface in an eruption, the pressure falls, water exsolves to form bubbles or voids, and the melt viscosity increases. This is responsible for the (by now) well-known large departures of the pressure gradient from hydrostatic. The exsolved gas is known to escape from the rising lava into the surrounding rock and since gas content is believed to have an important role in the transitions between eruption regimes, it is of some interest to understand the mechanism of escape. A key difficulty is to explain how small gas bubbles could migrate laterally across the conduit, in view of their small size and the large liquid viscosity. It is shown here that the increase of melt viscosity with height produces a horizontal pressure gradient directed so that the pressure at the centre of the conduit is larger than at the walls. Thus the fluid in the centre is less viscous and more dense; this configuration is (very probably) dynamically unstable. The resulting instability will have the effect of mixing the lava, thus tending to transport gas bubbles from the centre to the (depleted) regions near the walls. 0 1998 Elsevier Science B.V. All rights reserved. Keywords:
melts: lava; flow mechanism; viscosity; density; gases
1. Introduction The rise of magma through a vertical conduit has
been the subject of much theoretical development recently, and we may cite the papers by Dobran [I], Ramos [2], and Jaupart [3] as typical. We refer more particularly to Jaupart and Allegre [4]. Here, a study was made of the rise of silicic magmas and a principal objective was to understand and explain the transitions between explosive and effusive eruption regimes. Attention was drawn to the important role of gas content, that is bubble volume fraction, in this process. It appears that gas is transferred from the lava to the country rock at the conduit wall and then 0012-821X/98/$19.00 0 1998 Elsevier Science B.V. All rights reserved. PI/ SOOl2-821X(98)00019-3
migrates via pores and cracks, and so it is important to understand this process too. In the paper just cited, it is assumed that the gas escape is controlled by conditions at the conduit wall, in fact, by the pressure drop between the lava and the country rock at some distance. This model therefore incorporates the assumption that there is essentially no barrier to gas transport within the liquid phase itself; the discussion in Jaupart and Allegre [4] leaves this problem open. At present, there appears to be no acceptable explanation of how gas might be transported rapidly across the conduit; it seems inconceivable that gas bubbles of millimetre dimensions could migrate significantly relative to
14
S.D.R. Wilson/Earth and Planetarv Science Letters 156 (1998) 13-18
the liquid, whose viscosity is at least lo5 Pa s, and the permeable foam model of Eichelberger et al. [5] seems unrealistic as pointed out by Jaupart and Allegre [4] _ The present paper offers a contribution to the solution of the problem of gas transport within the conduit. Attention is drawn to the possibility that the basic vertical flow of lava may be subject to a dynamic instability, which, if present, would lead to large scale (i.e. on the scale of the conduit width) overturning and mixing, thus transporting liquid and gas bubbles together from the centre to the walls. In this introduction, we describe the mechanism in qualitative terms. The detailed calculations, for a somewhat idealised model, are given in Section 2 followed by numerical examples. A key feature of silicic magmas is the well-known extreme sensitivity of the value of the viscosity of the melt to changes in the water content. Roughly speaking, the viscosity increases by a factor 10 (or more) for each 1% decrease in water mass fraction. As the lava rises, the pressure falls, exsolved water forms bubbles of water vapour, and the water content falls. Thus the viscosity of the lava depends on the pressure, effectively, and so on the depth. The rapid increase in viscosity as the pressure falls towards the exit value is responsible for the (by now) well-known large departures of the pressure gradient from the hydrostatic value, described in [l--3]. However, there is a second effect, which is small by comparison, but nonetheless important. The viscosity variation with depth creates viscous stresses which set up a horizontal pressure gradient (i.e. transverse to the conduit) directed so that the pressure in the centre is greater than at the walls. (This effect is reinforced by a second, smaller, effect associated with the overall volume expansion of the lava.) This means, in turn, that the fluid in the centre has lower viscosity and higher density than the fluid at the walls. This configuration, or at any rate one very similar to it, is known to be dynamically unstable. There is a very large body of literature on this subject because of the industrial importance of the so-called ‘pipe-lining effect’, and it has recently been surveyed by Joseph et al. [6]. Its possible role in lava flows has been pointed out by Carrigan and Eichelberger [7], who discussed the problem of the flow in a conduit
of a combination of two varieties of lava (mafic and silicic) distinguished by composition and therefore by viscosity. In the present work, we consider only a single material, the point being, of course, that the unstable vertical stratification of viscosity and density is produced automatically by the flow itself. We refer particularly to a paper by Renardy [8], who studied a system consisting of two distinct liquids symmetrically disposed in three layers in sandwich fashion. Renardy showed that when the viscosity and density contrast are as in the present case, the flow is always unstable. The most unstable disturbance mode is the ‘snake’ mode in which the core fluid moves as a whole from one wall to the other, in zig-zag fashion, with amplitude which increases exponentially with time. The growth rates are probably large enough to be significant in the present case. It is worth noting that nearly all of the many articles on this subject are confined to calculation of the neutral curves, that is the boundary in parameter space which marks the transition from stability to instability. Renardy [8], however, provides calculation of the growth rates in the unstable region which is particularly useful here.
2. Theory We propose to regard the mixture of liquid and gas bubbles as a single continuum, with suitably averaged properties, modelled as a (compressible) Newtonian liquid. (This is implied in the theories developed in [ 141.) Thus, no account is taken of the presence of crystals. The constitutive equation is O;j
=
-pSij
+
tij
(1)
where aij is the mechanical stress tensor, p is the thermodynamic pressure and tij is the deviatoric stress tensor, given by T;j =
(K
-
$/A)ekkhij
+
2@ij.
(2)
Here eii is the rate of strain tensor and p and are the shear and bulk viscosity coefficients. (For a derivation of this equation see Hunter [9], for example.) The pressure p is linked to the local temperature and density by an equation of state. We can derive K
S.D.R. Wilson/Earth
and Planetary
this from a combination of an equation expressing the law of mass conservation for water with a solubility law. We are going to assume, for simplicity, that water is the only volatile present and that bubbles consist of water vapour which can be regarded as a perfect gas. Consider a volume of the mixture comprising a volume V, of liquid of density pu and a volume V, of gas of density ps. We assume that as the mixture expands, Vt and pe remain constant. Let (T denote the mass fraction of water in the liquid phase. Then we have pr Vpa + pg V, = constant = mass of water
(3)
(4)
P
where p is the average density of the mixture. The perfect gas law is
in the usual notation. Combining these we find
M
( >
Pea+
or
Pe-P
---p
= constant.
P
o = @‘I2
(7)
(with k = 4.11 x lO-‘j), quoted by Jaupart and Allegre [4]. Eliminating u gives an equation connecting p. p and T, which is the required equation of state: JPrPP)
constant. P
the forms of the components of tij are (9)
(10) (11) The equations of motion, neglecting inertial effects, are aP
-ax =
a
-PS
+ ~“‘X
(8)
Inspection of the magnitudes in this equation shows that the second term is negligible compared with the first, unless p is extremely close to pt. Thus, for void fractions of 2% or more, we can delete the second term and this will be done in most of what follows. (However the accurate expression is also needed, as will be seen.) Next, we turn to the equations of motion. Let us consider the two-dimensional problem of flow between parallel vertical walls, with x as the vertical coordinate, measured upwards, and y horizontal, and U, u as the corresponding velocity components. Then
a
+ Yj-j?
(12) (13)
Eq. 13 shows the origin of the transverse pressure gradient @lay. For the stress component rXvvaries with x, because p does; and the stress component rvYis not zero because au/& is not zero - an effect of the general expansion of the liquid. The system is completed by the overall continuity (mass conservation) equation, assuming steady flow: (14)
(6)
A suitable solubility law is
1.5
13-18
a ap a - -t,” + -z,,,. ay -. G-ax -.
and also V P--t-P _Vf
Science Letters 156 (1998)
We now use Eq. 13 to find a formula for the pressure difference between the centre of the channel and the wall. This can be done conveniently if we accept an answer in terms of the local values of the pressure p and the vertical pressure gradient ap/ax. The pressure p can be regarded as a rough measure of depth. Of course the model can be used to find p as a function of x, but this is rather long and complicated and inessential to the present purpose. Let us assume (and a detailed analysis confirms this) that the transverse velocity component IJ can be set equal to zero as a first approximation. Then, using the simplified forms of r_rvand rvv, we can integrate Eq. 13 with respect to y, giving the pressure contrast Ap as
(15) where urnaxis the velocity on the centre line of the channel (about 3/2 times the average). We have also, from Eq. 14, PU max =
constant
(16)
16
S.D.R. Wilson/Earth
and Planetary
and the simplified version of Eq. 8, namely
P(Pe- P)
(17)
and as we see in Section 3 these are sufficient to express Ap in terms of p and LIplax.
3. Numerical examples First, we turn to the dependence of viscosity on water content. As an example, we take the data for rhyolite given in Lange [lo]. The data points up to 0 = 4% can be reasonably well fitted by an exponential of the form =
plexp(-ba)
(18)
with pt = 1O’o.9= 7.9 x lo9 Pa s and b = 357. Using the solubility law (Eq. 7) gives
P
(19)
Here p.e is the shear viscosity of the liquid phase. The viscosity p which appears in Eq. 15 is, of course, the shear viscosity of the mixture. To keep matters simple we shall set /.L= /..Q,thus disregarding, for the moment, the effect of bubbles on p. We also need an estimate of the bulk viscosity K and since this is highly sensitive to the bubble volume fraction at small voidage we need to account for this. An approximate formula is given by Prud’Homme and Bird [ 111 and can be written K
=
P $LcLc----. Pe -
(20)
P
(The singularity when p = Pe is a false alarm; the stresses produced are finite. See Prud’home and Bird
[ill). In order to use the equation of state (Eq. 17) we need to evaluate the constant on the right and this depends on the original state of the lava. Let us suppose that the original water mass fraction (i.e. at great depth) is 3%. Using Eq. 8, we find that when the pressure is reduced sufficiently to create bubbles occupying 5% of the volume (i.e. p = 0.95~~) the pressure is 464 bar and the water content is 2.8%. Thereafter, as the fluid rises, we can
0.05 = constant = 464 x lo5 X -
0.95
= 2.44 x lo6 Pa.
(21)
(We have set pe = 2400 kg/m”, T = 1120 K in accordance with the data from Lange [lo] and assume that conditions are isothermal.) Eq. 21 can now be used to calculate the density in terms of p. We are now able to express Ap, from Eq. 15, in terms of p, p and dp/ti. Note that dp -=-dx
dp dp
dp dx and, from Eqs. 16 and 17 du - max = ---umax dP dx PdJ
= _U,,,~_(Pe - P) dp PeP dx’
(22)
(23)
Then Ap =
p,e = plexp(-kbp’/*).
13-18
use Eq. 17 which now reads P(Pe - P)
= constant
P
Fe
Science Letters 156 (1998)
dp - (K + ;,)@&@ dp
S&,&X. (24) I Note that dp/dp is negative and so is dpldx, so that Ap is the sum of two positive terms. Any values of p and dp/dx thought to be of interest can now be inserted into Eq. 24, which is completely explicit, apart from the necessity to use Eq. 21 to calculate p in terms of p. As an example, suppose we take p = 50 bar = 5 x lo6 Pa and dp/dx = 3 times the hydrostatic value = 3pg; and put u,,, = 3 mm SC’, corresponding to an average rise velocity of 2 mm s-‘. Estimates of these magnitudes are necessary rough and indirect; these values are in line with a number of sources, for example Stasiuk et al. [12], Eichelberger et al. [5], DeSilva et al [13] and Jaupart [3]. Then we find Ap = 2.3 bar approximately, i.e. a pressure contrast of only about 4%, but, using Eq. 19, creating a viscosity contrast of 7%, easily enough to produce a dynamic instability. To check this last point, we need to refer to the results of Renardy [S] in some detail. Renardy uses the label m to denote the viscosity ratio, which thus takes the value 1.08 in the example just worked out. G is the label given to the ratio of the pressure gradient to the hydrostatic value, taking the value 3 here. Renardy gives, in Fig. 2, the growth rate of long waves in the case m = 1.1, G = 3 (and sev-
S.D.R. Wilson/Earth and Planetay Science Letters 156 (1998) 13-18
era1 other cases). (The numerical coincidence here was deliberately contrived.) The growth rate is very small, about 10v7 in the chosen dimensionless units; but Renardy’s Fig. 5 shows that the growth rate of short waves is much larger. Only the case m = 2 is shown; the growth rates are larger than the long wave values by a factor of about 104. Assuming that the case m = 1.1 behaves in roughly the same way, we obtain a maximum growth rate in the range lo-“- 10e4. The unit of time in these computations is ,/@ where C is the channel (conduit) half width. Rutting f? = lOm, the conclusion is that the time scale of the growth of the instability is about 103lo4 s, say 1 h. The most unstable waves have a wavelength of the same order as the conduit width. The depth of the bubbly region can be calculated to be on the order of hundreds of metres, using Eq. 21 and assuming that the mixture remains cohesive up to about 50% voidage. This gives a residence time, at 3 mm/s, of about 40-50 h, much longer than the growth time. Note that a lower value of urnax will result in slower growth, but longer residence time, the two effects partially offsetting each other. (A unit of time based on viscosity, e.g. pf?/p, might seem more natural. This differs from Renardy’s actual choice by a factor of the Reynolds number. It would be natural, for R < I, to rescale the time as indicated; then the dimensionless growth rates would differ by the corresponding factor. But the physical time scale would, of course, come out the same.) These numbers are quite typical, as a little trial and error will show. Very similar results are obtained with umaX= 20 mm-‘, p = 100 bar, for example. The choices made were, as noted, aimed at enabling a reasonable comparison with Renardy’s computations. Furthermore, the pressure gradient may well exceed the lithostatic value by more than a factor 3. The general conclusion is that the instability is most likely to occur in the upper sections of the conduit where there are marked departures from hydrostatic equilibrium.
4. Concluding
remarks
Attention has been drawn to a possible hydrodynamic instability of vertical conduit flow of silicic
17
magmas. This results from a horizontal pressure gradient, stemming mainly from the variation of the shear viscosity with height. The mechanism does not depend on Newtonian behaviour of the melt; that assumption serves only to make possible simple explicit calculations. A somewhat crude model for the dependence of the melt viscosity pe on pressure has been developed (Eq. 19). Again, this was used simply to make the calculations explicit; any model similarly sensitive would give similar results. It would be of considerable interest to refine Eq. 19 so as to include the direct dependence of the liquid phase viscosity on pressure (i.e. at fixed composition), an effect discussed by Lange [lo]. Also of possible importance is the effect of void fraction on the viscosity of the two-phase mixture. Here, there is conflicting evidence. Jaupart and Allegre [4] and Dobran [ 11, among others, make use of models according to which the viscosity increases with void fraction, presumably resulting from ‘bubble interference’. To the present author, this seems unlikely to hold in the relevant parameter regime, where surface and inertial effects are completely negligible. Experiments on porous glasses by Rahaman et al. [ 141 and Ducamp and Raj [ 151 show a decrease in viscosity with void fraction, and these appear to be more relevant. The application of Renardy’s analysis [8] to the present case perhaps needs some caution. Renardy studied a system of two contrasting immiscible liquids, whereas we have one stratified liquid. The difference lies not in interfacial effects, which do not matter much in the cases cited, but in the fact that in the two-liquid case, the shear viscosity is a material constant convected with the liquid, whereas here it depends on the local pressure. The stability of a fluid whose viscosity varies continuously (e.g. as a result of temperature variations) has also been studied; see, for example, Wall and Wilson [16]. But only the neutral curves are given, not the growth rates in the unstable region. An extension of Renardy’s work to something like the present case would not be practicable because in addition to the feature just mentioned, the basic flow exhibits a rapidly varying pressure gradient, and the normal mode decomposition, leading to ordinary differential equations for the disturbance variables, is not possible. [FA]
18
S.D.R. Wilson/Earth
and Planetarv
References [I] F. Dobran, Non-equilibrium
[2]
[3] [4]
[5]
[6] [7] [S] [9]
flow in volcanic conduits and applications to the eruption of Mt. St. Helens on May 18 1980 and Vesuvius in AD79, J. Volcanol. Geotherm. Res. 49 (1992) 285-311. J.1. Ramos, One-dimensional, time-dependent, homogeneous 2-phase flow in volcanic conduits, Int. J. Numer. Methods Fluids 21 (3) (1995) 253-278. C. Jaupart, Physical models of volcanic eruptions, Chem. Geol. 128 (1996) 217-227. C. Jaupart, C.J. Allegre, Gas content, eruption rate and instabilities of eruption regime in silicic volcanoes, Earth Planet. Sci. Lett. 102 (1991) 413-429. J.C. Eichelberger, C.R. Carrigan, H.R. Westrich, R.H. Price, Non-explosive silicic volcanism, Nature 323 (1986) 598602. D.D. Joseph, R. Bai, K.P. Chen, Y.Y. Renardy, Core-annular flows, Ann. Rev. Fluid Mech. 29 (1997) 65-90. C.R. Carrigan, J.C. Eichelberger, Zoning of magmas by viscosity in volcanic conduits, Nature 343 (1990) 248-25 I. Y. Renardy, Viscosity and density stratification in vertical Poiseuille flow, Phys. Fluids 30 (1987) 1638-1648. S.C. Hunter, Mechanics of Continuous Media, Ellis Hor-
Science Letters
156 (1998)
13-18
wood, 1976, Chichester, pp. 135 and 149. [IO] R.A. Lange, The effect of H20, CO1 and F on the density and viscosity of silicate melts, in: M.R. &roll, J.R. Holloway (Eds.), Volatiles in Magma& Rev. Mineral. 30, 1994. [I I] R.K. Prud’Homme, R.B. Bird. The dilatational properties of suspension of gas bubbles in incompressible Newtonian and non-Newtonian fluids, J. Non-Newt. Fluid Mech. 3 (1977/S) 261-280. [12] M.V. Stasiuk. C. Jaupart, R.S.J. Sparks, On the variations of flow rate in non-explosive lava eruptions, Earth Planet. Sci. Lett. 114 (1993) 505-516. [ 131 S.L. DeSilva, S. Self, PW. Francis, R.E. Drake, C. Ramirez, The Chao dacite and other young lavas of the altiplanopuna volcanic complex, J. Geophys. Res. B 99 (1994) 17805-17825. [ 141 M.N. Rahaman, L.C. DeJonghe, G.W. Scherer, R.J. Brook, Creep and densification during sintering of powder compact, J. Am. Ceram. Sot. 70 (1987) 766-774. [ 151 V.C. Ducamp, R. Raj, Shear and densification of powder compact, I. Am. Ceram. Sot. 72 (1989) 798-804. [ 161 D.P. Wall, SK. Wilson, The linear stability of channel flow of a fluid with temperature-dependent viscosity, J. Fluid Mech. 323 (1996) 107-132.
Earth and Planetary
Science Letters
I56 ( 1998)
13-18
A mechanism for the lateral transport of gas bubbles in silicic lava rising in a vertical conduit S.D.R. Wilson Department
of Mathematics.
University
Manchester
Received
13 June 1997: revised version received
Ml3
of Manchester:
0,rford
Road.
9PL. UK
10 December
1997: accepted
26 December
1997
Abstract It is well-known that the shear viscosity of silicate-rich melts is very sensitive to water content. Thus, as such, a material rises towards the earth’s surface in an eruption, the pressure falls, water exsolves to form bubbles or voids, and the melt viscosity increases. This is responsible for the (by now) well-known large departures of the pressure gradient from hydrostatic. The exsolved gas is known to escape from the rising lava into the surrounding rock and since gas content is believed to have an important role in the transitions between eruption regimes, it is of some interest to understand the mechanism of escape. A key difficulty is to explain how small gas bubbles could migrate laterally across the conduit, in view of their small size and the large liquid viscosity. It is shown here that the increase of melt viscosity with height produces a horizontal pressure gradient directed so that the pressure at the centre of the conduit is larger than at the walls. Thus the fluid in the centre is less viscous and more dense; this configuration is (very probably) dynamically unstable. The resulting instability will have the effect of mixing the lava, thus tending to transport gas bubbles from the centre to the (depleted) regions near the walls. 0 1998 Elsevier Science B.V. All rights reserved. Keywords:
melts: lava; flow mechanism; viscosity; density; gases
1. Introduction The rise of magma through a vertical conduit has
been the subject of much theoretical development recently, and we may cite the papers by Dobran [I], Ramos [2], and Jaupart [3] as typical. We refer more particularly to Jaupart and Allegre [4]. Here, a study was made of the rise of silicic magmas and a principal objective was to understand and explain the transitions between explosive and effusive eruption regimes. Attention was drawn to the important role of gas content, that is bubble volume fraction, in this process. It appears that gas is transferred from the lava to the country rock at the conduit wall and then 0012-821X/98/$19.00 0 1998 Elsevier Science B.V. All rights reserved. PI/ SOOl2-821X(98)00019-3
migrates via pores and cracks, and so it is important to understand this process too. In the paper just cited, it is assumed that the gas escape is controlled by conditions at the conduit wall, in fact, by the pressure drop between the lava and the country rock at some distance. This model therefore incorporates the assumption that there is essentially no barrier to gas transport within the liquid phase itself; the discussion in Jaupart and Allegre [4] leaves this problem open. At present, there appears to be no acceptable explanation of how gas might be transported rapidly across the conduit; it seems inconceivable that gas bubbles of millimetre dimensions could migrate significantly relative to
14
S.D.R. Wilson/Earth and Planetarv Science Letters 156 (1998) 13-18
the liquid, whose viscosity is at least lo5 Pa s, and the permeable foam model of Eichelberger et al. [5] seems unrealistic as pointed out by Jaupart and Allegre [4] _ The present paper offers a contribution to the solution of the problem of gas transport within the conduit. Attention is drawn to the possibility that the basic vertical flow of lava may be subject to a dynamic instability, which, if present, would lead to large scale (i.e. on the scale of the conduit width) overturning and mixing, thus transporting liquid and gas bubbles together from the centre to the walls. In this introduction, we describe the mechanism in qualitative terms. The detailed calculations, for a somewhat idealised model, are given in Section 2 followed by numerical examples. A key feature of silicic magmas is the well-known extreme sensitivity of the value of the viscosity of the melt to changes in the water content. Roughly speaking, the viscosity increases by a factor 10 (or more) for each 1% decrease in water mass fraction. As the lava rises, the pressure falls, exsolved water forms bubbles of water vapour, and the water content falls. Thus the viscosity of the lava depends on the pressure, effectively, and so on the depth. The rapid increase in viscosity as the pressure falls towards the exit value is responsible for the (by now) well-known large departures of the pressure gradient from the hydrostatic value, described in [l--3]. However, there is a second effect, which is small by comparison, but nonetheless important. The viscosity variation with depth creates viscous stresses which set up a horizontal pressure gradient (i.e. transverse to the conduit) directed so that the pressure in the centre is greater than at the walls. (This effect is reinforced by a second, smaller, effect associated with the overall volume expansion of the lava.) This means, in turn, that the fluid in the centre has lower viscosity and higher density than the fluid at the walls. This configuration, or at any rate one very similar to it, is known to be dynamically unstable. There is a very large body of literature on this subject because of the industrial importance of the so-called ‘pipe-lining effect’, and it has recently been surveyed by Joseph et al. [6]. Its possible role in lava flows has been pointed out by Carrigan and Eichelberger [7], who discussed the problem of the flow in a conduit
of a combination of two varieties of lava (mafic and silicic) distinguished by composition and therefore by viscosity. In the present work, we consider only a single material, the point being, of course, that the unstable vertical stratification of viscosity and density is produced automatically by the flow itself. We refer particularly to a paper by Renardy [8], who studied a system consisting of two distinct liquids symmetrically disposed in three layers in sandwich fashion. Renardy showed that when the viscosity and density contrast are as in the present case, the flow is always unstable. The most unstable disturbance mode is the ‘snake’ mode in which the core fluid moves as a whole from one wall to the other, in zig-zag fashion, with amplitude which increases exponentially with time. The growth rates are probably large enough to be significant in the present case. It is worth noting that nearly all of the many articles on this subject are confined to calculation of the neutral curves, that is the boundary in parameter space which marks the transition from stability to instability. Renardy [8], however, provides calculation of the growth rates in the unstable region which is particularly useful here.
2. Theory We propose to regard the mixture of liquid and gas bubbles as a single continuum, with suitably averaged properties, modelled as a (compressible) Newtonian liquid. (This is implied in the theories developed in [ 141.) Thus, no account is taken of the presence of crystals. The constitutive equation is O;j
=
-pSij
+
tij
(1)
where aij is the mechanical stress tensor, p is the thermodynamic pressure and tij is the deviatoric stress tensor, given by T;j =
(K
-
$/A)ekkhij
+
2@ij.
(2)
Here eii is the rate of strain tensor and p and are the shear and bulk viscosity coefficients. (For a derivation of this equation see Hunter [9], for example.) The pressure p is linked to the local temperature and density by an equation of state. We can derive K
S.D.R. Wilson/Earth
and Planetary
this from a combination of an equation expressing the law of mass conservation for water with a solubility law. We are going to assume, for simplicity, that water is the only volatile present and that bubbles consist of water vapour which can be regarded as a perfect gas. Consider a volume of the mixture comprising a volume V, of liquid of density pu and a volume V, of gas of density ps. We assume that as the mixture expands, Vt and pe remain constant. Let (T denote the mass fraction of water in the liquid phase. Then we have pr Vpa + pg V, = constant = mass of water
(3)
(4)
P
where p is the average density of the mixture. The perfect gas law is
in the usual notation. Combining these we find
M
( >
Pea+
or
Pe-P
---p
= constant.
P
o = @‘I2
(7)
(with k = 4.11 x lO-‘j), quoted by Jaupart and Allegre [4]. Eliminating u gives an equation connecting p. p and T, which is the required equation of state: JPrPP)
constant. P
the forms of the components of tij are (9)
(10) (11) The equations of motion, neglecting inertial effects, are aP
-ax =
a
-PS
+ ~“‘X
(8)
Inspection of the magnitudes in this equation shows that the second term is negligible compared with the first, unless p is extremely close to pt. Thus, for void fractions of 2% or more, we can delete the second term and this will be done in most of what follows. (However the accurate expression is also needed, as will be seen.) Next, we turn to the equations of motion. Let us consider the two-dimensional problem of flow between parallel vertical walls, with x as the vertical coordinate, measured upwards, and y horizontal, and U, u as the corresponding velocity components. Then
a
+ Yj-j?
(12) (13)
Eq. 13 shows the origin of the transverse pressure gradient @lay. For the stress component rXvvaries with x, because p does; and the stress component rvYis not zero because au/& is not zero - an effect of the general expansion of the liquid. The system is completed by the overall continuity (mass conservation) equation, assuming steady flow: (14)
(6)
A suitable solubility law is
1.5
13-18
a ap a - -t,” + -z,,,. ay -. G-ax -.
and also V P--t-P _Vf
Science Letters 156 (1998)
We now use Eq. 13 to find a formula for the pressure difference between the centre of the channel and the wall. This can be done conveniently if we accept an answer in terms of the local values of the pressure p and the vertical pressure gradient ap/ax. The pressure p can be regarded as a rough measure of depth. Of course the model can be used to find p as a function of x, but this is rather long and complicated and inessential to the present purpose. Let us assume (and a detailed analysis confirms this) that the transverse velocity component IJ can be set equal to zero as a first approximation. Then, using the simplified forms of r_rvand rvv, we can integrate Eq. 13 with respect to y, giving the pressure contrast Ap as
(15) where urnaxis the velocity on the centre line of the channel (about 3/2 times the average). We have also, from Eq. 14, PU max =
constant
(16)
16
S.D.R. Wilson/Earth
and Planetary
and the simplified version of Eq. 8, namely
P(Pe- P)
(17)
and as we see in Section 3 these are sufficient to express Ap in terms of p and LIplax.
3. Numerical examples First, we turn to the dependence of viscosity on water content. As an example, we take the data for rhyolite given in Lange [lo]. The data points up to 0 = 4% can be reasonably well fitted by an exponential of the form =
plexp(-ba)
(18)
with pt = 1O’o.9= 7.9 x lo9 Pa s and b = 357. Using the solubility law (Eq. 7) gives
P
(19)
Here p.e is the shear viscosity of the liquid phase. The viscosity p which appears in Eq. 15 is, of course, the shear viscosity of the mixture. To keep matters simple we shall set /.L= /..Q,thus disregarding, for the moment, the effect of bubbles on p. We also need an estimate of the bulk viscosity K and since this is highly sensitive to the bubble volume fraction at small voidage we need to account for this. An approximate formula is given by Prud’Homme and Bird [ 111 and can be written K
=
P $LcLc----. Pe -
(20)
P
(The singularity when p = Pe is a false alarm; the stresses produced are finite. See Prud’home and Bird
[ill). In order to use the equation of state (Eq. 17) we need to evaluate the constant on the right and this depends on the original state of the lava. Let us suppose that the original water mass fraction (i.e. at great depth) is 3%. Using Eq. 8, we find that when the pressure is reduced sufficiently to create bubbles occupying 5% of the volume (i.e. p = 0.95~~) the pressure is 464 bar and the water content is 2.8%. Thereafter, as the fluid rises, we can
0.05 = constant = 464 x lo5 X -
0.95
= 2.44 x lo6 Pa.
(21)
(We have set pe = 2400 kg/m”, T = 1120 K in accordance with the data from Lange [lo] and assume that conditions are isothermal.) Eq. 21 can now be used to calculate the density in terms of p. We are now able to express Ap, from Eq. 15, in terms of p, p and dp/ti. Note that dp -=-dx
dp dp
dp dx and, from Eqs. 16 and 17 du - max = ---umax dP dx PdJ
= _U,,,~_(Pe - P) dp PeP dx’
(22)
(23)
Then Ap =
p,e = plexp(-kbp’/*).
13-18
use Eq. 17 which now reads P(Pe - P)
= constant
P
Fe
Science Letters 156 (1998)
dp - (K + ;,)@&@ dp
S&,&X. (24) I Note that dp/dp is negative and so is dpldx, so that Ap is the sum of two positive terms. Any values of p and dp/dx thought to be of interest can now be inserted into Eq. 24, which is completely explicit, apart from the necessity to use Eq. 21 to calculate p in terms of p. As an example, suppose we take p = 50 bar = 5 x lo6 Pa and dp/dx = 3 times the hydrostatic value = 3pg; and put u,,, = 3 mm SC’, corresponding to an average rise velocity of 2 mm s-‘. Estimates of these magnitudes are necessary rough and indirect; these values are in line with a number of sources, for example Stasiuk et al. [12], Eichelberger et al. [5], DeSilva et al [13] and Jaupart [3]. Then we find Ap = 2.3 bar approximately, i.e. a pressure contrast of only about 4%, but, using Eq. 19, creating a viscosity contrast of 7%, easily enough to produce a dynamic instability. To check this last point, we need to refer to the results of Renardy [S] in some detail. Renardy uses the label m to denote the viscosity ratio, which thus takes the value 1.08 in the example just worked out. G is the label given to the ratio of the pressure gradient to the hydrostatic value, taking the value 3 here. Renardy gives, in Fig. 2, the growth rate of long waves in the case m = 1.1, G = 3 (and sev-
S.D.R. Wilson/Earth and Planetay Science Letters 156 (1998) 13-18
era1 other cases). (The numerical coincidence here was deliberately contrived.) The growth rate is very small, about 10v7 in the chosen dimensionless units; but Renardy’s Fig. 5 shows that the growth rate of short waves is much larger. Only the case m = 2 is shown; the growth rates are larger than the long wave values by a factor of about 104. Assuming that the case m = 1.1 behaves in roughly the same way, we obtain a maximum growth rate in the range lo-“- 10e4. The unit of time in these computations is ,/@ where C is the channel (conduit) half width. Rutting f? = lOm, the conclusion is that the time scale of the growth of the instability is about 103lo4 s, say 1 h. The most unstable waves have a wavelength of the same order as the conduit width. The depth of the bubbly region can be calculated to be on the order of hundreds of metres, using Eq. 21 and assuming that the mixture remains cohesive up to about 50% voidage. This gives a residence time, at 3 mm/s, of about 40-50 h, much longer than the growth time. Note that a lower value of urnax will result in slower growth, but longer residence time, the two effects partially offsetting each other. (A unit of time based on viscosity, e.g. pf?/p, might seem more natural. This differs from Renardy’s actual choice by a factor of the Reynolds number. It would be natural, for R < I, to rescale the time as indicated; then the dimensionless growth rates would differ by the corresponding factor. But the physical time scale would, of course, come out the same.) These numbers are quite typical, as a little trial and error will show. Very similar results are obtained with umaX= 20 mm-‘, p = 100 bar, for example. The choices made were, as noted, aimed at enabling a reasonable comparison with Renardy’s computations. Furthermore, the pressure gradient may well exceed the lithostatic value by more than a factor 3. The general conclusion is that the instability is most likely to occur in the upper sections of the conduit where there are marked departures from hydrostatic equilibrium.
4. Concluding
remarks
Attention has been drawn to a possible hydrodynamic instability of vertical conduit flow of silicic
17
magmas. This results from a horizontal pressure gradient, stemming mainly from the variation of the shear viscosity with height. The mechanism does not depend on Newtonian behaviour of the melt; that assumption serves only to make possible simple explicit calculations. A somewhat crude model for the dependence of the melt viscosity pe on pressure has been developed (Eq. 19). Again, this was used simply to make the calculations explicit; any model similarly sensitive would give similar results. It would be of considerable interest to refine Eq. 19 so as to include the direct dependence of the liquid phase viscosity on pressure (i.e. at fixed composition), an effect discussed by Lange [lo]. Also of possible importance is the effect of void fraction on the viscosity of the two-phase mixture. Here, there is conflicting evidence. Jaupart and Allegre [4] and Dobran [ 11, among others, make use of models according to which the viscosity increases with void fraction, presumably resulting from ‘bubble interference’. To the present author, this seems unlikely to hold in the relevant parameter regime, where surface and inertial effects are completely negligible. Experiments on porous glasses by Rahaman et al. [ 141 and Ducamp and Raj [ 151 show a decrease in viscosity with void fraction, and these appear to be more relevant. The application of Renardy’s analysis [8] to the present case perhaps needs some caution. Renardy studied a system of two contrasting immiscible liquids, whereas we have one stratified liquid. The difference lies not in interfacial effects, which do not matter much in the cases cited, but in the fact that in the two-liquid case, the shear viscosity is a material constant convected with the liquid, whereas here it depends on the local pressure. The stability of a fluid whose viscosity varies continuously (e.g. as a result of temperature variations) has also been studied; see, for example, Wall and Wilson [16]. But only the neutral curves are given, not the growth rates in the unstable region. An extension of Renardy’s work to something like the present case would not be practicable because in addition to the feature just mentioned, the basic flow exhibits a rapidly varying pressure gradient, and the normal mode decomposition, leading to ordinary differential equations for the disturbance variables, is not possible. [FA]
18
S.D.R. Wilson/Earth
and Planetarv
References [I] F. Dobran, Non-equilibrium
[2]
[3] [4]
[5]
[6] [7] [S] [9]
flow in volcanic conduits and applications to the eruption of Mt. St. Helens on May 18 1980 and Vesuvius in AD79, J. Volcanol. Geotherm. Res. 49 (1992) 285-311. J.1. Ramos, One-dimensional, time-dependent, homogeneous 2-phase flow in volcanic conduits, Int. J. Numer. Methods Fluids 21 (3) (1995) 253-278. C. Jaupart, Physical models of volcanic eruptions, Chem. Geol. 128 (1996) 217-227. C. Jaupart, C.J. Allegre, Gas content, eruption rate and instabilities of eruption regime in silicic volcanoes, Earth Planet. Sci. Lett. 102 (1991) 413-429. J.C. Eichelberger, C.R. Carrigan, H.R. Westrich, R.H. Price, Non-explosive silicic volcanism, Nature 323 (1986) 598602. D.D. Joseph, R. Bai, K.P. Chen, Y.Y. Renardy, Core-annular flows, Ann. Rev. Fluid Mech. 29 (1997) 65-90. C.R. Carrigan, J.C. Eichelberger, Zoning of magmas by viscosity in volcanic conduits, Nature 343 (1990) 248-25 I. Y. Renardy, Viscosity and density stratification in vertical Poiseuille flow, Phys. Fluids 30 (1987) 1638-1648. S.C. Hunter, Mechanics of Continuous Media, Ellis Hor-
Science Letters
156 (1998)
13-18
wood, 1976, Chichester, pp. 135 and 149. [IO] R.A. Lange, The effect of H20, CO1 and F on the density and viscosity of silicate melts, in: M.R. &roll, J.R. Holloway (Eds.), Volatiles in Magma& Rev. Mineral. 30, 1994. [I I] R.K. Prud’Homme, R.B. Bird. The dilatational properties of suspension of gas bubbles in incompressible Newtonian and non-Newtonian fluids, J. Non-Newt. Fluid Mech. 3 (1977/S) 261-280. [12] M.V. Stasiuk. C. Jaupart, R.S.J. Sparks, On the variations of flow rate in non-explosive lava eruptions, Earth Planet. Sci. Lett. 114 (1993) 505-516. [ 131 S.L. DeSilva, S. Self, PW. Francis, R.E. Drake, C. Ramirez, The Chao dacite and other young lavas of the altiplanopuna volcanic complex, J. Geophys. Res. B 99 (1994) 17805-17825. [ 141 M.N. Rahaman, L.C. DeJonghe, G.W. Scherer, R.J. Brook, Creep and densification during sintering of powder compact, J. Am. Ceram. Sot. 70 (1987) 766-774. [ 151 V.C. Ducamp, R. Raj, Shear and densification of powder compact, I. Am. Ceram. Sot. 72 (1989) 798-804. [ 161 D.P. Wall, SK. Wilson, The linear stability of channel flow of a fluid with temperature-dependent viscosity, J. Fluid Mech. 323 (1996) 107-132.