Chapter 14 Continuous and campaign-style gravimetric investigations on Montserrat 2006 to 2009

Abstract Gravimetric time series can provide vital clues about subsurface dynamics associated with active volcanism. Here, we report on continuous and campaign-style gravimetric observations on Montserrat between 2006 and 2009. More than 240 days of continuous gravimetric records enabled us to derive a first local joint solid Earth tides and ocean loading model for Montserrat, and we report the tidal harmonics for 14 major wave groups. Compared to global predictions, the new model (MTY11) achieves an up to one order of magnitude better precision over diurnal and semi-diurnal frequencies. We anticipate that the model will help reduce the effects of tidal perturbations on other geodetic time series recorded on Montserrat. Abrupt variations in gravity accompanied Vulcanian explosions and probably reflect the response of a shallow aquifer to stress changes during pressurization and depressurization of the subvolcanic plumbing system. Campaign data enabled the quantification of mass variations during a cycle of activity including dome formation and repose. Both forward and inverse modelling of the spatio-temporal time series indicates that the source of the recorded gravity variations is situated beneath central Montserrat. Our favourite interpretation of the campaign data is that the gravity variations reflect volcano-tectonic interaction beneath the Centre Hills of Montserrat that are triggered by changes in the active magmatic system of Soufrière Hills Volcano (SHV). We also discuss our findings on subsurface mass variations in relation to annual precipitation records and active dome formation. Both continuous and discrete gravimetric observations indicate coupling between the dominant magmatic sources responsible for the ongoing eruption at SHV and shallow-seated local sources such as aquifers and fluid-saturated fault-damage zones. Our investigations demonstrate the value of including gravimetric observations over a wide frequency range for volcanic system characterization in a volcanic island arc setting. Supplementary material: Details on the inversion routine of the explorative source model GROWTH 2.0 and the resulting images from its application to time-lapse gravity data are available at http://www.geolsoc.org.uk/SUP18701.

Repeat gravity surveys in volcanic environments permit detection of mass changes associated with volcanic activity (Rymer 1994). Relating mass changes with volume changes from simultaneously occurring ground deformation enable an assessment of the nature of subsurface processes via the determination of source density. Gravimetric time series can, thus, provide vital clues about subsurface dynamics and are suited to discriminate between different processes such as those of magmatic or hydrothermal origin . The majority of volcano gravimetric studies have been performed at restless caldera volcanoes or persistently active volcanoes (Carbone et al. 2003a;Rymer et al. 2005;. Andesitic arc volcanoes undergoing an eruptive period have hitherto not been a main target, although major insights have been gained on subsurface dynamics; for instance, at Merapi (Jousset et al. 2000) and Mayon (Jentzsch et al. 2001). Substantial financial costs, limited accessibility, poor signal-to-noise ratio and issues of operator safety are perhaps chief amongst the deterrents to establish routine gravimetric monitoring under eruptive conditions. However, the current eruptive episode of Soufrière Hills Volcano (SHV) provides an outstanding opportunity to quantify spatio-temporal mass variations during a complex episode of andesitic volcanism on Montserrat with cycles of dome growth and destruction since 1995. Unfortunately, routine gravimetric investigations have not been included for most of this activity.
Owing to the relatively small size of the island of Montserrat and safety considerations within the exclusion zone, installation and operation of a gravimetric network faces a number of challenges and requires the design of survey strategies to overcome a number of challenges.
First, arc volcanic islands by virtue of their tectonogeographical setting undergo significant tidal loading. The effect of ocean loading in addition to solid Earth tides (SETs) needs be accounted for in detail in order to obtain a suitable detection level of changes in the gravimetric acceleration of a few mGal (1 mGal ¼ 10 8 m s 22 ) across a network of benchmarks and to improve the signal-to-noise ratio. While generic models for SET and ocean loading tides (OLT) (e.g. Dehant 1987;Ray 1999;Lyard et al. 2006) may be employed to predict tidal variations, on Montserrat, for example, a pronounced pattern of residual local tidal perturbation is apparent, which remains unaccounted for by generic predictions (Fig. 14.1). This observation is in line with residual gravimetric records indicating a shortfall of global model predictions in various areas on Earth (Baker & Bos 2003) including some arc islands (Jousset et al. 2000). Long-period geodetic signals on an island such as Montserrat are, thus, significantly affected by ocean loading, which requires calibration of a precision tidal model for data reduction.
Second, the presence of an actively erupting volcano on a small island puts restrictions on network design in terms of accessibility to proximal sites during initial installation, as well as during repeat occupations.
Third, subsurface changes across the network may also affect the reference benchmark to which all gravimetric data are usually related to calculate residual changes. In particular, changes in the mass distribution at greater depth over several years probably affect the gravitational field across the entire island, including the reference, and makes their detection rather difficult.
Fourth, volcanic arcs are highly dynamic settings with the potential for concurrent volcanic, magmatic and tectonic processes. Last, but not least, owing to its location in a zone of tropical cyclone activity and associated substantial variations in annual precipitation levels, gravity perturbations due to long-and shortterm hydrological cycles need to be accounted for.
The purpose of this chapter is threefold: (i) to document recent gravimetric efforts on Montserrat in the light of the abovementioned intricacies; (ii) to show that some intricacies provide outstanding research; and (iii) to demonstrate the value of the measurements for the characterization of subvolcanic system dynamics over a broad frequency range and, in particular, the coupling between different subsurface systems between 2006 and 2009.  The 120 h-long residual gravity time series obtained at station SDV after reduction for SETs (Dehant-Defraigne-Wahr (DDW) model: Dehant 1987), ocean loading (GOT00.2 model: Ray 1999) and linear drift. Note the pronounced diurnal and semi-diurnal nature of residual variations with peak-to-peak amplitudes of up to 130 nm s 22 indicating significant imprecision in generic tidal reduction. gravity data to a tidal recording station (Gravity Consult GmbH). To avoid aliasing, data were processed through an acausal low-pass filter with 0.05 mHz cut-off frequency and averaged over 10 s intervals. The resultant 0.1 Hz gravity data were time stamped using a time transfer GPS receiver and stored on a hard-drive. A temperature logger (Dickson SM325; sampling frequency of 2 mHz), a Vasaila Digital Barometer PTB210 coupled to the tidal recording station (sampling frequency of 0.1 mHz), and a TOPCON Hiperpro dual-frequency GNSS (Global Navigation Satellite System) receiver/antenna operating at 0.5 Hz were co-located with the gravimeter (Fig. 14.2).

Methodology
The gravimeter drift was non-linear over the period of entire deployment, consistent with maturing of the instrument with time. At the end of the deployment, B28 gave an instrumental drift of approximately 10 nm s 22 hPa 21 . Individual gravity readings are precise to within +10 nm s 22 in continuous mode. Raw data were processed and analysed in Matlab and TSOFT (Van Camp & Vauterin 2005).
Tidal model calibration and data reduction. We focus here on the analysis of the continuous gravimetric (cgrav) data to derive a tidal model. The tidal analysis was performed using software package TSOFT (Van Camp & Vauterin 2005) via a remove-restore technique on 1 min binned and averaged data from the original 0.1 Hz cgrav record to obtain a 0.0167 Hz time series. In a first step, we reduced the raw gravity time series for contributions from SETs using the tidal constituents reported in Dehant (1987), which yielded a first de-tided residual time series (TS dt ). This signal was then processed for gaps, steps and spikes, as well as for other instrumental disturbances to obtain time series TS cdt1 . The resulting time series was loaded in the ETERNA tidal analysis software (Wenzel 1995) to obtain tidal parameters and the local air pressure admittance simultaneously by data fitting using a least-square adjustment. To further improve the precision of the derived constituents, the remove -restore process was repeated on TS cdt1 now using the previously estimated tidal parameters. This joint removal of both solid Earth and ocean tidal loading signals facilitates additional detection and correction of other potential disturbances to obtain TS cdt2 . Finally, a new tidal analysis was performed on TS cdt2 providing refined tidal parameters reported in Table 14.1. We derive a local air pressure admittance of 1.77 nm s 22 hPa 21 (where 1 hPa ¼ 100 Pa) from the dataset, which deviates in value and sign from data reported in the literature (Merriam 1992). We attribute this to the presence of substantial breaks and other disturbances in the time series. For the purpose of data reduction, we propagate the generic value of -3 nm s 22 hPa 21 .

Results
Residual gravimetric record at tidal frequencies. Figure 14.3 shows an 18 day de-tided gravimetric record from station SDV using: (i) SET model DDW (Dehant 1987) and OLT model FES2004 (Lyard et al. 2006) (Fig. 14.3a); and (ii) the new precision model (Fig. 14.3b; we propose the acronym MTY11). RMSE (root mean square errors) of the 0.0167 Hz time series for 5 min data bins are shown in Figure 14.3c for the July 2008 record. We calculate RMSE ¼ s/ ffi ffi t p , where t is the bin size (5 min in our case) and s is the standard deviation within each data bin. RMSE values are on average 2 nm s 22 . Data reduction using MTY11 results in a significant reduction in residual amplitude by a factor of about 5 -10 at semi-diurnal and diurnal frequencies. A power-spectrum analysis gives no pronounced peaks for the residual time series during periods of volcanic quiescence and is consistent with a random walk process of the residuals .
Example of residual gravimetric record at high frequency associated with a Vulcanian explosion. A gravity perturbation accompanied the 29 July 2008 Vulcanian explosion (VE) at SHV. The event was preceded by constant gravity signals. About 3 min after the onset of the explosion, the record shows a 4 min-long transient that peaks at þ40 nm s 22 at 03:41 UTC. The gravity signal then returns to close to average background levels approximately 10 min later. Figure 14.4 shows the gravity time series, as well as a simultaneously recorded strain signature at Air Studios (AIRS; location shown in Fig. 14.2). The gravity record peaks as the strain signal indicates maximum decompression (240 nanostrain) of the shallow plumbing system upon conduit evacuation. Estimates of erupted dense rock equivalent (DRE) volume during this event are consistent with a model of evacuation of the entire shallow cylindrical conduit to a depth of around 1 km. Gottsmann et al. (2011) forward modelled the ground deformation as a function of distance assuming conduit depressurization of 10 MPa during eruption. The predicted field of ground subsidence falls off sharply beyond 1 km from the conduit to 0 at the SDV site. To induce the transient at SDV as a result of given conduit depressurization, either a Young's modulus of encasing rocks of the order 10 MPa or a conduit with a radius of more than 100 m must be invoked. Both appear highly implausible based on the evaluation of long-term geodetic observations (Hautmann et al. 2010;Odbert et al. 2014) and the observations of dome spine diameters at SHV. Ground deformation induced by conduit depressurization is therefore dismissed as a source of the gravity transient. Our preferred model to explain the positive gravity transient is the interaction between subvolcanic system dynamics and a local aquifer. A confined aquifer has been encountered during the installation of a borehole strainmeter at OLVN (see Odbert et al. 2014) and it is plausible that the gravimetric record indicates changes in this aquifer as a result of the eruption. Plausible scenarios include changes in the hydraulic conditions of the aquifer induced as a result of seismic excitation (Brodsky et al. 2003;Gottsmann et al. 2007) or dynamic strains on conduit decompression. Both mechanisms may induce abrupt changes in the hydraulic conductivity of the aquifer, which translate into pore pressure changes.
While our investigation of the interaction between local aquifers and magmatic/volcanic activity is still ongoing, a first-order approximation to explain the gravity changes (Dg) requires a 0.25 m rise and subsequent fall of the hydraulic head (Dh wt ) in unsaturated rock with a change in density, r, of 1000 kg m 23 based on a simple approximation of the groundwater systems as an infinite sheet, and an effective aquifer open porosity (w) of 0.4 for alluvium (as a proxy for the volcanic debris deposits: Lin & Lee 2007): where G is the gravitational constant. Albeit a rough first approximation, these water-table variations do not seem implausible (see Brodsky et al. 2003 for the case of well-level changes caused by distant earthquakes). At this stage, our hypothesis is that a combination of poroelastic effects, changes in the hydraulic conductivity and mass variations are the cause of the gravity transient. Further investigations are currently underway.

Discussion
The tidal model accounts for the effects of SET and ocean loading as well as the elastic crustal response, where the latter could be affected by the volcanic setting. The average residual diurnal and semidiurnal amplitude is about +10 nm s 22 , which represents the background gravity variation at such frequencies. We attribute the residuals to instrumental and environmental effects. Environmental factors include thermal expansion of the building, minor ground tilt and delayed aquifer responses due to tidal loading, all of which are common contributors to data uncertainty from spring gravimeters. Our available data do not allow us to quantify residual volcanic loading effects associated with semi-diurnal or diurnal frequencies over the observation period as they manifest an integral part of the time series and are accounted for by the de-tiding procedures. However, the model provides an acceptable level of certainty for the reduction of dynamic gravimetric investigations and potentially for other geodetic time series (such as strain, tilt and GPS) affected by tidal perturbations on Montserrat.

Time-lapse observations
Methodology Network set-up and data collection. Spatio-temporal gravity data were recorded across a network of benchmarks that covers the northern part of Montserrat island (Fig. 14.2). After having been installed and first deployed in June-July 2006, the network was reoccupied in January -February 2007 (first observation interval, 7 months) and August -September 2008 (second observation interval, 19 months). The original network comprised 10 sites with a distance to adjacent benchmarks of 2 -3 km, and a minimum distance of 4 km between STGH and the active vent. The northernmost benchmark (SLVH) was primarily considered as the base station to which all other measurements are referenced. After the loss of SLVH due to construction works in 2008, however, the neighbouring benchmark BASE in the NE of the network was taken as reference to calculate residual gravity changes to other sites. The gravity data were sampled via optical reading using a LaCoste & Romberg gravity meter (G-667). In each survey the gravity measurements were repeated up to 10 times at the individual benchmarks. In order to account for elevation changes, static GNSS occupations were performed using Leica 500, Trimble 5700 and Topcon HiperPro dual-frequency receivers operating at each site for 1 h at 0.5 Hz recording frequency.
Reduction for instrumental drift and tides. The sampled gravity data were reduced for SETs (Dehant 1987) and instrumental drift. The collection of continuous gravity data simultaneously with the first time-lapse survey campaign in 2006 on Montserrat revealed the gravity field to be significantly influenced by ocean loading that cannot be entirely eliminated using global loading models. Owing to the lack of a precise tidal model for Montserrat or, indeed, any other nearby islands of the West Indies at the time the dynamic gravity campaigns were conducted, a survey strategy was developed and applied that allowed for the reduction of gravity data, in particular semi-diurnal variations, such as ocean loading effects. The base station was reoccupied every 2 h, and interpolation every minute of the resulting differences permitted readings from other benchmark sites taken in between the base station readings to be corrected. Following this measurement protocol, the precision of multiple readings was within +10-100 nm s 22 for most of the benchmarks. An exception was found for JBYH, where repeated readings resulted in larger differences +200 nm s 22 (presumably due to anthropogenic disturbances). Owing to the resulting low signal-to-noise ratio, the data from JBYH were excluded from further analysis.
We reprocessed the raw time-lapse gravity data with the new precision model MTY11 in order to compare the quality and appropriateness of the above-described survey strategy. The direct comparison of the de-tided time series from both methodologies revealed that, within the uncertainties, both datasets match. This validates the initial survey strategy to yield highly precise data. This strategy can thus be recommended for future surveys in other volcanic settings that are affected by short-term gravity oscillations which cannot be fully eliminated with generic ocean loading models.
Correction for free-air effect. In order to account for the free-air effect -that is, to reduce the gravity data for changes in elevation caused by processes in the subsurface -GPS data were collected simultaneously with the gravity data along the network. Our postprocessed data, however, revealed vertical benchmark displacements only within the precision of GPS survey observations (i.e.,2 cm). In contrast to static campaign occupations, continuously recorded GPS (cGPS) data have a markedly higher precision. We therefore applied cGPS data that are routinely recorded on Montserrat (MVO and CALIPSO sites) in order to infer on small-scaled vertical surface displacements in our gravity survey area (Hautmann et al. 2010;Mattioli et al. 2010). The data revealed differences in vertical surface deformation in northern Montserrat of 0.5 and 1.5 cm during the first and second observation periods, respectively. As the positions of the cGPS network sites do not coincide with our gravity benchmarks, we cannot precisely constrain the individual uplift/subsidence amplitude at each of the individual gravity sites. We hence propagate in the residual gravity data the associated uncertainty due to vertical displacement of the benchmarks. The theoretical value for the free-air gradient, which represents the change in gravity with an elevation change, is 230.1 nm s 22 m 21 . According to the uncertainties in vertical surface deformation in each of the observation periods (0.5 and 1.5 cm, respectively), we propagate an error of +20 and +50 nm s 22 to data from the first and second monitoring interval, respectively.
Influence of mass redistribution at SHV. The first monitoring interval (June/July 2006 -January/February 2007) covers a period of volcanic deflation (effusive phase 3) at SHV, with a pressure decrease in the subsurface as magma ascends from the reservoir and extrudes at the surface. In contrast, the second observation interval (January/February 2007-August/September 2008) represents the end of Phase 3 (20 April 2007) and mainly a phase of volcanic quiescence (pause 3) that correlates with surface inflation due to internal pressure build-up as result of fresh magma intruding into the reservoir. Changes in the surface gravity field derived from magma extrusion in times of volcanic deflation result from: (i) mass loss at depth due to magma withdrawal out of the reservoir; and (ii) mass gain at the volcanic vent due to dome growth. Both effects are predicted to generate a negative gravity anomaly across the survey network. In contrast, intrusion of fresh material into the magmatic system (volcanic inflation) leads to a mass gain at depth and, thus, larger gravitational attraction. The effect of a point mass m on vertical gravitational attraction, Dg, can be quantified via: where G is 6.67 Â 10 211 m 3 kg 21 s 22 , z is the source depth (positive downwards) and r the distance to the mass. Routinely conducted monitoring of the active dome carried out by the Montserrat Volcano Observatory (MVO) provide constraints on the extruded magma volumes, which are estimated at 140 Â 10 6 m 3 DRE for the time between the first and second gravity surveys. In addition, cGPS data analysis from both survey intervals (Hautmann et al. 2010;MVO pers. comm.) reveals the surface deformation to be generated by decompression/compression of the mid-crustal magma chamber. Results from surface deformation and strain data modelling (Hautmann et al. 2010(Hautmann et al. , 2013 constrain the mid-crustal reservoir to be vertically elongated with its centre at 13 km below the surface. The source volume estimate is 25 km 3 . For the first monitoring interval, we hence assume m ¼ 23.5 Â 10 11 kg and z ¼ 13 000 m (to account for magma withdrawal out of the magma chamber), as well as m ¼ 3.5 Â 10 11 kg and z ¼ 2500 m (to allow for dome loading). Predicted gravity changes (Fig. 14.5) relative to BASE show that variations due to deep reservoir extrusion between June -July 2006 and January -February 2007 are too small to be resolved in a microgravity survey. The limiting factors in order to detect the changes generated by the deep reservoir are the large signal wavelength and the comparably small size of the island/network. In contrast, dome loading markedly influences the surveyed gravity field, with a gravity contrast of between 2300 nm s 22 at STGH and 10 nm s 22 at BASE inferred for the first observation interval. As the surface deformation rates are similar (although opposite) for the first and second survey interval, we assume the intrusion/extrusion rate of magmatic material into/out of the deep reservoir to be the same, but accounted for the longer survey interval between the second and the third campaign. We thus quantified the effect of deep magma intrusion by assuming m ¼ 8 Â 10 11 kg and z ¼ 13 000 m. The results (Fig. 14.5) show that the gravity changes in the survey area vary between 80 and 240 nm s 22 , and are thus slightly smaller than the uncertainties in the recorded data.
Note that our calculations are based on mass changes inferred from visual observations of dome growth at the active vent. Additional dynamics, such as density increase due to degassing or mass gain due to influx of fresh magma into the deep reservoir, are not accounted for as quantification of these processes bears large uncertainties, while at the same time their influence on the gravity field is only minor (e.g. assuming a constant influx rate of 2 m 3 s 21 into the deep reservoir, as suggested by Voight et al. 2010, affects predictions on gravity changes by c. 10 nm s 22 ).
For each benchmark, we applied the inferred correction for the observed mass redistribution at SHV to the recorded gravity data, in order to derive a residual gravity signal that documents only the hitherto undefined subsurface mass variations on Montserrat. The residual time-lapse gravity data thus reflect de-tided and de-drifted values corrected for observed mass redistribution at SHV, and errors accounting for measurement uncertainties and potential free-air effects.

Results
Residual gravity data from the first monitoring interval display a negative gravity anomaly at a minimum of 2350 nm s 22 in NW Montserrat (site STPT), together with positive gravity changes of 150 and 440 nm s 22 at the southernmost sites GRBD and STGH, respectively. In contrast, residual gravity data from the second observation period document positive gravity variations at all stations, and a maximum of 700 nm s 22 at MVO2 (see also  Table 14.2). The spatial distribution of the gravity changes from both the first and second observation interval is coherent over the survey area. The anomalies are not centred in the SHV region, but in the Centre Hills, with the 2006-2007 gravity anomaly located about 2 km north of the focus of the 2007-2008 anomaly. Based on the half-width method (Skeels 1963) of the gravity anomalies projected on a NE -SW profile, the causative sources (assumed in a simple horizontal cylindrical geometry) are both estimated to be at a maximum depth of 3.3 km. It is, however, important to note that any assumptions on maximum source depths reflect only rough approximations as the signal cannot be captured in its entire wavelength, because of the small size of the island and, hence, the survey area. A Gaussian integration of the net mass changes yields 22.3 Â 10 10 kg for the first survey period and þ4.4 Â 10 10 kg for the second, which is about an order of magnitude smaller than what was observed at the same time as mass transfers at SHV.

Data modelling
To obtain a model of the sources of the recorded anomalies, the time-lapse gravity data were analysed using different modelling approaches. Owing to the absence of significant ground deformation over the survey period, causative source volume changes cannot be quantified. As such, the gravity residuals must be regarded to primarily reflect mass changes. Gravity data display a coherent distribution of residuals from all sites, with exception of STGH. This suggests that STGH is possibly affected by a different gravity source, which could either be the upper plumbing system of SHV or the upper Belham Valley. The lack of additional data, however, limits our understanding of processes that trigger the gravity variations recorded at STGH. Residuals from STGH are, therefore, excluded from the data modelling presented below. Data from the first survey interval were offset by a value Dg 1off ¼ 2160 nm s 22 (Dg 1off is the offset value applied to the gravity data from monitoring interval I), as analysing the data with negative and positive residuals would require the solutions to have more than one source, which appears unwarranted given the coherent distribution of the residuals.
Finite source models. As a first approach to assess characteristics of the gravity source, we employ a data inversion based on the simple assumption of a spherical source. We then run forward models that assume a source geometry of a horizontal cylinder, a horizontal sill and a horizontal dyke. For details on the applied methods, the reader is referred to Battaglia et al. (2008). In accordance to the preferred direction of major tectonic features on Montserrat and the surrounding seafloor (faults, alignment of older lava domes), we orientated the elongated sources in a NW -SE direction. Given the spatial distribution of the recorded changes, we centred the sources at the maxima of the anomalies (i.e., first anomaly: N1854900, E584450 and second anomaly: N1852500, E583250). Model predictions are similar for both observation periods, yielding source depths of 2200 m for a spherical and a cylindrical source, 800 m for a sill (with a best-fit lateral extension of 4400 m), and a vertical extension of between 1200 and 2000 m depth for a dyke source (see Table 14.3 for a summary of model-inferred source characteristics). The inferred source sizes -that is, the radius or thickness -depend linearly on the assumed density contrasts, with small density contrasts yielding larger source sizes and vice versa. In particular, for a density contrast of 400 kg m 23 (corresponding to the density difference between solid rock (2700 kg m 23 ) and magma (2300 kg m 23 )), the inferred radii of the spherical and cylindrical sources are 280 and 80 m, respectively, for the first monitoring interval, while the inferred thickness of the dyke and the sill is 17 and 3.5 m, respectively. For a density contrast of 1700 kg m 23 (corresponding to the difference between solid rock and aqueous fluids), model-derived source volumes are consistently smaller by a factor of 4.25. In accordance to the differences in signal amplitudes from both observation periods, the model-derived source volumes The errors in the data account for measurement uncertainties and potential free-air effects. The data show a coherent distribution of the gravity changes, with the exception of the records from STGH. It is, thus, assumed that the STGH site is additionally influenced by magma dynamics in the shallow plumbing system of SHV. Data from STGH are, therefore, excluded from the analysis of the gravity anomalies that are associated with the dynamics in the western Centre Hills. for the second observation period are about a factor of 1.3 larger compared to the volumes inferred for the first monitoring interval.
Assuming a source length of 5 km for the sill, dyke and cylinder sources, the model-inferred mass changes fit the integrationderived changes. The same can be confirmed for a spherical source. For the first monitoring interval, the average fit quality of model to data is within 50 nm s 22 for a sphere, and within 40 nm s 22 for a dyke, sill and cylindrical source geometries (Fig. 14.7). The fit quality slightly worsens for the second monitoring interval, however, the average misfit is still smaller than the data uncertainty. A chi-square (x 2 ) test revealed no statistically significant differences (p , 0.01) between data and model predictions from the various approaches. In particular, for the first observation interval, x 2 is 0.89 for a sphere, 0.53 for a dyke, 0.69 for a sill and 0.65 for a cylindrical source. Values of x 2 are similar also for the second observation interval. Although this result might slightly hint towards a dyke as preferred source geometry, the quality of fit is similar for all source models and fit residuals from all approaches are within the error of the measurements, thus, not allowing us to favour one model over the other.
Explorative source models. In order to determine a free geometry rather than a finite geometry for the causative source(s), we performed a non-linear inversion using the inversion package GROWTH 2.0 (Camacho et al. 2011). The inversion is based on a 3D aggregation of M parallelepiped cells, which are filled, in a growth process, by means of prescribed positive and/or negative density contrasts. This methodology provides, via an automatic approach, a free 3D geometry of the anomalous body, which matches as much as possible the observed gravity anomaly.
Inverting the gravity residuals, we obtain models that show cylindrical to slightly flattened (i.e., dyke-like) structures that are WNW -ESE directed, dipping steeply towards the ESE. The inferred sources are situated in the western Centre Hills, matching the locations of source centres for the finite source models. The inferred mass loss during the first survey interval is 1.4 Â 10 10 kg, while the inferred mass gain during the second observation period is 5.5 Â 10 10 kg. The weighted standard deviation of the inversion residuals is within 16 nm s 22 , representing an excellent quality of fit-to-the-data.
The model results support our a priori assumptions on source orientation and location. The results of the explorative approach indicate that the time-lapse gravity anomalies are caused by cylindrical to dyke-like sources, although the inversion fails to resolve sources at very shallow depths (i.e. sources with a sill geometry). We cannot, therefore, unambiguously deduce dykes to be the causative source geometry behind the gravity variations. In order to constrain on the gravity sources, we thus need to make assumptions on the nature of the observed mass changes. These are discussed in the following section.

Discussion and interpretation
In order to assess whether the long-period gravity changes are affected by seasonal changes in the cumulative rainfall and hence, in the groundwater level, we estimated the aquifer head changes that are required in order to generate the observed gravity variations. The Bouguer sheet attraction corresponds to 419 (nm s 22 ) m 21 Â w, with w being the open porosity of the host rock. Assuming an average open porosity of 0.2 for Montserrat (Moreau-Fournier 2008), the recorded gravity changes from the first observation interval would be induced by a groundwater level drop of 5.9 m, while the data from the second monitoring period would correspond to a rise of the aquifer head of 8.3 m. A comparison with precipitation data from Montserrat (Montserrat Water Authority pers. comm.; Jones et al. 2010) shows, however, that groundwater level changes generated by seasonal variations in the cumulative rainfall would generate a gravity signal that is: (i) inverse to the observed gravity changes; and (ii) about an order of magnitude smaller than the recorded data. We thus disregard seasonal climatic changes as cause for gravity variations on Montserrat.
As the recorded data are corrected for the gravitational effect of mass redistribution in the SHV active magmatic system, the gravity signals cannot be explained with magma movement due to eruptive activity. In addition, the centre of the source responsible for the documented gravity variations must lie beneath the Centre Hills rather than beneath SHV. The Centre Hills are, however, an extinct volcanic complex with last known activity at 400 ka (Harford et al. 2002), making current magmatic activity in this region highly unlikely. In addition, the generation of the recorded gravity changes correlates with the redistribution of 75 Â 10 6 m 3 magma at shallow depths. This would simultaneously produce near-field-surface deformation of the order of tens of centimetres if one forward models a simple elastic crustal response to magma intrusion in a spherical source (Mogi 1958). Ground deformation at such amplitude was not recorded in the Centre Hills during the surveys.
The polarity of the gravity anomalies varies with the activity of SHV, as the negative gravity changes found in the first observation interval were recorded during volcanic deflation (effusive phase 3), while the positive gravity changes documented in the second dataset relate with a period of volcanic inflation (pause 3). Thus, it can be assumed that dynamics in the volcanic system of SHV trigger secondary processes that are documented in the long-period gravity data. The distribution of the observed gravity variations clearly indicate a non-isotropic structure of the Centre Hills subsurface. The gravity anomalies are centred along a NW -SEdirected alignment, which corresponds to the prevailing direction of volcano-tectonic features that are documented on-island and offshore of Montserrat (e.g. the shallow plumbing system at SHV: Mattioli et al. 1998;Linde et al. 2010; submarine fault systems such as the Bouillante -Montserrat Graben: Feuillet et al. 2002; alignment of earlier lava domes at SHV, and the direction of faults and older horst structures: Harford et al. 2002). In particular, the gravity anomaly of the second observation interval is centred in the close vicinity of the Belham Valley Fault, which separates the volcanic complexes of SHV in the western part of the island from Centre Hills. Although the location of the gravity anomaly reported in the data from the first survey period is hitherto not known as a fracture zone, there are several indications for the existence of a fault in this area: (i) the anomaly is centred beneath the Soldier Ghaut, which is a markedly steep and prominent valley on Montserrat; (ii) mapping of roadcuts in this area revealed a number of small-scale faults and hydrothermal veins that are directed north -south to NW -SE with a dip angle of between 608 and 908 ( Fig. 14.8); (iii) seismic streamer data report young, faulted, submarine sediments to the west of Montserrat (Kenedi et al. 2010); and, (iv) the lateral extension of the Soldier Ghaut towards the NW coincides with an offset of the shelf edge. We hence assume that the locations of both gravity anomalies are associated with an active fault/fault zone, which we refer to in the following as Belham Valley Fault and Centre Hills Fault (Fig. 14.8).
Long-term gravimetric observations from Etna Volcano (Carbone et al. 2003b(Carbone et al. , 2009 revealed that changes in the local stress field can induce a dynamic response of dilation and compression of microfractures in pre-existing zones of structural weakness. For the case of Montserrat, the negative gravity anomaly from the first observation period could, thus, correspond to fracture dilation due to stress relaxation during deflation. Consequently, the positive gravity anomaly from the second monitoring interval could be associated with fracture closure and, hence, medium compression triggered by stress increase due to the onset of a new inflation period. The resulting mass changes deduced for Etna Volcano are of the same order of magnitude as that inferred from the time-lapse gravity records from Montserrat. As to enhanced fracture permeability, fault damage zones can provide traces for stress-induced fluid migration (Tamagawa & Pollard 2008). Groundwater level changes along the Belham Valley Fault and the Centre Hills Fault in response to stress cycling at SHV might thus be an alternative or additional mechanism to generate the observed gravity variations. Water-level fluctuations of up to 7 m triggered by changes in volcanic stressing are also reported from other active volcanoes (e.g. Kilauea: Hurwitz & Johnston 2003;Mount Usu: Shibata et al. 2008). The same scale of aquifer head oscillations would also produce a gravity signal similar to that we recorded at SHV (see the calculations given above). The hydrological system of Montserrat is not yet explored in great detail, although it is known that the Belham Valley hosts two stacked aquifers separated by a clay aquiclude (Young 2002;Davies & Peart 2003). Recently, changes in the Belham Valley aquifers revealed that these systems can evolve rapidly, which might be caused by volcanic activity at SHV. Between September 2010 and February 2011, a groundwater production well located in the lower Belham Valley turned into an artesian well. As the borehole was surrounded by a concrete building in order to protect it from lahars, the change to an artesian well was not recognized immediately and the precise onset of overflow is, therefore, unknown. Regardless of this, the event implies a rise of several metres of the aquifer head in the subsurface of the Belham Valley prior to the overflow of the borehole. It is thus possible that the positive gravity changes recorded between 2007 and 2008 (i.e. the second observation interval) reflect the onset of groundwater-level rise.
We thus propose that the recorded long-term gravity variations reflect subsurface dynamics in central Montserrat that are triggered by changes in the SHV active magmatic system. We put forward two mechanisms to explain, either separately or in combination, the recorded gravity signals: (i) fracture dilation/compression in fault/fault zones as a tectonic response to perturbation in the stress field; and/or (ii) fluctuations in the hydrological system with groundwater flow channelized along fracture zones. While the scenario (i) would correspond to a dyke source geometry, (ii) would relate to a sill geometry. Models to fit the recorded data, however, did not allow us to discriminate between either source geometry as either source fits the data equally well.

Implications and outlook
Time-lapse gravity variations on Montserrat recorded between 2006 and 2008 are best explained by subsurface mass changes beneath the Centre Hills triggered by changes in the eruptive activity at SHV. In order to better assess and quantify the processes related to the mass changes in the subsurface, additional data acquisition is required. The collection of new time-lapse gravity data is currently underway, wherefore the benchmark network of the new gravity surveys has been adjusted based on the findings of the study presented herein. The adjustments include: (i) an increase in the number of survey benchmarks and an overall densification of the network in the Centre Hills in order to better constrain the minimum depth of the gravity anomalies; (ii) an extension of the network towards the north and the south of the island in order to better image the wavelength of the signal; and (iii) the installation of observation sites around SHV, in order to investigate dynamics that primarily relate to magmatic processes in the shallow volcanic plumbing system (intrusion, effusion, degassing). In particular, the installation of observation sites around SHV will help to investigate and understand the gravity signals recorded at STGH, which are either affected by dynamic processes in the shallow SHV plumbing system or the upper Belham Valley. Additional data from southern Montserrat, however, are crucial in order to distinguish and quantify the effect of volcano v. aquifer dynamics on the gravity field in this area.
Absolute gravity measurements could provide constraints on the signal amplitude independently from the wavelength of the signal. The employment of an absolute gravity meter, however, is very cost-intensive and, therefore, currently not feasible.
For the first time, the newly recorded datasets will not only be evaluated and compared with activity monitoring data from SHV, but also with respect to groundwater-level data from the island, which will be collected in order to explore possible hydrological responses triggered by changes in magmatic activity and to quantify its contribution to recorded geophysical signals.
Preparations to install a network of continuously recording gravimeters on Montserrat are underway. This installation aims to quantify spatio-temporal subsurface mass and density variations over periods of seconds to a few months, and, therefore, to bridge the gap between time-lapse gravimetric observations and continuous routine monitoring efforts by providing insights on subsurface processes at higher frequencies than achievable by conventional micro-gravity surveys. Initial results of exploiting the cgrav time series in terms of eruption mechanics indicate a substantial potential to investigate the effusive to explosive transition in domeforming eruptions via multi-component geophysical observations (Gottsmann et al. 2011). Further gravimetric studies over a wide frequency range at active island arc volcanoes will, undoubtedly, provide new and exciting insights on the inner workings of such volcanoes.