Lava flow shapes and dimensions as reflections of magma system conditions

ELSEVIER Journal of Volcanology and Geothermal Research 78 (1997) 3 l-50

Lava flow shapes and dimensions as reflections of magma system conditions

Mark V. Stasiuk *, Claude Jaupart fnstitut de Physique du Globe de Paris et Universitt! Paris 7, 4 place Jussieu, 75252 Paris, France

Received 1 July 1996; accepted 16 December 1996

Abstract

The physical conditions producing variation in the shapes of lava extrusions, from thin flows to domes to spines, are investigated with focus on the influence of the lava thickness. The accumulation of lava above the vent acts to increase the

vent pressure and reduce the eruption rate. In turn, the eruption rate determines how the volume and thickness of erupted lava change with time. The spreading of lava and the flow of magma from a deep reservoir through a conduit system are therefore coupled. A fluid dynamical model is developed to evaluate the influence of magma chamber and conduit dimensions, chamber overpressure, magma viscosity and dome bulk viscosity on the eruption behaviour and dimensions of the lava flow. The model assumes laminar flow of Newtonian, constant-density magma from an overpressured reservoir with elastic walls, up a vertical cylindrical conduit, and onto a planar horizontal surface, on which grows a radially symmetrical lava flow of constant bulk viscosity greater than the magma viscosity. Two dimensionless numbers characterize the resulting behaviour: the ratio of a scale erupted lava volume to the total volume available for eruption, and the ratio of the pressure due to a scale thickness of lava to the magma chamber overpressure. The dimensionless numbers are used to map the flow behaviours and reveal two extreme extrusion regimes separated by a narrow transition zone. One extreme corresponds to thin and voluminous flows, and occurs for large, deep, strongly overpressured magma reservoirs, narrow conduits and low viscosity contrasts between the magma and lava. The other extreme corresponds to small volume domes and spines, and is produced by shallow, weakly overpressured chambers, wide conduits and high viscosity contrasts. Near the transitional domain, small changes of parameters are sufficient to cause major changes of eruption regime. The model provides explanations for the variation in gross morphology of lava extrusions, transitions between flow, dome and spine eruptions and for long-lived episodic extrusion of lava. Model predictions are consistent with features of the eruptions of Montagne Pelte (1902-1903), Mount St. Helens (1980-19861, and Montserrat (1995-1997 period). The model is used to estimate the magnitude of overpressure driving eruptions and suggests that magma buoyancy is important in driving many terrestrial lava dome eruptions. The model further suggests that thick extrusions on Venus were driven primarily by buoyancy sourced at great depth.

Keywords: lava; dome; spine; conduit; magma chamber; eruption rate; overpressure

1. Introduction * Corresponding author. Present address: Environmental Sci-

ence Division, Lancaster University, Lancaster LA1 4YQ, UK.

Phone: +44-1524-593970; Fax: +44-1524-593985; E-mail:

[email protected]

Lava flows of similar composition exhibit a wide range in shape. Rhyolite, for example, occurs as

0377-0273/97/$17.00 0 1997 Elsevier Science B.V. All rights reserved. PIf SO377-0273(97)00002-4

32 M.V. Stasiuk, C. Jaupart/ Journal of Volcanology and Geothemzal Research 78 (1997) 31-50

flows with lengths of tens of kilometres and thicknesses of a few tens of metres (Hausback, 1987; Bonnichsen and Kauffman, 1987; Manley, 19961, domes with diameters of a few kilometres and thicknesses of a few hundred metres, or spines with widths of only tens of metres. A fundamental problem is to understand what controls the differences in morphology, and whether it contains information about the parent volcanic systems. Numerous factors influence the gross morphology of lava flows, such as rheology, total erupted volume, eruption rate, cooling rate and slope. A factor which is less well known, however, is the thickness of lava built up directly above the vent. Where the lava is thick, its effect may be large because it exerts considerable back pressure; for a lava dome 200 m thick, for example, this pressure is about 5 MPa and comparable to the values of chamber overpressure which can drive magma to the surface (see Tait et al., 1989). In this paper we investigate the influence of lava thickness on the eruption behaviour and gross morphology of lava extrusions.

An important aspect of the role of lava thickness is that it cannot be investigated independently of morphology and eruption behaviour because of dynamic coupling. The eruption rate and lava thickness are coupled and time varying: the thickness increases with time and with eruption rate (Huppert et al., 19821, while the eruption rate decreases due to chamber decompression (Scandone, 1979; Wadge, 1981) and the pressure exerted by the thickening lava. In this case the eruption rate varies in a way which cannot be assumed, but must be solved for along with the thickness. Observations show that decreasing eruption rate can be correlated with increasing thickness of lava (Huppert et al., 1982; Jaupart and Allegre, 1991; Stasiuk et al., 1993a,b). The increasing lava thickness may ultimately cause the eruption to stop, thus altering the total volume erupted. Lava morphology is also involved since it can be correlated with eruption rate and volume erupted (Huppert et al., 1982; Jaupart, 1991; Fink and Griffiths, 1990, 1992; Fink and Bridges, 19951, two variables which are influenced by the flow thickness. A practical aspect of the coupling is that it connects the dynamics of the surface activity to numerous features of the magma plumbing system. This means that lava morphology and behaviour

contain information on the structure of the volcanic system at depth, information that would not be available if the eruption rate was imposed externally as in a laboratory experiment. An understanding of the coupling is profitable in cases where one only knows the flow dimensions, as on other distant planets like Venus.

We examine the problem through a simplified physical model which emphasizes the new aspects introduced. In the first part of the paper the model is developed and numerical solutions derived. In the last part the results are considered in the light of observations of volcanic eruptions. Implications of the analysis are discussed in relation to general features of lava extrusions on Montagne PelCe, Montserrat, Mount St. Helens and Venus.

2. Theoretical framework

2.1. Equations and assumptions

We consider the case of simple lava flows fed from a reservoir located at depth L. Magma of constant dynamic viscosity p,,, moves to the surface through a cylindrical conduit of constant radius R,

and erupts onto a horizontal planar surface (Fig. 1). Thus, the lava flow builds up over and covers the

Cdduit

I

Fig. 1. Idealized geometry of a volcanic system. Magma flows

from chamber, up cylindrical conduit to flat surface.

M.V. Stasiuk, C. Jaupart/ Journal of Volcanology and Geothermal Research 78 (1997) 31-50 33

vent. In lava eruptions the erupted material contains relatively small quantities of gas. As the gas compo- nent is compressible and has solubility increasing with pressure, in a closed system the gas volume fraction below the vent must be less than at the surface. Calculations show that for subaerial eruptions it will become significant only within tens of metres of the surface (Jaupart and Tait, 1990; Stasiuk et al., 1993b). For submarine eruptions, wherein the ‘surface’ pressure is that at the sea floor, the erupting lavas have even smaller vesicularities. Therefore, in order to calculate the eruption rate, the ascending magma can be treated as bubble free and of constant density p,. The rheology of bubble-free magma is approximated well as Newtonian if it has less than about 40% crystals and for natural deformation rates (Pinkerton and Stevenson, 1992; Iejeune and Richet, 19951, and we shall restrict our study to such cases. We assume that flow in the conduit is in the laminar regime, which is valid for most non-explosive eruptions on Earth (Jaupart and Tait, 1990). In addition, typical ascent times for magma are much shorter than eruption durations, and hence flow in the conduit is in a quasi-steady state. These assumptions are discussed in detail by Stasiuk et al. (1993b), who tested them against data from a number of lava eruptions. Flow in the conduit is assumed to be fully developed, valid for high-viscosity fluids and conduits much longer than they are wide. In the laminar regime, for a vertical conduit much smaller in horizontal than vertical dimensions, the equations of motion reduce to:

(2.la)

where P is the pressure in the conduit, r is radius measured from the conduit axis, z is the vertical coordinate measured from the magma chamber, U is the flow velocity and g is the acceleration of gravity. Fq. (2.la) shows that pressure depends only on z. Integrating Eq. (2.lb) twice with respect to r for constant ,u, and p, determines the radial profile of velocity. The volume flux is then equal to:

Q = iRcU2=rdr (2.2a)

where R, is the conduit radius. The integral is given by:

(2.2b)

The equation shows separately the effects of the driving term and the conduit parameters. For constant density and viscosity of magma in the conduit and constant conduit radius, the equation shows that the driving pressure at any time must vary linearly with depth. Hence the volume flow rate can be written as:

-p,g 1 (2.2c)

where Q(t) is written to emphasize that the eruption rate may change slowly with time; V(t) is the volume of lava and L is the conduit length (or chamber depth). Flow is driven by the pressure difference between the chamber and the vent. P(0) is the time-varying pressure at the base of the conduit, or the top of the magma chamber, and P(L) is the time-varying pressure at the top of the conduit, or the base of the lava. Eq. (2.2~) may be derived for other conduit shapes, for which it has the same functional form but a geometrical factor slightly different from the value rRz/8 for a cylinder. (Stasiuk et al., 1993b3).

The pressures can be written as:

P(L) =Pa+PmgW) (2.3a)

P(0) =P,+p,gL+AP(t) (2.3b)

where h,(t) is the thickness of lava over the vent, A P(t) is the overpressure in the magma chamber, Pa is the atmospheric pressure and p, the country rock density. For lava flows erupted on the sea floor, the constant g in Eqs. (2.2b), (2.2c), (2.3al and (2.3b) must be replaced with a reduced gravity. For the moment, we assume that the pressure P(L) is negligibly influenced by overpressures developing within a solid carapace encasing the lava flow. The influence of such overpressures on the model results will be discussed later.

In order to calculate how the chamber pressure varies with time, one needs to specify how the chamber walls respond to volume changes. Close to

34 M.V. Stasiuk, C. Jaupati/Joumal of Volcanology and Geothermal Research 78 (1997) 31-50

the chamber, country rock is heated to high tempera-

tures and hence may behave viscously. Heat balance

considerations show that the volume of heated coun-

try rock must be significantly smaller than the vol-

ume of magma (Huppert and Sparks, 1988). By contrast, deformation involves a much larger volume

than the chamber, i.e., a volume which is at rela-

tively low temperature. This explains in part why

deformation appears to be rapid and reversible (i.e.,

inflation followed by deflation) on the time scale of

the eruptive process. In keeping with these various

observations, we assume elastic behaviour for the

chamber walls. Under such conditions, a loss of

volume from the reservoir produces a proportional

decrease in its pressure (Blake, 1981; Tait et al.,

1989). Thus, one may write:

AP(t) = AP,,, - gV(t) (2.4) "m

where A P,,, is the chamber overpressure at the start of eruption, and the rate of pressure change with

volume is the constant ffl/V, where f is a geomet-

rical factor, p is the bulk modulus of the country

rocks and V, is the magma chamber volume. f is 0.5 for a spherical chamber and changes by about a

factor of two for a range of realistic shapes (i.e., a cylinder or ellipsoid). Substituting Eqs. (2.3a), (2.3b)

and (2.4) into Eq. (2.2~) and assuming that pm = p,

= p gives:

dV Q(r)=%

(2.5a)

When the magma and rock densities are not equal,

this equation has an additional constant term which can be conveniently lumped with the chamber over-

pressure, giving:

(2Sb)

where Ap is the density difference p, - pm. These two equations are essentially the same, except that

Eq. (2.5b) effectively has an elevated chamber over-

pressure due to magma buoyancy. Thus, for simplic-

ity, we will use Eq. (2.5a) for our model calculations

and discuss the effects of buoyancy later. Eq. (2.5a)

shows how the eruption rate depends on the chamber

and conduit dimensions and geometries, the chamber

overpressure, volume erupted, magma properties and lava thickness. Coupling in the problem enters mainly

through the term for the lava height, which reduces

the pressure difference causing flow. In order to calculate the height of the lava, it is

necessary to have a quantitative description for the

lava spreading. Spreading is controlled by the erup-

tion rate, the rheological properties of lava and their

dependence on cooling, and the substrate topogra-

phy. Here, we restrict our attention to a flat topography, i.e., a plain or a crater floor, which is valid for

many cases of interest. The apparent viscosity of

lava flows commonly increases over that of the

supplying magma with time (e.g., Walker, 1967;

Huppert et al., 1982; Fink and Zimbelman, 1989; Naranjo et al., 1992; Borgia and Linneman, 1989),

eventually tending to infinity as flow stops. The causes of this are cooling, producing crystallization

and viscosity increases at the margins (Griffiths and

Fink, 1993; Crisp et al., 1994; Stasiuk et al., 1993a),

and degassing from the flow interior (Sparks and Pinkerton, 1978; Lipman and Banks, 1987). Cooling

can produce a viscous, plastic or brittle surface layer

(Griffiths and Fink, 1993). Fink and Griffiths (1992)

have examined experimentally the morphologic and

dynamic consequences of a brittle crust. Its influence is strongest during the formation of compound flows,

whose irregular form and motion are produced by the repeated formation and rupture of the crust on

the flow surface. Where flows are not compound but spread in a smooth, geometrically regular way, their

motion appears more consistent with control by ei- ther a plastic or viscous surface layer.

At present there is no general quantitative model which accounts for the various physical effects of cooling on spreading. For the purposes of modelling gross morphology and general eruption behaviour, previous workers (e.g., Huppert et al., 1982) have shown that simple lava flow units can be usefully treated as isoviscous gravity currents with bulk viscosities incorporating the various influences of cooling. Although a crude approximation of a complex

M.V. Stasiuk, C. Jaupart/Journal of Volcanology and Geothermal Research 78 (1997) 31-50 35

physical problem, this has been found to work rea- sonably well and represents a considerable simplification. Using the concept of a bulk viscosity enables us to focus directly on the problem of dynamic coupling between the lava thickness and the eruption rate. The observed bulk viscosity of a lava flow is typically greater than the viscosity of the magma in the subsurface, and increases with time during an eruption. Predictions of bulk viscosity of lava flows with time are not currently available, and some functional dependence must be assumed. Calcula- tions of bulk viscosity with time have been achieved with a set of simplified equations for momentum and energy for conditions of constant eruption rate, New- tonian temperature-dependent viscosity and radial temperature variations in the flow (Bercovici, 1994). However, to do so here would obscure the physics of interest and, given the assumptions, still would not guarantee realistic bulk viscosity values. Sakimoto and Zuber (1995) took prescribed time functions for the bulk viscosity of lavas based on observations, and solved for the spreading characteristics, but they did not investigate how the time-dependent bulk viscosity might be coupled to the eruption rate through the flow thickness. For our present purposes, the main effect is that a larger bulk viscosity leads to a greater flow thickness. Thus, we choose simply to specify the bulk viscosity of the flow as an independent value which is larger than the magma viscosity. One may then consider the influence of changes in the lava bulk viscosity with time, using the same model framework, by considering a sequence of the constant viscosity results. This approach is discussed further in Section 4.3 and Appendix A, where it is shown that increasing bulk viscosity does not funda- mentally alter the results of our relatively simple model, and in fact the physical effects described by the model will be accentuated. The influence of non-Newtonian rheology is discussed briefly in Sec- tion 5.

Huppert ( 1982) determined similarity solutions for the spreading of viscous-gravity currents fed by an eruption rate function (power law) prescribed independently. This cannot be done in the present problem since the eruption rate is a variable to be determined. In addition the general similarity solution does not estimate the flow central thickness, clearly important here, and it treats only point sources

of fluid volume. To overcome these difficulties, we develop an approximate equation for the spreading of viscous-gravity currents which is valid for almost any supply rate function, and hence which can be used in problems where the supply rate function is part of the solution. This solution is derived and its predictions are shown to compare favourably with the height and radius of laboratory flows in Ap- pendix A. The central height of the flow is given by the following equation:

-l/4

h(O,r) = V(t) pg /‘V( t’)4 dt’ + AR,” Gw 0 1

(2.6)

where p, is the bulk dynamic viscosity, and the constant parameter A = 12.5.

Since we model lavas as viscous flows, they continue to spread after the eruption rate has dropped to zero. Observations reveal that simple lava flows halt approximately when extrusion ends (Borgia et al., 1983; Guest et al., 1987; Borgia and Linneman, 1989), although it can be difficult to define precisely when an eruption stops because the extrusion rate decreases gradually to negligible values (e.g., Hup- pert et al., 1982). Late-stage adjustments of lava flow shapes by spreading and sagging at constant volume are small (Huppert et al., 1982; Swanson et al., 19871, presumably because cooling and solidifi- cation rapidly gain dominance over spreading pro- cesses. Once the eruption rate becomes small and the flow has nearly constant volume, its emplacement is essentially complete. In our model we observe the expected result that, when the supply rate becomes small, the flow behaviour becomes indistinguishable from that of spreading at constant volume. Thus the end of the eruption is taken to occur when the eruption rate drops below a threshold value, and at this point the lava emplacement is considered complete. Calculations show that, provided the selected threshold value is small compared to the initial eruption rate, its precise value is irrelevant (see Section 5). In the results below, the threshold value is fixed at 0.01 of the initial eruption rate.

The system of Eq. (2.5a) and Eq. (2.6) along with the conditions, V(0) = 0, h,(O) = 0, Q(O) =

36 M. V. Stasiuk, C. Jaupart / Journal of Volcanology and Geothermal Research 78 (1997) 31-50

rR~AP,/8p,,,L and Q<t,> = O.OlQ(O) where t, is

the time at the end of the eruption, are sufficient to

specify the behaviour.

2.2. Model scaling

The behaviour of the coupled system is investi-

gated in terms of dimensionless groups of variables

in order to provide the most compact description of

the dynamics. We proceed by scaling Eq. (2.5a) and Eq. (2.6), where scales will be denoted by a sub-

script ‘s’. For a pressure scale, we take the initial

chamber overpressure:

P,=AP,,,

For the volume flux Q(t), we choose the

eruption rate:

2.7a)

initial

2.7b)

The scale height of the lava is the viscous current

scale (Huppert, 1982):

(2.7~)

and the lava flow radius scales on the conduit radius:

R, = R, (2.7d)

Time and volume scales can be defined from V, =

Q,T, = H, Rt, yielding:

and

V,~H,Rf=R:( -i’l’(%,“’ (2.7f)

Dimensionless values of the variables will be denoted with overbars, such that t = T,i, rN = R,f, h, = H,h, V= V,v and Q = Q,a. Substituting these into Eq. (2.5a) and Eq. (2.6) yields the dimensionless eruption rate as:

and lava flow thickness as:

1 -l/4

(2.9)

where the dimensionless groups V and p are given by:

5= (3 ApJFL)l:l)( y4( :)“‘=; (2.10a)

(2.10b)

The numbers V and jj have simple physicalinterpre- tations. V is a ratio of a lava volume to the maximum volume which can be erupted V, = V, A P,/fp

(determined from Eq. (2.4) for zero final overpres-

sure). Thus, V is a ‘volume number’ and measures the sensitivity of the volcanic system to the volume

erupted. A volcano with V z=- 1 runs out of magma

before a large dome can grow. Only for V -=K 1 can

an extensive flow be emplaced. From Eq. (2.10a),

the conditions favouring such a situation are a large,

deep, strongly overpressured magma reservoir, a narrow conduit and a low viscosity contrast between

magma and lava.

p is a ‘pressure number’, a ratio of the weight of the column of lava above the vent to the initial

chamber overpressure. Small values of p imply that

the lava has little influence on the eruption rate and permits extensive flows to be emplaced. For large values of p the eruption rate will be strongly re-

duced by the weight of lava. Conditions leading to

large p are a shallow, weakly overpressured magma chamber, wide conduit and high viscosity contrast.

Note that these conditions are not required simulta-

neously. These conditions favour a relatively high initial eruption rate into a thick, viscous extrusion: the lava will initially build up rapidly in height, its weight soon matching the small chamber overpressure and halting further extrusion before the lava can spread. Apart from the volume and geometry of the magma chamber and the country rock elastic properties, which do not enter the pressure number, the conditions for strong coupling are the opposite of

M. V. Stasiuk, C. Jaupart / Journal of Volcanology and Geothermal Research 78 (1997) 31-50 31

those for emplacing voluminous flows. This implies that extrusions of small volume, such as domes and spines, are strongly coupled to the subsurface feeding system.

3. Results

3.1. Behaviour in time

Eq. (2.8) and Eq. (2.9) were integrated numeri- cally with fourth- and fifth-order Runge-Kutta methods for volume and pressure numbers ranging from 10e9 to lo-*. Initial values were obtained from the behaviour at small time, when Q = 1, so that ? z t and Eq. (2.9) gives h = A- ‘14t.

Solutions are shown in Fig. 2 for a volume number V = 0 and selected values of the pressure number j. For no coupling and constant supply rate (5 = p = 0) the numerical solution matches the analytical solution derived by integrating Eq. (2.9). The numerical solutions show an initial rapid increase in flow height and achieve a constant value for i > 10, although the time scale for growth and maximum height shrink with increasing p, evident from Eq. (B.l). Fig. 2 also shows that the solution B.l, for large p, approximates well the numerical results for ji 2 2.

0- I

0 2 4 6 6 10

t/T,

Fig. 2. Dimensionless central height of lava flow as function of

dimensionless time, for case where volume number C = 0. Solid curves = numerical solutions; dashed curve = Eq. (B.l).

1

0.5

I

0 0 5 10 15 20

t/T,

Fig. 3. Dimensionless height versus time for j5 = 0. Solid curves

= numerical solutions; dashed curve = Eq. (B.2). Eq. (B.2)

matches for all nonzero values of V but is only shown for the case

V = 0.1. The dot terminating the curve for V = 1 indicates the end

of eruption.

For pressure number jj = 0 the solutions for different values of Y are shown in Fig. 3, along with Eq. (B.2). The numerical and analytical solutions agree well, the solution B.2 matching for all nonzero values of 5. The flow initially increases in thickness rapidly but thereafter thins. For increasing V the maximum dimensionless height decreases due to the rapid decrease of chamber pressure.

In many physical systems, dynamic coupling can produce naturaloscillations as a result of inertial effects. For this problem, however, no oscillations in the lava thickness occur since the system is in the regime of creeping flow and effectively overdamped.

3.2. Eruption products

The time-dependent behaviour contains the most information on the dynamics of volcanic systems, and it is unfortunate that very few eruptions have been monitored with an appropriate degree of detail to permit comparison with the model results. One exception is the 1979 eruption of the Soufribre de St. Vincent, previously analysed (Huppert et al., 1982; Stasiuk et al., 1993b). Where possible, we shall discuss qualitatively salient points in the observations of time-dependent behaviour, but we focus on measurements and observations at the end of eruptions. For example, in Fig. 4a the final erupted

38 M.V. Stasiuk, C. Jaupart/ Journal of Volcanology and Geothermal Research 78 (1997) 31-50

volume of lava (dimensionless) is shown as a function of p for four values of i. For jj < 0.1 the erupted volume is independent of j and given by V- ’ (see Appendix B). This case is dominated by theplumbing system and is little affected by the

1 O-3 1 o-2 10-t 100 10'

pressure numberp

I I I

lo-3 lo-2 1D' 100

pressure numberp

10'

surface activity. For j 2 2, however, the volume erupted is independent of V and given by p- ’ (Appendix B). In this case the behaviour is dominated by the surface activity and independent of the chamber size and shape and the country rock elastic

10.1 -- ’ 1 I

1 o-3 1 o-2 IO-’ 100 10’

pressure number p

100

v=l

101

IO-’ - 102

103

10.2 L--_._..._._._-l_ IO-3 1 o-2 1 o-1 100 10’

pressure n urn ber p

pressure number p

Fig. 4. (a) Final dimensionless volume as function of pressure ratio p for various values of T. (b) Dimensionless eruption duration. (c) Final

dimensionless flow radius. (d) Final (unscaled) aspect ratio. (e) Final dimensionless chamber pressure.

M.V. Stasiuk, C. Jaupati/ Journal of Volcanology and Geothermal Research 78 (1997) 31-50 39

properties. The behaviour can thus be separated into two regimes, large and small volume eruptions, which are connected by a sharp transition for p - 1. In the vicinity of this value of p, the final volume of all the flows decreases dramatically with increasing j and the coupling between the surface activity and all features of the subsurface plumbing system is strong.

Similar features are evident for other aspects of the eruption behaviour, shown in Fig. 4b-e, in which the results are given as functions of p for the same values of 5 as in Fig. 4a. In general the large volume regime, for small p, is also characterized by eruptions of long duration, producing flows with great lateral extents and small aspect ratios, for which the overpressure in the chamber becomes completely exhausted. The large pressure number regime ex- hibits the opposite features. There are some interest- ing details in the behaviour of these variables. For example, there is visible in Fig. 4b, near the regime transition, an increase in the duration before the dramatic fall-off. This is a result of the increasing importance of coupling between the surface flow and the deeper system. Such behaviour is dominated by decreasing pressure in the chamber until the pressure is sufficiently small for the flow thickness to act on the eruption rate. With increasing coupling (increasing p) the eruption is stopped by the lava weight at earlier times. Fig. 4d shows that the values of aspect ratios for the two extrusion regimes are joined by an oscillation rather than a step. With increasing ii, the aspect ratios decrease slowly before increasing steeply to a peak value of about 0.5 for p w 1. Thereafter the aspect ratios fall with increasing p. The initial decrease is produced by the same long- time effect discussed for Fig. 4b, in which the eruption duration increases and permits the flow to spread. The peak in the curves occurs at the value of p at which the final radii are equal to the conduit radius. The fall after the peak is due solely to a decreasing final flow thickness, for which little lava can be extruded before the eruption is stopped by its weight.

The final dimensionless chamber overpressure is only weakly dependent on V for most reasonable values of this number (Fig. 4e). Of particular interest is that the overpressure at the end of the eruption increases with increasing p. For p _ 1, the pressure in the magma chamber at the end of the eruption is indistinguishable from that at the start. For V 2 1, the

overpressure falls as a result of rapid depressuriza- tion with lost volume, corresponding in general to very small magma chambers.

In summary, the results show that there are two eruption regimes separated by a sharp transition and dependent mainly on the pressure number. Small values of p, corresponding to deep, strongly overpressured chambers, narrow conduits and low viscosity contrasts, produce voluminous, long-lived eruptions and extensive, low aspect ratio lava flows. Large values of j correspond to the opposite conditions, and lead to thick flows which balance the weak overpressure in the chamber after a short time. Thus, large values of p favour the short-lived eruption of small volume domes and spines. The key result is the sharp transition between the large and small volume regimes. Dome eruptions are predicted in a small region of the parameter domain, and hence require specific conditions.

4. Comparisons with nature

In this section we make comparisons between the basic model predictions and general features of lava extrusions. We do not attempt detailed studies of specific eruptions because the model is simplified. Relevant extensions are suggested in Section 5. Our main intentions are to verify that the model predictions are of a relevant order of magnitude, and to comment on the importance of flow thickness in various observed eruptions.

4.1. Eruption scale values

For general illustration we focus on representative compositions of basalt, andesite and rhyolite. Ranges of the parameters, and corresponding median values to be used in a scaling analysis, are shown in Table 1. For basalt, the maximum lava bulk viscosity is based on field estimates by Fink and Zimbelman (1989). For andesite a maximum is taken from observations made at Lonquimay volcano, Chile (Naranjo et al., 1992). For rhyolite there are no appropriate data: a maximum lava bulk viscosity is taken from analyses of the spreading behaviour of the 1979 Soufriere de St. Vincent dome (Huppert et al., 1982). We consider only sub-aerial eruptions of simple flow


Table 1

Parameters used in calculations

Parameters d Basalt a Andesite b Rhyolite ’

(lloo°C) (1ooo”c) (900°C)

p (kgmm3) 2730 2580 2330

pCLm (Pas) 102 104 109

wI (Pas) 102-104 104-109 109-10’5

cL,/cL, l-102 l-105 l-106

R, (m) 0.1-2 OS- 10 l-50

L Cm) 500-10000 500-10000 500-loo00

V, Cm) 1os-10’2 1os-10’2 1os-10’2

AP, (Pa) 0.5-10x 106 0.5-10x 106 0.5-10x 106

f 0.5 0.5 0.5

p (Pa) 10’0 10’0 10’0

g (ms-‘1 9.81 9.81 9.81

v 10-7 10-5 0.005

L, (m3 s-l)

0.005 0.05 0.3

10 60 0.5

H, (m) 0.8 10 60

T, (s) 0.1 4 105

V, cm31 0.8 250 4x 104

H, /Rs 0.8 2 2.5

a wt.% oxides: SiO, 49; TiO, 2; Al,O, 16; MgO 7; CaO 10; Fe0

11; Na,O 3; K,O 1.

b wt.% oxides: SiO, 58.5; TiO, 1.5; A120, 16.3; MgO 2.4; CaO

5.5; MnO 0.2; Fe0 6.8; Na,O 5; K,O 0.9.

’ wt.% oxides: SiO, 76.5; TiO, 0.1; Al,O, 12.6; MgO 0.1; CaO

0.3; MnO 0.1; Fe0 1.1; Na,O 4.3; K,O 4.6.

d Density and viscosity (magma) at given temperatures calculated

using the methods of Bottinga et al. (1982) and Shaw (19721,

respectively. Bulk modulus of country rocks is a representative

value from Touloukian et al. (1981). Lava viscosity from literature

references are discussed in text.

units, but note that enhanced cooling under water is

likely to increase the lava bulk viscosity. Ranges in

conduit radii are based on commonly observed dyke dimensions and published data (e.g., Wada, 1994).

The conduit lengths, or magma chamber depths, are based on geophysical inferences (e.g., Iyer et al.,

1990; Endo et al., 1990; Murray, 1990; Ishihara,

1990; Ewart et al., 1990; Lees, 1992; Moran, 1994). The chamber volumes have a lower limit of a sphere of 500 m diameter, chosen arbitrarily, and an upper limit of the volumes of the largest known pyroclastic eruptions (e.g., Wilson, 1986).

Lava dimensions and eruption durations can be predicted for the values of V and 3 in Table 1 by scaling the dimensionless results. For all three compositions, predicted median volumes are of order lo7

m3. For basalt and andesite, predicted durations are of the order of hundreds of days, flow extents are kilometres and thicknesses are a few metres. For

rhyolite, median durations are years to tens of years, forming flows with extents less than a kilometre and

thicknesses of about 100 m. These values compare

well with nature. For example the dome eruptions at

Santiaguito and Bezymianny took place over several

tens of years while those at Mount Lamington and

Mount St. Helens occurred over several years (Swan-

son et al., 1987).

4.2. Flow aspect ratios

The coupling between lava weight and the subsur-

face volcanic system can result in strong variations in the final dimensions reached by lava flows. A

simple shape index is the aspect ratio, which is a

function of the volume and pressure numbers. The parameter values of Table 1 can be used to deter-

-- mine a field for each magma composition in ( p, v) space (Fig. 5a-c). Within the fields each point corre-

sponds to a dimensionless value of the final aspect ratio. Scaling these values with maximum values of

H,/R, results in contour plots of maximum aspect

ratios h/r of lava flows. For purposes of description we separate the extrusions into three classes. Spines

are extrusions with a height greater than radius, or

h/r > 1. We define domes in the same way as Moriya (19781, whose data are summarized with

values from Mount St. Helens by Swanson and

Holcomb (19901, wherein domes have 1 > h/r > 0.2.

Flows have h/r < 0.2. Basalts are predicted to form flows with small

aspect ratios, reaching a maximum of about 0.01 (Fig. 5a), consistent with observations of subaerial

flows. For submarine flows, cooling rates are greater and should lead to greater values of the viscosity contrast and hence of pressure number. Subaerial

andesite lavas (Fig. 5b) are expected to form flows, domes and spines, although the bulk of the field corresponds to flows. Again, this is generally consistent with observations. Andesites occur commonly as lava flows but sometimes form domes, such as at Mount Lamington and Bezymianny (Swanson et al., 1987); we know of no andesitic lava spines but such features erode rapidly. The rhyolite field (Fig. 5~) shows a wide range of aspect ratios, but the majority

M.V. Stasiuk, C. Jaupart/ Journal of Volcanology and Geothermal Research 78 (1997) 31-50 41

I

BASALT

Fig. 5. (a) Contour map of maximum (scaled) aspect ratios of

basalt flows versus 5 and jj, in field defined by ranges in

parameter values in Table 1. (bl Contour map of andesite aspect

ratios. (c) Contour map of rhyolite aspect ratios. In all cases,

arrow shows expected variation with time during eruption due to

increasing lava viscosity. Stippled regions in (bl and (c) delineate

conditions to produce domes.

of the range corresponds to domes and spines. Low aspect ratio rhyolite flows are relatively rare on the Earth, to the extent that their interpretation as lavas rather than rheomorphic ignimbrites has been contro- versial (Bonnichsen and Kauffman, 1987; Manley, 1996). It is worth noting in particular the sensitivity of flow dimensions to the values of jj and V, and hence to eruption conditions, visible on the andesite and rhyolite fields of Fig. 5.

4.3. Eruption transitions

The only parameter likely to change during eruptions is the bulk viscosity of the lava. As discussed earlier, this is normally observed to increase with time. It is known that a cooling viscous flow fed at a constant rate can be described by a constant bulk viscosity whose value is greater for lower eruption rate (Stasiuk et al., 1993a,b; Bercovici, 1994). Thus, one approach is to approximate the case of a lava flow fed by steadily decreasing eruption rate by a sequence of constant rate episodes, and this can explain the temporal behaviour of the bulk viscosity of lava flows. Our relatively simple model is not designed to take account of such changes explicitly, but we may show that the basic results remain valid and in fact are reinforced when the lava viscosity increases with decreasing eruption rate. In the con- text of the model, the eruption rate is most sensitive to a thick, strongly coupled flow with high bulk viscosity relative to the magma. In nature, such conditions are expected to lead to increasing bulk viscosity, which would in turn further increase the flow thickness and decrease more the eruption rate. This feedback can be illustrated using the results of our model shown in Fig. 5. From Eqs. (2.10a) and (2.10b) any increase in the viscosity contrast will produce equal changes in both j and 5. Example trends on the fields of aspect ratios in Fig. 5 show how the predicted final aspect ratio will evolve during an eruption. For subaerial basalt the effect is small because the trend may be nearly parallel to the contours. In nature, transitions to compound flow behaviour could increase the aspect ratio. For andesite and rhyolite the effects of increasing viscosity are likely to be great, potentially driving an eruption

42 M.V. Stasiuk, C. Jaupart/Joumal of Volcanology and Geothermal Research 78 (1997) 31-50

from producing a flow toward a dome or spine. In particular, the change from dome to spine is very sharp and takes place over a change of lava viscosity of only a factor of two. The sensitive feedback between bulk viscosity, eruption rate and flow thickness for conditions of strong coupling implies that there can be rapid transitions from dome to spine growth, accompanied by equally rapid variations in eruption rate and bulk behaviour of the extrusion. Such eruption transitions have been documented. During the 1902- 1903 eruption of Mt. Pelee (Lacroix, 1904), a lava dome grew from approximately May to November 1902, at which point a spine grew from the top of the dome and continued for the next year. We note that the exponential spine growth described by Eq. (B.2) works well for the Mt. Pelee spine (Jaupart and Allegre, 1991). Some of the eruptive episodes of Mount St. Helens between 1980 and 1986 also exhibited transitions from dome (or lobe) extrusion to spine growth (Swanson et al., 1987). Typically, the spines emerged from the sur- faces of extrusive lobes. A spine about 10 m high grew on 12 April 1981, the final day of a three-day eruptive episode, and another grew to a height of about 60 m from 23 February 1983 to 21 April 1983. Similarly, the dome complex growing in the Soufriere Hills of Montseirat from December 1995 to Septem- ber 1996 exhibited repeated growth of spines.

4.4. Overpressures

In principle, for a given lava flow, one may use the model to derive constraints on the characteristics

Table 2 Lava extrusions and implied overpressures in the subsurface

of the volcanic system at depth, such as conduit length, chamber volume and overpressure. In prac- tice, estimates of the chamber size and depth are very sensitive to poorly constrained variables such as the conduit radius. For example, we have performed calculations with parameters available in the literature for the 1980-1986 dome building phases of Mount St. Helens, and found that they are consistent with a chamber lying at a depth of several kilometres but the specific value may change by more than one order of magnitude if the average conduit radius is varied within a reasonable range. The model, however, can provide reliable estimates of the overpressures which drive eruptions. Within our physical framework, the eruption stops as a result of the pressure at depth being balanced by the weight of lava, and hence the final lava thickness over the vent gives a direct estimate of the overpressure remaining in the plumbing system. For domes and spines fed from typical magma chambers (i.e., V-=K l), overpressures at the beginning and end of extrusion are almost equal (Fig. 4e).

Height data from a number of terrestrial lava extrusions are shown in Table 2, along with the implied initial chamber overpressures. The overpressures were calculated from A P,,, = pgh, where h is the lava height over the vent and p is lava density, and where we have assumed that all the lava densities are approximately 2000 kg/m3. This density corresponds to a dacite with about 20% by volume gas bubbles. For Mount St. Helens from 1980-1986, the dacite lava dome was erupted episodically and so there are two categories: one for the total height

Extrusion

Obsidian Dome

Mount St. Helens 1980-1986 (total)

incremental Montagne Pelee 1902-1903

Soufribre de St. Vincent 1979

Mt. Lamington 1951-1953

various (142) in Japan maximum

top 5%

top 11%

Height over vent

Cm) - 100

261

20-40 -300

434

- 550

760

>500

>400

Overpressure

(MPa) -2

-5

0.4-0.8 -6

-9

- 11

- 15

> 10

>8

Reference

Westrich et al. (1988)

Swanson and Holcomb (1990)

Lacroix ( 1904)

Huppert et al. (1982)

Taylor (1958)

Moriya (1978)

M. V. Stasiuk, C. Jaupart/ Joumal of Volcanology and Geothermal Research 78 (1997) 31-50 43

achieved by the dome over the six years of growth and an incremental category for the increases in height for each episode. The group of 142 domes in Japan, summarized by Moriya (19781, are indicated in terms of the maximum dome height measured, the 5% of domes with the greatest heights and the highest 11%.

Tait et al. (1989) discuss likely values of the chamber overpressure required to initiate eruptions by tensile fracturing of the country rock, and quote a range of 0.5-8 MPa. Calculated values for a number of cases are listed in Table 2 and most fall within this range. Some larger values of overpressure are found in several cases (Table 2). One explanation is that they are a result of buoyancy, i.e., a contrast in density between country rock and magma in the conduit. This effect is shown explicitly by Eq. (2Sb) to be identical to an effective value of overpressure which may be higher than the maximum initial overpressure values discussed so far. The data of Table 2 suggest that magma buoyancy may be an important and common driving force.

A large number of spectacular lava domes have been found on Venus (Pavri et al., 1992). These domes are large volume (25-3400 km31 extrusions with great thicknesses (average about 700 m> and lateral extents (average radius about 12 km). By the definitions used earlier these extrusions are lava flows, but we use the accepted term. The conditions and flow dimensions on Venus offer a good example of how the model results can be used when only simple observations have been made. The elevated surface temperature (500°C) implies a relatively low value of the viscosity contrast, favouring extensive flows and low values of both the pressure and volume numbers. Witbin the framework of our model, this conclusion is also consistent with the combina- tion of large volumes and low aspect ratios. There- fore, these eruptions were negligibly affected by the weight of the lava at the surface. Yet, the great thickness of the domes guarantees large pressures on the vent. The requirement of small jJ along with great thickness implies very high values of overpressure. For example, for a (probably unrealistic) minimum lava density of order 1000 kg/m3, a maximum 13 of about 0.1 gives a minimum (average)overpressure of about 100 MPa (Eq. (2.1Ob)). Such a large value of overpressure is not consistent with the

mechanical properties of rocks, and hence must be regarded as an effective value which includes the effects of magma buoyancy (see Eq. (2.5a) and Eq. (2.5b) above). This suggests that dome eruptions on Venus are dominated by buoyancy. The driving pressure resulting from buoyancy is proportional to the magma source depth and the density contrast between magma and crust. If the density contrasts on Venus are similar to those on Earth, then the implica- tion for Venus is that the plumbing system feeding the dome magmas extends to great depths. One possibility is that the dome magmas themselves arise from great depths and do not involve crustal magma chambers, and hence may not be silicic in composition. Another possibility, consistent with the interpretation of the domes as silicic in composition, is that they arise from shallow, evolved chambers during replenishment by large volumes of primitive magma in turn sourced at great depth.

For strong coupling conditions, i.e., in dome and spine eruptions, the final magma chamber overpressure remains close to its initial value. This implies that the chamber is not exhausted of magma. A small increase in the pressure difference between the chamber and surface could provoke renewed activity. This could happen as a result of an increase of chamber overpressure or the removal of the extrusion, for example by gravitational collapse. These concepts are consistent with the observed long-lived, high-frequency episodic behaviour of dome and spine eruptions, well-documented examples of which are the emplacement of the composite lava dome on Mount St. Helens from 1980 to 1986 (Swanson and Holcomb, 1990) and the 1995-1996 activity on Montserrat.

5. Discussion

The present model takes a simple approach by leaving out features like explicit time-dependent effects of cooling, and topography. Inclusion of these aspects would improve the numeric details, but is unlikely to substantially change the conclusions. The effects of cooling are likely to increase with time during any particular eruption. Provided this occurs slowly and can be approximated as an increase in lava bulk viscosity, these time-dependent effects can


10-z i 1 1

1 o-2 IO.1 100

pressure number p

IO’

Fig. 6. Final dimensionless height versus jj for V = 0.001 and different values of final eruption rate. Variation of final flux by one order of magnitude does not change the essential behaviour: a sharp increase in flow thickness with increasing jj in vicinity of ji= 1.

be approximated by a succession of quasi-steady state solutions, as discussed in Section 4.3. Another aspect not explicitly treated is the halting of the flow front, a part of the emplacement history which we approximated simply by selecting a small final eruption rate @t,> = 0.01 at which the emplacement and eruption were taken as complete. The choice is not critical to our conclusions since, for simple flow units, it has a weak effect on the solutions. This is demonstrated by calculations made for three values of the final eruption rate, shown in Fig. 6. For the case of compound flow development, where individ- ual flow lobes stop while the eruption continues, the lava flow cannot be modelled as a viscous-gravity current but the conceptual framework of our model still applies. The numeric solutions could be redone using approximate expressions for flows with strong surface crust (Griffiths and Fink, 1993), but would produce qualitatively similar results to cases with large viscosity contrasts since a strong crust reduces lateral spreading and favours thick accumulations of lava. A surface crust may further increase vent pressure because the lava may become overpressured (Griffiths and Fink, 1993). This would enhance the coupling, allowing even a small extrusion to influence the dynamics. A lava flow with strongly non- Newtonian rheology, for example possessing a large yield stress, would again enhance the coupling, in the same way as a surface crust.

Our simplification of flat topography is reasonable for silicic lavas on very shallow slopes, on

caldera and wide crater floors, parts of the sea floor, and the surface of Venus. Steep slopes could play an important role in some situations but cannot be studied within this framework. The effects of coupling would be strongly reduced for vents located on steep volcano flanks because the lava may not accu- mulate over the vent. On the other hand, coupling could be increased where lava ponds over the vent, such as within pyroclastic cones or summit craters.

Simplifications aside, the model provides insight into the large-scale features of the volcanic system which feeds lava flows. Consider for example evolved lavas with large viscosities. High aspect ratio, small volume spines and episodically emplaced domes are expected from shallow, small chambers. These conditions correspond well to the small, shallow chambers in convergent continental settings such as the Cascade Range of western North America (e.g., Iyer et al., 1990). Large volume, extensive rhyolite flows are expected from a deep and large source. These conditions may correspond to the Snake River plain volcanic province (Bonnichsen and Kauffman, 1987). Similarly, the morphology of underwater basaltic flows could be sensitive moni- tors of the upper mantle magma extraction systems below mid-ocean ridges.

The model results show how the various parameters of a volcano system are coupled. This may shed some light on the debate about which parameter, total erupted volume or average eruption rate, controls the length and morphology of lava flows (e.g., Malin, 1980). For basic compositions there is now convincing field (Rowland and Walker, 1990) and laboratory (Fink and Griffiths, 1990, 1992) evidence that eruption rate is a main control on the transition to compound flow behaviour for given cooling conditions. For simple flows of broad compositional range the situation is less clear, probably due to the large number of subtly linked variables in the problem, and hence attributing lava lengths to a single variable is likely a severe oversimplification. For example, in our model the lateral extent of a lava flow correlates with its volume, since to produce large volumes requires thin flows, but the correlation reveals an effect rather than a cause. The total erupted volume is a function of many parameters, and it is how these balance, in groups such as the pressure and volume numbers, which dictates the dynamics.

M.V. Stasiuk, C. .Jaupart/ Journal of Volcanology and Geothermal Research 78 (1997) 31-50 45

Our model predicts a similar correlation between auerage eruption rate and flow extent (evident from Table l), since higher average eruption rates correspond to eruptions which are never ‘plugged’ by thick accumulations of lava over the vent and hence small values of the pressure number, accompanied by large flow extents. In nature, the eruption rate has a further significance through its interplay with cooling, since a volume of lava erupted over a shorter period has felt less the influence of cooling and spreads farther. Griffiths and Fink (1993) and Fink and Bridges (1995) have established that the morphology and aspect ratio of flows depends on the interplay between cooling and both eruption rate magnitude and variation, a result supported by experiments with fluids which solidify upon cooling (Fink and Griffiths, 1990, 1992) and those which become more viscous (Stasiuk et al., 1993a,b). The importance of eruption rate adds weight to our model results, where it is found to vary significantly and to be dependent, under certain conditions, on the surface activity.

6. Conclusions

A simplified fluid mechanical model of lava eruptions couples the growth and spreading of lava over a vent to the flow of magma from an overpressured chamber and up a cylindrical conduit. The model reveals the physical conditions which are likely to produce extensive lava flows, domical extrusions, or high lava spines. A scaling analysis shows that the eruptive behaviour depends on two dimensionless numbers: a volume number V representing a ratio of a scale volume to the total volume available for eruption, and a pressure number j5 representing a ratio of a scale lava weight to the initial magma chamber overpressure. Analytical results show that there are two lava eruption regimes separated by a sharp transition and dependent mainly on the value of the pressure number. These results show that small differences in the pressure number may lead to large variations in flow shape. Large, deep and strongly overpressured magma chambers erupting into narrow conduits produce voluminous, long-lived eruptions and low aspect ratio flows. The opposite conditions lead to relatively high initial eruption

rates of viscous lava, leading to thick domes and spines which act to stop the eruption before the lava can spread far. In (5, j5) space, an intermediate domain of small extent separates the flow and spine regimes, and delineates the conditions for the formation of domes. Small changes in the eruption conditions, such as the contrast in viscosity between the magma at depth and the lava, can produce sudden transitions from dome extrusion to spine growth which are consistent with a variety of observations of dome eruptions. The model results can be used to show that magma buoyancy plays an important part in driving many terrestrial dome eruptions. The results also suggest that the lava domes on Venus were driven primarily by buoyancy sourced at great depths. Two possible explanations are that the domes are not silicic in composition, or were erupted when their crustal source reservoirs were impinged upon by primitive magma arising from great depth.

Acknowledgements

This work was initiated while MVS was at the University of Bristol with a Commonwealth Scholar- ship and under the supervision of R.S.J. Sparks, and completed in France with support from a NSERC PDF.

Appendix A. Integral method solution for radial viscous-gravity currents

A.]. Theory

For the case of an isoviscous, Newtonian current with negligible surface tension and inertia which flows across a horizontal surface, two integral expressions of continuity can be given. First, the rate of volume change is:

r,(t) rh( r,t) dr (A4

where h(r, t) is the height of the current surface as a function of the coordinates radius and time, and TN(t) is the radial position of the flow front. Altema- tively, this volume change may be written for a cylinder of radius a(t). Liquid flows out of this

46 M.V. Stasiuk, C. Jaupart/Joumal of Volcanology and Geothenal Research 78 (1997) 31-50

cylinder through its vertical edges and contributes to an increase of its height. For convenience the radius a(t) is a constant fraction E of the radial extent so that a(t) = erN(t). The volume change may be written as follows:

(A4

where U(r, z, t) is the vertical profile of horizontal velocity and z. is the vertical coordinate. For thin currents, and away from the flow outlet, this velocity is largely horizontal and one may use the lubrication approximation (e.g., Batchelor, 1967). The vertical profile of horizontal velocity for a viscous-gravity current was given by Huppert (1982) as:

U(r,z,t) = +zz{2h(r,t) -z} (A.3) I

where p is the current density, ,FL, is dynamic viscosity (both assumed to be large compared to values for the ambient fluid) and g is gravity. At long time a(t) will be large relative to the outlet radius, guaran- teeing that the velocity profile Eq. (A.31 will be accurate. At early times a(t) is comparable to the outlet radius, but it will be shown that using Eq. (A.31 leads to a reasonable description of the flow height for all time and hence is sufficient for our purposes.

We impose a dimensionless form for the profile shape. Huppert (1982) found that the first term of an expansion about the flow front is:

QJ) 1 _ r Y -= h,(t) [ 1 rt-J(t) (A4

with y = l/3. This describes well the shapes of flows in the laboratory, except close to r = Okee

below). Using the velocity profile Eq. (A.31 and the shape Eq. (A.41 in Eq. (A.l) and Eq. (A.2), integrating and combining to eliminate h,(t) yields:

(ASa)

where

@= y(y + l)4(y+ 2)4 E(l - e)3Y-2

87r3 (1+2e) ( ASb)

and

OR= (2e-4&2ye- y2e2- ye2)

(1- l )(l +2e) (A&)

Eq. (ASa) can be solved using the integrating factor Vv4, provided that 0 = - 8. For h,(O) = 0, r,(O) =

R, the radial extent of the current is:

@Pg

1

l/8

54(t) = / ‘V( to4 dt’ + R;

ww 0 (A4

and solving for the height from the volume, derived by integrating over the shape Eq. (A.41, gives:

h (t) = b+ l>(r+2) v(t) @P&T 0 2rr [ 3t-N t )

1 -l/4

x I

i’( tf)4 dt’ + R; 0

(A.7)

In order to determine 0 from Eq. (A.5b), it is sufficient to define one of y or E. A single shape leads to accurate results over a range of flow conditions. Laboratory experiments (see below) suggest a value of about y = l/3, which implies E = 0.894 and 0 = 0.381. A comparison between our solution for the flow radius for R, = 0 and V(t) = Qta and the similarity solution shows a difference of less than 2% for y = l/3 for all time and any volume supply intherangeOsa<l.

In laboratory flows which are erupted at a constant rate, there is a central bulge above the feeding pipe. This bulge typically has a radius about three times that of the feeding pipe and its height represents approximately 25% of the total flow height. It may be approximated by a Gaussian curve which is added to the shape Eq. (A.41, yielding:

h( r,t) 1 _ r ’ -= h,(t) [ 1 5df)

r +ii, l-- I 1 e -(r/R,)%

IN(t)

(A.8)

where h,(t) is the flow central height minus the bulge height and k, = 0.25 and ii, = 0.2 are dimensionless parameters representing the height and width of the central bulge. It is possible to use the shape A.8 to recast the problem into the same one as that

M.V. Stasiuk, C. Jaupart/ Joumal of Volcanology and Geothermal Research 78 (1997) 3X-50 41

d 0.8 . c 0.6

0.4

031291.1 (3IJ.2) 031291.2 (31.1.6) 031291.2 (43.1.6)

0 I I

0 0.2 0.4 0.6 0.8 1

r I IN

Fig. 7. Dimensionless flow shapes from the syrup experiments with superimposed profile shape functions from Eq. (A.4) (fine) and Eq. (AIt) (bold). The value of the ratio R, /RN chosen for the shape Eq. (A.8) is an average of the values for the various data profiles. Legend gives the experiment number followed by the profile radius and height in cm.

for the flow without a bulge. This is done by defin- ing a ‘fictive’ flow whose volume is identical to that of the real flow (with the bulge) but which has the shape of a flow without a bulge. Thus the fictive flow has a form given by Eq. (A.4) but with t&) replaced by the fictive radius_?,(t) and h(r,r) replaced by the fictive height h(r,t). By comparing volume expressions, the fictive radius is given by:

&(t) =P-&)(l +z(t))“* (A.9a)

where

(A.9b)

and /I(t) = &/r,(r). The fictive flow satisfies the same conditions and equations of continuity Eq. (A.l) and Eq. (A.2) as the real flow, and has the same velocity profile Eq. (A.31 and shape Bq. (A.4). Hence the expressions for FN(t> and h,(t) are given by Eq. (A.@ and Eq. (A.7), respectively, except that the initial (fictive) extent is now given by Eqs. (A.9a) and (A.9b) at t = 0 rather than R,. The real flow radius is given by the transcendental Eqs. (A.9a) and (A.9b) and the real height by Eq. (A.8) for r = 0, given explicitly in Bq. (2.6).

A.2. Laboratory experiments

The integral method solution provides an accurate description of flow radius, as already stated through the comparison with the similarity solution of Hup- pert (1982). However, flow height data have not previously been published and since this is an important element of the coupled flow problem experiments were done to test the estimates from the integral method solution. The experimental fluid was Lyle’s Golden Syrup, a glucose syrup with Newto- nian rheology, and the apparatus is that described in Stasiuk et al. (1993a). The syrup was pumped at constant flow rate onto a levelled perspex surface, into air, where it spread axisymmetrically. The flow

Table 3 Experimental parameters for radial syrup flows

Experiment number: 310193 031291.1 031291.2 3081 041291 280891

Supply rate Q (units of lOF6 m3 s- ‘) Viscosity r.r, (Pas) Density p (kg me3 1 Pipe radius R, (m) Expt duration (s) Final height (m) Final radius (m> Height scale a H, (m) Time scale b T. (s)

1.78 0.87 2.69 91.0 68.7 68.0

1438 1438 1438 0.02 0.01 0.01

3180 1800 2170 0.017 0.012 0.016 0.465 0.27 0.46 0.010 0.008 0.011 2.33 0.927 0.397

0.8 8.41 8.33 61.3 46.0 78.5

1438 1438 1438 0.02 0.02 0.02

6730 1100 1382 0.012 0.019 0.022 0.52 0.53 0.55 0.008 0.013 0.015 3.84 0.61 0.705

a From F$ (2.7~). b From Eq. (2.7e).

48 M. V. Stasiuk, C. Jaupart/ Journal of Volcanology and Geothermal Research 78 (1997) 31-50

I ’ F

0.1 I- 0.1 1 10 100 1000 IO4

fls Fig. 8. Dimensionless central height of experimental currents and theoretical predictions as functions of the dimensionless time. The dashed line is an estimate from the similarity solution, determined by integrating the shape EXq. (A.4) for volume and substituting the similarity solution for the flow radius and constant supply rate; it starts discontinuously at zero time. Also shown are the integral method solution without central bulge (fine) and with central bulge (hold). Central height has ken normalized by the scale of Eq. (2.7~1, and time by the scale of Eq. (2.7e).

shapes were determined at various times by measur- ing the thickness at a number of radii using probes lowered from the tank lid (error + 2 mm). In three experiments a scale was mounted centrally over the feeder pipe, and measurements had an estimated error of + 1 mm. The main experimental parameters are shown in Table 3.

Some flow profile shapes from the experiments are compared with the shapes Eq. (A.4) and Eq. (A.81 in Fig. 7, where the radius and height have been normalized by rN, the front radius, and h,, the flow height in the absence of the bulge. Eq. (A.41 provides a good fit to the data away from the feeding pipe, and Eq. (A.8) gives a good fit for all positions.

To allow simple comparisons between the experiments and theory, the height and time scales of Eq. (2.9) were used to normalize the measured and predicted values. The values of the height and time scales for the experiments are shown in Table 3. The measured ( normalized) flow centre heights and the curves representing the theoretical predictions of Eq. (A.7) and that derived using the shape with a Gauss- ian bulge are shown in Fig. 8. The data show that the central height of the flows increased rapidly in the approximate interval 0 < t < 10, and was nearly constant for t > 100. The integral method solution which includes a central bulge agrees with the data within

error for t I 2 and 2 > 100. For dimensionless time 2 <i < 100 there is a greater discrepancy but the agreement is still quite good, the difference between experiment and theory is maximum at small time (25% at i G 3) and decreases with increasing time. For our eruption model the integral method provides a reasonable continuous description of the flow thickness.

A.3. Cooled flows

Cooling acts to introduce changes of viscosity and to generate a brittle carapace. The transition between brittle and viscous behaviour for the cooled outer shell of a lava flow is not understood, and we restrict our discussion to viscous effects. Cooling acts to modify the shape of the lava flow and to modify its spreading rate. As shown by Stasiuk et al. (1993a,b) and Bercovici (19941, a cooled flow develops an almost vertical leading edge, and has small thickness variations. The shape may be described by function (A4) with a small value of exponent y (Stasiuk et al., 1993a,b). For our present purposes, the change of shape is not important because it does not affect significantly the volume balance. However, changes of bulk viscosity are important because they affect the rate of spreading and hence, for a given volume, change the aspect ratio of the flow. These effects can be understood by considering time changes of bulk viscosity. As already discussed in the main text, they are not likely to affect the various spreading regimes found.

Appendix B. Asymptotic cases

For certain conditions Eq. (2.8) and Eq. (2.9) have solutions convenient for comparison with the numerical results. For sufficiently large V and/or jj the eruption stops before spreading occurs, corresponding to the growth of spines, Integrating Eq. (A.8) for the flow volume yields V = k?z for ? = 1, where k = 1.8. Using this result in Eq. (2.9) and solving for h gives:

;= &[I -exp{-(v+P/k)i}] (B.1)

M. V. Stasiuk, C. Jaupafl/ Journal of Volcanology and Geothermal Research 78 (1997) 31-50 49

As j + 0, Eq. (2.8) becomes a = 1 - 57 and leads to the exponential volume relationship v = (1 - e-“j/F derived by Scandone (19791, Wadge (1981) and Stasiuk et al. (1993b). Using this volume in Eq. (2.9) yields the flow height for an exponen- tially decreasing eruption rate as:

h z (1 - ewEi)

v

I

52 + de-C’ _ 3e-2ZT + oe-38i _ &-4Bi _ g

jj”( 1 - e-“)

1 -l/4

SA (B.2)

Although the flow height will be described by Eq. (B.2) for short and intermediate times, the flux and height at long time are influenced by the height since eventually j% N 1 - TV, as demonstrated by the numerical solutions.

References

Batchelor, C&K., 1967. An Introduction to Fluid Dynamics. Cam-

bridge University Press, Cambridge, 615 pp.

Bercovici, D., 1994. A theoretical model of cooling viscous

gravity currents with temperature-dependent rheology. Geo-

phys. Res. L&t., 21: 1117-1180.

Blake, S., 1981. Volcanism and the dynamics of open magma

chambers. Nature, 289: 783-785.

Botichsen, B. and Kauffman, D.F., 1987. Physical features of

rhyolite lava flows in the Snake River Plain volcanic province,

southwestern Idaho. Geol. Sot. Am., Spec. Pap., 212: 119-

14.5.

Borgia, A. and Linneman, S.R., 1989. On the mechanisms of lava

flow emplacement and volcano growth: Arenal, Costa Rica.

In: J.H. Fink (Editor), Lava Flows and Domes: Emplacement

Mechanisms and Hazard Implications. IAVCEI Proc. Vol-

canal., 2: 208-243.

Borgia, A., Linneman, S., Spencer, D., Morales, L.D. and Andre,

J.B., 1983. Dynamics of lava flow fronts, Arenal Volcano,

Costa Rica. J. Volcanol. Geotherm. Res., 19: 303-329.

Bottinga, Y., Weill, D. and Richet, P., 1982. Density calculations

for silicate melts. I. Revised methods for aluminosilicate com-

positions. Geochim. Cosmovhim. Acta, 46: 909-919.

Crisp, J., Cashman, K.V., Bonim, J.A., Hougen, S.B. and Pieri,

DC., 1994. Crystallization history of the 1984 Mauna Loa

lava flow. J. Geophys. Res., 99: 7177-7198.

Endo, E.T., Dzurisin, D. and Swanson, D.A., 1990. Geophysical

and observational constraints for ascent rates of dacitic magma

at Mount Saint Helens. In: M.P. Ryan (Editor), Magma Trans-

port and Storage. Wiley, Chichester, pp. 317-334.

Ewart, J.A., Voight, B. and Bjomsson, A., 1990. Dynamics of

Krafla caldera, north Iceland: 1975-1985. In: M.P. Ryan

(Editor), Magma Transport and Storage. Wiley, Chichester,

pp. 225-276.

Fink, J.H. and Bridges, N.T., 1995. Effects of eruption history and

cooling rate on lava dome growth. Bull. Volcanol., 57: 229-

239.

Fink, J.H. and Griffiths, R.W., 1990. Radial spreading of

viscous-gravity currents with solidifying crust. J. Fluid Mech.,

221: 485-509.

Fink, J.H. and Griff~ths, R.W., 1992. A laboratory analogue study

of the surface morphology of lava flows extruded from point

and line sources. J. Volcanol. Geotherm. Res., 54: 19-32.

Fink, J.H. and Zimbelman, J., 1989. Longitudinal variations in

rheological properties of lavas: Puu 00 basalt flows, Kilauea

Volcano, Hawaii. In: J.H. Fink (Editor), Lava Flows and

Domes: Emplacement Mechanisms and Hazard Implications.

IAVCEI Proc. Volcanol., 2: 157-173.

Griffiths, R.W. and Fink, J.H., 1993. Effects of surface cooling on

the spreading of lava flows and domes. J. Fluid Mech., 252:

667-702.

Guest, J.E., Kilbum, C.R.J., Pinkerton, H. and Duncan, A.M.,

1987. The evolution of lava flow-fields: observations of the

1981 and 1983 eruptions of Mount Etna, Sicily. Buil. Vol-

canal., 49: 527-540.

Hausback, B.P., 1987. An extensive, hot, vapor-charged rhyo-

dacite flow, Baja California, Mexico. Geol. Sot. Am. Spec.

Pap., 212: 111-118.

Huppert, H.E., 1982. The propagation of two-dimensional and

axisymmetric viscous gravity currents over a rigid horizontal

surface. J. Fluid Mech., 121: 43-58.

Huppett, H.E., Shepherd, J.B., Sigurdsson, H. and Sparks, R.S.J.,

1982. On lava dome growth with application to the 1979 lava

extrusion of the Soufribre of St. Vincent. J. Volcanol.

Geotherm. Res., 14: 199-222.

Huppert, H.E. and Sparks, R.S.J., 1988. The generation of granitic

magmas by intrusion of basalt into continental crust. J. Petrol.,

29: 599-624.

Ishihara, K., 1990. Pressure sources and induced ground deforma-

tion associated with explosive eruptions at an andesitic vol-

cano: Sakurajima volcano, Japan. In: M.P. Ryan (Editor),

Magma Transport and Storage. Wiley, Chichester, pp 335-356.

Iyer, H.M., Evans, J.R., Dawson, P.B., Stat&r, D.A. and Achauer,

U., 1990. Differences in magma storage in different volcanic

environments as revealed by seismic tomography: silicic vol-

canic centres and subduction-related volcanoes. In: M.P. Ryan

(Editor), Magma Transport and Storage. Wiley, Chicester, pp.

293-316.

Jaupart, C., 1991. Effects of compressibility on the flow of lava.

Bull. Volcanol., 54: 1-9.

Jaupart, C. and Allegre, C.J., 1991. Gas content, eruption rate and

instabilities of eruption regime in silicic volcanoes. Earth

Planet. Sci. Lett., 102: 413-429.

Jaupart, C. and Tait, S., 1990. Dynamics of eruptive phenomena.

50 M. V. Stasiuk, C. Jaupart/ Journal of Volcanology and Geothermal Research 78 (1997) 31-50

In: J. Nicholls and J.K. Russell (Editors), Modem Methods in

Igneous Petrology. Rev. Mineral., 24: 213-238.

Lacroix, A., 1904. La Montagne PelCe et ses Eruptions. Masson et

Cie, Paris, 662 pp.

Lees, J.M., 1992. The magma system of Mount St. Helens:

non-linear high-resolution P-wave tomography. J. Volcanol.

Geotherrn. Res., 53: 103-116.

Lejeune, A. and Richet, P., 1995. Rheology of crystal-bearing

silicate melts: an experimental study at high viscosities. J.

Geophys. Res., 100: 4215-4229.

Lipman, P.W. and Banks, N.G., 1987. Aa flow dynamics, Mauna

Loa 1984. U.S. Geol. Surv., Prof. Pap., 1350: 1527-1568.

Malin, M.C., 1980. Lengths of Hawaiian lava flows. Geology, 8:

306-308.

Manley, C.R., 1996. Physical volcanology of a voluminous rhyo-

lite lava flow: The Badlands lava, Owyhee Plateau southwest-

em Idaho. J. Volcanol. Geotherm. Res., 71: 129-153.

Moran, SC., 1994. Seismicity at Mount St. Helens, 1987-1992:

evidence for repressurization of an active magmatic system. J.

Geophys. Res., 99: 4341-4354.

Moriya, I., 1978. Morphology of lava domes. Bull. Dep. Geogr.

Komazawa Univ. Jpn, 14: 55-69 (in Japanese).

Murray, J.B., 1990. High-level magma transport at Mount Etna

volcano, as deduced from ground deformation measurements.

In: M.P. Ryan (Editor), Magma Transport and Storage. Wiley,

Chichester, pp. 357-384.

Naranjo, J.A., Sparks, R.S.J., Stasiuk, M.V., Moreno, H. and

Ablay, G.J., 1992. Morphological, structural and textural vari-

ations in the 1988-1990 andesite lava of Lonquimay Volcano,

Chile (38” S). Geol. Mag., 129: 657-678.

Pavri, B., Head, J.W.I., Klose, K.B. and Wilson, L., 1992. Steep-

sided domes on Venus: characteristics, geologic setting, and

eruption conditions from Magellan data. J. Geophys. Res., 97:

13,445-13,478.

Pinkerton, H. and Stevenson, R.J., 1992. Methods of determining

the rheological properties of magmas at sub-liquidus tempera-

tures. J. Volcanol. Geotherm. Res., 53: 47-66.

Rowland, S.K. and Walker, G.P.L., 1990. Pahoehoe and aa in

Hawaii: Volumetric flow rate controls the lava structure. Bull.

Volcanol., 52: 615-628.

Sakimoto, S.E.H. and Zuber, M.T., 1995. The spreading of vari-

able-viscosity axisymmetric radial gravity currents: applica-

tions to the emplacement of Venusian ‘pancake’ domes. J.

Fluid Mech., 301: 65-77.

Scandone, R., 1979. Effusion rate and energy balance of Paricutin

eruption., 1943-1952) Michoacan, Mexico. J. Volcanol.

Geothenn. Res., 6: 49-56.

Shaw, H.R., 1972. Viscosities of magmatic silicate liquids: an

empirical method of prediction. Am. J. Sci., 272: 870-893.

Sparks, R.S.J. and Pinkerton, H., 1978. Effect of degassing on

rheology of basaltic lava. Nature, 276: 385-386.

Stasiuk, M.V., Jaupart, C. and Sparks, R.S.J., 1993a. Influence of

cooling on lava flow dynamics. Geology, 21: 335-338.

Stasiuk, M.V., Jaupart, C. and Sparks, R.S.J., 1993b. On the

variations of flow rate in non-explosive lava eruptions. Earth

Planet. Sci. Lett., 114: 505-516.

Swanson, D.A., Dzurisin, D., Holcomb, R.T., lwatsubo, E.Y.,

Chadwick Jr., W.W., Casadevall, T.J., Ewert, J.W. and He-

liker, C.C., 1987. Growth of the lava dome at Mount St.

Helens, Washington (USA), 1981-1983. Geol. Sot. Am.,

Spec. Pap., 212: 1-16.

Swanson, D.A. and Holcomb, R.T., 1990. Regularities in growth

of the Mount St. Helens Dacite Dome 1980-1986. In: J.H.

Fink (Editor), Lava Flows and Domes: Emplacement Mecha-

nisms and Hazard Implications. IAVCEI Proc. Volcanol., 2:

3-24.

Tait, S., Jaupart, C. and Vergniolle, S., 1989. Pressure, gas

content and eruption periodicity of a shallow crystallising

magma chamber. Earth Planet. Sci. Lett., 92: 107-123.

Taylor, G.A.M., 1958. The 1951 eruption of Mount Lamington,

Papua. Aust. Bur. Min. Resour. Bull., 38: 117.

Touloukian, Y.S., Judd, W.R. and Roy, R.F., 1981. Physical

Properties of Rocks and Minerals. McGraw-Hill, New York,

NY.

Wada, Y., 1994. On the relationship between die width and

magma viscosity. J. Geophys. Res., 99: 17,743-17,755.

Wadge, G., 1981. The variation of magma discharge during

basaltic eruptions. J. Volcanol. Geotherm. Res., 11: 139-168.

Walker, G.P.L., 1967. Thickness and viscosity of Etnean lavas.

Nature, 213: 484-485.

Westrich, H.R., Stockman, H.W. and Eichelberger, J.C., 1988.

Degassing of rhyolitic magma during ascent and emplacement.

J. geophys. res., 93: 6503-6511.

Wilson, C.J.N., 1986. Pyroclastic flows and ignimbrites. Sci.

Prog. Oxf., 70: 171-207.

Date post:	11-Nov-2023
Category:	Documents
Upload:	ipgp
View:	1 times
Download:	0 times

Lava flow shapes and dimensions as reflections of magma system conditions

Documents