Abstract

Dome growth at the Soufriere Hills volcano (1996 to 1998) was frequently accompanied by repetitive cycles of earthquakes, ground deformation, degassing, and explosive eruptions. The cycles reflected unsteady conduit flow of volatile-charged magma resulting from gas exsolution, rheological stiffening, and pressurization. The cycles, over hours to days, initiated when degassed stiff magma retarded flow in the upper conduit. Conduit pressure built with gas exsolution, causing shallow seismicity and edifice inflation. Magma and gas were then expelled and the edifice deflated. The repeat time-scale is controlled by magma ascent rates, degassing, and microlite crystallization kinetics. Cyclic behavior allows short-term forecasting of timing, and of eruption style related to explosivity potential.

Growth of the andesite lava dome at Soufriere Hills volcano, Montserrat, has been unsteady and frequently accompanied by cyclic patterns of ground deformation, seismicity, pyroclastic flow generation, and explosive eruptions, overprinted on a background of continuous lava extrusion. We describe the cycles, interpret the causes of cyclic behavior, and develop models for deformation and magma flow. Recognition of cyclic behavior enables short-term forecasting and improved management of volcanic crises. The cycles provide insights into eruption dynamics at andesite volcanoes, showing that degassing, rheological stiffening of the magma, and pressurization in the uppermost parts of volcanic conduits are intimately coupled and control many of the geophysical and dynamical phenomena observed.

The eruption of Soufriere Hills began on 18 July 1995 (1). Earthquake swarms and phreatic explosions preceded the extrusion of an andesite lava dome in mid-November 1995. Dome growth was accompanied by rockfalls and dome collapses producing pyroclastic flows. Some substantial dome collapses were followed by explosive eruptions: a subplinian eruption on 17 September 1996 (2), and sequences of vulcanian explosions on 3 to 12 August and September to October 1997. Activity has been monitored by telemetered short-period and broadband seismic networks (3–5) (Fig. 1), electronic tilt stations as close as 150 m from the crater rim (6), global positioning system (GPS) and electronic distance measurement (EDM) methods (7, 8), correlation spectrometer measurements of SO2 flux (9), visual observations, and surveys of dome and pyroclastic deposit volumes (10).

Seismic cycles have been recognized since 20 July 1996. Before this date, individual hybrid earthquakes had typically been repetitive, of similar size, and regular spacing, but the earthquake swarms showed no clear regularity (5). After this date, the earthquakes ceased to be repetitive, but the events were clustered in swarms. The period in November and early December 1996 was particularly instructive and was characterized by major deformation of the steep southern flank of the volcano known as the Galways Wall. Large fractures and seismically triggered landslides indicated that the Galways Wall was under severe stress. Edifice landslides could be distinguished seismically from dome rockfalls, and deformation of the Galways Wall correlated with hybrid swarms. Dome growth was more active in the aseismic periods.

Tilt cycles were first recognized in January 1997, when low-amplitude (∼2 μrad) cyclic inflation and deflation with 6- to 8-hour periods were observed (6). Tilt cycles represent deformation of the volcanic edifice. More distinct rhythmic patterns (amplitude 15 to 20 μrad, period 12 to 18 hours) were observed in May 1997 (Fig. 2). Deflations are generally more rapid than inflations. In many individual cycles the edifice returned to nearly its initial position, showing that the deformation was largely recoverable. Longer term changes in tilt, however, suggested systematic permanent deformation, or changes in lava loading or magma storage.

Correlation of cyclic CP2 tilt time series (radial component, 067 azimuth) with triggered earthquakes, hybrid earthquakes, LP earthquakes, VT earthquakes, and seismic amplitude RSAM for 18 to 22 May 1997 (see text). RSAM provides a simple real-time quantitative measure of seismicity but does not discriminate between event types; large spikes reflect pyroclastic flows (11). Triggered earthquakes are the numbers recorded at a given station during 10-min periods with data presented as triggers per hour; they are a measure in real-time of hybrid earthquake swarms.

Seismic swarms were coeval with the deformation cycles (Figs. 2 to 4). Seismicity included volcano-tectonic (VT) earthquakes, monochromatic long-period (LP) earthquakes, hybrid earthquakes with impulsive high-frequency starts and long-period codas, rockfall and pyroclastic flow signals, explosion signals, and tremor (5,11). The hybrid earthquakes were most common, characterized by sharp onsets, lack of identifiable S phase, shallow foci (<2 km), and occurrence in swarms (>10 events per hour). Swarms started and stopped abruptly and alternated with periods that were either aseismic or characterized by signals related to rockfalls and pyroclastic flows generated by dome collapses. The hybrid swarms began when inflations were underway and stopped near the peak of inflation or the start of deflation. Dome instability and sometimes a few LP earthquakes generally began at the cycle peak (Fig. 2). Hundreds of cycles occurred during 1996 and 1997. Cycle period (trough to trough) varied from 3 to 30 hours.

Correlation of CP3 radial tilt with triggers, hybrids, LPs, VTs, and seismic amplitude, RSAM for 20 to 27 June 1997. Sharp increases in tilt amplitude, hybrid earthquakes, and RSAM occurred on and after 22 June. LP earthquakes to 22 June are at background level, but VTs representing deep rock fracture on 22 June are notable. A major dome collapse on 25 June generated pyroclastic flows that destroyed several villages and killed 19 persons. Tilt cycles continued after the collapse, with hybrid seismicity reduced. Spikes on tilt record are noise.

Tilt at CP3 and seismicity for 30 July to 6 August. Strong amplitude tilt rises above background level, in correspondence with hybrid earthquake swarms. Sporadic LP activity is above background. Y-axis tilt data are used because radial tilt goes off-scale; the tiltmeter was destroyed on 5 August, after the onset of a series of major vulcanian explosions.

The cycles were visually correlated with major eruptive phenomena. Dome growth, sometimes with spine formation, and rockfall and pyroclastic flow activity became more pronounced after hybrid swarms at the peak of a cycle and during deflation. Dome growth was accommodated by aseismic slip on ductile shear faults with areas of tens to hundreds of square meters with grooved surfaces that bounded one or more sides of an active lobe or large spine. Degassing phenomena were enhanced near cycle peaks and continued into deflations. Degassing was most vigorous when closely spaced hybrid seismicity coalesced into continuous tremor. Nonexplosive continuous ash emission often occurred for tens of minutes during deflations and in some cases correlated with pulsed venting of gas and ash from cracks on the dome surface. Another manifestation of degassing was a monochromatic (1 Hz) signal that accompanied and followed vulcanian explosions.

From late May through mid-June 1997 cycle amplitudes as measured by tilt had declined from 10 to 25 μrad to 5 to 10 μrad, and periods increased to >20 hours; hybrid seismicity also declined (Figs. 2 and3). On 22 June a burst of VT earthquakes occurred and from then on tilt amplitudes sharply increased to ∼30 μrad, periods shortened to 8 hours, and hybrid swarm seismicity resumed (4, 6) (Fig. 3). On 25 June tilt inflation started at about 0900 and peaked at about 1200. The accompanying hybrid swarm began at 1050 and merged into tremor at about 1245, when rapid deflation commenced. Ten minutes later 5 × 106 m3 of the dome collapsed, producing violent pyroclastic flows that moved 5.5 km down the northeast flank and caused 19 deaths. Removal of about the upper 100 m of the 250-m-thick dome on 25 June caused a pronounced deflation, but no major changes in periodicity of the cycles, suggesting that cyclicity was not governed fully by surficial controls such as dome overburden pressure. However, the cycle shape changed as inflations sharpened, and seismicity declined (Fig. 3).

During early August strong tilt cycles and hybrid swarms developed with 10- to 14-hour periodicity (Fig. 4). On 3 August a relatively large pyroclastic flow, anticipated from the tilt cycle and previous activity, swept into the island capital Plymouth. On the next day the first of a series of vulcanian explosions occurred. These eruptions lasted a few minutes at high intensity, then decayed over tens of minutes. They produced fountain-collapse pyroclastic flows and vertical columns rising 7 to 12 km, scattering coarse ash and pumice over Montserrat. The explosions initiated at 10- to 12-hour intervals near the peak of the tilt cycles, accompanied by sharp deflations indicating a nearly instantaneous elastic relaxation of pressure. The tiltmeter systems were destroyed by explosions on 5 August. A second cyclic explosion period began on 22 September and lasted until 21 October, with 75 explosions and related pyroclastic flows. The period between explosions ranged from 7 to 12 hours, with the mean of about 9.5 hours.

Two processes controlled the cyclic activity: (i) degassing and crystallization in the upper conduit, which created a viscoplastic magma plug that inhibited conduit flow, and (ii) pressurization of magma and gas under the plug, which eventually caused extrusion of magma or resulted in explosive activity. Gas exsolution from ascending magma and consequent microlite crystallization (12) increased magma viscosity. We estimate that the undegassed magma had a viscosity of about 106 Pa s, and the completely degassed dome lava had a viscosity of about 1014 Pa s (13) with the material non-Newtonian and with yield strength. Stable heights of lava spines ∼100 m high suggest that the shear strength of extruded hot dome lava was ≥1 MPa (14). Degassing thus formed a rheologically stiffened plug that could resist extrusion and prevent steady conduit flow. We emphasize that cooling had little role in this stiffening process. Conductive cooling penetrated too slowly to be effective on the short time-scales considered here; latent heat was liberated in degassing-promoted crystallization, dome lava was incandescent under its thin outer carapace, and temperatures >700°C have been recorded in pyroclastic flows from dome lava.

Overpressure from the magma chamber was transferred to magma under the plug, and additional high-level pressurization of magma was produced by gas exsolution and microlite crystallization (15). The buildup of pressure in the upper conduit caused edifice inflation, and the amount of inflation (tilt) can be used as a pressure gauge. In May, two tiltmeters were operating at different distances (630 and 770 m) from the dome center (and assumed conduit epicenter) (Fig. 1), and we use the ratio of tilt inflations at the two stations (1.28 average) to model depth of pressurization. For a strain nucleus source embedded in a homogeneous elastic half-space (16), depth is 770 m. However, the pressure source was almost certainly not spherical, and a better model is a finite line source in an elastic half-space (16,17), representing a pressurized vertical cylindrical conduit. The resulting displacements U at the surfacez = 0 are given by(1)where c1 and c2are the depths from the origin of coordinates at the surface to the top (z = c1,r = 0) and bottom of the line source at (z = c2, r = 0),Ri = (r2 +ci2)1/2, are the distances from any point at the surface to the top and bottom of the line source, r̂ and ẑ are unit vectors in the r and z coordinate directions, andf is a constant.

This case produces a range of possible solutions, with the minimum depth to the top of the pressurizing part of the conduit (interpreted as the base of the plug) of 400 m. The value is not sensitive to shifts <200 m of the assumed vent position. The reference instrument altitude (that is, the assumed half-space surface) is 890 m above sea level, so 400 m implies a depth of about 270 m below the base of the dome at 760 m above sea level. The base of the line source (c2) is limited by the depth to the magma chamber, ∼5 km on the basis of seismic (3) and phase equilibria constraints (18). However, if the top depth is >500 m, then base depth is <1400 m, so there is interplay between conduit top and base positions with acceptable solutions. The constant f may be approximated as,where a is the radius of the conduit and Δp is the pressure change within the conduit, μ is the rigidity modulus, and ν is Poisson's ratio. For the conduit top = 400 m, f = 11.6 m2. For the average tilt of 17.5 μrad at CP3 in May 1997, an edifice rigidity of 0.9 GPa (on the basis of compression tests on altered Soufriere Hills tuff), ν = 0.25, anda = 15 m (discussed below), the pressure change is 60 MPa. This is an upper bound resulting from geometric departures from the half-space model and other considerations (19) and is a factor ∼2 greater than the upper bound of 27 MPa derived from ballistic launch velocities (180 m/s) in the 17 September 1996 explosion (2).

The data imply that the hybrid earthquake swarms that initiated during inflations began once a given pressurization threshold was exceeded, and declined when pressure dropped below this threshold. When peak pressures recalibrated for the vulcanian explosions are used (19), the average tilt at hybrid swarm onset in May 1997 (∼5 μrad at CP3) used in Eq. 1 suggests a seismic threshold overpressure of roughly 3 to 8 MPa; these values are similar to the tensile strength (∼4 MPa) of intact Soufriere Hills andesite. Thus, a plausible explanation for the hybrid earthquakes is hydrofracturing and escape of pressurized gas from the upper conduit (<2 km). This explanation is not inconsistent with seismic waveform analysis (4), and we suggest that multiple hydrofracturing events occurred as the conduit gas pressure built toward the critical level required for extrusion or explosion.

Eventually the overpressure increased sufficiently to drive the degassed plug from the conduit, overcoming dome overburden, plug weight, and shear resistance; boundary shear resistance was roughly 0.4 to 1 Pa (20). Extrusion and escape of pressurized gas resulted in a relaxation of pressure in the upper conduit, causing edifice deflation. A fresh batch of magma then ascended to replace the extruded magma, and a new cycle of degassing and microlite nucleation occurred. A new plug formed to retard conduit flow, pressure rebuilt, and the inflation cycle repeated. Petrologic observations support the above interpretations (12). Finally, the pressure buildup at high conduit level in some cases resulted in explosive destruction of the plug, causing vulcanian explosions. The potential for explosivity was indicated by the twin pressure gauges of tilt amplitude and hybrid swarm intensity. Thus, these indicators can often be used to anticipate the possibility of vulcanian explosions, or the potential for autoexplosivity of collapsing pressurized lava domes to generate unusually energetic pyroclastic flows.

The cyclic processes are constrained by the measured extrusion rateQ. Since June 1997, taking the range of Q as 5 to 8 m3/s dense-rock equivalent (DRE) (10), a nominal 11.5-hour cycle yielded about 210,000 to 330,000 m3(DRE) of extruded magma. Equivalent conduit radius was about 15 m, a size compatible with observations of early vents and spine dimensions; thus, the conduit was evacuated in each cycle to a depth >290 to 470 m, with actual depth influenced by vesiculation porosity.

The conduit pressure gradient is constrained by conduit dimensions, extrusion rates, Q, and viscosity:(2)where U is a function of conduit shape (U = 8/πa4 for a cylindrical tube of radiusa), P(z) is the excess pressure above magmastatic at depth z, Vois outlet viscosity at z = 0, and band n are constants (15). The minimum conduit outlet pressure is constrained by the pressure at the base of the 250-m-thick dome, about 5 MPa. The “effective average viscosity” of the magma column is estimated from Eq. 2 as 5 × 106to 1 × 107 Pa s with the method described in (15). Using the ranges of Q above and radius 15 m yields the average pressure gradient 1300 to 4000 Pa/m from top to bottom of the conduit. With depth of magma chamber about 5 km (3, 18), and taking the (minimum) outlet pressure at the base of the dome as 5 MPa, magma chamber overpressure is ∼11 to 25 MPa (21). Comparison of these results with those presented above suggests that shallow conduit overpressures may equal or briefly exceed magma chamber overpressures, thus supporting the proposition (15) that the pressure gradient in the conduit of an andesitic volcano is highly nonlinear.

The long-term average velocity of ascent is Q/area and thus 26 to 41 m/hour; the maximum (centerline) velocity is twice this value, 52 to 82 m/hour, for cylindrical viscous flow. These estimates are consistent with the range of ascent speeds 4 to >42 m/hour estimated from application of experimental data to hornblende reaction rim breakdown; the latter value, representative for the period of cyclic behavior after July 1996, is a minimum rate associated with no observable reaction rim (22). As noted previously, short-term values during deflations may be several times greater than the average value. The correlation between increased ascent rate and increased rate of dome growth over the course of the eruption indicates that the increased extrusion rate resulted from accelerated flow velocity, and not simply increased conduit size (22). Likewise, the correlation between accelerated ascent rate and cyclic behavior after 20 July 1996 is consistent with a reduced rate of volatile loss during ascent at rates >1 km/day; thus, volatile-rich magma reached the upper conduit in < 5 days, whereupon shallow degassing altered critical magma properties in the upper conduit to the extent necessary to trigger unstable flow behavior.

Steady flow can become unstable when flow resistance (or pressure drop) decreases while velocity increases (23). The condition of instability can be simply expressed as dF/dw >k, where F is the resisting force of the upper conduit plug, k reflects the elastic properties of the system (compressibility of conduit magma and elastic conduit wall deformation) being loaded at constant rate (long-term average magma flux at magma chamber outlet), and w is slip displacement. A similar expression can be developed for fluid systems with spatially and time-varying viscosity. The significance of this general approach is its emphasis on stability in explaining cyclic extrusion, and the simplest solution yields expressions for slip durationt′ = π(m/k)1/2, where mis mass, and for duration of reloading t" with the full cycle period thus t′ + t"(24). Although obviously too simple (constant load rate, instantaneous drop in strength coefficients) (25), these expressions provide a first-order stick-slip description of the cyclicity at Montserrat and are qualitatively consistent with strain-weakening and severe viscosity reduction observed near-failure of Soufriere Hills lava in laboratory creep experiments at 900°C (13).

The observations and analyses presented above establish clear links between the geophysical parameters that are monitored, the behavior of the volcano, and the dynamics of extrusion of volatile-bearing magma. Edifice deformation is shown to be related to conduit pressurization with a nonlinear pressure gradient, as influenced by magma ascent rates, degassing, and microlite crystallization. Pressure buildup is indicated by tilt amplitude and by intensity of hybrid seismicity, and both parameters can be monitored in real-time to provide measures of explosivity potential. Overall this understanding advances our ability to interpret monitoring data in terms of the causative physical processes and assists us to forecast the timing and, to a useful extent, the eruptive style. These results improve our ability to mitigate the dangerous effects of andesite volcanism.

↵* To whom correspondence should be addressed at the Department of Geosciences, Penn State University, University Park, PA 16802, USA. E-mail: voight{at}ems.psu.edu

↵† Also affiliated with U.S. Geological Survey Volcano Hazards Program, Cascades Volcano Observatory, Vancouver, WA, and on Montserrat, with British Geological Survey.

Types of seismicity are discussed by A. D. Miller (5) and B. A. Chouet [Nature 380, 309 (1996)]. Here we correlate cyclic tilt time series with triggered earthquakes, hybrid earthquakes, LP earthquakes, VT earthquakes, and seismic amplitude. The latter is recorded by real-time seismic amplitude measurement (RSAM), described by E. T. Endo and T. L. Murray [Bull. Volcanol. 53, 533 (1991)], to yield a simple quantitative measure of overall seismicity. RSAM measures the average absolute amplitude of seismic signals by using a sampling rate of about 60 samples per second. It does not discriminate between various event types, apart from spikes that reflect large pyroclastic flows. Triggered earthquakes are the number of earthquakes automatically recorded at a given station during a 10-min period (data are presented in triggers per hour) (6) and are mainly a measure of hybrid earthquake swarms.

Magma properties vary from those of a hydrated crystal-rich magma with a viscosity of about 106 Pa s to those of a degassed crystalline lava (5 to 10 weight % melt) with viscosity of about 1014 Pa s and yield strength. We have estimated the viscosity of undegassed magma for a crystal content of 65%, temperature of 850°C, and rhyolitic melt phase containing 5% dissolved water using experimental results from K.-U. Hess and D. Dingwell [Am. Mineral. 81, 1297 (1996)]. We have estimated the viscosity of fully degassed and crystalline dome lava from experimental creep tests at 900°C carried out at University College, London, by A. M. Lejeune. These same experiments indicate severe strain weakening with ductile shear zones forming as sample failure is approached. The viscosity just before failure reduces to about 1011 Pa s.

Tilt is likely exaggerated relative to theoretical model values because of topographic and homogeneity departures from the half-space model; as a result the calculated pressure is too large for a given set of conditions by a factor ≥2. Uncertainty of conduit location of 100 m yields a pressure variation <20%. Pressure is sensitive to assumed radius, thus assuming that a ≥20 m can yield “acceptable” pressures ≤30 MPa for the given amplitudes; however, this size is incompatible with other evidence, for example mineralogical constraints on magma ascent rate, and conduit flow modeling. The medium is inhomogeneous, with μlava < μedifice. Full resolution of the issue requires three-dimensional elastic-plastic modeling. Seismic threshold overpressures in the text are calculated by recalibrating the maximum amplitude pressure to 10 to 27 MPa, as deduced by R. E. A. Robertson et al. (2), and applying similar methods for the August explosions.

Shear resistance was calculated by assuming an extrusion overpressure of 10 to 27 MPa and setting it equal to shear strength times boundary area of plug. Plug depth was assumed as ∼270 m. Weight of extruded magma is accounted for by definition of overpressure.

Feedback between pressurization and flow from the chamber result in decline of flux as magma degasses and stiffens, and gas pressure builds, followed by accelerated flux after pressure is relieved. For one set of constraints, the system may be stable and independent of velocity perturbation; for another, it may be stable under quasi-static loading but unstable if subjected to a large pressurization or velocity jump [

We thank our colleagues at Montserrat Volcano Observatory for assistance. Supported by the Department for International Development (UK), the British Geological Survey (BGS), the Seismic Research Unit of the University of the West Indies, and the U.S. Geological Survey (USGS). B.V. acknowledges support from NSF, BGS, and USGS, and R.S.J.S. from the Natural Environmental Research Council, and the Leverhulme Trust.