The mechanics of harmonic gas pressurization and failure of lava domes

Summary

Solutions are developed for the harmonic interior pressurization of idealized hemispherical lava domes to determine the evolution of instability with time. Integrated gas pressure distributions are evaluated to determine uplift on inclined planes of potential detachment that traverse the dome, and to evaluate the limit‐equilibrium stability of blocks isolated by these planes. Average pressure distributions show an initial transient rise to peak dynamic gas pressures that corresponds to a harmonic steady state. Progress towards potential instability mirrors this rise, and oscillates with the harmonic steady state. The timing of pressurization of the dome lavas is modulated by a diffusive time. The magnitude of this pressurization is modulated by dome geometry and the magnitude and non‐dimensional frequency of the interior inflation. Two distinct pressurization behaviours are apparent: the first, at high relative inflation frequencies, is where the dynamic steady state corresponds to static pressurization at the mean interior pressure; and the second, at low relative inflation frequencies, is where the dynamic state synchronously responds to the oscillating interior‐pressure boundary condition, and peak uplift is proportional to the mean pressure plus the maximum harmonic overpressure. These two conditions are separated by a frequency shift, or change in fluid diffusivity, of only two orders of magnitude, and are significant because the destabilizing uplift force may increase by a factor of two, where all other parameters, except excitation frequency or diffusivity, remain fixed. These results are applied to stability evaluations of the lava dome at the Soufriere Hills volcano, to examine failure mechanisms potentially consistent with the observed delay in collapse. Linear models incorporating time‐invariant strength and transport parameters are only able to represent collapse for interior gas pressures at the upper limit of observed magnitudes. Conversely, behaviours incorporating fatigue‐weakening, or hysteretic permeability enhancement, of the dome rocks are able to approach failure for modest pressurization magnitudes, to match observed delays in collapse and potentially to explain why some domes collapse explosively whereas others do not.

1 Introduction

Explosive eruptions of the Soufriere Hills volcano, Montserrat, BWI, directly followed major collapses of the lava dome on September 17 1996, and June 25, August 3 and September 21 1997 (Robertson 1998; Cole 1998; Young 1998), indicating the presence of volatile‐rich magma high in the conduit and suggesting a possible role of gas pressurization in dome instability (Voight & Elsworth 2000). Tilt deformation of the edifice prior to these four collapse events confirmed shallow pressurization and indicated a coincidence of collapse with the timing of peak cyclic pressurization (Voight 1998a). These instances, and a number of analogous cases elsewhere (Newhall & Melson 1987; Sato 1992; Miller 1994; Voight 1998b; Elsworth & Voight 1995; Voight & Elsworth 1997), imply two types of gravitational lava dome failure—one with little or no gas overpressure, and the other with significant gas overpressure (Massol & Jaupart 1999) that influences the failure process (Newhall & Voight 1997). This mode of pressure‐induced failure is consistent with volatile‐rich lava recorded in the crystallinity (Cashman 1992) and vesicularity of fresh lava (Anderson & Fink 1990), and hydrous phenocrysts lacking reaction rims (Rutherford & Hill 1993; Devine 1998). Dome collapse is significant as it may trigger seismic events (Chouet 1996) and pyroclastic flows (Fink & Kieffer 1993; Rose 1976) and result in the explosive unloading of the gas‐pressurized conduit (Robertson 1998; Fisher & Heiken 1982; Alidibirov & Dingwell 1996).

A rigorous analysis of the forces acting within a harmonically gas‐pressurized lava dome is completed to define the relative roles of cohesive and frictional strengths, fluid transport parameters, dome geometry and pressurization magnitudes and frequencies on the potential mode of failure. This paper develops diffusion models to calculate the distribution of gas overpressures in an idealized lava dome, embeds these overpressures into stability analyses and demonstrates that gas pressurization, at observed frequencies and magnitudes, can promote deep‐seated instability and that this approach can explain the timing of collapse. The analysis is applied to dome collapses on Montserrat, and is used to explain mechanistically the potential transition from quiescent to explosive collapse.

2 Pressurization Mechanics

The idealized case of a hemispherical lava dome of radius b, resting above a pressurized conduit, is first considered, as illustrated in Fig. 1. To aid the mathematical development, the connection between dome and conduit is represented as a hemispherical cavity of radius a. In modelling gas overpressure as a function of position and time, general transient solutions are developed to the harmonic pressurization of the dome core. The radial diffusion of gas overpressure, p, within the dome is controlled by the diffusion equation

where κ is the gas or liquid diffusivity of a slightly compressible fluid, r is the radial ordinate and t is time (Carslaw & Jaeger 1959).

Density and buoyant gas transport effects are ignored, as are the dependence of permeability on gas saturation and stress level and the anticipated minor role of ‘undrained’ gas pressure changes resulting from the inflation of the gas cap. These effects are anticipated to be of second order to the role of interior gas pressurization at levels approaching lithostatic stress. The models more closely approximate the case of endogenous dome growth than that of exogenous dome growth associated with lobes that break out to the dome surface, with associated release of gas.

The steady solution for pressure diffusive flow is controlled by the dome geometry, and transient and harmonic behaviour is additionally modulated by the fluid diffusion coefficient, κ, and the frequency of the pressurization signal, ω. Dome destabilization and collapse is ultimately controlled by the distributed form of this diffusing pressure signal, as embodied in gas pressures acting on a potentially unstable block, as illustrated in Fig. 1(a).

Solutions are developed in the following for the transient harmonic non‐inertial gas pressurization of a hemispherical dome and applied to define the effect on time‐dependent instability. Of principal interest is the distribution of gas pressures on planes of potential detachment that may develop within the lava dome and that may lead ultimately to its failure.

2.1 Instantaneous inflation

The first step in developing the harmonic solution is to define the transient changes in interior gas pressure that result from the Heaviside application of initial pressure, p_i(r=a), at the dome core at time t=0⁺. This internal pressurization is superposed on the existing gas pressure profile within the dome, and is typically orders of magnitude larger than the non‐diffusive component (Voight 1999). For gas pressure diffusion with radius r, the solution of

for the boundary conditions

yields the following distribution of pressures:

for a pressure of p_a=p_i=constant at the dome core. This solution contains two primary terms: a first term, (b−r)/(b−a), representing the steady pressure profile and a second term,

diminishing with time that represents progress towards the steady state. This solution is the basic component to be used in defining harmonic behaviour.

2.2 Harmonic inflation

Where interior harmonic pressurization of amplitude dp_i and frequency ω is overprinted on the Heaviside pressurization, the interior pressure changes as p_a=p_i+dp_i cos(2πωt). The resulting radially decaying pressure profile will develop a harmonic overprint as illustrated in Fig. 1(a). This radial distribution of pressure may be determined by applying Duhamel's theorem (e.g. Carslaw & Jaeger 1959),

to represent behaviour of the linearized system by superposition. For this harmonic excitation, the forcing function is φ(τ)=dp_i cos(2πωt), incorporating both the amplitude, dp_i, and the forcing frequency, ω, of the interior pressurization. Substituting eq. (4) for F[t, (t−τ)] in eq. (5) and using the forcing φ(τ)=dp_i cos(2πωt) yields the harmonic solution for null mean pressurization, p_i=0. Further adding the Heaviside solution of eq. (4) results in a full solution for the instantaneous application of core pressure p_a=p_i+dp_i cos(2πωt) against a null background pressure,

with the component ℱ_D defined as

in terms of the ancillary variables

and χ=[1−(a/b)]². The accompanying coefficients, π_D, χ, ω_D and t_D, are non‐dimensional. The most significant in a physical sense are those representing dimensionless inflation frequency, ω_D, and diffusive time, t_D. Again, eq. (6) comprises two components, the first representing steady behaviour [(b−r)/(b−a)] and the second,

the combined transient and harmonic effects, embodied in the term ℱ_D[ω_D, t_D, π_D, dp_i/p_i, a/b] of eq. (7). The square‐bracketed terms represent the functional dependence of ℱ_D.

In eq. (7), the first term represents the transient behaviour and the second the harmonic response. As t_D→∞, the first exponential term collapses, leaving the steady harmonic portion. In this dynamic steady state there are two end‐member behaviours of note. The first is at high inflation frequencies where ω_D→∞, π_D→0 and ℱ_D→0; the resulting radial pressure distribution approaches that for inflation at a constant mean pressure of p_a=p_i, and inflation frequency has no influence. The second is at very low inflation frequencies, relative to the gas diffusivity, where ω_D→0 and consequently π_D→1 and the interior pressurization becomes synchronous with the core pressurization as (dp_i/p_i) cos(2πω_Dt_D). In this mode, at very low frequencies the pressure profile will be at all times steady, and will monotonically approach the instantaneous core pressure with decreasing radius; for ω_D→0 the resulting inflation‐synchronous pressure distribution is given everywhere as a multiple of the steady pressure distribution as

2.3 Block uplift force distribution

In determining the influence of gas pressure diffusion on stability, it is the integrated effect of induced overpressures, acting on a critical surface, that is of interest. Dome gas pressures may concurrently produce downslope block disturbing forces and a corresponding decrease in strength on critical failure surfaces. Correspondingly, it is the effect of gas pressures, integrated over a potential detachment surface, that provides fluid uplift and a corresponding reduction in frictional strength. Eq. (6) may be integrated over a plane that cuts the hemispherical dome geometry, as illustrated in Fig. 1(a), to yield an equivalent uplift force, P. Noting that the components of ℱ_D, viz.ω_D, t_D, dp_i/p_i, π_D and a/b are all independent of the space dimension, r, enables the straightforward integration of eq. (6),

where the relations

and

may be extracted from the geometry described in Figs 1(a) and (b), where the symmetry of the solution with ξ is noted. Eq. (11) is evaluated numerically to define the relative effects of harmonic pressurization of the dome core, and ultimately to define this effect on stability.

3 Pressure Response

The pressure response is examined in this section to define limiting behaviours, and in Section 4 to define how these pressures may mediate dome collapse. These results are subsequently applied to explain the dome‐collapse chronology observed at Soufriere Hills volcano in Section 5.

The uplift force acting on an inclined surface provides a useful index of instability, and is a key component in evaluating the evolution of instability. The full transient form of the uplift, embodied in eq. (11) degenerates to a harmonic steady state at late times (t_D>0.5), and has straightforward limiting magnitudes at high and low pressurization frequencies. Consequently, significantly simplified relations may be used to define the destabilizing magnitude of uplift if both the duration of the transient inflation period is evaluated and a narrowed range of limiting excitation frequencies is determined. These ranges are defined as follows.

3.1 Parameter ranges for Soufriere Hills volcano

As identified in eqs (11) and (7), the magnitude of uplift is controlled by dome geometry, the inclination of the potential detachment surface and the crucial parameters of forcing frequency, ω, and the relevant transport parameter of fluid diffusivity, κ.

Broad ranges of these parameters are available for collapses of the dome at Soufriere Hills Volcano on Montserrat, BWI (Voight & Elsworth 2000). Conduit radius is estimated at 15 m (Voight 1999) and, combined with an observed pre‐collapse external radius of about 200 m, yields a ratio b/a of close to 10. Tilt measurements (Voight 1998a) indicate a steady inflation frequency of the order of 2 cycles day⁻¹ that probably results from oscillatory losses of volatiles in the melt to the conduit wall rocks, with a corresponding increase in melt viscosity (Wylie 1999). Magma properties were reviewed by Dingwell (1997). Magnitudes of fluid diffusivities may be estimated from permeabilities, k, fluid viscosities, μ, and elastic parameters of modulus, E, and Poisson ratio, ν, as

(Rice & Cleary 1976). For a Poisson ratio of ν=0.25, a modulus of the dome materials with between 10 and 0 per cent melt of 2.8–77 GPa (Dingwell 1997) and a viscosity of steam at 800 °C of 10⁻⁴ Pa s results in diffusivities of the order of κ≈k/10¹⁹ m² day⁻¹ for k in metres. For unmeasured permeabilities in the anticipated range 10⁻¹⁸≥k≥10⁻¹³ m² (Voight & Elsworth 2000), this results in diffusivities in the range 10⁵≥κ≥10¹ m² day⁻¹. Together with an observed con‐duit pressure cycling of 2 per day (Voight 1999; Wylie  1999; Denlinger & Hoblitt 1999), these parameters give non‐dimensional inflation frequencies of the order 10⁴≥ω_D≥10⁰.

3.2 Time to steady states

The time taken to reach the dynamic steady state is controlled by the fluid diffusivity and geometry through the non‐dimensional diffusive time, t_D, of eq. (10). The onset and termination of the transient effect on the aggregate uplift force may be indexed relative to the time taken to reach 5 and 95 per cent of the steady pressure distribution through the diffusive times t⁵_D and t_D⁹⁵. These respectively index the onset of potential inflation‐induced instability, its duration, and the utility of using simplified expressions for steady behaviour in the evaluation of stability. Results are documented in Fig. 2 for a variety of dome dimensions. Pressurization times are shortest for potential detachment surfaces passing closest to the dome core, where aggregate uplift forces are correspondingly largest. For relatively large domes straddling relatively small‐diameter conduits (b/a≥10), the duration of the transient component is of the order t_D⁹⁵≈0.3. For the Soufriere Hills volcano, for the estimated range of fluid diffusivities, t⁹⁵ is of the order of 1–10 000 day. The breadth of this estimate reflects the relative uncertainty in determining the fluid transport properties of the dome lavas.

3.3 Steady harmonic response

Subject to the constraint that the transient portion of the behaviour is essentially complete, degenerate forms of eq. (11) provide a convenient means of evaluating uplift and defining its influence on dome stability.

3.3.1 Limiting behaviours: ω_D→∞ or ω_D→0

The radial symmetry embodied in the problem description of eq. (11) enables closed‐form solution for limiting frequencies of the steady non‐harmonic case where t_D→∞ and alternatively ω_D→∞ or ω_D→0. For the first of these cases (ω_D→∞), the steady uplift magnitude, P, is defined in terms of the mean inflation pressure, p_i, the dome geometry, b/a, and the inclination of the failure plane, α, only as

Conversely, at very low frequencies, where the pressure distribution reaches a succession of steady states with the interior forcing pressure, the uplift force is obtained by multiplying eq. (15) by the factor

described previously in Section 2.2. The maximum uplift intensifier is of the order of 1+dp_i/p_i, which for an expected inflation cycling between null pressure and a peak amplitude of double the mean represents a doubling of the uplift. The peak uplift force predicted by eq. (15) is shown in Fig. 3 for the limiting cases of ω_D→∞ and ω_D→0, where the non‐dimensional parameter [P/(p_ib²)](b−a)/a is used. As noted previously, the limiting uplift magnitude for low non‐dimensional frequencies is double that for high frequencies. The question remains as to what actual ranges of non‐dimensional frequency, interior to the full range between zero and infinity, represent effective limits to this behaviour.

3.3.2 Interior to the limits ∞≫ω_D≫0

The simplified expressions representing steady harmonic response, defined in eq. (15), are essentially correct for a range of non‐dimensional frequencies less broad than the range 0→−∞ (as is apparent in Fig. 4). For realistic dome dimensions, peak harmonic uplift magnitudes are shown to prevail for 10²≥ω_D≥10⁻¹ for relatively large domes (b/a≥10), and to move slightly upscale, 10³≥ω_D≥10⁰, for relatively small domes (b/a≤2), as illustrated in Fig. 4. For the known geometry and inflation history of Soufriere Hills volcano, non‐dimensional inflation frequencies are estimated to be of the order of 10⁰–10⁴, suggesting that the ω_D→0 solution is applicable only over the lower part of the frequency range. Correspondingly, the anticipated uplift force may approach that corresponding to double the mean pressure, as suggested by the peak multiplier [1+(dp_i/p_i) cos(2πω_Dt_D)] applied to eq. (15).

3.3.3 Harmonic rise times

For very high‐ and very low‐frequency cycling of pressures within the dome core, relative to the reciprocal of fluid diffusivity, the steady harmonic solutions adequately describe the onset of pressurization and the development of average uplift forces. The rise time to steady behaviour is adequately described using the non‐harmonic form of eq. (11) or Fig. 2. For intermediate frequencies, viz. 10³≥ω_D≥10⁻¹, depending on the geometric ratio b/a and inclination angle α, pressure rise has an overprinted harmonic, as illustrated in Fig. 5 for dome dimensions b/a from 2 to 100. These transient results are asymptotic to the same peak steady magnitudes as reported in Fig. 4, but the overprinted harmonic in the initial rise time is significant in representing the cycling towards instability. At non‐dimensional frequencies, ω_D, of the order of 10⁰–10² and centred on 10¹, peak pressures build cyclically in the initial rise‐time period for t_D≤0.3. This opens the potential, at a relatively narrow range of material and excitation parameters, for progress towards failure building non‐monotonically and not occurring on the first pressurization cycle. This behaviour is consistent with that observed for the 1996 and 1997 dome collapses of the Soufriere Hills volcano, Montserrat, where collapse occurred only after a number of cycles (Voight 1999). Evolution towards instability and the potential for ultimate failure may be evaluated by considering the limiting equilibrium of the dome and associated potential collapse structures. This analysis follows.

4 Collapse Mechanics

The full influence of gas pressurization on dome instability may be evaluated if the calculated pressure distribution is applied to an assumed failure geometry. A potentially unstable, isolated block is defined by an arbitrary basal sliding plane passing through the hemispherical dome inclined at an angle α (Fig. 1a). Gas uplift pressures act on this basal surface. The block is assumed to detach along a steep rear surface, on which gas pressures can act in a direction that tends to destabilize the block. The theoretical stability for a block with this geometry is less than that for a block defined by a single basal plane that passes completely through the hemispherical dome. The specific geometry of this rear block surface is arbitrary, but for analytical convenience it is assumed that the unstable block is defined by the basal plane, and by two planes that form perpendicular bisectors with the basal plane. The two bisector planes are separated by the sector angle ξ, and thus the idealized lava dome failure is a ‘sector failure’ (Fig. 1b). The assumed geometry is consistent with the natural development of radial fractures in growing lava domes.

4.1 Factor of safety

For the block defined in Fig. 1(a), the ratio of resisting to disturbing forces defines stability in terms of a ‘factor of safety,’F_s, as

where W is the block weight, c and φ represent cohesive strength and frictional angle, A is the area of the basal failure plane, δ is a horizontal pseudo‐static seismic acceleration expressed as a fraction of the acceleration of gravity g, and F*_m is the downslope force due either to ‘magmastatic’ pressures, F^m_m, or to diffusing gas pressures, F^g_m, acting on the rear of a block of included sector angle ξ. Failure along the planar basal surface is assumed in order to simplify mathematical treatment for this ‘feasibility’ analysis. The relation F_s≤1 implies failure.

4.2 Geometric components

The weight, W, of the block both drives failure and provides frictional resistance on the mobile sector of the circular failure plane. The downslope disturbing force is further resisted by cohesion applied over the area of the failure plane, A. For a failure plane inclined at angle α to the horizontal (Fig. 1), the geometric relationships of failure plane area, A, block volume, V, and weight, W, for a dome of unit weight, γ_r, are

and

where these quantities are evaluated for the entire spherical cap and may be substituted directly in the stability relation of eq. (16).

4.3 Gas and magma pressures

In addition to the uplift pressures, defined previously, that contribute to the reduction of effective stresses and strength on the failure plane, two additional sources of disturbing forces may alternately act. Destabilizing pressures on the block rear may result from the radial outwardly diminishing distribution of gas pressures or the downward‐increasing magmastatic pressures. The magmastatic force per unit width across the block rear is evaluated from the unit weight of magma, γ_m, as

with y²_max=(b²−x²), to yield the total magmastatic force

Conversely, if radially diminishing gas pressure is the only disturbing component present, the resulting gas force is

with r²=x²+y². This solution closely follows that for uplift pressurization of the circular detachment plane, and ℱ_D[ω_D, t_D, π_D, dp_i/p_i, a/b] is as defined in eq. (7).

These components provide the final unknowns in defining the prevailing mechanisms and timing of dome collapse.

4.4 Non‐dimensional parameters

Evolution towards collapse may be indexed directly through the factor of safety, F_s. The factor of safety may be defined by a minimum number of non‐dimensional parameter groups that represent important component responses of the system. Noting the functional dependences of eqs (11), (17), (18), (20) and (21) as they contribute to eq. (16) enables a minimum set of non‐dimensional groups to be determined. These are

representing on the right‐hand side, respectively, components of strength, gas or magma pressurization, geometry and seismic loading. Representation in this form enables broad classes of instability behaviour to be defined.

5 Collapse Behaviour for Soufriere Hills Volcano

Similar to the poor constraint on the fluid transport parameters that are representative of cooling dome lavas, the strength characteristics are only slightly better defined. For the Soufriere Hills volcano on Montserrat, observed tiltmeter deformations and ballistic throwout suggest cavity pressures in the range 5–10 MPa (Voight 1999), and the observed slow ejection of a 100 m spine from the dome top is consistent with c–φ pairs in the range 0.9 MPa–15° to 0.5 MPa–45° (Voight & Elsworth 2000). These data are combined with the observed external radius of the dome of the order of 200–250 m and a conduit radius of about 15 m. To simplify the interpretation of the results, a ratio of b/a=10 is used as broadly representative of the geometry.

For an unpressured homogeneous dome, the stability of the dome is illustrated in Fig. 6(a) for a variety of friction angle–cohesion pairs. Where the dome strength is increasingly derived from frictional resistance, the inclination of the critical failure plane increases. Cooled and fractured lavas on the exterior of the dome with high frictional strengths may fail on steep detachment interfaces that pass through the dome toe, but will be unlikely to cut close to the dome core; as a consequence, unroofing of the cap will be unlikely. When cohesion comprises a progressively greater component of the mobilized strength, perhaps representative of partly molten materials closer to the core, the critical failure surface becomes shallower but is also innately more stable, as is apparent in Fig. 6(a). Thus, interior pressurization, of sufficient magnitude, may be a key component for failure to progress close to the core in the hot and cohesive materials that comprise the interior. The appropriate overpressures that may result in collapse, and the form of any resulting collapse, depend jointly on the magnitude, the periodicity and the transport properties of the dome lavas. These behaviours may be examined both in terms of a harmonic steady state, to evaluate the potential for collapse, and as a transient harmonic response, to evaluate the timing of this collapse.

5.1 Maximum steady harmonic response

The stability of an unpressured dome, of geometry similar to that of the Soufriere Hills lava dome, is metastable for dome material cohesion of 0.5 MPa and a friction angle of 25°. Interior pressurization will decrease stability, as noted for the relevant geometry of Fig. 6(b). Ultimate stability is represented for the maximum and minimum effects of pressurization at a mean pressurization magnitude of 5 MPa and at the harmonic steady state (Fig. 6b). The maximum stability under pressurization is for high‐frequency inflation with ω_D>10² and for a collapse surface inclined at about 35° to the horizontal. The minimum ultimate stability results at low dimensionless frequencies with ω_D<10⁻¹ and exhibits a propensity for increasingly deep failures (on shallowly inclined detachment planes) to develop that are able to tap the high‐pressure magnitudes close to the dome core and potentially unroof the conduit cap. Superimposing a pseudo‐static seismic acceleration of 0.1 g, representative of the seismic activity at Soufriere Hills volcano, further destabilizes the cap.

As the cohesive strength of the dome is increased, the unpressured dome transits from metastable to stable. For the c–φ pair 2.5 MPa–15°, the effect of interior pressurization is insufficient to bring the dome to failure, as illustrated in Fig. 6(c), where F_s>1. Even the overprinted effects of slow harmonic pressurization (ω_D<10⁻¹) and moderate‐intensity pseudo‐static earthquake loading (δ=0.1 g) are insufficient to bring the dome to failure.

Where the effective diameter of the conduit is increased, as illustrated in Fig. 6(d), the stability is further reduced from the case of a smaller conduit, as the effect of the internal pressure is effectively transferred further into the dome rocks. This situation may also represent the case where the interior of the dome is ruptured, or of higher than average permeability, with a corresponding reduction in the steady‐state stability.

It is clear from these scoping studies that a relatively narrow range of strength parameters constrain the evolution of an unruptured dome, with the cohesive parameter c/γ_rb tan φ controlling this behaviour. For relatively large cohesive strengths, the dome is increasingly more stable, and rupture, if it is able to develop, will follow increasingly deep trajectories. Thus, deep‐seated failure of the dome cap, with the attendant potential for explosive unroofing of the conduit (Fig. 1a), requires both sufficiently high cohesion to avoid shallow ravelling failure during dome growth and sufficient late‐stage overpressurization to cause deep fracturing of the core. If the dome rocks are insufficiently cohesive, resulting from viscous or block‐interlock effects, high‐angle collapses of the outer shell will occur episodically with dome growth, dissipating core pressures in a relatively benign manner. Conversely, cohesive materials may build large‐diameter, steep‐sided domes that remain stable under null pressurization, but may fail on deeply incised detachment planes with appropriate magnitude, duration and frequency of the applied interior pressurization. With this collapse potential defined, the remaining issues relate to the required duration and frequency of pressurization and the expected timing of collapse.

5.2 Transient harmonic response

As anticipated from the previous discussion of integrated uplift pressure as a suitable, but insufficient, index for collapse, the progress to failure may be followed by using eqs (16) or (22) directly. Results are shown for the transient stability of the pressured dome for the two dimensionless geometries corresponding to b/a values of 10 and 5 in Fig. 7. A nominal inflation frequency of 2 cycles day⁻¹ (Voight 1999; Wylie 1999) is used for a moderately cohesive dome of mid‐range cohesion and friction angle representative of the metastable cases of steady harmonic instability of Figs 6(a) and (b). As anticipated, for both geometries the approach to instability occurs over the initial rise‐time interval of 0<t_D<0.3. Each of the curves in Figs 7(a) and (b) represents an order‐of‐magnitude change in fluid diffusivity, κ(ω_D∝1/κ), but with the same excitation of the 5 MPa interior pressurization at 2 cycles day⁻¹ and for constant strength parameters. At low fluid diffusivities, corresponding to low permeabilities and high non‐dimensional frequencies, ω_D, the harmonic overprint of the interior pressurization is completely masked in the failure response, i.e. for ω_D>10². Progress to instability is essentially monotonic, and will occur if the mean interior pressurization is sufficiently high to overcome the strength of the dome. Conceivably, even for this monotonic progress to failure, a discernible signal of harmonic dome tilt and deformation may be apparent and measurable, although the effect on instability is small. Conversely, at intermediate fluid diffusivities and permeabilities and non‐dimensional frequencies (10²>ω_D>10⁻¹), the role of the harmonic overprint is crucial, as progress to potential failure is no longer monotonic but for a linear system is still constrained to within the rise time t_D<0.3. Failure is assumed where the factor of safety trace first crosses F_s=1, with the potential for a variable number of cycles to have been completed within this rise time, as is apparent for ω_D=10¹ in Fig. 7(a). As non‐dimensional frequency drops, corresponding to high fluid diffusivities and permeabilities and low non‐dimensional pressurization frequencies, the peak instability is met on either the first or the second cycle of pressurization, synchronous with the interior forcing, and results in failure conditioned by the cycle time alone, independent of the aforementioned rise time, t_D<0.3. In this situation, failure can only occur on subsequent cycles either if there is a fatigue‐related loss of strength, a loss of strength by some other means, i.e. injection of weaker material or a change in transport characteristics of the dome materials, or if source overpressurization in the dome core builds up.

For plausible dimensions and strengths of the Soufriere Hills volcano in Montserrat, mean core pressurization at 5 MPa with a harmonic overpressure of 5 MPa is shown to be sufficient to induce failure (Fig. 7a), even for a relatively shallow inclination of the potential detachment surface at 35° above the horizontal. It has already been shown that the critical rupture surface may be more shallowly inclined, so these transient harmonic results underestimate the full potential for instability. Where the core radius is further extended, perhaps representing the presence of a soft magmatic core interior, then stability is further reduced. For b/a=5, as shown in Fig. 7(b), the degree of instability is compounded by the effect of interior pressures permeating further out towards the dome periphery.

6 The 1996 and 1997 dome failures of the Soufriere Hills volcano

Dome failures at the Soufriere Hills volcano in 1996 and 1997 followed periods of harmonic tilt deformation and RSAM activity, generally inferred as proxies for gas overpressurization of the underlying conduit cap (Voight 1998a, 1999). These cycles of activity were remarkably harmonic, with a mean frequency of 2.6 day⁻¹, with failure in one case occurring on the ninth cycle after 3.5 day. Tilt measurements and the observed expulsion of a lava spine to 100 m above the dome top constrain maximum interior gas overpressures to the range 5–10 MPa.

These well‐constrained observations allow two different hypotheses for collapse to be tested. The first is that failure occurred within the initial rise time of the interior pressurization as a peak pore pressure built at a critical location within the dome carapace; in this case the harmonic overprint of pressurization is merely incidental and contributes only marginally to the ultimate failure. These failures result from pressurization at high non‐dimensional frequencies, ω_D>10², and may be the most quiescent, as failure is reached over a number of repeated cycles, enabling failure to occur in the external skin of the dome. The second is that failure occurred beyond the rise time, after the establishment of a harmonic steady‐state pressurization within the dome materials, and the advance to failure was prompted by weakening of the dome rocks, either by fatigue, or an expanding volume of weak, less‐crystalline magma in the dome core, or loss of interlock between blocks, or rapid hydrothermal weakening on critical detachment surfaces, or by block detachment at the dome toe (Voight & Elsworth 2000), or from the inertial forces of a well‐timed earthquake.

For failure within the rise time of pressurization, the transport parameters of the dome lavas may be narrowly constrained. For example, for collapse on the ninth pressurization cycle, after 3.5 day, for a 200 m external‐radius dome and noting that during the rise time t_D⁹⁵<0.3 (Fig. 2), an upper limit of fluid diffusivity of κ<3.4×10³ m² day⁻¹ is given. For interior pressurization between 5 and 10 MPa, the frequency of pressurization is 2.6 cycles day⁻¹ (Voight 1998a) or ω_D>30. The progress towards failure for dome inflation at 2.6 cycles day⁻¹ and for κ=3.4×10³ m² day⁻¹ is shown in Fig. 8(a) for contrasting overpressures and average dome strengths. For the plausible c–φ pair of 0.5 MPa–25° [c/γ_rb tan φ=0.2] (Fig. 6b), inflation to a mean pressure of 5 MPa with a harmonic overprint of ±5 MPa is insufficient to cause failure on a detachment surface inclined at 35°. Raising the interior pressure to 10±5 or 10±0 MPa is capable of bringing the dome to failure, but the second of these is incompatible with cyclic pressurization inferred from the tilt measurements. In the first case, the peak overpressure is 15 MPa, also near the extreme limit of plausible pressurization magnitudes. Notably for this threshold fluid diffusivity and observed pressurization frequency, the rise‐time behaviour is given by the mean interior pressure alone (i.e. with ω_D≈10²).

Alternatively, reduced average strengths of the dome lavas may be considered in an attempt to bring the structure closer to failure. If cohesive strength is successively reduced synchronously with an increase in frictional resistance to 0.33 MPa–35° [c/γ_rb tan φ=0.1] and 0.1 MPa–45° [c/γ_rb tan φ=0.02], the dome remains stable on the prescribed 35° detachment surface (Fig. 8a). It should be noted, however, that the dome will be unstable on steeper detachments, even in the absence of pressurization, as illustrated in Fig. 6(a), and a block release at the toe is a possible mechanism for successively deeper, dome‐transecting collapses.

A second hypothesis is that the failure condition is not met within the rise time and failure therefore results from progressive fatigue weakening of the dome, progressive ‘interior‐weakening’ of the dome via cumulative injection of weak lava, a second long‐period (monotonic inflation over a period of days) rise in pressurization superimposed on the twice‐daily inflation or the transport parameters of the dome lavas changing with any of these fatigue‐weakening or long‐period pressurization modes. Some of these mechanisms may act simultaneously and cumulatively to push the dome to failure. For all these prior failure modes, the fluid diffusivity of the dome lavas is no longer constrained by the rise time, and stability behaviour for increases in diffusivity by one and two orders of magnitude are shown in Fig. 8(b) for c–φ values of 0.5 MPa–25°. Relatively low‐order inflation pressures of 5±5 MPa are sufficient (for κ≥3.4×10⁵) to bring the dome to failure, and here the magnitude of the overprinted harmonic component of ±5 MPa is crucial. Notably, for κ≤3.4×10⁵, ω_D≥0.3 and induced uplift pressures are close to the maximum steady harmonic response of ω_D<10⁻¹ of Fig. 6(b). Importantly, for failure in these modes the dome is close to failure on a broad range of low‐angle detachment surfaces, evidenced by the flat slope of the F_S versus α trend in Fig. 6(b) for ω_D<10⁻¹, representing post‐rise‐time failure. This, in contrast to generally higher‐angle initial failures that might occur within the initial rise‐time period, may denote a fundamental difference in failure mode. Intrarise‐time failures may be relatively quiescent, as steep‐angled shallow detachment encroaches insignificantly towards the dome core, whereas late‐time dome‐weakened failures may more directly transect and unroof the dome core, resulting in explosive disintegration of the dome and catastrophic defragmentation of lavas in the underlying conduit.

Most notably, this transition from relatively quiescent progressive collapse to energetic, explosive rupture is conditioned by narrow ranges of strength, overpressure magnitudes, inflation frequencies and transport parameters, and may partly explain why only some domes rupture explosively. Alternative dome models, including those represented by a soft cohesive core and a hard frictional outer shell, could also explain the transition to deep transection and unroofing. These models are similar, in the limit, to the more collapse‐susceptible thin‐shelled domes described, e.g. for b/a=5 in Fig. 7(b).

7 Conclusions

Our analysis has considered the potential role of harmonic interior gas pressurization as a critical cause of the delayed failure of lava domes. In describing delayed failure and the range of observed collapse modes—from relatively quiescent gravity collapse to explosive failure—two general hypotheses appear feasible.

The first, that failure occurs as gas pressures build to a critical level in the dome lavas, requires that the fluid diffusivity is sufficiently low that maximum average overpressure is within the rise‐time period. For the 1996 and 1997 Soufriere Hills volcano dome failures this collapse mode seems only marginally feasible for two reasons. The first is that relatively high interior sustained overpressures are required for failure, of the order of 10 MPa, with an overprinted harmonic component of ±5 MPa, which are at the upper limit of plausible magnitudes (Voight 1999). The second is that for a broad range of c–φ strength magnitudes and subject to expected pressurization magnitudes, only failure of the exterior of the dome carapace will occur (Fig. 6a). If the dome lavas are highly cohesive, low‐angle rupture of the dome is possible and failure is insensitive to gas overpressures, as is apparent in Figs 6(a) and (c). This mechanism favours a dome with a hard exterior shell and a soft cohesive core. If the dome lavas are predominantly frictional, the dome will be potentially more highly affected by gas pressurization, but the incursion of a detachment surface to unroof the dome core is progressively less likely. It is conceivable that the core may be unroofed by the retrogressive failure sequence of Fig. 1(c), with a relatively low friction angle, of the order of φ≤25°, required for transection close to the core (Fig. 6a). For φ≈25°, large gas overpressures are required for failure (Fig. 8a), and this makes failure by this mode less likely. However, the case of φ≪25° is not tested here. Observed failures have been explosive and are also retrogressive.

An alternative hypothesis is that failure occurs after the initial transient period, when a harmonic steady state has been reached. Reconciling how collapse may occur after the first or second cycle of pressurization suggests that this is only possible if structural change of the dome materials is feasible; fatigue weakening, hydrothermal weakening, core weakening by repeated injection of weak lava, or increases in dome permeability are all feasible means of progressively approaching this metastable state. In this mode, the expected gas pressure diffusion response, and final resulting collapse, is not constrained by an upper limit of fluid diffusivity, derived from the time to the steady state (e.g. t⁹⁵). For progressively larger fluid diffusivities, steady peak uplift pressures approach the harmonic upper limit, p_i+dp_i cos 2πωt, synchronously with the exterior pressurization. Thus a mean pressure of 5 MPa with a harmonic overpressure component of ±5 MPa approaches the high‐frequency response for a mean pressure of 10 MPa with any magnitude of superimposed harmonic component. Correspondingly, the dome approaches metastability for much lower average pressures for high fluid diffusivities than it does for low diffusivities.

The critical issue in defining the preference to failure by this low‐frequency mode during the harmonic steady state period is the potential that significant changes in the average mechanical and transport properties of the dome may develop within the time frame of a few days. Strengths may be reduced by fatigue weakening as a result of the cyclic pressurization through plausible mechanisms of thermal cracking, minor hydraulic fracturing, loss of block interlock with arching deformations, hydrothermal weakening, as well as repeated injection of weak volatile‐rich magma into the core. Although timescales are probably too short for significant hydrothermal alteration, all the remaining modes remain feasible. In the absence of long‐period growth in core pressurization, the observed growth of a permanent set in recorded tilt measurements (Voight 1999) underscores, in particular, the potential for loss of interlock. The remaining weakening modes, viz. hydraulic and thermal fracturing, block rotations and loss of interlock, are likely to have the greatest reduction in cohesion, with a smaller effect on the frictional coefficient. These changes may be sufficient to bring the dome structure to a condition of incipient failure.

Accompanying the arching of the dome, loss of interlock will also probably result in changes in fluid diffusivities. Fluid diffusivity is directly related to permeability and deformation modulus as

(Rice & Cleary 1976). Expansive deformations will create new fractures and dilate existing ones, potentially increasing the permeability, k, and reducing the effective modulus, E. Neither modulus or permeability will return to initial magnitudes upon cyclic deflation, and the cubic dependence of bulk permeability on change in fracture aperture will probably outstrip the corresponding reduction in modulus. The timescales for this mechanically mediated process are congruent with the inflation process, and a significant increase in fluid diffusivity could develop during the enhanced inflation period of a few days that is associated with cyclic pressurization. Indeed, the rate of diffusivity increase is controlled jointly by the mechanical properties of the dome lavas and the severity of the gas pressurization, independent of any other timescale. A readily attainable increase in fluid diffusivity of two orders of magnitude may double the disturbing fluid pressurization force if the prior pressurization effect was at or near the threshold minimum effect (ω_D=10²). This potential doubling of the disturbing uplift pressure not only decreases the stability of the edifice but promotes failure on progressively shallower detachment surfaces, with the potential to more nearly transect the dome core. A similar response will result if the interior conduit is widened with the injection of a soft core (Fig. 6d), effectively transmitting core pressures further into the dome rocks. However, a mechanism for the transient expansion of the ‘idealized cavity’ at the dome core by a factor of perhaps two is not apparent.

These discussions note two potential modes of failure. The first, involving an injection of soft material into the core materials, reduces strength with its progression and enables core pressures to be transmitted further towards the outer hard shell of the dome. This dome may fail as a result of the overall reduction in strength, defined in terms of total stresses applied to the system. The second requires a progressive and achievable increase in the fluid diffusivity of the dome materials to bring the dome to failure and ultimately to result in its transection. This mechanism requires no change in strength of the dome with time, and may occur even after the dynamic steady state has been reached. Although in reality weakening of the dome materials and their increase in diffusivity may occur concurrently and additively, the plausible change in diffusivity is potentially most damaging, since only small extensional strains are necessary to amplify fluid diffusivities vastly. This is only true for a specific range of diffusivity, inflation frequency and dome radius triplets and perhaps helps to explain why only some domes exhibit delayed failure and why only some exhibit explosive collapse.

Acknowledgments

This work has been supported by NSF under grants EAR‐9628413, EAR‐9614622 and EAR‐9909673. The constructive comments of S. Vinciguerra and K. Wohletz have added to the final paper and are appreciated.

References