- Email: [email protected]
Effects of flow density on the dynamics of dilute pyroclastic density currents
Available online at www.sciencedirect.com R
Journal of Volcanology and Geothermal Research 132 (2004) 269^281 www.elsevier.com/locate/jvolgeores
E¡ects of £ow density on the dynamics of dilute pyroclastic density currents Sarah E. Nield a; , Andrew W. Woods b a b
Earth Sciences, Bristol University, Wills Memorial Building, Queen’s Road, Bristol BS8 1RJ, UK BP Institute for Multiphase Flow, Cambridge University, Madingley Rise, Cambridge CB3 0EZ, UK Received 30 June 2003; accepted 8 August 2003
Abstract Dilute turbulent pyroclastic density currents may have densities as large as 10^100 kg m33 . Laboratory experiments using dense gases show that for such currents, the speed of the head of the current depends on the density difference between the current and the surrounding air in a complex way. Here, a model of a dilute pyroclastic density current generated by a short-lived volcanic eruption is developed to account for this effect. The model predicts that the head of the current travels more slowly than in existing models, leading to shorter run-out distances and smaller dynamic flow pressures. The effects are most significant for dense pyroclastic density currents such as those produced by dome collapse events. We apply the new model to the eruptions of Mount St. Helens 1980, Taupo 1800 years ago and Soufrie're Hills 1997. 3 2003 Elsevier B.V. All rights reserved. Keywords: pyroclastic £ow; gravity current; mathematical model; sedimentation; non-Boussinesq
1. Introduction The dynamics of pyroclastic density currents is fundamental to understanding the origin of volcanic £ow deposits and for hazard assessment. Density di¡erences between dilute suspensiondriven pyroclastic £ows, which run along the ground, and the surrounding air play a key role in determining the £ow behaviour, since it is this density di¡erence that drives the £ow. Here we focus on the dynamics of £ows pro-
* Corresponding author. Fax: +0044-117-9253385. E-mail addresses: [email protected] (S.E. Nield), [email protected] (A.W. Woods).
duced by relatively short-lived violent explosions. A number of well-documented pyroclastic density currents have been produced by such explosions. For example the eruption of Taupo 1800 years ago (Wilson and Walker, 1985; Wilson, 1985; Dade and Huppert, 1996) produced a pyroclastic density current which travelled 80 km in approximately 10 minutes at speeds of 100^200 m s31 , while the lateral blast produced during the eruption of Mount St. Helens on 18 May 1980 travelled 20 km at speeds of 100^150 m s31 (Kie¡er, 1981; Pallister et al., 1992). The dome explosion on 26 December 1997 of Soufrie're Hills Volcano, Montserrat travelled over 4 km at speeds of 80^ 100 m s31 (Sparks et al., 2002). In models of short-lived dilute £ows, one of the
0377-0273 / 03 / $ ^ see front matter 3 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0377-0273(03)00314-7
VOLGEO 2719 6-4-04
270
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
key inputs which determine the speed of the current is the so-called Froude number condition at the head of the current (Huppert and Simpson, 1980). For small di¡erences in density between the ambient and current £uids, the Froude number is obtained from the ratio between excess hydrostatic pressure and stagnation pressure, and is thus de¢ned as: uF ffi Fr ¼ pffiffiffiffiffiffi g0 h
ð1Þ
where gP = gvb/bref is reduced gravity, bc and ba denote current and ambient density with vb = bc 3ba , g is the gravitational constant, uF denotes frontal velocity of the current and h denotes the height of the current. This ratio would be constant for similar £ows. For the case of a lock exchange £ow release mechanism, pffiffiffi Benjamin (1968) proved analytically that Fr ¼ 2, however this has been found to overpredict the frontal velocity due to e¡ects such as dissipation of mechanical energy by turbulence in the head and viscous drag. In a series of ¢nite-release gravity current laboratory experiments involving small vb, Huppert and Simpson (1980) measured the Froude number value to be 1.2. The established models for the Froude number at the head of a gravity current are based on the assumption that the current satis¢es the Boussinesq approximation, which requires the density of the current to be similar to that of the ambient £uid. However, owing to the presence of dense particles in dilute pyroclastic density currents, the density of the £ow may be much greater than that of the surrounding environment. Gro«belbauer et al. (1993, hereinafter referred to as GFB) investigated the motion of such non-Boussinesq gravity currents using a range of gases whose density ratio varied from 1 to 20. In this case, as vb increases, de¢nition (1) becomes inappropriate since for values of vb s bref , the reduced gravity exceeds gravitational acceleration. GFB therefore rede¢ned the Froude number as: uF Frð b Þ ¼ pffiffiffiffiffi gh
ð2Þ
where the Froude number varies with b*, which is de¢ned as (Yih, 1965):
b ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Mv b M : bc þ ba
ð3Þ
They ran lock-exchange experiments using different ambient and current densities, producing currents of two di¡erent regimes : currents which were dense relative to the ambient, which ran along the bottom of the tank; and currents which were less dense than the ambient £uid, which ran under the top of the tank. Hereinafter we use the terms ‘‘dense current’’ and ‘‘light current’’ respectively to describe these two regimes. The data from these experiments established that the Froude number diverges substantially from the classical Boussinesq value (Huppert and Simpson, 1980), as a function of density di¡erence between the current and ambient £uids. Motivated by these results and their potential importance in calculating the speed of a pyroclastic density current, we now develop a revised model of the propagation of such currents, which includes the parameterised version of Froude number for non-Boussinesq currents. To set the scene, ¢rst we investigate the e¡ect on the dynamics of purely gaseous gravity currents. We then extend the model to account for the e¡ects of particles being transported by and sedimented from the £ow. We use this to predict the propagation speed, the dynamic pressure and the associated deposit shape of dilute pyroclastic density currents. Finally we apply our non-Boussinesq model to the previously documented eruptions of Montserrat, Taupo and Mount St. Helens, and compare results with those predicted by the Boussinesq model.
2. Non-Boussinesq e¡ects on a homogeneous current Fig. 1 illustrates the experimental data of GFB for the value of the Froude number as a function of b* for the case of in¢nite ambient depth. After running lock exchange experiments varying vb for two depths of ambient gas, they extrapolated their data to obtain values of b* for in¢nite ambient depth. They found that the £ow speed varies with b* in a di¡erent manner for dense currents com-
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
271
and for light currents:
pffiffiffi Fr1 ð b Þ ¼ 1:2 2 b 30:5 b :
ð7Þ
In order to illustrate the e¡ect of large di¡erences in density on the propagation speed of the current as measured experimentally, (Eqs. 2, 6 and 7), with the speed of propagation of a current of density close to the host £uid, it is convenient to express the Boussinesq frontal velocity condition as a function of b*, g and h. In this case the Froude number in Eq. 2 may be represented by: Fig. 1. Froude number as a function of density variable b*, for non-Boussinesq light and dense currents, compared with the classical Boussinesq model (Eqs. 8 and 4). For the nonBoussinesq cases, experimental data points (Gro«belbauer et al., 1993) are shown, with empirical laws Eqs. 6 and 7. The Boussinesq model strictly only applies to currents of density which is very close to that of the host £uid, that is, for small b*.
pared to light currents. Although the main interest in this paper is in the motion of pyroclastic density currents, in this section, for completeness we illustrate the di¡erences between the propagation speed of very dense and very light currents. We have developed an analytic expression to represent this data and this is also shown in Fig. 1. The lower part of the curves corresponds to small density di¡erences for which the classical Boussinesq model applies. The non-Boussinesq and Boussinesq Froude numbers therefore tend asymptotically to the same value as b*C0. Indeed combining the Boussinesq limit de¢ned as (Huppert and Simpson, 1980): uF Fr ¼ 1:2 ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffi vb gh b ref
ð4Þ
with the non-Boussinesq Froude number condition (2), we ¢nd that as b*C0, the frontal velocity has the value: pffiffiffi pffiffiffiffiffi pffiffiffiffiffiffiffi uF ¼ 1:2 g0 h ¼ 1:2 2 b gh
ð5Þ
We thus ¢nd that the GFB data for dense currents is well-represented by the empirical law: pffiffiffi Frd ð b Þ ¼ 1:2 2 b þ 2:2
b 23 b
sffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffi Mv b M 2 FrB ð b Þ ¼ 1:2 ¼ 1:2 b b ref 13 b 2
by combining Eqs. 2, 3 and 4. The two non-Boussinesq empirical laws (Eqs. 6 and 7) are plotted in Fig. 1 with the Boussinesq empirical law (Eq. 8) for comparison. For both light and dense currents, the conventional Boussinesq model predicts greater values of Froude number and hence current speed than has been found by experimental methods for non-Boussinesq currents. We would expect the frontal velocity of non-Boussinesq currents to be overpredicted when using the Boussinesq approximation for values of b* s 0.2. 2.1. The one-dimensional box model for channelized currents To illustrate the importance of the Froude number parameterisation in gravity current model predictions, we incorporate this new law for frontal velocity into the one-dimensional box model, where variables are assumed to be uniform with respect to height due to turbulent mixing within the current (Huppert and Simpson, 1980). The governing equations are the conservation of volume (Eq. 9), the new Froude number condition (Eq. 10) and a kinematic condition at the front (Eq. 11): V ¼ qt Q ¼ Lh pffiffiffiffiffi uF ¼ Frð b Þ gh
6:6 ;
ð6Þ
ð8Þ
uF ¼
dL dt
VOLGEO 2719 6-4-04
ð9Þ ð10Þ ð11Þ
272
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
where L and h denote the length and height of the current, q the volume of £ow from the source with q constant, and V denotes the volume per unit width. Fr(b*) is de¢ned by the empirical laws (Eqs. 6 and 7). In the case of an instantaneous ¢nite release of £uid, Q = 0; for constant £ux, Q = 1. Eliminating h and solving this system gives the distance propagated as a function of Froude number and time: L¼
3Fr Q þ2
2=3
ðgqÞ1=3 tð Q þ2Þ=3 :
ð12Þ
In the case of an instantaneous release of gas, volume is independent of time so we have Q = 0 and hence L is proportional to t2=3 . By introducing non-dimensional variables: pffiffiffi 3 g ^ ¼ pffiffiffi 1 ^t ¼ L L and t; 2 q1=4 qFrB ð b Þ2=3
ð13Þ
where the length is scaled by the Boussinesq Froude number, we obtain the dimensionless governing equation: ^¼ L
FrnB ð b Þ FrB ð b Þ
2=3 ^t2=3 :
ð14Þ
In the case of a constant mass £ux from the source, the total volume V increases linearly with time t so Q = 1 and hence the current length, L, also increases linearly with time. The dimensionless propagation length of a current in the case of a ¢nite instantaneous release is plotted as a function of dimensionless time in Fig. 2 using the non-Boussinesq model, for both (a) relatively dense gas and (b) relatively light gas. Curves are shown for di¡erent values of b*. The Boussinesq solution is plotted (dashed curve) for comparison. The Boussinesq model predicts higher propagation speeds for all density di¡erences, hence the non-Boussinesq and Boussinesq predictions of the current position diverge with time. This overprediction is greater for larger density di¡erences. For a given value of b*, the divergence is greater for a relatively light current, as expected from Fig. 1.
Fig. 2. The distance propagated by channelized currents of di¡erent densities as a function of dimensionless time for an instantaneous ¢nite release of £uid, with length scaled by the Boussinesq Froude number. (A) corresponds to dense gas, (B) light gas. Values of the densities are given by b* = 0.1, 0.25, 0.5, 0.75 and 0.9, as shown on the ¢gures.
2.2. The axisymmetric box model For an axisymmetric current, we describe the frontal velocity in terms of the radial position of the nose of the current which we refer to as the radius of spread, r, given by: dr 1 ¼ Frð b Þ dt r
rffiffiffiffiffi gq Q =2 t : Z
ð15Þ
The curvature is small, owing to the lengthscale of the £ows, therefore locally at the head the £ow is two-dimensional (Huppert and Simpson, 1980, Bonnecaze et al., 1995), and we may use the same
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
value of Froude number as in the channellized case. Eq. 15 may be solved to obtain the radius of spread as a function of time and Froude number: r¼
4 Frð b Þ Q þ2
rffiffiffiffiffi1=2 gq tð Q þ2Þ=4 : Z
ð16Þ
In the case of instantaneous release r is proportional to t1=2 , and for a source with constant volume £ux r is proportional to t3=4 . In Fig. 3 we illustrate the distance travelled as predicted by the non-Boussinesq model as a fraction of that predicted by the Boussinesq model. The variation of this ratio is shown as a function of b* for both a channelized and a radially spreading current. For larger density di¡erences and thus larger b*, the Boussinesq and non-Boussinesq solutions diverge. For example, for b* = 0.5, the di¡erence in the predicted propagation distance is 5^10%, whereas for b* = 0.8 di¡erences of up to 40% arise. This leads to a di¡erence in the predicted £ow speed of 40% and a di¡erence in the prediction of dynamic pressure of over 60% (Eq. 26). Larger di¡erences between the non-Boussinesq and Boussinesq solutions are predicted for channelised £ow than axisymmetric £ow. Although the motion of short-lived pyroclastic density currents is strongly controlled by sedimentation, before extending the model to account for this, it is of interest to consider typical values of
Fig. 3. The ratio between the non-Boussinesq and Boussinesq predictions of frontal position as a function of the density variable, for light and dense currents. Dashed lines correspond to axisymmetric £ow; solid lines to channellized £ow.
273
b* to gain insight into the possible magnitude of these non-Boussinesq e¡ects. For typical dilute pyroclastic density currents, in which the volume of gas exceeds 80^90% volume, the bulk density will be smaller than 100 kg m33 so the density variable lies in the range 0.2 6 b* 6 0.97. GFB obtained data for values of b* up to 0.95, therefore we extrapolate their data for currents of density 100 kg m33 . However, such extrapolation is only required during initial spreading since the density of the current evolves quickly as particles are sedimented (Section 3). A value of b* as high as 0.97 suggests that the £ow velocity may be reduced by a factor as large as 2 (Figs. 1 and 3). Such a reduction in £ow speed has important implications for the prediction of total run-out distance, deposit geometry and also the dynamic pressure associated with the £ow. This forms the subject of the remainder of the paper.
3. Dilute pyroclastic density currents 3.1. Theoretical model We now extend the model (Eqs. 3 and 15) to examine the propagation of radially spreading sedimenting pyroclastic density currents along near-horizontal terrain following the instantaneous release of material, for example, after the collapse of an eruption column or lava dome. We allow the density to evolve as the current propagates owing to sedimentation of particles from the £ow. This model represents a development of the work of Dade and Huppert (1996) and Woods (2000), who examined the propagation of Boussinesq pyroclastic density currents in which they implicitly assumed that density di¡erences were small. For the more concentrated pyroclastic density currents, the initial density may be such that b* exceeds the experimental range of GFB, we therefore extrapolate their data. However, such extrapolation is only required for a short time during the initial propagation since the density of the current decreases with sedimentation and hence time.
VOLGEO 2719 6-4-04
274
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
Key predictions of the model will be the £ow speed and hence dynamic pressure, and the runout distance of the current when there has been su⁄cient sedimentation that the £ow becomes less dense than air (Sparks et al., 1993; Dade and Huppert, 1996; Bursik and Woods, 1996) and thus lifts o¡ the ground. We consider currents which are only a few hundred metres deep so that the pressure in the current is close to atmospheric, hence internal strati¢cation of the £ow due to compression of the gas may be neglected (Woods, 2000). For simplicity, to demonstrate the non-Boussinesq e¡ects on behaviour, we consider a current with particles whose fall speed may be characterised by a constant mean value. However, the model may be readily extended to account for multiple particle sizes (Bursik and Woods, 1996). We examine the £ow for cases with di¡erent gas mass fractions, eruption masses and mean settling velocities of particles. The bulk density of the current is a function of the densities of the particles, bs , and the interstitial gas, bg , and of the mass fraction n of the interstitial gas. We approximate the current density using typical values of the parameters (Woods, 2000):
bc ¼
n 13n þ bg bs
31 W
bg : n
ð17Þ
This approximation is valid for n s O(10 ), which relates to s 80^90% volume. Here, bg = p/ Rg Tg , where Rg and Tg denote the gas constant and temperature of the erupting gas, and p the pressure in the current. n evolves as the particles sediment from the £ow. The density variable (Eq. 3) may be expressed in the form: 33
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M13V nM ; b ¼ 1þVn
Woods, 2000). Therefore we may neglect air entrainment as long as particle fallout has a greater e¡ect on the £ow density than air entrainment (Sparks et al., 1993). We consider the radial spreading of a gas and particle mixture. In addition to Eq. 18, the equations governing the £ow include the conservation of mass: V g ¼ Z r2 h;
the Froude number condition at the nose of the current: uF ¼
pffiffiffiffiffi dr ¼ Frð b Þ gh; dt
ð20Þ
where Fr(b*) is de¢ned by the empirical laws Eqs. 6 and 7 and an additional equation to describe the evolution of the density of the current due to sedimentation of particles. This was de¢ned by Einstein (1968), and expressed as a function of gas mass fraction (Dade and Huppert, 1996; Woods, 2000) has the form: dn vs nð13nÞ ¼ dt h
ð21Þ
where vs denotes the mean settling velocity of the particles. Solving for the radius in terms of the gas mass fraction, we obtain the expression: "
pffiffiffi 4V g 3=2 g r ¼ rðn0 Þ þ Z 3=2 vs 4
R
n
n0
Frðn0 Þ dn0 n0 ð13n0 Þ
#1=4
ð22Þ
where r(n0 ) and n0 denote the radius of the £ow and bulk gas mass fraction immediately after the blast and n denotes the gas mass fraction when the front has reached radius r. Also, we obtain an expression relating the time after release to the radial extent of the £ow:
ð18Þ t¼
where V = Rg Tg /Ra Ta denotes the density ratio between the current and ambient. We neglect the e¡ects of entrainment of air, and also heat transfer to the ground over the relatively short timescale of the £ow. Therefore we take V to be constant. Neglecting air entrainment is an approximation, but still allows a good model for sedimentation from aqueous solutions in the laboratory (Bonnecaze et al., 1995; Gladstone and
ð19Þ
rffiffiffiffiffiffiffiffiffiR rðtÞ 0 Z r dr0 gV g rð0Þ Frðr0 Þ
ð23Þ
Note that the Froude number is a function of radius since the current density and hence b* is dependent upon radius. 3.2. Flow dynamics We now use this model to examine the dependence of the non-Boussinesq and Boussinesq pre-
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
275
dictions of speed and dynamic pressure on several key parameters, including the eruption mass mer = 108 ^1013 kg, the settling velocity vs = 0.1^10 ms31 and the mass of air entrained into the eruption column. The interstitial gas in the current is typically made up of water vapour originating from the magma and air which may mix into the £ow during fountain collapse. We represent the bulk gas constant by: Rg ¼
nv Rv þ ðn0 3nv ÞRa n0
ð24Þ
where Rv = 462 J kg31 K31 and Ra = 287 J kg31 K31 denote the gas constant for water vapour and air respectively. In the case of a collapsing fountain, the gas in the current is a mixture of air entrained during the collapse and magmatic gas. In this case we assume a constant initial gas mass fraction of vapour, nv = 0.03, and an initial bulk gas mass fraction in the range n0 = 0.03^0.2. In the case of dome collapse, we take n0 to have values in the range 0.001^0.01, where the gas composed only of air for gravitational dome collapse and only of magmatic vapour for explosive dome collapse (Woods et al., 2002). The temperature of the current is found by applying the conservation of thermal energy to the erupted mass and, in the case of a collapsing fountain, to any air entrained into the £ow. In making these calculations, we assume that the particles are small enough to attain thermal equilibrium with the interstitial gas over a timescale which is much shorter than that of the £ow. This leads to the expression : Tc ¼
mer T er ½ner C v þ ð13ner ÞC s þ ma C a T a mer ½ner C v þ ð13ner ÞC s þ ma C a
ð25Þ
for the temperature of the current, where subscripts er , a , v , s and c denote eruption, air, vapour, solid and bulk current, and C denotes speci¢c heat, where the speci¢c heat for vapour at eruption temperature, ash solids and air are approximately 2000, 1000 and 993 J kg31 K31 respectively. The temperature of the erupted mass is typically 1000 K, and we take surrounding air temperature to be 273 K.
Fig. 4. (A) the ratio between non-Boussinesq and Boussinesq radius of spread as a function of time, for various initial mass fractions. n0 = 0.001, 0.005 represent the case of dome collapse, with nv = 0; the other curves represent collapsing eruption fountains with an initial water vapour mass fraction of 0.03 and initial bulk vapour and air gas densities as marked. (B) the evolution of non-Boussinesq (solid) and Boussinesq (dashed) velocity and density ratio as a function of radial spread, for a current with initial gas mass fraction of 0.1, eruption mass of 109 kg and mean particle settling velocity of 0.1 m s31 . The density di¡erence between the ambient and current drives the £ow and hence determines frontal velocity.
3.3. Propagation speed Fig. 4A compares the non-Boussinesq prediction for distance propagated as a function of time with that predicted by the Boussinesq model
VOLGEO 2719 6-4-04
276
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
(Dade and Huppert, 1996). The distance propagated is up to 35% less than the Boussinesq prediction for the case of a collapsing fountain (n0 v 0.01), and up to 65% in the case of dome collapse (n0 = 0.001, 0.005). The di¡erences are largest for currents with smaller initial volumes of air in the £ow, since the current density is larger for smaller gas content and thus the value of the densimetric Froude number which governs the £ow dynamics is relatively smaller than the value of the Boussinesq Froude number. Fig. 4 also shows that the di¡erences in the model predictions are largest at smaller times and shorter distances from the source, before substantial sedimentation has occurred. Here the current is most dense, leading to the largest di¡erences between Boussinesq and non-Boussinesq Froude numbers and hence frontal velocity. In the case of an initial gas mass fraction of 0.1, for non-Boussinesq propagation speeds of 80 m s31 , the Boussinesq model overpredicts the £ow speed by about 30^40% (Fig. 4B) leading to comparable di¡erences in the predicted run-out distance. 3.4. Lift-o¡ conditions Lift-o¡ occurs when the density of the current evolves to the density of the surrounding air. This is governed by: (1) the rate of sedimentation of dense particles; and (2) the amount of gas initially mixed into the £ow. In the present model, neutral buoyancy occurs when b* = 0, and the current then lifts o¡. This requires that n = 1/V = bg /ba . Fig. 5A shows that there is a critical value of the initial gas mass fraction, ncrit W0.1, such that lift-o¡ length is maximised. For any given lift-o¡ radius smaller than this maximum there are two possible initial gas mass fractions, related to different lift-o¡ times. In the case of a gas mass fraction smaller than the critical value, less air is entrained into the eruption column. This leads to a larger concentration of particles in the £ow, which causes faster sedimentation, but also a hotter current due to reduced cooling by the air. For a gas mass fraction larger than ncrit , the reduced concentration of particles causes the £ow density to decrease towards the neutral buoyancy con-
Fig. 5. (A) Lift-o¡ radius of the current as a function of initial gas mass fraction. There is a critical value of the initial gas mass fraction, at ncrit = 0.1, for which runout distance is maximised. For smaller n0 there is a higher concentration of particles which causes faster sedimentation, whereas for larger n0 the current has a smaller initial bulk density and hence lifts o¡ after less sedimentation. Eruption mass is ¢xed at 109 kg. (B) Non-Boussinesq (solid) and Boussinesq (dashed) lift-o¡ radius as a function of eruption mass. In this ¢gure, the volatile mixture for each curve (labelled with n0 ) has composition na = 0.12 nv = 0.03, na = 0 nv = 0.03 and na = 0.001 nv = 0. Mean particle settling velocity is ¢xed at 0.1 m s31 .
dition following a smaller amount of sedimentation. In comparison with the predictions of the Boussinesq model, the non-Boussinesq model predicts a shorter runout distance because the associated slower £ow has more time to sediment. The di¡erence between non-Boussinesq and Boussinesq lifto¡ radii is largest for small gas mass fractions and
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
277
hence denser £ows, such as those associated with dome collapse (see Fig. 5A,B). It is worth noting from Fig. 5B that any attempt to estimate eruption mass based on runout distance and mean particle size is subject to uncertainty by over an order of magnitude owing to uncertainty in the initial gas content. 3.5. Deposit geometry The reduction in £ow speed as predicted by the non-Boussinesq model leads to more sedimentation at each point in the current. As a result, the currents form deposits with aspect ratio that is slightly higher than predicted by the Boussinesq model, while the total volume of sediment deposited remains the same. Fig. 6 illustrates this e¡ect. Curves are shown for Boussinesq (dashed) and non-Boussinesq (solid) deposit shapes. The di¡erence using the Boussinesq model is largest in the case of high mean settling velocity. This e¡ect would be larger for channelized £ows, in which sediment is deposited over a much smaller area. The ¢gure also illustrates the very strong impact of mean particle size on deposit geometry which, in fact, dwarfs any di¡erences between the Boussinesq and non-Boussinesq models.
Fig. 7. Evolution and magnitude of dynamic pressure with distance from the source, using the Boussinesq (dashed) and non-Boussinesq (solid) models. For initial bulk gas mass fractions of 0.15 and 0.03, the Boussinesq model predicts pressure up to three times as high as the non-Boussinesq model.
3.6. Dynamic pressure For hazard analysis of pyroclastic density currents, the impact of the current on obstacles in the path of the £ow may be quanti¢ed in terms of the dynamic pressure, pD , of the £ow (Valentine, 1998). This depends on the density bc and velocity uF of the £ow, according to: 1 pD ¼ b c u2F : 2
ð26Þ
In Fig. 7 we compare the prediction of dynamic pressure as calculated using the Boussinesq and non-Boussinesq models. The ¢gure illustrates that for a range of gas mass fractions, the dynamic pressure is overpredicted by the Boussinesq model by as much as a factor of 3. The di¡erences in predictions are greatest within a proximal distance from the source.
4. Discussion
Fig. 6. Non-Boussinesq (solid) and Boussinesq (dashed) sediment deposit pro¢les for currents composed of particles with di¡erent mean settling velocities. Eruption mass is 109 kg and initial gas mass fraction 0.1. Settling velocities are in m s31 .
We have developed a model of a sedimenting dilute pyroclastic density current, accounting for e¡ects of the large density di¡erence between the £ow and the overlying air on the propagation speed of the £ow, applicable to currents with densities up to 100 kg m33 . The model predictions
VOLGEO 2719 6-4-04
278
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
imply slower velocities than for the classic Boussinesq model. This leads to the prediction of greater sedimentation at each point in the £ow and hence shorter run-out distances. Also, the dynamic pressure of the £ow is smaller, impacting models of damage associated with pyroclastic £ows (Valentine, 1998). The model calculations indicate that non-Boussinesq e¡ects increase with the density di¡erence between the current and the ambient. We therefore expect the e¡ects to be more signi¢cant in denser pyroclastic density currents such as those produced by explosive dome collapse, for which the initial gas content is small. For example, recent models of the explosive dome collapse on 26th December 1997 at Soufrie're Hills, Montserrat suggest that a volume of about 5U107 m3 of the andesite dome was disrupted, and the associated blast cloud formed a dilute pyroclastic density current. The deposits from this were laden with ¢ne particles (Sparks et al., 2002) thus we may use the model to illustrate the propagation given these initial conditions. The initial density has been estimated to have been of order 100 kg m33 (Woods et al., 2002). The current spread from the dome towards the coastline, approximately 4 km from the dome, with speeds which have been estimated to be of order 80^100 m s31 (Sparks et al., 2002). Based on the deposit geometry, the £ow may be modelled as a radially spreading current, con¢ned by the topography to spread over an angle of approximately 70‡. In Fig. 8A, we illustrate the predictions of the non-Boussinesq model for the velocity of the head of this current as a function of radial distance from the dome, where, for simplicity, we assume that the material in the head of the current was released rapidly from the fragmenting dome compared to the travel time of the £ow. This is a somewhat simpli¢ed end-member calculation, given that the theoretical predictions by Woods et al. (2001) suggest that the dome break-up took a few tens of seconds, but serves to illustrate the importance of non-Boussinesq e¡ects. Curves are shown for mean particle settling speed of 1 m s31 , with a total mass of 1011 kg, initial radius 500 m and initial density 100 kg m33 corresponding to initial
Fig. 8. (A) Non-Boussinesq (solid) and Boussinesq (dashed) velocity as a function of frontal position for the 1997 eruption at Soufrie're Hills, Montserrat. Curves correspond to initial densities of 100 kg m33 . (B) corresponding dynamic pressure as a function of frontal position. Eruption mass is 1011 kg, initial radius 500 m and mean particle settling speed 1 m s31 .
gas mass fraction of 2U1033 , the gas here taken to be magmatic vapour. The curves show that a mean settling velocity of about 1 m s31 is appropriate, using the non-Boussinesq model. For comparison, the dashed lines illustrate the predictions of the Boussinesq model for the same initial conditions. It is seen that the Boussinesq model predicts velocities for such a dense current of about an order of magnitude greater than the non-Boussinesq model. For such large velocities the £ow is expected to be compressible and not fully describable by the model.
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
Fig. 8B compares the dynamic pressure associated with the head of the current for both Boussinesq and non-Boussinesq models. The nonBoussinesq model predicts pressure to be smaller by several orders of magnitude. In the case of more explosive events of higher gas content, the di¡erences between the Boussinesq and non-Boussinesq models are expected to be smaller (Section 3), since the density of the £ow decreases as gas content increases. For example, during the lateral blast event at Mount St. Helens on 18 May 1980, about 2.5U1011 kg of magma was erupted from the shallow cryptodome with a gas mass fraction of about 0.04. This trav-
279
elled about 15^20 km before lifting o¡. Fig. 9A,B shows that the Boussinesq model predicts a current speed about 60% greater than the non-Boussinesq model, and dynamic pressure to be about 150% greater. During the Taupo event (Dade and Huppert, 1996), around 4U1012 kg was erupted with a gas mass fraction of about 0.01 and spread 80 km radially. Again, the Boussinesq model predicts speeds greater by approximately 60%, as shown in Fig. 9C. Channellized currents, such as those which initially begin as a radially spreading current but then develop into a channellized current owing
Fig. 9. Non-Boussinesq (solid) and Boussinesq (dashed) velocity as a function of frontal position and corresponding dynamic pressure as a function of frontal position for Mount St. Helens 1980 (A,B) and for the Taupo eruption 1800 years ago (C,D). Mean particle settling speed is 1 m s31 in both cases.
VOLGEO 2719 6-4-04
280
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
to the topography, for example, running into a valley, would be overpredicted by a greater amount than a purely axisymmetric current. In this paper we have examined the e¡ects of adapting a model to include a non-Boussinesq parameter, which has important e¡ects on the £ow. The model gives good insight into the dynamics and deposit of pyroclastic density £ows, however it would be desirable to incorporate other e¡ects such as the slope of the topography and grain-size distribution.
Acknowledgements ‘‘This research supported in part by the Engineering and Physical Sciences Research Council (Nield) and by the BP Institute, Cambridge.’’
Appendix Notation
Frl FrnB g gP h L LΠma mer n n0 ner nv p pD q r Ra Rg Rg Rv t
“t Ta Tc Ter Tg uF V Vg vs Q vb
V ba bc bg bref bs b*
dimensionless time ambient air temperature current temperature eruption temperature erupting gas temperature frontal velocity of the current volume per unit width 2-dimensional gas volume mean particle settling velocity power of time, relating to volume £ux density di¡erence between the current and ambient £uids density ratio between the volcanic gas and ambient ambient density current density interstitial gas density reference density solids density density parameter
References
The following symbols are used in this paper: Ca Cs Cv Fr FrB Frd
Appendix Notation (Continued).
speci¢c heat of air speci¢c heat of the volcanic solids speci¢c heat of vapour Froude number densimetric Froude number for a Boussinesq current densimetric Froude number for a dense (heavy) current densimetric Froude number for a light current densimetric Froude number for a non-Boussinesq current gravitational constant reduced gravity height of the density current current length dimensionless current length mass of entrained air eruption mass gas mass fraction initial gas mass fraction gas mass constant of the erupted gas gas mass fraction of vapour pressure dynamic pressure volume of £ow radius of the axisymmetric current gas constant of the ambient (air) gas constant of the erupting gas bulk gas constant gas constant for vapour time
Benjamin, T.B., 1968. Gravity currents and related phenomena. J. Fluid Mech. 31, 209^249. Bonnecaze, R.T., Hallworth, M.A., Huppert, H.E., Lister, J.R., 1995. Axisymmetric particle-driven gravity currents. J. Fluid Mech. 416, 187^195. Bursik, M.I., Woods, A.W., 1996. The dynamics and thermodynamics of large ash £ows. Bull. Volcanol. 58, 175^193. Dade, W.B., Huppert, H.E., 1996. Emplacement of the Taupo ignimbrite by a dilute turbulent £ow. Nature 381, 509^512. Einstein, H.A., 1968. Deposition of suspended particles in a gravel bed. J. Hydraul. Div. Proc. Am. Soc. Civil Eng. 94, 1197^1205. Gladstone, C., Woods, A.W., 2000. On the application of box models to particle-driven gravity currents. J. Fluid Mech. 294, 93^121. Gro«belbauer, H.P., FannelTp, T.K., Britter, R.E., 1993. The propagation of intrusion fronts of high density ratios. J. Fluid Mech. 250, 669^687. Huppert, H.E., Simpson, J.E., 1980. The slumping of gravity currents. J. Fluid Mech. 9, 785^799. Kie¡er, S.W., 1981. Fluid dynamics of the May 18 blast at Mount St. Helens. In: Lipman, P., Mullineaux, D. (Eds.), The 1980 Eruptions of Mount St. Helens, Washington. Geol. Survey Prof. Pap. 1250, 379^400. Pallister, J.S., Hoblitt, R.P., Crandell, D.R., Mullineaux, D.R., 1992. Mount St. Helens a decade after the 1980 eruptions magmatic models, chemical cycles and a revised hazards assessment. Bull. Volcanol. 54, 126^146. Sparks, R.S.J., Barclay, J., Calder, E.S., Herd, R.A., Komorowski, J.-C., Luckett, R., Norton, G.E., Ritchie, L., Voight, B., Woods, A.W., 2002. Generation of a debris avalanche
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281 and violent pyroclastic density current on 26 December (Boxing Day) at Soufrie're Hills Volcano, Montserrat. In: Druitt, T.H., Kokelaar, B.P. (Eds.), The Eruption of Soufrie're Hills Volcano, Montserrat, from 1995 to 1999. Memoir of the Geological Society of London. Sparks, R.S.J., Bonnecaze, R.T., Huppert, H.E., Lister, J.R., Hallworth, M.A., Mader, H., Phillips, J., 1993. Sedimentladen gravity currents with reversing buoyancy. Earth Planet. Sci. Lett. 114, 243^257. Valentine, G.A., 1998. Damage to structures by pyroclastic £ows and surges, inferred from nuclear weapons e¡ects. J. Volcanol. Geotherm. Res. 87, 117^140. Wilson, C.J.N., 1985. The Taupo eruption, New Zealand. 2. the Taupo ignimbrite. Philos. Trans. R. Soc. Lond. A 314, 229^310.
281
Wilson, C.J.N., Walker, G.P.L., 1985. The Taupo eruption, New Zealand. 1. general-aspects. Philos. Trans. R. Soc. Lond. A 314, 199^228. Woods, A.W., 2000. Dynamics of hazardous volcanic £ows. Philos. Trans. R. Soc. London A 358, 1705^1724. Woods, A.W., Sparks, R.S.J., Ritchie, L.I., Batey, J., Gladstone, C., Bursik, M.I., 2002. The explosive decompression of a pressurized volcanic dome: the 26 December 1997 collapse and explosion of Soufrie're Hills Volcano, Montserrat. In: Druitt, T.H. and Kokelaar, B.P. (Eds.), The eruption of Soufrie're Hills Volcano, Montserrat, from 1995 to 1999. Geological Society, London, Memoirs, pp. 21, 409^434. Yih, C.-S., 1965. Dynamics of Non-Homogeneous Fluids. Macmillan.
VOLGEO 2719 6-4-04
Journal of Volcanology and Geothermal Research 132 (2004) 269^281 www.elsevier.com/locate/jvolgeores
E¡ects of £ow density on the dynamics of dilute pyroclastic density currents Sarah E. Nield a; , Andrew W. Woods b a b
Earth Sciences, Bristol University, Wills Memorial Building, Queen’s Road, Bristol BS8 1RJ, UK BP Institute for Multiphase Flow, Cambridge University, Madingley Rise, Cambridge CB3 0EZ, UK Received 30 June 2003; accepted 8 August 2003
Abstract Dilute turbulent pyroclastic density currents may have densities as large as 10^100 kg m33 . Laboratory experiments using dense gases show that for such currents, the speed of the head of the current depends on the density difference between the current and the surrounding air in a complex way. Here, a model of a dilute pyroclastic density current generated by a short-lived volcanic eruption is developed to account for this effect. The model predicts that the head of the current travels more slowly than in existing models, leading to shorter run-out distances and smaller dynamic flow pressures. The effects are most significant for dense pyroclastic density currents such as those produced by dome collapse events. We apply the new model to the eruptions of Mount St. Helens 1980, Taupo 1800 years ago and Soufrie're Hills 1997. 3 2003 Elsevier B.V. All rights reserved. Keywords: pyroclastic £ow; gravity current; mathematical model; sedimentation; non-Boussinesq
1. Introduction The dynamics of pyroclastic density currents is fundamental to understanding the origin of volcanic £ow deposits and for hazard assessment. Density di¡erences between dilute suspensiondriven pyroclastic £ows, which run along the ground, and the surrounding air play a key role in determining the £ow behaviour, since it is this density di¡erence that drives the £ow. Here we focus on the dynamics of £ows pro-
* Corresponding author. Fax: +0044-117-9253385. E-mail addresses: [email protected] (S.E. Nield), [email protected] (A.W. Woods).
duced by relatively short-lived violent explosions. A number of well-documented pyroclastic density currents have been produced by such explosions. For example the eruption of Taupo 1800 years ago (Wilson and Walker, 1985; Wilson, 1985; Dade and Huppert, 1996) produced a pyroclastic density current which travelled 80 km in approximately 10 minutes at speeds of 100^200 m s31 , while the lateral blast produced during the eruption of Mount St. Helens on 18 May 1980 travelled 20 km at speeds of 100^150 m s31 (Kie¡er, 1981; Pallister et al., 1992). The dome explosion on 26 December 1997 of Soufrie're Hills Volcano, Montserrat travelled over 4 km at speeds of 80^ 100 m s31 (Sparks et al., 2002). In models of short-lived dilute £ows, one of the
0377-0273 / 03 / $ ^ see front matter 3 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0377-0273(03)00314-7
VOLGEO 2719 6-4-04
270
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
key inputs which determine the speed of the current is the so-called Froude number condition at the head of the current (Huppert and Simpson, 1980). For small di¡erences in density between the ambient and current £uids, the Froude number is obtained from the ratio between excess hydrostatic pressure and stagnation pressure, and is thus de¢ned as: uF ffi Fr ¼ pffiffiffiffiffiffi g0 h
ð1Þ
where gP = gvb/bref is reduced gravity, bc and ba denote current and ambient density with vb = bc 3ba , g is the gravitational constant, uF denotes frontal velocity of the current and h denotes the height of the current. This ratio would be constant for similar £ows. For the case of a lock exchange £ow release mechanism, pffiffiffi Benjamin (1968) proved analytically that Fr ¼ 2, however this has been found to overpredict the frontal velocity due to e¡ects such as dissipation of mechanical energy by turbulence in the head and viscous drag. In a series of ¢nite-release gravity current laboratory experiments involving small vb, Huppert and Simpson (1980) measured the Froude number value to be 1.2. The established models for the Froude number at the head of a gravity current are based on the assumption that the current satis¢es the Boussinesq approximation, which requires the density of the current to be similar to that of the ambient £uid. However, owing to the presence of dense particles in dilute pyroclastic density currents, the density of the £ow may be much greater than that of the surrounding environment. Gro«belbauer et al. (1993, hereinafter referred to as GFB) investigated the motion of such non-Boussinesq gravity currents using a range of gases whose density ratio varied from 1 to 20. In this case, as vb increases, de¢nition (1) becomes inappropriate since for values of vb s bref , the reduced gravity exceeds gravitational acceleration. GFB therefore rede¢ned the Froude number as: uF Frð b Þ ¼ pffiffiffiffiffi gh
ð2Þ
where the Froude number varies with b*, which is de¢ned as (Yih, 1965):
b ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Mv b M : bc þ ba
ð3Þ
They ran lock-exchange experiments using different ambient and current densities, producing currents of two di¡erent regimes : currents which were dense relative to the ambient, which ran along the bottom of the tank; and currents which were less dense than the ambient £uid, which ran under the top of the tank. Hereinafter we use the terms ‘‘dense current’’ and ‘‘light current’’ respectively to describe these two regimes. The data from these experiments established that the Froude number diverges substantially from the classical Boussinesq value (Huppert and Simpson, 1980), as a function of density di¡erence between the current and ambient £uids. Motivated by these results and their potential importance in calculating the speed of a pyroclastic density current, we now develop a revised model of the propagation of such currents, which includes the parameterised version of Froude number for non-Boussinesq currents. To set the scene, ¢rst we investigate the e¡ect on the dynamics of purely gaseous gravity currents. We then extend the model to account for the e¡ects of particles being transported by and sedimented from the £ow. We use this to predict the propagation speed, the dynamic pressure and the associated deposit shape of dilute pyroclastic density currents. Finally we apply our non-Boussinesq model to the previously documented eruptions of Montserrat, Taupo and Mount St. Helens, and compare results with those predicted by the Boussinesq model.
2. Non-Boussinesq e¡ects on a homogeneous current Fig. 1 illustrates the experimental data of GFB for the value of the Froude number as a function of b* for the case of in¢nite ambient depth. After running lock exchange experiments varying vb for two depths of ambient gas, they extrapolated their data to obtain values of b* for in¢nite ambient depth. They found that the £ow speed varies with b* in a di¡erent manner for dense currents com-
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
271
and for light currents:
pffiffiffi Fr1 ð b Þ ¼ 1:2 2 b 30:5 b :
ð7Þ
In order to illustrate the e¡ect of large di¡erences in density on the propagation speed of the current as measured experimentally, (Eqs. 2, 6 and 7), with the speed of propagation of a current of density close to the host £uid, it is convenient to express the Boussinesq frontal velocity condition as a function of b*, g and h. In this case the Froude number in Eq. 2 may be represented by: Fig. 1. Froude number as a function of density variable b*, for non-Boussinesq light and dense currents, compared with the classical Boussinesq model (Eqs. 8 and 4). For the nonBoussinesq cases, experimental data points (Gro«belbauer et al., 1993) are shown, with empirical laws Eqs. 6 and 7. The Boussinesq model strictly only applies to currents of density which is very close to that of the host £uid, that is, for small b*.
pared to light currents. Although the main interest in this paper is in the motion of pyroclastic density currents, in this section, for completeness we illustrate the di¡erences between the propagation speed of very dense and very light currents. We have developed an analytic expression to represent this data and this is also shown in Fig. 1. The lower part of the curves corresponds to small density di¡erences for which the classical Boussinesq model applies. The non-Boussinesq and Boussinesq Froude numbers therefore tend asymptotically to the same value as b*C0. Indeed combining the Boussinesq limit de¢ned as (Huppert and Simpson, 1980): uF Fr ¼ 1:2 ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffi vb gh b ref
ð4Þ
with the non-Boussinesq Froude number condition (2), we ¢nd that as b*C0, the frontal velocity has the value: pffiffiffi pffiffiffiffiffi pffiffiffiffiffiffiffi uF ¼ 1:2 g0 h ¼ 1:2 2 b gh
ð5Þ
We thus ¢nd that the GFB data for dense currents is well-represented by the empirical law: pffiffiffi Frd ð b Þ ¼ 1:2 2 b þ 2:2
b 23 b
sffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffi Mv b M 2 FrB ð b Þ ¼ 1:2 ¼ 1:2 b b ref 13 b 2
by combining Eqs. 2, 3 and 4. The two non-Boussinesq empirical laws (Eqs. 6 and 7) are plotted in Fig. 1 with the Boussinesq empirical law (Eq. 8) for comparison. For both light and dense currents, the conventional Boussinesq model predicts greater values of Froude number and hence current speed than has been found by experimental methods for non-Boussinesq currents. We would expect the frontal velocity of non-Boussinesq currents to be overpredicted when using the Boussinesq approximation for values of b* s 0.2. 2.1. The one-dimensional box model for channelized currents To illustrate the importance of the Froude number parameterisation in gravity current model predictions, we incorporate this new law for frontal velocity into the one-dimensional box model, where variables are assumed to be uniform with respect to height due to turbulent mixing within the current (Huppert and Simpson, 1980). The governing equations are the conservation of volume (Eq. 9), the new Froude number condition (Eq. 10) and a kinematic condition at the front (Eq. 11): V ¼ qt Q ¼ Lh pffiffiffiffiffi uF ¼ Frð b Þ gh
6:6 ;
ð6Þ
ð8Þ
uF ¼
dL dt
VOLGEO 2719 6-4-04
ð9Þ ð10Þ ð11Þ
272
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
where L and h denote the length and height of the current, q the volume of £ow from the source with q constant, and V denotes the volume per unit width. Fr(b*) is de¢ned by the empirical laws (Eqs. 6 and 7). In the case of an instantaneous ¢nite release of £uid, Q = 0; for constant £ux, Q = 1. Eliminating h and solving this system gives the distance propagated as a function of Froude number and time: L¼
3Fr Q þ2
2=3
ðgqÞ1=3 tð Q þ2Þ=3 :
ð12Þ
In the case of an instantaneous release of gas, volume is independent of time so we have Q = 0 and hence L is proportional to t2=3 . By introducing non-dimensional variables: pffiffiffi 3 g ^ ¼ pffiffiffi 1 ^t ¼ L L and t; 2 q1=4 qFrB ð b Þ2=3
ð13Þ
where the length is scaled by the Boussinesq Froude number, we obtain the dimensionless governing equation: ^¼ L
FrnB ð b Þ FrB ð b Þ
2=3 ^t2=3 :
ð14Þ
In the case of a constant mass £ux from the source, the total volume V increases linearly with time t so Q = 1 and hence the current length, L, also increases linearly with time. The dimensionless propagation length of a current in the case of a ¢nite instantaneous release is plotted as a function of dimensionless time in Fig. 2 using the non-Boussinesq model, for both (a) relatively dense gas and (b) relatively light gas. Curves are shown for di¡erent values of b*. The Boussinesq solution is plotted (dashed curve) for comparison. The Boussinesq model predicts higher propagation speeds for all density di¡erences, hence the non-Boussinesq and Boussinesq predictions of the current position diverge with time. This overprediction is greater for larger density di¡erences. For a given value of b*, the divergence is greater for a relatively light current, as expected from Fig. 1.
Fig. 2. The distance propagated by channelized currents of di¡erent densities as a function of dimensionless time for an instantaneous ¢nite release of £uid, with length scaled by the Boussinesq Froude number. (A) corresponds to dense gas, (B) light gas. Values of the densities are given by b* = 0.1, 0.25, 0.5, 0.75 and 0.9, as shown on the ¢gures.
2.2. The axisymmetric box model For an axisymmetric current, we describe the frontal velocity in terms of the radial position of the nose of the current which we refer to as the radius of spread, r, given by: dr 1 ¼ Frð b Þ dt r
rffiffiffiffiffi gq Q =2 t : Z
ð15Þ
The curvature is small, owing to the lengthscale of the £ows, therefore locally at the head the £ow is two-dimensional (Huppert and Simpson, 1980, Bonnecaze et al., 1995), and we may use the same
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
value of Froude number as in the channellized case. Eq. 15 may be solved to obtain the radius of spread as a function of time and Froude number: r¼
4 Frð b Þ Q þ2
rffiffiffiffiffi1=2 gq tð Q þ2Þ=4 : Z
ð16Þ
In the case of instantaneous release r is proportional to t1=2 , and for a source with constant volume £ux r is proportional to t3=4 . In Fig. 3 we illustrate the distance travelled as predicted by the non-Boussinesq model as a fraction of that predicted by the Boussinesq model. The variation of this ratio is shown as a function of b* for both a channelized and a radially spreading current. For larger density di¡erences and thus larger b*, the Boussinesq and non-Boussinesq solutions diverge. For example, for b* = 0.5, the di¡erence in the predicted propagation distance is 5^10%, whereas for b* = 0.8 di¡erences of up to 40% arise. This leads to a di¡erence in the predicted £ow speed of 40% and a di¡erence in the prediction of dynamic pressure of over 60% (Eq. 26). Larger di¡erences between the non-Boussinesq and Boussinesq solutions are predicted for channelised £ow than axisymmetric £ow. Although the motion of short-lived pyroclastic density currents is strongly controlled by sedimentation, before extending the model to account for this, it is of interest to consider typical values of
Fig. 3. The ratio between the non-Boussinesq and Boussinesq predictions of frontal position as a function of the density variable, for light and dense currents. Dashed lines correspond to axisymmetric £ow; solid lines to channellized £ow.
273
b* to gain insight into the possible magnitude of these non-Boussinesq e¡ects. For typical dilute pyroclastic density currents, in which the volume of gas exceeds 80^90% volume, the bulk density will be smaller than 100 kg m33 so the density variable lies in the range 0.2 6 b* 6 0.97. GFB obtained data for values of b* up to 0.95, therefore we extrapolate their data for currents of density 100 kg m33 . However, such extrapolation is only required during initial spreading since the density of the current evolves quickly as particles are sedimented (Section 3). A value of b* as high as 0.97 suggests that the £ow velocity may be reduced by a factor as large as 2 (Figs. 1 and 3). Such a reduction in £ow speed has important implications for the prediction of total run-out distance, deposit geometry and also the dynamic pressure associated with the £ow. This forms the subject of the remainder of the paper.
3. Dilute pyroclastic density currents 3.1. Theoretical model We now extend the model (Eqs. 3 and 15) to examine the propagation of radially spreading sedimenting pyroclastic density currents along near-horizontal terrain following the instantaneous release of material, for example, after the collapse of an eruption column or lava dome. We allow the density to evolve as the current propagates owing to sedimentation of particles from the £ow. This model represents a development of the work of Dade and Huppert (1996) and Woods (2000), who examined the propagation of Boussinesq pyroclastic density currents in which they implicitly assumed that density di¡erences were small. For the more concentrated pyroclastic density currents, the initial density may be such that b* exceeds the experimental range of GFB, we therefore extrapolate their data. However, such extrapolation is only required for a short time during the initial propagation since the density of the current decreases with sedimentation and hence time.
VOLGEO 2719 6-4-04
274
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
Key predictions of the model will be the £ow speed and hence dynamic pressure, and the runout distance of the current when there has been su⁄cient sedimentation that the £ow becomes less dense than air (Sparks et al., 1993; Dade and Huppert, 1996; Bursik and Woods, 1996) and thus lifts o¡ the ground. We consider currents which are only a few hundred metres deep so that the pressure in the current is close to atmospheric, hence internal strati¢cation of the £ow due to compression of the gas may be neglected (Woods, 2000). For simplicity, to demonstrate the non-Boussinesq e¡ects on behaviour, we consider a current with particles whose fall speed may be characterised by a constant mean value. However, the model may be readily extended to account for multiple particle sizes (Bursik and Woods, 1996). We examine the £ow for cases with di¡erent gas mass fractions, eruption masses and mean settling velocities of particles. The bulk density of the current is a function of the densities of the particles, bs , and the interstitial gas, bg , and of the mass fraction n of the interstitial gas. We approximate the current density using typical values of the parameters (Woods, 2000):
bc ¼
n 13n þ bg bs
31 W
bg : n
ð17Þ
This approximation is valid for n s O(10 ), which relates to s 80^90% volume. Here, bg = p/ Rg Tg , where Rg and Tg denote the gas constant and temperature of the erupting gas, and p the pressure in the current. n evolves as the particles sediment from the £ow. The density variable (Eq. 3) may be expressed in the form: 33
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M13V nM ; b ¼ 1þVn
Woods, 2000). Therefore we may neglect air entrainment as long as particle fallout has a greater e¡ect on the £ow density than air entrainment (Sparks et al., 1993). We consider the radial spreading of a gas and particle mixture. In addition to Eq. 18, the equations governing the £ow include the conservation of mass: V g ¼ Z r2 h;
the Froude number condition at the nose of the current: uF ¼
pffiffiffiffiffi dr ¼ Frð b Þ gh; dt
ð20Þ
where Fr(b*) is de¢ned by the empirical laws Eqs. 6 and 7 and an additional equation to describe the evolution of the density of the current due to sedimentation of particles. This was de¢ned by Einstein (1968), and expressed as a function of gas mass fraction (Dade and Huppert, 1996; Woods, 2000) has the form: dn vs nð13nÞ ¼ dt h
ð21Þ
where vs denotes the mean settling velocity of the particles. Solving for the radius in terms of the gas mass fraction, we obtain the expression: "
pffiffiffi 4V g 3=2 g r ¼ rðn0 Þ þ Z 3=2 vs 4
R
n
n0
Frðn0 Þ dn0 n0 ð13n0 Þ
#1=4
ð22Þ
where r(n0 ) and n0 denote the radius of the £ow and bulk gas mass fraction immediately after the blast and n denotes the gas mass fraction when the front has reached radius r. Also, we obtain an expression relating the time after release to the radial extent of the £ow:
ð18Þ t¼
where V = Rg Tg /Ra Ta denotes the density ratio between the current and ambient. We neglect the e¡ects of entrainment of air, and also heat transfer to the ground over the relatively short timescale of the £ow. Therefore we take V to be constant. Neglecting air entrainment is an approximation, but still allows a good model for sedimentation from aqueous solutions in the laboratory (Bonnecaze et al., 1995; Gladstone and
ð19Þ
rffiffiffiffiffiffiffiffiffiR rðtÞ 0 Z r dr0 gV g rð0Þ Frðr0 Þ
ð23Þ
Note that the Froude number is a function of radius since the current density and hence b* is dependent upon radius. 3.2. Flow dynamics We now use this model to examine the dependence of the non-Boussinesq and Boussinesq pre-
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
275
dictions of speed and dynamic pressure on several key parameters, including the eruption mass mer = 108 ^1013 kg, the settling velocity vs = 0.1^10 ms31 and the mass of air entrained into the eruption column. The interstitial gas in the current is typically made up of water vapour originating from the magma and air which may mix into the £ow during fountain collapse. We represent the bulk gas constant by: Rg ¼
nv Rv þ ðn0 3nv ÞRa n0
ð24Þ
where Rv = 462 J kg31 K31 and Ra = 287 J kg31 K31 denote the gas constant for water vapour and air respectively. In the case of a collapsing fountain, the gas in the current is a mixture of air entrained during the collapse and magmatic gas. In this case we assume a constant initial gas mass fraction of vapour, nv = 0.03, and an initial bulk gas mass fraction in the range n0 = 0.03^0.2. In the case of dome collapse, we take n0 to have values in the range 0.001^0.01, where the gas composed only of air for gravitational dome collapse and only of magmatic vapour for explosive dome collapse (Woods et al., 2002). The temperature of the current is found by applying the conservation of thermal energy to the erupted mass and, in the case of a collapsing fountain, to any air entrained into the £ow. In making these calculations, we assume that the particles are small enough to attain thermal equilibrium with the interstitial gas over a timescale which is much shorter than that of the £ow. This leads to the expression : Tc ¼
mer T er ½ner C v þ ð13ner ÞC s þ ma C a T a mer ½ner C v þ ð13ner ÞC s þ ma C a
ð25Þ
for the temperature of the current, where subscripts er , a , v , s and c denote eruption, air, vapour, solid and bulk current, and C denotes speci¢c heat, where the speci¢c heat for vapour at eruption temperature, ash solids and air are approximately 2000, 1000 and 993 J kg31 K31 respectively. The temperature of the erupted mass is typically 1000 K, and we take surrounding air temperature to be 273 K.
Fig. 4. (A) the ratio between non-Boussinesq and Boussinesq radius of spread as a function of time, for various initial mass fractions. n0 = 0.001, 0.005 represent the case of dome collapse, with nv = 0; the other curves represent collapsing eruption fountains with an initial water vapour mass fraction of 0.03 and initial bulk vapour and air gas densities as marked. (B) the evolution of non-Boussinesq (solid) and Boussinesq (dashed) velocity and density ratio as a function of radial spread, for a current with initial gas mass fraction of 0.1, eruption mass of 109 kg and mean particle settling velocity of 0.1 m s31 . The density di¡erence between the ambient and current drives the £ow and hence determines frontal velocity.
3.3. Propagation speed Fig. 4A compares the non-Boussinesq prediction for distance propagated as a function of time with that predicted by the Boussinesq model
VOLGEO 2719 6-4-04
276
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
(Dade and Huppert, 1996). The distance propagated is up to 35% less than the Boussinesq prediction for the case of a collapsing fountain (n0 v 0.01), and up to 65% in the case of dome collapse (n0 = 0.001, 0.005). The di¡erences are largest for currents with smaller initial volumes of air in the £ow, since the current density is larger for smaller gas content and thus the value of the densimetric Froude number which governs the £ow dynamics is relatively smaller than the value of the Boussinesq Froude number. Fig. 4 also shows that the di¡erences in the model predictions are largest at smaller times and shorter distances from the source, before substantial sedimentation has occurred. Here the current is most dense, leading to the largest di¡erences between Boussinesq and non-Boussinesq Froude numbers and hence frontal velocity. In the case of an initial gas mass fraction of 0.1, for non-Boussinesq propagation speeds of 80 m s31 , the Boussinesq model overpredicts the £ow speed by about 30^40% (Fig. 4B) leading to comparable di¡erences in the predicted run-out distance. 3.4. Lift-o¡ conditions Lift-o¡ occurs when the density of the current evolves to the density of the surrounding air. This is governed by: (1) the rate of sedimentation of dense particles; and (2) the amount of gas initially mixed into the £ow. In the present model, neutral buoyancy occurs when b* = 0, and the current then lifts o¡. This requires that n = 1/V = bg /ba . Fig. 5A shows that there is a critical value of the initial gas mass fraction, ncrit W0.1, such that lift-o¡ length is maximised. For any given lift-o¡ radius smaller than this maximum there are two possible initial gas mass fractions, related to different lift-o¡ times. In the case of a gas mass fraction smaller than the critical value, less air is entrained into the eruption column. This leads to a larger concentration of particles in the £ow, which causes faster sedimentation, but also a hotter current due to reduced cooling by the air. For a gas mass fraction larger than ncrit , the reduced concentration of particles causes the £ow density to decrease towards the neutral buoyancy con-
Fig. 5. (A) Lift-o¡ radius of the current as a function of initial gas mass fraction. There is a critical value of the initial gas mass fraction, at ncrit = 0.1, for which runout distance is maximised. For smaller n0 there is a higher concentration of particles which causes faster sedimentation, whereas for larger n0 the current has a smaller initial bulk density and hence lifts o¡ after less sedimentation. Eruption mass is ¢xed at 109 kg. (B) Non-Boussinesq (solid) and Boussinesq (dashed) lift-o¡ radius as a function of eruption mass. In this ¢gure, the volatile mixture for each curve (labelled with n0 ) has composition na = 0.12 nv = 0.03, na = 0 nv = 0.03 and na = 0.001 nv = 0. Mean particle settling velocity is ¢xed at 0.1 m s31 .
dition following a smaller amount of sedimentation. In comparison with the predictions of the Boussinesq model, the non-Boussinesq model predicts a shorter runout distance because the associated slower £ow has more time to sediment. The di¡erence between non-Boussinesq and Boussinesq lifto¡ radii is largest for small gas mass fractions and
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
277
hence denser £ows, such as those associated with dome collapse (see Fig. 5A,B). It is worth noting from Fig. 5B that any attempt to estimate eruption mass based on runout distance and mean particle size is subject to uncertainty by over an order of magnitude owing to uncertainty in the initial gas content. 3.5. Deposit geometry The reduction in £ow speed as predicted by the non-Boussinesq model leads to more sedimentation at each point in the current. As a result, the currents form deposits with aspect ratio that is slightly higher than predicted by the Boussinesq model, while the total volume of sediment deposited remains the same. Fig. 6 illustrates this e¡ect. Curves are shown for Boussinesq (dashed) and non-Boussinesq (solid) deposit shapes. The di¡erence using the Boussinesq model is largest in the case of high mean settling velocity. This e¡ect would be larger for channelized £ows, in which sediment is deposited over a much smaller area. The ¢gure also illustrates the very strong impact of mean particle size on deposit geometry which, in fact, dwarfs any di¡erences between the Boussinesq and non-Boussinesq models.
Fig. 7. Evolution and magnitude of dynamic pressure with distance from the source, using the Boussinesq (dashed) and non-Boussinesq (solid) models. For initial bulk gas mass fractions of 0.15 and 0.03, the Boussinesq model predicts pressure up to three times as high as the non-Boussinesq model.
3.6. Dynamic pressure For hazard analysis of pyroclastic density currents, the impact of the current on obstacles in the path of the £ow may be quanti¢ed in terms of the dynamic pressure, pD , of the £ow (Valentine, 1998). This depends on the density bc and velocity uF of the £ow, according to: 1 pD ¼ b c u2F : 2
ð26Þ
In Fig. 7 we compare the prediction of dynamic pressure as calculated using the Boussinesq and non-Boussinesq models. The ¢gure illustrates that for a range of gas mass fractions, the dynamic pressure is overpredicted by the Boussinesq model by as much as a factor of 3. The di¡erences in predictions are greatest within a proximal distance from the source.
4. Discussion
Fig. 6. Non-Boussinesq (solid) and Boussinesq (dashed) sediment deposit pro¢les for currents composed of particles with di¡erent mean settling velocities. Eruption mass is 109 kg and initial gas mass fraction 0.1. Settling velocities are in m s31 .
We have developed a model of a sedimenting dilute pyroclastic density current, accounting for e¡ects of the large density di¡erence between the £ow and the overlying air on the propagation speed of the £ow, applicable to currents with densities up to 100 kg m33 . The model predictions
VOLGEO 2719 6-4-04
278
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
imply slower velocities than for the classic Boussinesq model. This leads to the prediction of greater sedimentation at each point in the £ow and hence shorter run-out distances. Also, the dynamic pressure of the £ow is smaller, impacting models of damage associated with pyroclastic £ows (Valentine, 1998). The model calculations indicate that non-Boussinesq e¡ects increase with the density di¡erence between the current and the ambient. We therefore expect the e¡ects to be more signi¢cant in denser pyroclastic density currents such as those produced by explosive dome collapse, for which the initial gas content is small. For example, recent models of the explosive dome collapse on 26th December 1997 at Soufrie're Hills, Montserrat suggest that a volume of about 5U107 m3 of the andesite dome was disrupted, and the associated blast cloud formed a dilute pyroclastic density current. The deposits from this were laden with ¢ne particles (Sparks et al., 2002) thus we may use the model to illustrate the propagation given these initial conditions. The initial density has been estimated to have been of order 100 kg m33 (Woods et al., 2002). The current spread from the dome towards the coastline, approximately 4 km from the dome, with speeds which have been estimated to be of order 80^100 m s31 (Sparks et al., 2002). Based on the deposit geometry, the £ow may be modelled as a radially spreading current, con¢ned by the topography to spread over an angle of approximately 70‡. In Fig. 8A, we illustrate the predictions of the non-Boussinesq model for the velocity of the head of this current as a function of radial distance from the dome, where, for simplicity, we assume that the material in the head of the current was released rapidly from the fragmenting dome compared to the travel time of the £ow. This is a somewhat simpli¢ed end-member calculation, given that the theoretical predictions by Woods et al. (2001) suggest that the dome break-up took a few tens of seconds, but serves to illustrate the importance of non-Boussinesq e¡ects. Curves are shown for mean particle settling speed of 1 m s31 , with a total mass of 1011 kg, initial radius 500 m and initial density 100 kg m33 corresponding to initial
Fig. 8. (A) Non-Boussinesq (solid) and Boussinesq (dashed) velocity as a function of frontal position for the 1997 eruption at Soufrie're Hills, Montserrat. Curves correspond to initial densities of 100 kg m33 . (B) corresponding dynamic pressure as a function of frontal position. Eruption mass is 1011 kg, initial radius 500 m and mean particle settling speed 1 m s31 .
gas mass fraction of 2U1033 , the gas here taken to be magmatic vapour. The curves show that a mean settling velocity of about 1 m s31 is appropriate, using the non-Boussinesq model. For comparison, the dashed lines illustrate the predictions of the Boussinesq model for the same initial conditions. It is seen that the Boussinesq model predicts velocities for such a dense current of about an order of magnitude greater than the non-Boussinesq model. For such large velocities the £ow is expected to be compressible and not fully describable by the model.
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
Fig. 8B compares the dynamic pressure associated with the head of the current for both Boussinesq and non-Boussinesq models. The nonBoussinesq model predicts pressure to be smaller by several orders of magnitude. In the case of more explosive events of higher gas content, the di¡erences between the Boussinesq and non-Boussinesq models are expected to be smaller (Section 3), since the density of the £ow decreases as gas content increases. For example, during the lateral blast event at Mount St. Helens on 18 May 1980, about 2.5U1011 kg of magma was erupted from the shallow cryptodome with a gas mass fraction of about 0.04. This trav-
279
elled about 15^20 km before lifting o¡. Fig. 9A,B shows that the Boussinesq model predicts a current speed about 60% greater than the non-Boussinesq model, and dynamic pressure to be about 150% greater. During the Taupo event (Dade and Huppert, 1996), around 4U1012 kg was erupted with a gas mass fraction of about 0.01 and spread 80 km radially. Again, the Boussinesq model predicts speeds greater by approximately 60%, as shown in Fig. 9C. Channellized currents, such as those which initially begin as a radially spreading current but then develop into a channellized current owing
Fig. 9. Non-Boussinesq (solid) and Boussinesq (dashed) velocity as a function of frontal position and corresponding dynamic pressure as a function of frontal position for Mount St. Helens 1980 (A,B) and for the Taupo eruption 1800 years ago (C,D). Mean particle settling speed is 1 m s31 in both cases.
VOLGEO 2719 6-4-04
280
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281
to the topography, for example, running into a valley, would be overpredicted by a greater amount than a purely axisymmetric current. In this paper we have examined the e¡ects of adapting a model to include a non-Boussinesq parameter, which has important e¡ects on the £ow. The model gives good insight into the dynamics and deposit of pyroclastic density £ows, however it would be desirable to incorporate other e¡ects such as the slope of the topography and grain-size distribution.
Acknowledgements ‘‘This research supported in part by the Engineering and Physical Sciences Research Council (Nield) and by the BP Institute, Cambridge.’’
Appendix Notation
Frl FrnB g gP h L LΠma mer n n0 ner nv p pD q r Ra Rg Rg Rv t
“t Ta Tc Ter Tg uF V Vg vs Q vb
V ba bc bg bref bs b*
dimensionless time ambient air temperature current temperature eruption temperature erupting gas temperature frontal velocity of the current volume per unit width 2-dimensional gas volume mean particle settling velocity power of time, relating to volume £ux density di¡erence between the current and ambient £uids density ratio between the volcanic gas and ambient ambient density current density interstitial gas density reference density solids density density parameter
References
The following symbols are used in this paper: Ca Cs Cv Fr FrB Frd
Appendix Notation (Continued).
speci¢c heat of air speci¢c heat of the volcanic solids speci¢c heat of vapour Froude number densimetric Froude number for a Boussinesq current densimetric Froude number for a dense (heavy) current densimetric Froude number for a light current densimetric Froude number for a non-Boussinesq current gravitational constant reduced gravity height of the density current current length dimensionless current length mass of entrained air eruption mass gas mass fraction initial gas mass fraction gas mass constant of the erupted gas gas mass fraction of vapour pressure dynamic pressure volume of £ow radius of the axisymmetric current gas constant of the ambient (air) gas constant of the erupting gas bulk gas constant gas constant for vapour time
Benjamin, T.B., 1968. Gravity currents and related phenomena. J. Fluid Mech. 31, 209^249. Bonnecaze, R.T., Hallworth, M.A., Huppert, H.E., Lister, J.R., 1995. Axisymmetric particle-driven gravity currents. J. Fluid Mech. 416, 187^195. Bursik, M.I., Woods, A.W., 1996. The dynamics and thermodynamics of large ash £ows. Bull. Volcanol. 58, 175^193. Dade, W.B., Huppert, H.E., 1996. Emplacement of the Taupo ignimbrite by a dilute turbulent £ow. Nature 381, 509^512. Einstein, H.A., 1968. Deposition of suspended particles in a gravel bed. J. Hydraul. Div. Proc. Am. Soc. Civil Eng. 94, 1197^1205. Gladstone, C., Woods, A.W., 2000. On the application of box models to particle-driven gravity currents. J. Fluid Mech. 294, 93^121. Gro«belbauer, H.P., FannelTp, T.K., Britter, R.E., 1993. The propagation of intrusion fronts of high density ratios. J. Fluid Mech. 250, 669^687. Huppert, H.E., Simpson, J.E., 1980. The slumping of gravity currents. J. Fluid Mech. 9, 785^799. Kie¡er, S.W., 1981. Fluid dynamics of the May 18 blast at Mount St. Helens. In: Lipman, P., Mullineaux, D. (Eds.), The 1980 Eruptions of Mount St. Helens, Washington. Geol. Survey Prof. Pap. 1250, 379^400. Pallister, J.S., Hoblitt, R.P., Crandell, D.R., Mullineaux, D.R., 1992. Mount St. Helens a decade after the 1980 eruptions magmatic models, chemical cycles and a revised hazards assessment. Bull. Volcanol. 54, 126^146. Sparks, R.S.J., Barclay, J., Calder, E.S., Herd, R.A., Komorowski, J.-C., Luckett, R., Norton, G.E., Ritchie, L., Voight, B., Woods, A.W., 2002. Generation of a debris avalanche
VOLGEO 2719 6-4-04
S.E. Nield, A.W. Woods / Journal of Volcanology and Geothermal Research 132 (2004) 269^281 and violent pyroclastic density current on 26 December (Boxing Day) at Soufrie're Hills Volcano, Montserrat. In: Druitt, T.H., Kokelaar, B.P. (Eds.), The Eruption of Soufrie're Hills Volcano, Montserrat, from 1995 to 1999. Memoir of the Geological Society of London. Sparks, R.S.J., Bonnecaze, R.T., Huppert, H.E., Lister, J.R., Hallworth, M.A., Mader, H., Phillips, J., 1993. Sedimentladen gravity currents with reversing buoyancy. Earth Planet. Sci. Lett. 114, 243^257. Valentine, G.A., 1998. Damage to structures by pyroclastic £ows and surges, inferred from nuclear weapons e¡ects. J. Volcanol. Geotherm. Res. 87, 117^140. Wilson, C.J.N., 1985. The Taupo eruption, New Zealand. 2. the Taupo ignimbrite. Philos. Trans. R. Soc. Lond. A 314, 229^310.
281
Wilson, C.J.N., Walker, G.P.L., 1985. The Taupo eruption, New Zealand. 1. general-aspects. Philos. Trans. R. Soc. Lond. A 314, 199^228. Woods, A.W., 2000. Dynamics of hazardous volcanic £ows. Philos. Trans. R. Soc. London A 358, 1705^1724. Woods, A.W., Sparks, R.S.J., Ritchie, L.I., Batey, J., Gladstone, C., Bursik, M.I., 2002. The explosive decompression of a pressurized volcanic dome: the 26 December 1997 collapse and explosion of Soufrie're Hills Volcano, Montserrat. In: Druitt, T.H. and Kokelaar, B.P. (Eds.), The eruption of Soufrie're Hills Volcano, Montserrat, from 1995 to 1999. Geological Society, London, Memoirs, pp. 21, 409^434. Yih, C.-S., 1965. Dynamics of Non-Homogeneous Fluids. Macmillan.
VOLGEO 2719 6-4-04