## Abstract

Volatiles dissolved in silicic magma at depth exsolve as the magma nears the surface and cause an increase in viscosity of the magma. A model of a volcanic conduit within an elastic medium and a viscosity dependent on the volatile content of the magma produces oscillatory magma flow for a critical range of steady input flow rates. Oscillatory flow is recognized as a fundamental mode of behavior at silicic volcanoes, and understanding it allows improved short-term forecasting of timing and eruption style.

Recent well-documented extrusions of silicic magma at Soufriere Hills volcano, Montserrat, were unsteady and accompanied by oscillating patterns of ground deformation and seismicity (1, 2). The tiltmeters recorded surface deformation (Fig. 1, top panel)—with periods ranging from 3 to 30 hours—that indicated pressure changes in magma at shallow depths (<1 km). Earthquake swarms occurred about the conduit when a critical magma pressure was exceeded (Fig. 1). Analogous oscillatory behavior has been reported at other volcanoes, such as Unzen, Japan, and Pinatubo, Philippines (3).

The oscillating periods of rapid effusion of magma or explosive activity at Montserrat and elsewhere indicate instability in the magma flow (1). Factors that influence magma flow include complex magma-generation and magma-pressurization processes, irregular conduit geometries, elastic and permeability property variations of the conduit wall rock, and complicated volatile-exsolution processes with feedback on crystallization and rheological changes. Our purpose is to identify the critical dynamic elements that could cause oscillatory magma flow, and we propose a model of oscillating flow that results primarily from viscosity changes induced by degassing and microlite crystallization processes in an elastic conduit. These processes have been documented in eruptions involving andesitic or dacitic magmas (1, 4–6). Although our dynamic model is intentionally simple and caution should be exercised in interpreting the qualitative results, behavior is consistent with observations at Montserrat in 1996–97. We exclude in our discussion mechanisms associated with oscillatory flows of basaltic magmas involving the independent rise of bubbles through fluidal magma to produce large venting events (7); the high viscosity of silicic magmas prohibits the high bubble rise velocities required for this mechanism (8).

Several dynamic models have been proposed for oscillating flows of silicic magma. Whitehead and Helfrich (9) demonstrated flow oscillations due to temperature change in a fluid whose viscosity is temperature-dependent. However, the Montserrat data had oscillations too rapid to be caused by cooling of magma, which is a slow process (1). A model by Ida (10) required rapidly changing conduit outlet dimensions by viscous flow of wall rock; however, this assumption is not acceptable given the time scale of the changes and the brittle nature of shallow wall rock suggested by shallow seismicity and elastic surface deformation (1). Woods and Koyaguchi (11) demonstrated for certain conditions a nonuniqueness of eruption rate for small chamber overpressures that could lead to sharp transitions between effusion and explosive eruptions. Although some aspects of this model seem attractive in relation to Montserrat for the period of September and October 1997, the paper did not treat directly the case of oscillating effusion that is the focus of our report. Finally, a stick-slip condition for compressible magma along the conduit wall has recently been proposed as a first-order mechanism for the oscillations observed at Montserrat (1, 12). The argument for stick slip can be presented in relation to plastic solids (for example, highly crystalline magma) or to fluid systems with spatially and time-varying viscosity as developed by this report.

We assumed that the conduit is a vertically oriented cylindrical tube. Most of the pressure drop along a conduit occurs in approximately the top half kilometer because silicic magma ascending toward the surface becomes heavily degassed and therefore has larger viscosity (13–15). With this in mind, we divided the conduit into two distinct parts. The upper portion is where the majority of the viscous flow resistance occurs, and the lower portion acts like a pressurized elastic short-period storage chamber. The model is idealized, but division of the conduit into these two portions with distinct mechanical functions simplifies mathematical development. Magma is fed into the lower conduit from a magma chamber whose volume is substantially larger than the conduit volume. We assumed a fixed flow-rate inlet condition for the magma entering the conduit from the magma chamber, although we recognize that actual magma input can be more complicated.

Magma is fed into the lower conduit from below by a constant flux*Q̄*
_{0} (overbar denotes dimensional variables), and if this exceeds outflow, the walls dilate. An overpressure driving this dilation is thus recorded by elastic stresses in the surrounding rock. Impelled by this pressure, magma then flows into the upper conduit with a time-dependent flux. We assume that the magma in the lower conduit experiences no degassing and has a dissolved volatile contentν̄_{0}, but on entering the upper conduit, the magma exsolves the dissolved volatiles exponentially with time. This assumption reflects the solubility decrease due to pressure decrease as the magma nears the surface (13) and the dynamics of the process of exsolution (14,16).

In the upper conduit, the dissolved volatile content of the magma obeys(1)where *t̄* is time,*Q̄* is flow rate through the upper conduit,*z̄* is the vertical distance measured upward from the bottom of the upper conduit, ν̄ is the dissolved volatile content, *r̄* is the radius of the upper conduit, and τ̄ is the volatile-exsolution decay time. We assumed that the dissolved volatile decay time is given by the time required for volatiles to diffuse within the melt to a bubble nucleation site. We further assumed that in the upper conduit the dissolved volatiles exsolve and either escape from the magma because of system permeability (17) or that their presence has little effect on the bulk viscosity compared with the effect of the remaining dissolved volatiles (18). For simplicity, pressure effects on volatile decay time are also neglected.

Assume the flow has a low Reynolds number and that the dissolved volatile distribution is constant across the radius of the conduit. The momentum balance is obtained by integrating the equation of motion across the conduit radius (19):(2)Here, *L̄* is the length of the upper conduit, μ̄ is the bulk viscosity of the magma, *p̄* is the pressure in the upper conduit, and *P̄* is the pressure in the lower conduit.

The pressure causing elastic dilation of the lower conduit is reflected by stresses in the lower conduit wall, calculated with elasticity theory (20) to give(3)
*Ē* andσ̄ are the Young's modulus and Poisson's ratio of the surrounding rock, respectively, and *H̄*is lower conduit length. This simplification allows body forces and flow resistance in the lower conduit to be ignored. Equations 1, 2, and3 form a closed model: We nondimensionalize (21) to obtain equations governing the dimensionless variables (overbars absent):(4)
(5)
(6)Two dimensionless numbers ɛ = [16μ̄_{0}
*L̄H̄*(1 +σ̄)]/*Ēr̄*
^{2}τ̄ and *Q*
_{0} =*Q̄*
_{0}
*τ̄/πr̄*
^{2}
*L̄*along with the viscosity dependence on volatiles determine the behavior of the system. The parameter ɛ represents the ratio of viscous pressure drop in the upper conduit to the elastic pressure change in the lower conduit as the upper conduit is filled during one time scale of exsolution. The parameter *Q*
_{0}represents the ratio of the time scale of exsolution to the residence time of magma in the upper conduit.

For any desired μ(ν), these equations are solved by integrating over time and depth. For additional simplification, we reduce them to ordinary differential equations. We integrate eq. 4 from*z* = 0 to *z* = 1 to obtain(7)where *V* = ∫_{0}
^{1}ν(*t,z*)*dz* is the total dissolved volatile content in the conduit.

To estimate the term in square brackets in Eq. 7, we assume that dissolved volatile content is given by its steady-state solution, which yields(8)and that the integrated viscous resistance (13) in the upper conduit, μ_{i}, is a function of the integrated dissolved volatile content in the upper conduit, so Eq. 5 becomes(9)
Equations 6, 8, and 9 form our model. We suppose that(10)so that the total integrated viscous resistance has an Arrhenius-type law. The parameter β, which is a third dimensionless number, represents the sensitivity of the change in viscosity to melt degassing (22). A more accurate parameterization for the viscosity that includes non-Arrhenius effects (8) could be used but would require solution of the more complicated Eqs. 4 to 6. Elimination of the pressure from Eq. 6 withEq. 9 yields
(11)Some of the values for the parameters are rather poorly constrained, but in general the values are guided by observations at Montserrat, with order of magnitude accuracy presumed (1,2). For calculation, we use the following values:μ̄_{0} = 1.5 × 10^{5} Pa·s,*Q̄*
_{0} = 10 m^{3}s^{−1}, *L̄ *= 400 m,*H̄* = 5 km, *r̄* = 10 m, *Ē* = 3 × 10^{10} Pa,τ̄ = 1 hour, β = 5.6, σ̄ = 0.2, and ν̄_{0} = 5 weight %. These values produce dimensionless parameter values of ɛ = 0.53 and*Q*
_{0} = 0.29.

The steady-state relation between *P̄* and*Q̄* for this set of parameters is found by combining the steady-state part of Eq. 8 with Eqs. 9 and 10. A single curve relating *P̄ *to*Q̄ *is the result, and this curve depends only on β (Fig. 2). The dimensionless values have been multiplied by the scaling parameters to revert to physical units. For β = 5.6 (solid curve), a region with negative slope exists, which indicates that oscillatory solutions are possible if*Q*
_{0} is in that region. For a reduced value of β = 3.8 (Fig. 2), the region of negative slope has almost disappeared and little oscillatory flow is expected.

A typical eruption cycle (Fig. 3) can be described as follows: Starting from a pressure minimum at about hour 4 (time zero is arbitrary), the flow rate from the upper conduit is less than the constant flow rate of 10 m s^{−1} into the lower conduit, so pressure begins to build up in the elastic lower conduit. At the same time, viscosity is increasing in the upper conduit because the total dissolved volatile content is decreasing. The initially low value of *Q̄* into the base of the upper conduit is inadequate to resupply dissolved volatiles to this part of the system as rapidly as exsolution is removing them. Initially, *Q̄ *is still decreasing and only begins to increase after hour 7 in response to pressure buildup. At a threshold denoted by pressure of about 9 MPa (Fig. 3), which is reached at about hour 8, *Q̄ *begins to more rapidly increase as the newly invigorated resupply of dissolved volatiles leads to a lower viscosity. At hour 11, when lower conduit pressure reaches a maximum, the acceleration increases even though driving pressure is temporarily constant because viscous resistance is decreasing with the increasing volatiles. During the increasingly rapid effusion of magma, the dissolved volatile content in the conduit continues to increase so that viscosity in the conduit decreases. As a consequence, the outflow exceeds the inflow in the lower conduit so that driving pressure drops after hour 12. At about hour 16, the flow cycle repeats as pressure reaches a minimum.

The system is sensitive to all three dimensionless numbers, and only restricted ranges of each number lead to cyclic behavior. A necessary condition for flow instability is that β > 3.4. Given that β is large enough, *Q*
_{0} must have values within a certain range (which depends on the value of β and ɛ) for instability. For example, *Q*
_{0} < 10 m^{3} s^{−1} in Fig. 2, to the left of the pressure maximum, leads to steady, slow effusion of viscous, volatile-depleted magma.

This feature is consistent with observations at Montserrat, when unstable flow was recognized only after July 1996 after a substantial increase in extrusion rates (1). Slower, generally steady extrusion, consistent with flow to the left of the peak on Fig. 2, had occurred earlier between November 1995 and July 1996. Likewise, the longer period patterns over 6 to 7 weeks at Montserrat, involving nearly steady extrusion after several weeks of strongly oscillating, high-pressure and high-amplitude flowage (1), is most simply explained by a periodic rebuilding of high input flow rates (associated with high magma-chamber pressurization) followed by their gradual decline (23).

Increasing *Q*
_{0} in the unstable range leads to higher frequencies of oscillation. As *Q*
_{0}approaches the pressure minimum of Fig. 2 from the left, oscillation amplitude approaches zero. Beyond the minimum, steady effusion of volatile-rich magma occurs.

If ɛ is below a critical value, no oscillations develop, because exsolution does not have enough time to occur. Thus, more viscous flows over deeper conduit systems with smaller conduit radius and faster dissolution processes are more prone to oscillate. Conduit wall rocks of lesser elastic modulus or conduit shapes conducive to compliant behavior such as blades or high–aspect ratio cross sections promote oscillations.

A decreased amount of dissolved volatiles,ν̄_{0}, in the magma chamber and lower conduit has two effects. First, it reduces the contrast of viscosity between the input magma and the degassed magma and thus indirectly reduces β. If the viscosity ratio becomes sufficiently small, then the flow will be stable and steady. Second, it increases the viscosity of the magma in the upper conduit. This increases the pressure drop along the upper conduit and indicates that larger pressures are required in the lower conduit to extrude the upper conduit parcel of degassed magma. Because more time is required to build up pressures to higher levels, the period of the oscillations increases. Although not explicitly modeled here, the rate of volatile loss from the magma to permeable conduit walls (17) can influence the above arguments; the influence of this factor is reduced with high rates of flow.