## Abstract

We use histories of magma efflux and surface deformation from a continuously operating global positioning system (cGPS) to quantitatively constrain magma transfer within the deep crustal plumbing of the Soufrière Hills Volcano (SHV). Displacement records reach a surface aperture of approximately 11 km and are continuous over three successive cycles of eruption followed by a pause spanning 1995–2008, and we focus on data of this time period. The assumed geometry and flow topology is for twin vertically stacked spherical chambers pierced by a vertical conduit that transmits magma from the deep crust to the surface. For a compressible magma column within an elastic crust we use mean deformation rates measured at between 6 and 13 cGPS stations for periods of effusion then repose and the time-history of magma efflux to define optimal chamber depths and basal magma input. The best fit for a constrained constant basal input to the system is obtained for chambers at 5 and 19 km, and a constant magma input rate of approximately 1.2 m^{3} s^{−1}. Eruptive then pause episodes are, respectively, characterized by synchronous deflation then inflation of both shallow and deep chambers. Throughout this period of three repeated episodes of effusion then repose, the total effusive volume (*c.* 0.95 km^{3} dense rock equivalent, DRE) has been sourced half from the lower chamber (*c.* 0.5 km^{3}) and half from below this chamber (*c.* 0.45 km^{3}). A consistent observation, repeated through three episodes, is that the eruption restarts as the shallow chamber regains its original volume following the pause and that eruption rearrests when the shallow chamber has deflated by a small but constant volume change (*c.* 16–22 Mm^{3}). This magmatic metering is consistent with a control on eruption periodicity that involves overpressured breaching of the shallow chamber followed by underpressured sealing. We contrast these observations with other contemporary models that have consistently placed an upper chamber at a depth of approximately 5–6 km, and deeper chambers at 12 km and deeper.

Melting and buoyant transport in the upper mantle and lower crust results in the complex architecture of connected magmatic systems present beneath many arc volcanoes. The deep supply of melt is driven by the subduction rate, and may be considered constant over very long timescales (>Ma). However, the resulting deep storage and chemical evolution of magma, and its transport to and within the shallow crust and to the surface, is episodic at timescales encompassing millennia and centuries to hours. A variety of mechanisms act at these various timescales; in general, longer periodicities imply controlling processes rooted at greater depth and involving reservoirs of increasing volumes. These deeper processes are, however, less well understood, and longer period processes are less resolved in the historical record.

To explore this, we develop a method to constrain magma efflux using wide-aperture geodetic data. These data are combined with the extrusion record and estimations of magma compressibility and host rock compressibility to recover the history of magma redistribution within the crust. Specifically, the aperture of the global positioning system (GPS) array (see Odbert *et al.* 2014*a* for details) is capable of capturing magmatic exchange to a depth comparable to the surface aperture of the geodetic measurement. We apply this method to early activity at Soufrière Hills Volcano (SHV) on Montserrat, West Indies (Mattioli *et al.* 1998; Druitt 2002; Mattioli & Herd 2003; Wadge *et al.* 2014), where eruptive periodicity is manifested as major migrations of the volcanic centre at intervals of several hundred ka (Zellmer *et al.* 2003), eruptions at 10^{3}–10^{4} year intervals (Roobol & Smith 1998) and recent seismic-intrusive events at approximately 30 year intervals (Young *et al.* 1998), with the most recent eruption comprising a >15 year series of 2–3 year eruptive episodes interspersed with approximately 2 year pauses (Wadge *et al.* 2014). We use these data to describe patterns of magma migration within the crust and to define controls of short-term eruptive periodicity on the order of years.

## Fluid–solid system coupling

Magma transfer within the deep subsurface plumbing of active volcanic systems may be followed by applying geodetic reconstructions to surface-measured time series of displacements. These reconstructions are improved as more independent measurements are applied to the system. For example, measuring surface displacements at wide aperture allow point- or distributed-source volume changes to be defined at prescribed depths. However, the additional requirement that these volume changes must also honour observed magmatic fluxes, where available, forces an added constraint that may result in a more truthful representation of the system. This is the approach taken in this work where the response of the fluid-transport and mechanical systems are jointly represented to enable simultaneous coinversion of the efflux and geodetic data. As shown in Figure 12.1, we consider an idealized magmatic system for SHV that is arguably typical of many other andesitic volcanoes worldwide. The system contains separate deep and shallow magmatic chambers that are arrayed vertically and connected by intervening conduits.

### Fluid-transport system

Vertical conduits connect the two reservoirs, and link the system to the deeper crust and mantle, and the surface. The magma charges the upper (I) and lower (II) reservoirs at mass rates (d/d*t*)(*ρ**V*)^{I} and (d/d*t*)(*ρ**V*)^{II}, respectively, via the adjoining conduits, where *ρ* is average magma density and *V*^{i} is the volume of chamber *i*. The conduits have little capacity for magma storage, relative to the reservoirs, and storage within the conduits is neglected. The linking conduits are of small diameter (order 30 m: e.g. Voight *et al.* 1999); at this dimension, the surface geodetic response is insensitive to their effect at the wide aperture (>11 km) considered here. Mass is conserved within the flowing system where the mass rate of change of chamber *i* may be linked to the magma mass influx () and efflux (), where is the volume rate, as: .

The independent effect of magma expansion by decompression () and the volume change of the chamber () apparent in the geodetic signal may be separated as:
(1)
The rate of change in magma density () is related to the bulk modulus of the magma (*K*_{M}) as:
(2)
through the rate of change in pressure (), and may be related to the volume change of a spherical chamber of radius *a* embedded in an elastic medium of shear modulus *G*_{R} as:
(3)
This is recovered by defining the radial displacement of the empty chamber as it is inflated from within by fluid pressure of magnitude *p*. Substituting the magma pressure to chamber volume relationship of equation (3) to be congruent with the density change in the magma of equation (2), and inserting the result into equation (1) yields:
(4)
Further defining outputs from the system in terms of volume rates, *q*, defined at mean chamber density yields:
(5)
with the influence of magma compressibility represented in the multiplier, *C*^{i}, and with volume fluxes defined in terms of dense rock equivalents (DREs). When the geodetically observed volume change of the chamber is equivalent to the volume change in the magma, when *C*^{i}→1, then the magma can be considered incompressible. Conversely, when *C*^{i}>1, the volume change of the magma is greater than the volume change in the chamber alone, owing to finite compressibility of the magma. Therefore, the influence of magma compressibility on the resulting efflux is indexed relative to the ratios of the relative stiffnesses of the rock and magma (*G*_{R}/*K*_{M}) as embodied in the parameter *C*^{i}.

Assuming that there is no storage in the linking conduits, then for the dual-chamber system considered here, , enabling the relationships: (6) to be combined as: (7) where is the observable surface efflux from the system recorded in DRE.

### Mechanical system

Displacement at the surface results from the superposed inflation, or deflation, of both reservoirs. For an elastic system, with an inflating source at I, the radial () and vertical () surface velocities measured at some location (: Fig. 12.1) are proportional to the reservoir inflation rate as or . The terms and are coefficients relating to the disposition of the measuring point () relative to the source. For an elastic system the effect of multiple reservoirs may be accommodated by the superposition of velocities. This superposition is exact where the sources are of infinitesimal size and within a homogenous half-space but is an approximation as the relative size of the chambers is significant relative to their separation and where heterogeneity is important (Pascal *et al.* 2011). For the two-reservoir system considered here, the resulting radial velocities measured at locations *i*=1, 2 are:
(8)
where the coefficients linking the surface velocities to inflation rates are recovered from the Mogi solutions for a strain nucleus as (Mogi 1958; McTigue 1987). This represents the effect of an inflating chamber collapsed to a point within an homogeneous elastic half-space bounded by a horizontal surface, and for a Poisson ratio of *v*=0.25. These assumptions have proved adequate for more restrictive analyses on the same system (Widiwijayanti *et al.* 2005), and are used directly here.

The resulting surface radial velocities may be determined for a second reservoir () by permuting the coordinates (). Alternately, if vertical displacement rates are available, then the equivalent of equation (8) may be constructed with and .

### Co-inversion

Velocities measured at two different radial distances (say, radial velocities and ) together with the surface efflux () enable equations (7) and (8) to be solved simultaneously for the two reservoir inflation rates and , and the supply from the deep crust–mantle flux . From these, the full suite of flux terms defining exchange between the various reservoirs may be determined. Importantly, measured velocities must be available for at least two different radial locations – radial and vertical velocities measured at a single location are not independent – although, vertical and radial velocities from different locations may be mixed. In this discussion, the depths and locations of the reservoirs are assumed known, but are based on inversions of the entire surface displacement field at specific epochs (Mattioli *et al.* 1998, 2010; Mattioli & Herd 2003). If additional reservoirs were to be present, then the number of displacement-measuring stations needed to solve for magma flux history is increased by one for each additional reservoir.

## Deep magma transfer at SHV 1995–2008

Continuous and highly resolved geodetic and efflux records are available for very few volcanoes. The recent eruption at SHV is an exception, providing a valuable window into deep processes contributing to strato-volcano behaviour.

### Eruptive chronology

The current eruption at SHV followed a pattern of seismic crises separated by about 30 years (Young *et al.* 1998). The most recent volcanoseismic crises, in the 1890s, 1930s and 1960s, are interpreted as aborted eruptions, with a seismic crisis in the 1990s developing into the ongoing eruption. In the most recent eruption, phreatic activity began in July 1995 following several years of seismic unrest. Lava was extruded continuously in Eruptive Episode 1 from November 1995 until about 10 March 1998, followed by a period of an eruptive pause with passive dome collapse ending in November 1999 (Druitt & Kokelaar 2002; Norton *et al.* 2002). This cycle, of an active lava dome growth episode followed by a pause, was repeated a second time, with Eruptive Episode 2 lasting between December 1999 and mid-July 2003, followed by a pause lasting until October 2005 (MVO 1995–2008 unpublished data; Wadge *et al.* 2010, 2014). Eruptive Episode 3 commenced in August 2005 and ended in March 2007 (Loughlin *et al.* 2010), followed by a pause lasting until July or August 2008. Thereafter, there were two additional, but shorter, episodes (July or August 2008–January 2009 and October 2009–February 2010) (Wadge *et al.* 2010), which suggest a fundamental change in periodicity and related flow dynamics of the system. In this work, we focus on the first three cycles of eruption followed by repose.

The restriction of volcano-seismicity since magma breakthrough in 1995 to a depth of <5 km, the stability of crystal phases at pressures of approximately 130 MPa and inversion of early GPS data (1995–1998) were consistent with a magma chamber top at a depth of roughly 5 km (Aspinall *et al.* 1998; Barclay *et al.* 1998; Mattioli *et al.* 1998; Voight *et al.* 1999; Miller *et al.* 2010), and this has been corroborated by recent seismic tomography (Shalev *et al.* 2010; Paulatto *et al.* 2012). Small enclaves of basalt mixed in the erupted andesite indicate a deeper supply of hot mafic magma (Murphy *et al.* 2000; Annen *et al.* 2006; Barclay *et al.* 2010), and some crystal phases suggest a connected deep chamber at >10 km (Devine *et al.* 2003). Geodetic data have provided further support for deep storage (Mattioli *et al.* 2010), whether interpreted as single chambers of various shapes (Mattioli *et al.* 2010; Voight *et al.* 2010), or interconnected chambers in shallow (*c.* 5 km) and relatively deep (>12 km) crust (Elsworth *et al.* 2008). In a later study (Foroozan *et al.* 2011), using inversions of four continuously active GPS stations, it was shown that the deep chamber undergoes volume changes an order of magnitude larger than that of the shallow one.

We have shown that the difference between surface displacement signatures of GPS stations located near or far from the volcano conduit cannot be explained with a single pressure source (Foroozan *et al.* 2010). The significant cumulative volume of the eruption (*c.* 1.0 km^{3} from 1995 to 2009) (Wadge *et al.* 2010), and its continuity and chemical consistency, coupled with detailed observations of co-eruptive ground displacements suggests a voluminous upper magma source of several km^{3} (Voight *et al.* 2006). When taken together, these observations constrain our proposed model of two vertically stacked magma chambers. In this work, we report inversions for data from between six and 13 GPS stations (compared to four previously: Elsworth *et al.* 2008) and for new efflux data (Wadge *et al.* 2010). We then use a compressible magma phase (Huppert & Woods 2002) (incompressible in Elsworth *et al.* 2008) to define metered volumes of the upper chamber that control the periodicity of eruption (Foroozan *et al.* 2011).

### Constitutive data

Where magma is compressible in comparison with the host rock containing the chamber, the geodetic signal represents only a portion of volume change of the magma inflow. In this case, the ratio of the moduli of the rock and magma controls the partitioning of the geodetic signal between the volume change in the chamber at constant density and the density change in the magma at constant volume. The geodesy records only the former signal (1 in equation 5) and is blind to the component related to the expansion of the magma (4*G*_{R}/3*K*_{M} in equation 5). Thus, behaviour is indexed by the coefficient *C*^{i}=1+4*G*_{R}/3*K*_{M}, enabling magma transport to be recovered from the geodetic data even if the magma is compressible with *K*_{M}<*G*_{R}.

#### Magma compressibility

The bulk modulus of magma (*K*_{M}) may be related to density as: *K*_{M} = *ρ*_{0} ∂*p*/∂*ρ*. Density, *ρ*, is dependent on the exsolved gas content (Huppert & Woods 2002) as:
(9)
where *n* is the exsolved volatile content following Henry's law, with *s*=4.1×10^{−6} Pa^{−1/2} for the case of water vapour (Burnham 1975), and *x* and *N* are crystal and volatile contents. The gas phase follows the ideal gas law: . A melt water content of *N=*5% for the SHV magma (Barclay *et al.* 1998; Murphy *et al.* 1998) and a volumetric crystal fraction of *x*=40% (Devine *et al.* 1998; Murphy *et al.* 1998) are assumed in these calculations. Based on petrology, the magma temperature (*T*) is assumed to be constant at 850 °C (Devine *et al.* 1998; Murphy *et al.* 1998). Pressure (*p*) is lithostatic assuming constant rock density of *ρ*=2600 km m^{−3}, with densities of crystals (*ρ*_{c}) and melt (*ρ*_{m}) considered constant at *ρ*_{c}=2600 km m^{−3} and *ρ*_{c}=2300 km m^{−3} (Melnik & Sparks 1999; Costa *et al.* 2007).

#### Rock compressibility

For the rock we use a variable shear modulus recovered from seismic tomography profiles (Shalev *et al.* 2010) approximated with a 1D profile (Foroozan *et al.* 2010). The shear modulus of the wallrock, is estimated on the basis of a 1D seismic velocity profile (Shalev *et al.* 2010) and decreased by a factor of 10 to account for possible large-strain effects (Bonaccorso *et al.* 2005) and modulus ‘softening’ related to thermal effects under SHV that are unrepresented in the extrapolated 1D seismic velocity data. The compressibility factor due to this modulus choice linearly affects the magma volumetrics of the shallow chamber but does not change our basic conclusions.

#### Composite compressibility

The degree to which the geodetic data directly reflect magma redistribution is indexed to *C*^{i}. As a result, for the case of a relatively deep magma chamber (depth >12 km) the assumption of incompressibility will be reasonable; however, compressibility must be considered for a shallower magma chamber. A compressibility factor *C*^{i} is used for the Mogi source (Elsworth *et al.* 2008) defined in equation (5). Thus, magnitudes of the compressibility factor are evaluated with depth but, for a shallow crustal modulus of G_{R} of approximately 5 MPa, are of the order of 4.5 at 4 km and <1.2 at >12 km.

### Analysis of continuously operating GPS (cGPS) data

The idealized geometry of the magmatic plumbing system is that of dual, vertically stacked, chambers pierced by a vertical conduit that transits from deep crust and mantle to the surface (Fig. 12.1) (Elsworth *et al.* 2008). This system is constrained by measured efflux rates and surface displacements that are co-inverted to define: (i) optimal locations of the chambers; and (ii) the decadal history of mass transfer within the magma column.

#### Data

Magma efflux rate and episodic removal by dome failure has been recorded continuously since the onset of the magmatic eruption in November 1995 (MVO 1995–2008 unpublished data; Devine *et al.* 1998; Sparks *et al.* 1998; Wadge *et al.* 2010). Similarly, long-term surface displacements have been recorded since 1995 by cGPS stations between about 1.6 and 9.6 km from the conduit, and at a full range of azimuths (Mattioli *et al.* 1998, 2010; Mattioli & Herd 2003; Odbert *et al.* 2014*b*), as shown in Figure 12.2 (left). These are appropriately corrected for plate motion in a standard manner (DeMets *et al.* 2000, 2007; Jansma & Mattioli 2005). Consistent with surface displacement above a point pressure source, vertical velocities decrease radially outwards from the chamber centre, and horizontal velocities increase from zero above the chamber, through a peak and then towards zero in the far-field (Fig. 12.2, right). Estimates of surface efflux from the conduit provide average efflux rates (Fig. 12.3) in each of the three periods of eruption and three periods of pause (Wadge *et al.* 2010). Average velocities in each of these periods are derived from the cGPS (Fig. 12.3) and used to first constrain optimal depths of the chambers, and then deconvolve the inflation and deflation histories of these chambers that thereby define magma transport in the upper crust.

#### Constraining chamber depths

Elsworth *et al.* (2008) used pairs of GPS stations to individually solve for chamber inflation/deflation at depths prescribed *a priori*. Foroozan *et al.* (2011) directly constrained the locations of chambers based on a best fit to the geodetic data. This uses a least-squares fit inversion that search the full parameter space of potential chamber depths for a suite of up to 13 GPS stations available during 1995–2008. As in Mattioli *et al.* (2010), given the uncertainties involved (Fig. 12.2, right), no statistic is robust enough to unambiguously constrain the depths of the chambers. With each depth combination of the Mogi pressure sources, we constrain the mass balance throughout the crustal column with the observed efflux. Thus, the magma influx rate to the deep chamber (together with the chamber volumetrics) is calculated from the volume change rates of the chambers, in turn evaluated from geodetic inversions constrained for magma compressibility. In a volcanic system, the basal influx may be reasonably considered as constant over the timescale considered here (Voight *et al.* 2010). The feasibility of this is examined by plotting the coefficient of variation of deep chamber influxes between six episodes of eruption or pause in the time span 1995–2008 (Fig. 12.4). It is significant to note that for this inversion the full suite of cGPS data are included and for all six epochs (1995–2008), further details of the inversion method are contained in Foroozan (2011).

The (shallow, deep) chamber depth combinations (3–7.5, 13–23 km) have coefficients of variation for the deep influx of less than 0.75 standard deviations between 0.41 and 1.19 m^{3} s^{−1}, and mean influx values of 0.88–1.69 m^{3} s^{−1}. The shallow chamber depth must be consistent with petrological evidence that implies long residence times and pressures of approximately 130 MPa (Barclay *et al.* 1998; Devine *et al.* 2003). Constraining the shallow chamber depth to 5 km, minimal coefficients of variation for the deep influx (Fig. 12.4) arise for a deep chamber at at a depth of about 15–21 km with a minimum at 19 km, resulting in a mean deep influx of 1.21 m^{3} s^{−1} with a standard deviation of 0.47 m^{3} s^{−1} during all the episodes of activation and pause. The standard variation of the deep influx remains under 0.60 m^{3} s^{−1} for the deep chamber at 17–22 km and the shallow chamber at 4–5.5 km. The geometry that resulted from the inversion performed for Mogi sources in a homogeneous medium can be transferred into the heterogeneous case of SHV using a linear transfer function resulting in deeper sources (by *c.* 1 km or more) (Foroozan *et al.* 2010).

Least-squares inversions are then completed for average surface velocities for chambers fixed at depths of 5 and 19 km. These inversions simultaneously satisfy: (1) observed magma efflux rates; (2) conservation of magma mass within the crustal column; (3) compressibility conditions of the magma with depth; and (4) constant deep influx (>30 km) into the deep chamber. Geodetic results are shown in Figure 12.2 (right). Items (3) and (4) give significant extra constraint over previous results (Elsworth *et al.* 2008), specifically with regard to identifying triggering mechanisms.

Based on our analyses, constraining basal influx results in a slightly degraded fit to the data relative to the unconstrained case, but the difference is not significant as the change in the coefficient of variation between the two sets are less than 5%. The benefit of adding this constraint, however, is that the resulting model satisfies the expectation that deep influx remains unperturbed by feedbacks from fluid transfer within the shallow crust in the decade-long timescale considered here (i.e. that the magma supply rate driven by the subducting slab is to a first approximation uniform over long periods).

#### Constraining magma transport

Rates of magma transport and their distribution relative to the two chambers are outcomes from the analysis. This enables cumulative volume changes of the deep and shallow magma chambers to be evaluated (Fig. 12.5). Eruptive and subsequent pause episodes are characterized by synchronous deflation then inflation of both the shallow and deep chambers. The results imply that the depletion rate of the deep chamber exceeds its inflation rate, and has depleted by about 500 Mm^{3} in volume since the initiation of activity in 1995. The total erupted material (DRE) for the first three cycles is 0.95 km^{3}, implying that about half of the erupted material has been supplied from sources below the deep chamber at a rate of approximately 1 m^{3} s^{−1}.

A consistent observation, repeated through three episodes, is that the eruptive episode starts approximately at the time that the shallow chamber has (re)gained its original volume after each preceding pause (Fig. 12.5). Eruption is triggered as either the volume of the shallow chamber or its related internal pressure reaches a threshold magnitude and either dilates or (re)fractures the chamber exit, and thus allows magma to enter a conduit that connects to the ground surface (Voight *et al.* 2010).

A complementary observation is that the eruptive phase transitions to a pause when the shallow chamber has lost 16–22 Mm^{3} of its volume during the eruption. If this chamber volume change (average of 19 Mm^{3}) is equivalent to a stress change limited by the hot tensile strength of intrusives, about 5–20 MPa (Pinel & Jaupart 2003; Tuffen & Dingwell 2004), then the size of the shallow magma chamber can be estimated. The tangential tensile stress around a pressurized sphere is half of its internal pressure, so the pressure (in excess of lithostatic) of the chamber at the onset of the eruption phase is 10–40 MPa. We assume a static shear modulus in the range of 3–10 GPa for the shallow crust surrounding the magma chamber. Moduli based on island-wide 1D tomography data to <5 km depth (Shalev *et al.* 2010) are not reliable for this purpose because these data do not record the low seismic velocities of high-temperature rock enclosing the SHV magma system. Knowing the average volume of lava erupted per cycle (19 Mm^{3}), the pressure loss can be computed as a function of chamber radius. Capping the pressure loss to 40 MPa (maximum initial excess pressure), the chamber size cannot be smaller than about 1.8 km^{3} (radius >0.75 km) for *G*=3 GPa, or 6.0 km^{3} (radius >1.1 km) for *G*=10 GPa. Pressure losses in the range 10–20 MPa yield radii of 0.95–1.2 km and volumes of 3.6–7.2 km^{3} for *G*=3 GPa, and radii of 1.42–1.79 km and volumes of 12–24 km^{3} for *G*=10 GPa. Assuming a radius of 1 km and a volume 4.2 km^{3} (Voight *et al.* 2006), the pressure drop at the end of an eruption episode will be approximately 17 MPa for *G*=3 GPa and 57 MPa for *G*=10 GPa, the latter seemingly excessive for a volume of 10 km^{3} (radius 1.34 km), and the pressure drop is 7.2 MPa for *G*=3 GPa and 24 MPa for *G*=10 GPa (these results are halved for 20 km^{3}). The above calculations refer to the case for 40% crystallinity. For 60% volume crystal content, the average volume loss is 14±3 Mm^{3}, yielding slightly reduced chamber dimensions. Interpreting these results, if we take pressure losses of 10–20 MPa as the most probable, we find that chamber volumes >10 km^{3} are feasible if *G* is approximately 10 GPa, whereas small chambers of <10 km^{3} require a very low magnitude of modulus.

## Discussion

Concurrent observations of surface efflux and surface displacement illuminate rates and patterns of magma transport at crustal depths below SHV. The analysis presented in this chapter uses a two-chamber model coupling fluid mass balance and surface-measured displacement with compressible magmatic constituents for the full six epochs of eruption–repose from 1995 to 2008. We simultaneously invert geodetic data and measured magma efflux rates with the additional constraint of constant basal influx to determine rates of magma redistribution in the crust. Results constrain basal supply to be approximately 1.2 m^{3} s^{−1} for the initial eruptive activity (1995–2008) at SHV. Volume changes within the deep and shallow chambers are synchronous but differ in modality as the larger deep chamber is net depleting with time while the smaller upper chamber cycles between upper and lower limiting volumes – interpreted as limiting overpressures that either initiate or staunch the eruption. If limit overpressures are 10–20 MPa, shallow chamber volumes >10 km^{3} are feasible if the shear modulus of the rock is ∼10 GPa; a small chamber <10 km^{3} would require a very low-modulus host. These observations provide rare documentation of the cyclic control of an eruption through volumetric changes and associated magma overpressure fluctuations.

However, these characterizations of the subsurface plumbing are merely one set of evaluations within a broad and evolving suite of complementary models. These models have considered the subsurface plumbing system at a variety of length and timescales, and are summarized in Figure 12.6. Similarly, these documented analyses are preceded by others and this discussion is not meant to be exhaustive – rather, it provides a summary of contemporary understanding of the system that is supplemented by recent additional evidence of its structure derived from seismic tomography (Paulatto *et al.* 2012).

These various approaches have codified the supposition (Voight *et al.* 1999; Elsworth *et al.* 2008) that, as the time-separation increases between eruptions followed by pauses, these events are sourced progressively deeper and from more voluminous magma reservoirs (Elsworth *et al.* 2008). Thus, short-period extrusions, 5–7 weeks in duration, source from the shallow subsurface (Widiwijayanti *et al.* 2005; Costa *et al.* 2007) and where periodicity may be suitably explained by a rheologically non-linear melt combined with the mechanical compliance and volumetric capacitance of a shallow dyke (Costa *et al.* 2007). This dyke is fed from below about 5 km at constant pressure from a voluminous magma chamber (Fig. 12.6). Indeed, this dyke may be an important feature and is potentially aligned along the surface expression of a fault that transects the upper structure of the volcano aligned NNW–SSE. However, despite its presence and with a potential along-strike width of 500 m (Costa *et al.* 2007), its effect on the deformations around the volcano diminishes quickly with increased radius – the far-field velocities are predominantly radial (Elsworth *et al.* 2008; Foroozan *et al.* 2010). This is a manifestation of the increasing effect of the deep magmatic system on ground surface velocities measured far from the volcano conduit (>5 km) and the relatively small magnitudes of overpressures that may be generated at shallow depth.

The observed GPS velocities are largely radial, and the intermediate- and far-field velocities together capture, respectively, the shallow and deep inflation response in the subsurface. Indeed, this was the principal motivation for the work discussed in this chapter (Elsworth *et al.* 2008; Foroozan *et al.* 2010), as the differential sampling of intermediate- and far-field velocities allows magma transfer to be defined within the shallow crust. The best available petrological and seismic data have been used to constrain the depths of shallow (5–6 km) and deep (12 km) magma chambers. Fitting surface velocities for a full suite of GPS measurements at a range of azimuths from the volcano – in the intermediate (*c.* 6 km) and far field (*c.* 9 km), and throughout the then full sequence of the three sequences of eruption then – allowed the magma exchange in the shallow crust to be determined. When attempted for an incompressible magma (Elsworth *et al.* 2008), this showed that the periodic efflux was entirely consistent with a near-constant magma input at the base of the system (into the lower chamber) and that the periodic forcing appeared to result from below the upper chamber. Later analysis (Foroozan *et al.* 2011) for the same geometry and geodetic time series, but with a compressible magma, suggested that the periodicity resulted from the metering of magma by the upper chamber – eruption resulting when the upper chamber had refilled to a prescribed volume and, hence, re-established the fracturing pressure needed to reopen the system to the surface. Thus, the predictions were mechanistic in that they defined plausible mechanisms that could result in the observed periodicity – compressible magma supplied at a constant rate into a constant-volume upper magma reservoir prescribes the time taken to reach a critical overpressure and, hence, reinitiation of eruption.

In reality, the geometry of the system is ill constrained owing to the nature of the transfer functions linking the Mogi-type strain nucleus to the surface deformation – a large change in the depth or magnitude of the strain nucleus has a proportionately smaller influence on the displacement field. Thus, multiple and ambiguous presumed chamber geometries can equally satisfy a single set of noisy displacement data (Foroozan *et al.* 2010), this is an intrinsic challenge in inverse analysis and a limitation in all fits to geodetic data.

Other approaches have considered the subsurface chambers to be arranged both as a single vertically aligned ellipsoid linking at a depth of between 6 and 14 km (Fig. 12.6) and containing a compressible magma (Voight *et al.* 2010), and as twin vertically stacked chambers containing an upper dyke (Hautmann *et al.* 2010). Each of these approaches have used selected epochs of eruption-followed-by-repose (just one set of eruption followed by repose) as a preferred dataset to reach a best fit with geodetic data. This is in contrast to other analyses (Elsworth *et al.* 2008; Foroozan *et al.* 2011) that have attempted to fit all three epochs of eruption–repose rather than just the central (2003–2005 then 2005–2007) set. The use of a restricted set of data is justified by the presumed improved quality-of-fit to displacement data within these epochs. However, as a result, there is no guarantee that the improved fits match displacement data in the other epochs.

In general, the ‘magma sponge’ model (Voight *et al.* 2010) may be considered as a variant of a two-chamber model linked by an infinitely conductive conduit, requiring that inflation overpressures are instantly equilibrated between the upper and lower chambers (Fig. 12.6). Thus, this differs from the twin-chamber models principally in that both the upper and lower chambers must be synchronously overpressured as opposed to the two-chamber models that may be asynchronously pressured. Fitting the data with a single-chamber model is more challenging as providing more degrees of freedom (dual chamber models) allows a better fit to data but a less faithful representation of Occam's razor.

In overview, the earliest of the dual-chamber models (Elsworth *et al.* 2008; Hautmann *et al.* 2010) and the magma sponge model (Voight *et al.* 2010) all give remarkably similar geometries for the subsurface plumbing (Fig. 12.6), despite rather different assumptions. The Elsworth *et al.* (2008) model is for incompressible magma and the other two are for compressible magma; the ‘magma sponge model’ (Voight *et al.* 2010) is for a single chamber and the other two are for dual chambers; the Hautmann *et al.* (2010) model includes a near-surface dyke and heterogeneous rock modulus with depth and the others do not. This, in part, reflects the desire to honour petrological and seismic observations in addition to the constraints applied by the geodetic data. With upper chambers placed at about 5–6 km, these geometries are all similar to the most recent characterizations recovered from seismic tomography (Paulatto *et al.* 2012) and joint inversion of thermal models.

A further attempt (Foroozan *et al.* 2011) to define geometry of the chambers has been undertaken where formal minimization is attempted between measured displacements for the second (1999–2003–2005) and third (2005–2007–2008) epochs of eruption-followed-by-pause. This attempt included a heterogeneous distribution of rock modulus with depth, and resulted in the prediction of a significantly deeper lower chamber at about 17 km while the upper chamber remains at around 5 km. Reanalysis using a more complete record of about 16 cGPS stations, for the full three epochs of eruption-followed-by-pause and where the basal input of magma is additionally constrained to be constant throughout the 13 year period (1995–2008), results in an even deeper lower chamber at around 19 km. The basal supply rate is constrained to be constant in the analysis, and this constant influx is independently evaluated to be approximately 1.2 m^{3} s^{−1} over the 13 year period of the analysis – similar in magnitude to the average effusion rate during this same period. This model matches the greatest number of displacement data and over the longest period of all prior models, and represents the magma as a compressible phase that honours both the history of the surface efflux and of the constrained constant input rate into the base of the system. While optimally satisfying this large array of constraints, the outputs of the model also suggest that the eruption periodicity is controlled by volumetric metering in the upper chamber (Foroozan *et al.* 2011), and thus give a robust and repeatable mechanism to explain the periodicity of the eruption sequence.

## Acknowledgments

This work is as a result of partial support from the National Science Foundation under awards NSF-CD-0607691, NSF-CD-0507334, NSF-CD-0607782 and NSF-CD-1063248, and the National Aeronautics and Space Administration award NNX11AC03G. The CALIPSO Facility was supported by NSF-IF-0523097, NSF-IF-0732728 and NSF-IF-1063246. This support is gratefully acknowledged. In addition, we also acknowledge the support of staff and colleagues at MVO, and others who have assisted in data collection. The conclusions reported here are those of the authors.

- © The Geological Society of London 2014