The Transcrustal Magma Reservoir Beneath Soufrière Hills Volcano, Montserrat: Insights From 3‐D Geodetic Inversions

We invert intraeruptive ground displacements recorded between 2003 and 2005 on Montserrat to shed light on the magmatic plumbing system of Soufrière Hills Volcano. Incorporating 3‐dimensional crustal mechanical and topographic data in a finite‐element model, we show that the recorded displacements are best explained by a southeastward dipping (plunge angle of 9.3°) vertically extended triaxial ellipsoidal pressure source with semiaxes lengths of 1.9 and 2.0 km horizontally, and 5.0 km vertically. The source is centered at 9.35 km depth below main sea level and embedded in independently imaged anomalously weak crustal rocks. The source orientation appears to be controlled by the local stress field at the intersection of two major WNW‐ESE and NW‐SE striking tectonic lineaments. We derive an average volumetric strain rate of 8.4 × 10−12 s−1 by transcrustal pressurization which may have contributed to flank instability and mass wasting events in the southern and eastern sectors of the island.


Introduction
Volcano deformation can be caused by a variety of subsurface processes including the transport of magma from depth to the shallower regions of the crust (Dvorak & Dzurisin, 1997), cooling and crystallization of a melt body (Parker et al., 2014), degassing through decompression or phases changes in a magma reservoir , and hydrothermal processes (Bartel et al., 2003). Advances in remote sensing techniques have allowed for global coverage of the world's deforming volcanoes  and provide a means to study the relationship between deforming volcanoes and eruption . Inversion of deformation data can provide estimates on the location, shape, and volume changes in magma reservoirs supplying volcanoes (Anderson & Segall, 2011).
Analytical geodetic models rely on simplified assumptions regarding the geometry of the source undergoing stress changes and/or the mechanical behavior of surrounding rocks (e.g., Anderson, 1936;Fialko et al., 2001;Mogi, 1958;Yang et al., 1988). However, geochemical, petrological, and geophysical data illustrate complex architectures of subvolcanic plumbing systems and crustal rocks Magee et al., 2018). These heterogeneities fundamentally influence the stress versus strain relationship and therefore require detailed analysis beyond the isotropic, homogenous, and elastic (IHE) half-space approximation embodied in traditional geodetic models. Numerical models provide an opportunity to account for complexities found in volcanic areas, such as significant topography, crustal mechanical heterogeneity, variations in subsurface temperatures, and, hence, inelastic rheological behavior of crustal material (Currenti & Williams, 2014;Del Negro et al., 2009;Gregg et al., 2013;Hickey et al., 2015;Zhan et al., 2019). Differences in derived source parameters from IHE half-space models versus 3-D numerical models can be significant with implications for hazard assessment and risk mitigation (e.g., Hickey et al., 2016).
The Soufrière Hills Volcano (SHV) on the island of Montserrat (British West Indies; Figure 1) in the north of the Lesser Antilles island arc in the Caribbean is a medium-sized andesitic, composite volcano (Young et al., 1998) and one example where crustal mechanical complexities influence the subsurface stress versus strain relationship Hautmann et al., 2010a;Odbert, Ryan, et al., 2014;Young & Gottsmann, 2015). Volcanism on Montserrat has migrated from North to South over 2.6 Ma and formed four volcanic complexes (from North: Silver Hills, Centre Hills, Soufrière Hills, and South Soufrière Hills) which align NNW-SSE ( Figure 1). The most recent eruptive cycle began in 1995 on SHV and has so far had five distinct eruptive phases, characterized by periods of lava dome formation and destruction and intermittent periods of repose (Odbert, Stewart, et al., 2014). Becoming one of the best studied eruptive sequences of andesite volcanism (see Druitt &Kokelaar, 2002, andWadge et al., 2014, for further details), magma extrusion ceased in 2010, and since then, the volcano has been in a period of repose (Scientific Advisory Committee, 2018). The protracted nature of volcanism has left distinct mechanical fingerprints (Figure 2b) in the upper crust beneath Montserrat including zones of mechanical weakness beneath the southern part of the island, dense cores of relict magmatic plumbing systems beneath the Silver and Centre Hills and superficial pyroclastic aprons at SHV (Paulatto et al., 2010(Paulatto et al., , 2012. It is therefore conceivable that the substantial 3-D mechanical heterogeneity in the middle and upper crust promotes strain partitioning (i.e., the heterogenous distribution of strain amplitude and type governed by rheological contrasts of crustal rocks)  Baird et al., 2015;Feuillet et al., 2010;Hautmann et al., 2014). White broken lines highlight key geomorphological features associated with flank collapses at SHV and SSH, including the prominent collapse scars of English's crater (EC) and from the 1997 Boxing Day collapse (BD) (after Le Friant et al., 2004). The model accounts for onshore and offshore topography (a) in addition to a mechanically heterogeneous crust (b). (a) shows the Cartesian coordinate system (x, y, z) of the model, boundary conditions (fixed, roller, free) as applied to the main modeling domain (shown in gray), and the infinite element domain (shown in blue). Starting in 1995, the most recent episode of eruptive activity from Soufrière Hills Volcano (SHV) was accompanied by island-wide ground deformation recorded at a network of GPS receivers (Odbert, Ryan, et al., 2014). For the current study, we use eastward (x direction) and northward (y direction) and vertical (z direction) displacements U recorded between 2003 and 2005 at nine stations (marked by black dots; see also Table 1). (b) The colormap shows the distribution of the static Young's modulus in the domain which has been derived from seismic tomography data (Paulatto et al., 2012) using Equation 2 and scaled to static values. Thermal and mechanical alteration by protracted magmatic activity may explain the crustal weakening beneath the southern part of Montserrat where the best fit deformation source is located.

Geophysical Research Letters
from subsurface stress changes in the magmatic plumbing system. These complexities have not been captured in earlier geodetic models. Here we present the first inversion of intraeruptive ground deformation data at SHV accounting for 3-D topography and 3-D mechanical heterogeneity of the crust beneath Montserrat.

Model Formulation
We numerically solve the elastic stress and resultant deformation field induced in a mechanically heterogenous material by a buried pressure source such that the divergence of the elastic stress tensor is zero. Computations are performed in COMSOL Multiphysics® 5.4 using a 3-D modeling domain in the structural mechanics module coupled with an optimization procedure to solve the inverse problem (Hickey et al., 2015(Hickey et al., , 2016. A typical model uses~47,000 tetrahedral mesh elements and solves for~167,000 degrees of freedom. The main modeling domain is surrounded by an infinite-element domain to prevent boundary effects ( Figure 2a). The topographic data to construct the digital elevation model of Montserrat are from pre-1995 SRTM data, and bathymetric data are from Le Friant et al. (2004).
The magma reservoir is represented by a triaxial ellipsoidal cavity with semiaxes a, b, and c such that The orientation of the ellipsoid in 3-D space is expressed by spherical coordinates using the following notations: The plunge of the ellipsoid is expressed by polar angle θ measured from the vertical, while ɸ is the clockwise azimuthal angle measured from the North. Rotation angle κ is about the ellipsoid's c axis and is measured counterclockwise from a to b. For κ = 0°, the ellipsoid's b axis is parallel to the xy plane.
Reservoir pressure change is applied as a boundary load on the ellipsoid. The influence of temperature-dependent crustal rheology and reservoir properties on the stress and strain relationship are explored by Gottsmann and Odbert (2014) for the 2-D case and not discussed further here. The final domain size for all models is 80 km × 80 km × 31 km, with a higher mesh density applied around the source, the free surface, and Global Positioning System (GPS) site locations. Benchmarking against analytical IHE half-space models was performed according to procedures outlined in Hickey and Gottsmann (2014).
We implement the 3-D variations in crustal rock mechanics using seismic compressional wave velocities (v p ) presented in Paulatto et al. (2012) to parameterize rock density ρ r and dynamic Young's modulus E d following (Brocher, 2005): where ν is Poisson's ratio. In the absence of shear wave velocity data, we select an average Poisson's ratio (ν = 0.27) appropriate for the crustal rocks of Montserrat (Sevilla et al., 2010). Static values E can be a factor of a few to up to an order of magnitude lower than the dynamic modulus (Cheng & Johnston, 1981;Gudmundsson, 1983;Heap et al., 2020); we employ a scaling factor of 0.5 to convert E d to its static value E s , such that E s = 0.5E d . The resulting distribution of E s across the modeling domain is shown in Figure 2b.

Inverse Model
We solved the inverse problem by applying an optimization routine (e.g., Hickey et al., 2016) to minimize an objective function J (an indicator of model misfit between the observed data D and modeled data M; supporting information Table S1) using a weighting, derived from the root mean square error (RMSE) σ of the GPS data. The algorithm searches for the optimum deformation source parameters (longitude, latitude, depth to source center, semimajor and semiminor axes, κ, ɸ, and θ) to best fit all three components of surface displacement (northward, eastward, and vertical) at each GPS site i using where N is the total number of GPS sites. The initial optimization followed a Monte Carlo approach by using random sampling with uniform distribution within a range of values of the design variables presented in Table 1. Data from earlier geodetic studies were used to inform the value range and geometry of the reservoir Hautmann et al., 2010aHautmann et al., , 2010bOdbert, Ryan, et al., 2014). Following the Monte Carlo optimization, a refinement routine using a Bound Optimization BY Quadratic Approximation algorithm was used to further minimize the objective function ( Figure S1). Bound Optimization BY Quadratic Approximation is a derivative-free bound-constrained optimization algorithm that solves a trust region subproblem to promote good linear independence in the interpolation conditions during an iterative minimization process (Powell, 2009). Our iterative and repetitive optimization approach aims to reduce the limits of design variables and builds nested parameter constraint grids to ensure the most robust final solution via the smallest misfit. To avoid finding solutions from local minima, we performed several iterations with different starting values of the design variables.

GPS Data
Ground deformation data on Montserrat have been collected since 1995 using a variety of methods (Odbert, Ryan, et al., 2014). Here we focus on the intraeruptive period of ground uplift recorded between July 2003 and August 2005. This episode marks a period of relative quiescence after Dome Formation Phase 2 and prior to Dome Formation Phase 3 and provides us with the opportunity to probe the recharge of the magmatic system. For the data inversion, we use displacements obtained from time-integrated and linearly averaged three-component (Wessel & Smith, 1998) (eastward Ux, northward Uy, and vertical Uz) daily positions of nine cGPS stations. Daily positions were calculated using the GAMIT/GLOBK code (Herring et al., 2010) (see also Text S1, Table S1, and Figure S3).

Results and Discussion
Figures 2c and 2d show predicted eastward (u), northward (v), and vertical displacements (w) from the optimal inversion results (Table 1) with respect to observed displacements recorded between 2003 and 2005. The optimal 3-D geometry of the deformation source is a triaxial ellipsoid centered at 9.35 km below main sea level with semiaxes (a, b, and c) lengths of 1.9, 2.0, and 5.0 km, respectively, undergoing a pressure change of 11 MPa, equivalent to a volume change of 0.043 km 3 . The volume of the ellipsoid is 80 km 3 , and it plunges at an angle of 9.3°toward the SE. Uncertainties in the observed displacement vectors have only a minor influence on the size and orientation of the ellipsoid (Table 1) with plunge angles ranging between  (Figure 2).
Published geodetic models of the upper-crustal magmatic plumbing system at SHV can be broadly subdivided into single or multiple source models. Single source models propose a vertically elongated structure and are best approximated by prolate ellipsoidal bodies Voight et al., 2010); multiple source models are a vertical superposition of sources, either spherical or a combination of spherical and ellipsoidal (Elsworth et al., 2008;Hautmann et al., 2014). Our results provide the first 3-D image of a single geodetic source at SHV driving the observed elastic ground deformation, taking into account 3-D variations in topography and crustal mechanics. Our model does not account for inelastic thermomechanical behavior at SHV, which for the 2-D case has been explored in Gottsmann and Odbert (2014). We expect a similar (by a factor of a few) decrease in required excess pressure for the thermomechanical conditions explored by that study.
The source volume change (ΔV s ) of 0.043 km 3 over the 2 year period (giving a volume strain of 5.3 × 10 −4 and average volume strain rate of 8.4 × 10 −12 s −1 ) is at the upper limit of inelastic volume changes (0.029 and 0.042 km 3 ) proposed by Gottsmann and Odbert (2014). Although volume change estimates from elastic models are often larger than those from viscoelastic models, our elastic models yield similar volume changes as models incorporating a pressure source in a viscoelastic medium due to the weak crustal rocks beneath SHV and 3-D subsurface strain partitioning not previously considered. However, the temporal evolution of purely elastic crustal deformation is different to deformation resulting from a time-dependent rheology; continuous ground deformation records should help discriminate mechanical from thermomechanical effects (Head et al., 2019). The interested reader is referred to Gottsmann and Odbert (2014) for an analysis of time-dependent crustal deformation at SHV.
A volume of 0.28 km 3 of lava (V ex ) extruded during Eruptive Phase 3 between August 2005 and April 2007. The ratio r = V ex /ΔV s = 1 + ß m /ß r (Rivalta & Segall, 2008) relates extruded magma volume to the source volume change. With a value of r = 6.6, our results imply a major role of magma compressibility ß m in modulating eruptive behavior at SHV. ß r is the reservoir compressibility and is given by 1.65(1 + ν)/E for the derived source geometry (Anderson & Segall, 2011). For values of E in the range of 30-40 GPa (Figure 2), we calculate ß m values between 3 and 4 * 10 −10 Pa −1 which match published values in the range of 4 * 10 −11 to 10 −9 Pa −1 (Christopher et al., 2015;Edmonds et al., 2014;. A vertically elongated source geometry is consistent with the model of transcrustal magmatic systems fueling many active volcanic systems , including SHV (Edmonds et al., 2016). Joint inversion of gravity and seismic data from SHV (Paulatto et al., 2019) revealed middle-and upper-crustal anomalous velocity and density structures consistent with the proposed source location, geometry, size, and volume, with the top of our modeled deformation source corresponding to the −400 m/s v p anomaly.
SHV is located at the intersection of two dominant and active tectonic systems (Figure 1) (Feuillet et al., 2011): the NW-SE striking Centre Hills fault segment as part of the extensional Bouillante-Montserrat fault system and the W-E to WNW-ESE striking Belham Valley, Richmond Hill, and Montserrat-Havers fault systems (extension and left-lateral slip) (Baird et al., 2015). It is conceivable that magma migration through the middle and upper crust utilizes structural weaknesses in the form of deep-seated fault zones. We find a close alignment of the modeled ellipsoid's a axis (strike of 123°) with the NW-SE strike of the Centre Hills Fault. Thẽ NW-SE direction agrees with the polarization orientation of the fast shear wave at SHV (Baird et al., 2015) and a stress reorientation of the background~E-W maximum horizontal stress along the left-lateral extensional Montserrat-Havers fault system to a NW-SE direction (Baird et al., 2015). The ellipsoid's b axis is oriented along the least compressive horizontal stress and may explain the slight elongation of the pressure source in the NE-SW direction (Figure 3)

Geophysical Research Letters
Geodetic models using volumetric strain data propose a dike-conduit system (Odbert, Ryan, et al., 2014) that connects the deeper part (considered here) of the magmatic plumbing system with the surface. The orientation of the dike is still under discussion, and while both NNW-SSE to NW-SE (Baird et al., 2015;Hautmann et al., 2010b;Linde et al., 2010) and NE-SW (Roman et al., 2011) alignments have been proposed, the orientation of the modeled ellipsoid supports an~NW-SE alignment with an opening along the NE-SW direction. In summary, we propose a magmatic plumbing system for SHV whose geometry and orientation are controlled by the local stress regime and consists of an inclined and vertically extended transcrustal magma reservoir that is connected to a shallow-seated dike-conduit (Figure 3).
The effects of deposit loading between 1995 and 2010 on surface deformation (Odbert et al., 2015) at SHV may explain lower than expected observed vertical displacements at GPS Stations HARR, SSOU, and WTYD that are located closest to the active volcanic center (Figures 2c and 2d and S2). Up to 0.6 m of subsidence are predicted over the 15 years of eruptive activity and volcanic sedimentation (Odbert et al., 2015) in the southern part of Montserrat. The deposition of volcanoclastics is heavily controlled by surface topography and geomorphological features whereby the highest accumulation of up to 381 m thickness is reported within English's crater (Odbert et al., 2015) (Figure 2) and a mean annual vertical accumulation rate of 0.4 m between 2002 and 2003 for the middle and lower Belham valley (Barclay et al., 2007). Surface loading contributes on average about −0.04 m/year in vertical displacement in areas proximal to the SHV dome and 0.01 to 0.02 m/year in surrounding areas (Odbert et al., 2015) and can explain the residuals between observed and modeled vertical displacements at Stations HARR, SSOU and WTYD as well as in the lahar-filled Belham valley (GPS Station OLVN; see Barclay et al., 2007; Figures 2c and S2).  Figure 1, and the center of the ellipsoid is located at 9.35 km below mean sea level. The hatched area marks an area which may be susceptible to flank instability due to the largest amplitude of surface deformation by the pressurization of the modeled transcrustal magma plumbing system. (b) Plan view sketch of the volcano-tectonic relationship of southern Montserrat. Stress rotation of a background~E-W directed maximum horizontal stress (S H ; dark green arrows) by a left-lateral and extensional (red arrows indicate kinematics) set of faults (abbreviations defined in Figure 1) results in a NW-SE directed S H at SHV (Baird et al., 2015). Semiaxis a of the best fit triaxial ellipsoidal source matches the strikes of S H and of the Centre Hill Fault (CHF). Semiaxis b is aligned with the direction of the minimum horizontal stress. (c) N-S cross-sectional sketch of the conceptual model of the transcrustal plumbing system of SHV. The main reservoir as modeled by the best fit pressure source is connected with the eruption center at SHV by a NNW-SSE to NW-SE striking dike-conduit system (Hautmann et al., 2009). 10.1029/2020GL089239

Geophysical Research Letters
Soufrière Hills and South Soufrière Hills volcanic complexes show abundant evidences for flank collapses and debris avalanche deposition along their southwestern and eastern shorelines; episodes of mass wasting correlate with episodes of heightened eruptive activity on Montserrat (Coussens et al., 2016). On shore, English's Crater is currently the most pronounced collapse feature at SHV (Figure 1), while submarine surveys highlight flank collapse structures from both onshore and submarine mass movements (Coussens et al., 2016;Feuillet et al., 2010Feuillet et al., , 2011Le Friant et al., 2010) in adjacent submarine sections to the east of SHV. Peak amplitudes of total ground displacement from the pressurization of the transcrustal plumbing system are predicted for the eastern and southeastern slopes of SHV including the offshore areas to the east and southeast of SHV (Figure 3). It is hence conceivable that flank collapses at SHV, including the formation of English's crater and the 1997 Boxing Day collapse (marked EC and BD, respectively, in Figure 1), have been influenced by subsurface and surface strain from the transcrustal magmatic pressurization. These findings may have implications for hazard assessment in relation to current and future flank instability on Montserrat and possible crossover to other volcanic islands of the Lesser Antilles (Deplus et al., 2001).

Conclusions
Our 3-D modeling reveals that surface displacement on Montserrat from magmatic forcing is profoundly controlled by (i) the heterogenous mechanical conditions of the middle and upper crust of Montserrat and (ii) surface topography. We image a southeastward dipping triaxial ellipsoidal pressure source embedded in mechanically weak crustal rocks to explain ground displacements during a period of intraeruptive unrest (2003)(2004)(2005) at SHV. The semiminor axes of the pressure source are aligned with the maximum and minimum horizontal stress directions (NW-SE and NE-SW), respectively, and indicate a control of the local stress conditions on the geometry and orientation of the magmatic plumbing system as part of an intersection between two major NW-SE and WNW-ESE striking tectonic lineaments. We neglect inelastic crustal rheology in the current study, which we expect to reduce required pressure change amplitudes while maintaining the predicted surface strains and introduce a temporal component to the modeled deformation. It remains to be seen how stress changes in the proposed transcrustal plumbing system can also explain other periods of intraeruptive ground deformation at SHV (Odbert, Ryan, et al., 2014). Our study demonstrates the significant influence of 3-D topography and heterogenous crustal mechanics on the inversion of geodetic data. We recommend incorporating multiparametric geophysical data in advanced modeling of volcano deformation studies to illuminate subsurface processes behind volcano deformation and to improve on results from traditional analytical deformation models.

Data Availability Statement
GPS data are available from the UNAVCO archive (https://www.unavco.org/data/data.html).