Topographical Effects on Volcano Deformation Signal Intensity: Implications for GPS Network Configuration

Volcano GPS networks can capture vital information during volcanic unrest to aid with hazard assessment and eruption forecasting, but can be hindered by their discrete point locations and possibly miss key spatial information. We show how numerical models can reveal controls on spatial deformation signal intensity compared against GPS network design. Using the GPS network at Soufrière Hills Volcano (SHV), Montserrat, and a range of models, we explore expected surface deformation patterns. Peak horizontal deformation is located offshore, highlighting the difficulties with geodetic monitoring on small ocean‐island volcanoes. Onshore areas where the deformation signal is expected to be high are also identified. At SHV, topography plays a greater role in altering the relative distribution of surface displacement patterns than subsurface heterogeneity. Our method, which can be adapted for other volcanoes, highlights spatial areas that can be targeted for effective GPS station placement to help improve deformation monitoring efficiency.


Introduction
Volcano deformation is a key observable during volcanic unrest and can ultimately help track the subsurface migration of magma and forecast its possible eruption (Biggs & Pritchard, 2017;Fernández et al., 2017;Hickey et al., 2017).Deformation of a volcano typically results from the accumulation of magma (e.g., Le Mevel et al., 2016), but can also be caused by magma migration (Sigmundsson et al., 2015) or thermodynamic processes within a magma reservoir (e.g., Caricchi et al., 2014).Non-magmatic mechanisms can also cause a volcano to deform (e.g., Rouwet et al., 2014;Schaefer et al., 2015).
A variety of geodetic methods can be used to record the spatial and temporal evolution of volcano deformation (Dzurisin, 2006).One of the most effective and widely used methods for monitoring volcano deformation is through Global Navigation Satellite System (GNSS) technology, particularly using the most developed and widely-recognised Global Positioning System (GPS) constellation (Janssen, 2007).We therefore use the term GPS throughout our work, but note that GNSS could be equally applied.GPS monitoring involves deploying multiple GPS receivers around the volcano to form a network (Figure 1).GPS stations can be established to record continuously, as a cGPS station, which excel at recording the temporal patterns of deformation, or may only be occupied episodically, as an eGPS station (Dzurisin, 2006).Both cGPS and eGPS can measure the three components of surface displacement (east-west, north-south, vertical), typically within <1 cm accuracy (Dzurisin, 2006).However, GPS stations can only record the deformation from their discrete point location within a network.Consequently, depending on their design and density, GPS networks could be missing insightful spatial information from a surface deformation pattern (e.g., Hickey et al., 2020).
Complexities in volcano deformation analysis can arise from the influence of topography, which can modulate the surface deformation pattern from a subsurface causative source (Cayol & Cornet, 1998;Williams & Wadge, 1998, 2000), or lead to inaccurate results if neglected when interpreting geodetic data (Hickey et al., 2016;Johnson et al., 2019;Ronchin et al., 2015;Tammaro et al., 2021).In particular, abrupt topographic changes (e.g., caldera walls, steep cliffs, landslide scars) can have drastic impacts on surface deformation patterns, even reversing local deformation gradients (Johnson et al., 2019;Ronchin et al., 2015).If GPS stations are located in, or close to, these areas they may be collecting inefficient data that is difficult to interpret or misleading for volcano deformation monitoring.On the contrary, a thorough understanding of the influence of topography on an expected surface deformation pattern could be exploited to improve both interpretation of geodetic data and the location of GPS sites within an existing network.Similarly, but to a lesser extent given our inability to accurately classify it, an understanding of the impact of subsurface mechanical heterogeneity could also help improve GPS network configuration.Subsurface mechanical heterogeneity strongly influences the partitioning of stress and strain within the crust in response to a deformation-inducing process, with cumulative effects for surface deformation (Gudmundsson, 2012;Hickey et al., 2013;Manconi et al., 2007;Masterlark, 2007), and significant impacts on inferred results when simulating geodetic data (Hickey et al., 2016).
The demonstrated impact of topography and subsurface mechanical heterogeneity in influencing surface deformation patterns, and the potential limitations of discrete GPS station locations, set the context for our study.We test whether we can use numerical deformation models to identify regions that have higher potential surface deformation signals, and explore the key controls on the spatial distribution of signal magnitude.Ultimately, this work aims to contribute to improved understanding and interpretation of volcano deformation with respect to GPS network configuration.
We use Soufrière Hills Volcano (SHV) as the case study for our analyses.SHV is an andesitic dome-building volcano in Montserrat (Figure 1) currently in an intra-eruptive unrest period (Alshembari et al., 2024;Hickey et al., 2022) within a prolonged eruptive episode that began in 1995 (Wadge et al., 2014).Deformation monitoring has been well-established on the island since 1995 and managed by the Montserrat Volcano Observatory (MVO), using tilt, strain, EDM and GPS (Odbert et al., 2014).In the tropical environment, InSAR has not been used reliably due to loss of coherence on steep and heavily vegetated slopes, a relatively small deformation signal, and substantial atmospheric errors (e.g., Camejo-Harry et al., 2023;Wadge et al., 2006).Current deformation measurements indicate the volcano has been inflating since the end of the latest eruption in 2010 (Figure 1 & Figure S1 in Supporting Information S1), and is believed to be the result of magma pressurization in a reservoir centered at ∼9.5 km below sea level (BSL) (Alshembari et al., 2024;Hickey et al., 2022;Johnson et al., 2023;Neuberg et al., 2022).SHV is a particularly good candidate for our study as the topography (Stinton, 2015) and subsurface heterogeneity (Paulatto et al., 2012(Paulatto et al., , 2019) ) are suitably well characterized, but the small island geography and harsh terrain limit the spatial distribution of GPS stations.GPS network design at SHV has mostly been based on distance from the vent, azimuthal coverage, and ease of access (Odbert et al., 2014).Additionally, assuming a typical volcano deformation pattern where the greatest vertical displacement is located closest to the surface projection of the deformation source (e.g., Mogi, 1958), the vertical measurements at SHV show lower magnitudes closest to the source than might be expected (Figure 1d).

Model Setups
Our numerical models are constructed and solved using Finite Element Analysis (FEA) with COMSOL Multiphysics v5.6 (Figure 2).We use a three-dimensional (3D) model setup to incorporate 10 m resolution topography data (Stinton, 2015) and 50 m resolution bathymetry data (Le Friant et al., 2004).The topography and bathymetry define the free surface of the model, the base of the model has a zero-displacement boundary condition and the lateral edges have a roller condition (Hickey et al., 2022;Hickey & Gottsmann, 2014).Surrounding the model geometry is an infinite element domain to prevent boundary conditions impacting the results of the model interior; we do not load the model gravitationally.The deformation source is a pressurized prolate ellipsoidal cavity, initially located offset from the vent (Figure 1) at a central depth of 9.5 km BSL, with a semi-major axis of 4 km and semi-minor axes of 2 km (Gottsmann et al., 2020;Johnson et al., 2023).Magma reservoir overpressure is simulated in the deformation source by applying a uniform pressure boundary load of 10 MPa.The region surrounding the deformation source and the surface of the model have the greatest mesh density, with a typical model consisting of ∼800,000 mesh elements.
In our numerical models employing a 3D heterogeneous subsurface, mechanical properties are constrained with a 3D P-wave (v P ) tomography model (Paulatto et al., 2012).Dynamic Young's Modulus (E D ) and density (ρ) are calculated from P-wave velocities (Brocher, 2005): where v is the Poisson's Ratio of 0.25.The dynamic Young's Modulus is scaled to a static Young Modulus (E) where E = E D /2, producing values of ∼1-100 GPa.E D is typically 2-13 times larger than E (Cheng & Johnston, 1981), and E is more appropriate for volcanic deformation timescales (Gudmundsson, 2011;Heap et al., 2020).In our models with a homogeneous subsurface a uniform value of E = 25 GPa is applied.As we are only interested in exploring the spatial aspects of the deformation field we only consider an elastic crustal rheology.
We also use an analytical model to calculate surface deformation for the case of a homogeneous subsurface with a flat, topography-free, half-space (Figure 2d).Our analytical models are setup with the same source parameters and calculated using the solutions of Yang et al. (1988) within the dMODELS package (Battaglia et al., 2013).

Experimental Approach
We model surface deformation using an array of forward models, encompassing uncertainties on all deformation source characteristics.Our simulations are conducted using four different model setups:  S1 in Supporting Information S1).For each model setup, and every solution from the source grid, the surface deformation calculations are summed and normalized to identify areas of greatest displacement in any given deformation component.
Uncertainties on the deformation source axes lengths are not considered as they are insignificantly small compared to the uncertainties on the source location and can, for the current study, be considered to be incorporated into the source location uncertainty within the grid.As we sum and normalize our solutions to produce deformation intensity 'heat-maps', the absolute value of source overpressure is irrelevant so no variation is applied to this model parameter.

GPS Data
To validate our models we use data for the 2010-2022 period from the SHV GPS network (Figure 1).GPS data are processed using GAMIT/GLOBK (Herring et al., 2010), and station coordinates are calculated in a Caribbean reference frame (Odbert et al., 2014) to remove the regional tectonic component from the displacements.The GPS measurements reveal a pattern of long-term inflation at the volcano over the 12-year period studied (Figure 1).

Vertical Deformation
The influence of topography is immediately apparent when inspecting model results (Figure 3 and Figure S2 in Supporting Information S1).Peak vertical deformation is shifted east of the source grid in models incorporating the topography of the island, caused by the morphology of the volcanic edifice.It is easier to vertically deform the eastern flank of the volcano since it has a lower relief and is therefore closer to the deformation source(s).In contrast, in the topography-free models peak vertical deformation is concentrated around the source grid.For the homogeneous topography-free model, peak vertical deformation around the source grids is circular and followed by concentric circles of decreasing vertical deformation intensity (Figure 3).However, for the heterogeneous topography-free model, peak vertical deformation around the source grid is asymmetrical (Figure S2 in Supporting Information S1), reflecting an additional but lesser role of crustal Young's Modulus in modulating surface uplift patterns.
The homogeneous topography-free model predicts that 10 GPS stations should be within a circular region of the greatest vertical surface deformation, centered on the source grid (Figure 3).This pattern is not reflected in the GPS measurements, which show a number of stations with lower (<12 cm) deformation within a 5 km radius of the source grid center (Figure 1).In contrast, the full complexity FEA model suggests there would only be two GPS stations within the region of peak vertical deformation intensity (Figure 3a; RCHY and LGRD).However, these two stations do not have the greatest recorded vertical deformation magnitudes in our dataset (Figure 1).A possible reason for this mismatch is that our models simulate an inferred major mid-crustal deformation source but do not account for any potential shallower deformation sources located beneath the volcanic edifice, such as hydrothermal processes and/or small magmatic intrusions (Johnson et al., 2023), or for deposit loading (Odbert et al., 2015), all of which can impact vertical deformation magnitudes but were beyond the scope of this study.Deposit loading and cooling and contraction of shallow magmatic intrusions would cause subsidence with small surface footprints (Johnson et al., 2023;Odbert et al., 2015), and may be contributing to the lower-than-expected recorded vertical inflation close to the edifice.
All model setups indicate that vertical deformation should be recorded across the entirety of the island, which is reflected in the GPS records.A broad deformation field is expected for a source ∼10 km below the surface, but our models indicate where surface deformation produced by that source should be greatest, and that it is most strongly influenced by topography.

Horizontal Deformation
In all model setups, peak horizontal deformation mostly occurs offshore (Figure 3 and Figure S2 in Supporting Information S1), highlighting the problems with monitoring deformation on a small ocean-island volcano.However, when inspecting the full-complexity FEA model (Figures 3a-3e), which we consider to be the closest to natural conditions, some areas of high-magnitude northward and westward deformation do occur on the island.Existing GPS sites TRNT, MVO1, and GERD are suggested to capture some of the highest magnitude northward deformation, while BNBY, OLVN, and AIRS are suggested to catch some of the highest magnitude westward deformation.Examining the deformation measurement magnitudes of these sites shows they record the highest values (Figures 1b and 1c).The close match between deformation observations and the predicted deformation intensities of our models (Figure 3) provides qualitative verification for our approach.
As the majority of horizontal deformation occurs offshore, bathymetry influences the shape of the deformation footprints, along with topography.For example, southward deformation intensity bands bend around the south side of the island, while the northward deformation has a highly offset peak intensity band that is also warped to the shape of the island (Figure 3).In both cases this is likely the combined result of the near-offshore seabed surface being closer to the deformation source(s), as well as the morphology of the topography and bathymetry affecting surface strain partitioning.The effect is more apparent off the south of the island as the insular shelf is shorter, the bathymetric slope is steeper (Le Friant et al., 2004), and the deformation source(s) are in closer proximity.
The analytical model produces deformation intensity predictions that are broadly similar to the full-complexity FEA model; the majority of the peak surface displacements are still offshore.However, the analytical model results do not show the influence of topography in modulating the deformation signal.Topography has subtle but important ramifications for the predicted surface deformation patterns, altering the on-island regions of higher deformation intensity, particularly for the north and west displacement components.For example, in the north component the analytical models suggest that GERD should have the largest deformation magnitude, but TRNT and MVO1 have the greatest magnitudes of north displacement in the GPS data, matching the predictions from the full-complexity FEA model.
Similar to the vertical deformation, subsurface heterogeneity has a lesser impact on the predicted horizontal surface deformation patterns than topography.For northward deformation, from a homogeneous and topography-free starting point, heterogeneity pushes the peak deformation region to the south while topography pushes it offshore.However, topography has a more significant impact on the mid-range intensities, with deformation intensity bands being shaped around the islands relief.Together, our results suggest that topography has greater control than heterogeneity on influencing the spatial distribution of surface deformation signals at SHV.

Total Deformation
The normalized combination of all surface deformation components highlights on-island areas where the total deformation signal should be high (Figure 4).Focusing on the full-complexity numerical model indicates that several GPS sites in the central to northern part of the island are already wellpositioned to receive a theoretically high signal from the modeled deformation source.However, as seen for the horizontal components in isolation, the peak signal is offshore.There remain some areas that could be targeted for future GPS site development if based on predicted deformation intensities alone, for example, in the central part of the island and further toward the southern shoreline.However, these suggestions should be considered only first-order, as other aspects need to be considered when siting a new GPS station for an island like Montserrat, particularly site access, ground conditions, radio links for data transmission, and vegetation.
The analytical model predictions have a roughly similar shape to the numerical models, but have four regions of equal peak surface deformation intensity (Figure 4 & Figure S3 in Supporting Information S1).Consequently, this would suggest that four regions of the island should record large magnitude deformation, but this is not reflected in the recorded data (Figure 1).Instead, the numerical model with subsurface heterogeneity and  surface topography predicts one major region of surface deformation on the island and the GPS sites contained within it are often those recording the highest magnitudes of displacement (Odbert et al., 2014; Figure 1).

Implications for Volcano Deformation Monitoring
Our approach has identified regions most likely to experience higher surface deformation signals from a given subsurface deformation source.By using SHV as our case study we were able to limit and constrain key modeling parameters for a well-studied magmatic system.Even when significantly increasing uncertainty estimates on model source parameters, and exploring larger ranges of source location such as may be the case in lesser-studied volcanoes, the major patterns identified remain the same (e.g., Figure S4 in Supporting Information S1).For SHV, the key control on the spatial distribution of the signal magnitude appears to be topography (e.g., Neuberg et al., 2022), with a lesser role played by heterogeneity in subsurface elastic properties.Our results imply the role of topography should be strongly considered when analyzing volcano deformation datasets.
The approach we have presented could be applied to other volcanoes where there is prior information on likely deformation source parameters, and as topographical datasets are globally available this key model component could be implemented to test its significance in different settings.Information about crustal heterogeneity is likely only available at volcanoes with substantial research and/or monitoring resources, but nonetheless experience here for SHV shows it plays a lesser role than topography.Using an analytical model approach will capture broadly similar patterns to numerical models, and could possibly be improved by incorporating first-order analytical topographical corrections (Williams & Wadge, 1998, 2000).Such analytical techniques may be the only feasible approach in many resource-limited environments where numerical modeling software is unaffordable and/or outside the expertise of research staff.However, the full role of topography in modifying surface deformation patterns, which is not obtainable through analytical modeling, should not be overlooked; it does produce modifications that change deformation predictions and would impact the use of the results in understanding and interpreting volcano deformation with respect to GPS network configuration (e.g., knowing what signal(s) may be missed by the network design).
For potential future design modifications of volcano GPS networks our method would allow volcano geodesists to focus on various deformation components.The "best" choice is likely to be volcano and resource dependent, but given horizontal GPS observations are typically accompanied by lower uncertainty than vertical measurements they could favor higher in any decision making.Consideration still needs to be given to network design aspects like access, power supply, communication, and cost, but our approach provides a novel technique to potentially enhance strategic station placement with regards to recording a higher relative deformation signal, and could help improve deformation monitoring efficiency.Furthermore, enhanced measurement of volcano deformation signals will also lead to improvements in modeling and analyzing the observations, producing more robust inferences about the underlying causative deformation source(s).Overall, we show that numerical models can be applied to improve understanding of topographical effects on volcano deformation signal intensity and potentially enhance GPS network configurations to maximize the "value" of a GPS site from a monitoring perspective.

Conclusions
Our approach used a range of complexity of models to explore cumulative surface deformation patterns expected from a subsurface distribution of literature-derived deformation source configurations.We focused on the GPS network at SHV, Montserrat, to demonstrate how sites within the network may or may not be ideally situated to profit from an expected deformation signal intensity.Peak vertical deformation is shifted east of the deformation source due to the impact of edifice topography.Models neglecting topography show peak vertical deformation located around the edifice and directly above the deformation source grid, but these patterns are not reflected in the measured GPS data.Predicted peak horizontal deformation intensities are located offshore for all model setups, but some regions of high-intensity north and west deformation signal are expected onisland and correlate with recorded GPS magnitudes.By comparing the different model setups, topography appears to have a greater influence on impacting the spatial patterns of the deformation signals than variability in subsurface crustal mechanics at SHV.For the offshore regions, bathymetry has a similar effect to topography.Generally, analytical models predict broadly similar patterns to the numerical models but inherently cannot capture the localized modifications caused by topography which would be important when using model results to assess volcano deformation with respect to GPS network configuration.The modeling approach we present allows for specific spatial areas to be targeted for GPS site development depending on their predicted deformation intensity, benefitting volcano monitoring strategies.Our techniques can be adopted for application at other volcanoes to potentially help strategically redesign volcano GPS networks to target a higher relative deformation signal, with knock-on effects for improved deformation monitoring efficiency in an often resourcelimited environment.

Figure 1 .
Figure 1.The 2010-2022 GPS network and surface displacement records at Soufrière Hills Volcano.(a) Spatial distribution of cGPS (red) and eGPS (blue) stations on Montserrat.The pink triangle is the vent, and the purple circle is the inferred magma reservoir deformation source location (Johnson et al., 2023).(b) Observed east-west deformation.(c) Observed north-south deformation.(d) & (e) Observed vertical deformation, relative to the distance from the vent (d) or deformation source (e).
(1) heterogeneous subsurface and surface topography; (2) homogeneous subsurface and surface topography; (3) heterogeneous subsurface and flat surface; and (4) homogeneous subsurface and flat surface (Figure 2).Model setups 1-3 are solved numerically using COMSOL Multiphysics, while model setup 4 is solved analytically.Within each model setup the deformation source location is varied within a 3D grid consisting of 125 points.The initial source location, derived from Johnson et al. (2023), is used as the center of the grid.The solution of Johnson et al. (2023) is only marginally different to that of Gottsmann et al. (2020), but as Johnson et al. (2023) do not provide uncertainty estimates for their deformation source parameters we use the uncertainties from Gottsmann et al. (2020), multiplied by two, to define the upper and lower limits of our parameter grid (Table

Figure 2 .
Figure 2. Model setups and boundary conditions.(a) FEA model with 3D heterogeneous subsurface and surface topography.(b) FEA model with 3D heterogeneous subsurface and no topography.(c) FEA model with homogeneous subsurface and surface topography.(d) Analytically-solved model with homogeneous subsurface and no topography.
Figure 3. Figure 3. Modeled deformation intensity heatmaps.(a-e) Results from FEA models incorporating 3D subsurface heterogeneity and surface topography for different surface displacement components.(f-j) Same as a-e but for a model setup with no topography and a homogeneous subsurface.Warmer colors indicate a greater likelihood of deformation.The pink triangle is the vent, and the purple circles are the surface projections of the deformation source locations.cGPS stations are indicated by black circles; eGPS stations are shown by white squares.

Figure 3 .
Figure 3. Figure 3. Modeled deformation intensity heatmaps.(a-e) Results from FEA models incorporating 3D subsurface heterogeneity and surface topography for different surface displacement components.(f-j) Same as a-e but for a model setup with no topography and a homogeneous subsurface.Warmer colors indicate a greater likelihood of deformation.The pink triangle is the vent, and the purple circles are the surface projections of the deformation source locations.cGPS stations are indicated by black circles; eGPS stations are shown by white squares.

Figure 4 .
Figure 4. Combined total surface deformation intensity heatmaps.(a) Results from FEA models incorporating 3D subsurface heterogeneity and surface topography for the combined results of all surface displacement components.(b) Same as (a) but for a model setup with no topography and a homogeneous subsurface.Warmer colors indicate a greater likelihood of deformation.The pink triangle is the vent, and the purple circles are the surface projections of the deformation source locations.cGPS stations are indicated by black circles; eGPS stations are shown by white squares.