Rapid localized flank inflation and implications for potential slope instability

High rates of volcano surface deformation can be indicative of a forthcoming eruption, but can also relate to slope instability and possible ﬂank collapse. Tungurahua volcano, Ecuador, has been persistently active since 1999 and has previously experienced catastrophic ﬂank failures. During the ongoing eruptive activity, signiﬁcant surface deformation has been observed, with the highest rates contained within the amphitheatre-shaped scar from the 3000-year-old failure on the west ﬂank. However, the cause of this asymmetric deformation and how it might relate to slope stability has not been assessed. Here, for the ﬁrst time, we present a range of models to test physical processes that might produce asymmetric deformation, which are then applied to slope stability. Our models are informed by InSAR measurements of a deformation episode in November 2015, which show a maximum displacement of ∼ 3.5 cm over a period of ∼ 3 weeks, during which time the volcano also experienced multiple explosions and heightened seismicity. Asymmetric ﬂank material properties, from the rebuilding of the cone, cannot explain the full magnitude and spatial footprint of the observed west ﬂank deformation. The inﬂation is inferred to be primarily caused by shallow, short-term, pre-eruptive magma storage that preferentially exploits the 3 ka ﬂank collapse surface. Shallow and rapid pressurization from this inclined deformation source can generate shear stress along the collapse surface, which increases with greater volumes of magma. This may contribute to slope instability during future unrest episodes and promote ﬂank failure, with general application to other volcanoes worldwide displaying asymmetric deformation patterns.


Introduction
Magma transport and storage have important implications for volcanic hazards , particularly for volcanoes with a history of flank collapse, where pre-existing structural discontinuities may contribute to inherent instability and act as preferential pathways for magma/fluid migration (Donnadieu et al., 2001;Schaefer et al., 2013;Tibaldi, 1996). Constraints on magma storage and transport can be obtained from analysis of deformation data during volcanic unrest episodes (e.g., Hickey et al., 2017;Lloyd et al., 2018;Sigmundsson et al., 2014), and geodetic data can also be used to monitor slope stability and iden-tify conditions which may promote flank failure (e.g., Schaefer et al., 2017;Tridon et al., 2016). Analyses of shallow magma transport and storage are shown to promote flank instability and cause localised surface deformation (Donnadieu et al., 2001;Schaefer et al., 2017;Siebert, 1984). We use Interferometric Synthetic Aperture Radar (InSAR) deformation data to investigate a rapid and localised, co-eruptive inflation episode at Tungurahua volcano, Ecuador, in November 2015. Tungurahua (5023 m asl, above sea level) is a steep stratovolcano in the Eastern Cordillera of the Ecuadorian Andes ( Fig. 1) and has been frequently active since 1999, displaying Strombolian to sub-Plinain explosive activity interspersed with periods of quiescence (Hidalgo et al., 2015;Mothes et al., 2015). The present edifice of Tungurahua (Tungurahua III), with over 3000 m of relief and slope angles greater than 30 • , is built upon the remnants of two former cones (Tungurahua I and II) that each subsequently collapsed (Hall et al., 1999). The most recent collapse occurred 3000 years ago and directed an 8 km 3 (Hall et al., 1999), and the blue diamond is the summit vent. Coloured, transparent shapes indicate the footprints of deformation episodes: -2008(green, Fournier et al., 2010-2008(red, Biggs et al., 2010-2011(pink, Morales-Rivera et al., 2016(yellow, Muller et al., 2018 and November 2015 (blue, this study). GPS and tilt stations are shown by white triangles, and red squares, respectively. The dashed black line is the projection of the cross-sections shown in west that travelled 15 km (into drainages in the north-west and south-west) and covered 80 km 2 , leaving a 110 • wide amphitheatre collapse scar (Hall et al., 1999). A lahar immediately followed the debris avalanche, likely due to dewatering of the avalanche deposits or melting of small summit glaciers. Scoria bombs with cooling fractures in the lahar deposits imply the initial collapse was triggered, or accompanied, by an eruption (Hall et al., 1999). Subsequent dacitic lava flows extruded from within the collapse scar may further allude to the intrusion of viscous magma, possibly as a cryptodome, as the trigger for the collapse event (Hall et al., 1999). Continued eruptive activity has since rebuilt the cone, and the amphitheatre has been backfilled with 3 km 3 of unconsolidated material (Molina et al., 2005), overlying the 3 ka unconsolidated avalanche debris.
Geophysical and geochemical monitoring of Tungurahua is conducted by the Instituto Geofísico de la Escuela Politécnica Nacional (IG-EPN), who maintain seismic, infrasound, GPS and tiltmeter networks, as well as EDM (electronic distance measurement) lines and gas flux observations with DOAS (differential optical absorption spectroscopy) and FTIR (Fourier transform infrared spectrometry) (Alvarado et al., 2018). SO 2 measurements generally correlate with eruptive activity; 95% of emissions in the period 2007 -2013 occurred during active phases, with the remaining 5% during quiescent periods (Hidalgo et al., 2015). Seismicity at Tungurahua is mostly dominated by long-period (LP) events (Bell et al., 2018;Palacios et al., 2015), which have been proposed to occur due to coupled magma and gas ascent (Bell et al., 2017), resonance of an ash-laden gas in the conduit (Molina et al., 2004), or variations in magma ascent rate and conduit stick-slip behaviour (Neuberg et al., 2018). Volcano-tectonic (VT) seismicity is rare, with small events sometimes occurring prior to increased rates of LP's before the start of new explosive activity (Bell et al., 2017), and most often detected beneath the west flank (Palacios et al., 2015).
Surface deformation at Tungurahua is observed over a range of length and timescales (Biggs et al., 2010;Champenois et al., 2014;Fournier et al., 2010;Morales-Rivera et al., 2016;Muller et al., 2018, Table 1). At the smallest spatial scale (< 2 km), high-amplitude near-vent tilt cycles have been detected and are thought to be caused by shear stress in the conduit due to viscous magma flow resistance (Marsden et al., 2019;Neuberg et al., 2018). Separate large-scale deformation has been detected using GPS and InSAR (Champenois et al., 2014;Muller et al., 2018). Between 2003and 2009, Champenois et al. (2014 observe a longterm 50 km wide deformation footprint with Envisat InSAR data; their preferred model for this period of uplift is a spherical source at a depth of 11.5 km bsl (below sea level) with a volume increase rate of 7 × 10 6 m 3 /yr. Using combined GPS and TerraSar-X InSAR data for 2011, Muller et al. (2018 observe edificewide uplift and model the deformation using a slightly dipping prolate shaped source at a depth of 7.4 km below average elevation (∼5 km bsl) with a volume increase rate of 2 × 10 6 m 3 /yr.
Interestingly, at the small-medium spatial scale, surface deformation at Tungurahua is concentrated within the 3 ka flank collapse scar ( Fig. 1), and between the continuous GPS network . Numerous InSAR studies have detected this uplift pattern over timescales of months to years (Table 1), which has been modelled and interpreted as shallow magma intrusion or storage beneath the west flank, often preceding or accompanying an eruption (Biggs et al., 2010;Fournier et al., 2010;Morales-Rivera et al., 2016;Muller et al., 2018). It is unknown how the inferred presence of shallow magma may relate to the stability of the west flank, especially since this is also the site of unconsolidated material and high-frequency seismicity, and magmatic activity is inferred to have contributed to a previous collapse.
In this study, we similarly detect a period of uplift on the west flank of Tungurahua, coincident with an eruptive phase, but provide significantly better time constraints than previously available. We present new models of the deformation using Finite Element Analysis (FEA) to test different physical mechanisms which may promote asymmetric (west-flank) displacement. Model results are used to constrain a source model and provide the first assessments of the implications of shallow pressurization on the evolution of slope stability.

Surface deformation
We use SAR (Synthetic Aperture Radar) data from two satellite systems, Sentinel-1 (S1) and COSMO-SkyMed (CSK), employing both ascending and descending viewing angles from April 2015 to February 2016 to produce 47 interferograms. S1 TOPS (Terrain Observation by Progressive Scans) data (ascending path 18 and descending path 142) were processed using GAMMA (Werner et al., 2000). Topographic phase contributions were removed using 1-arc-second SRTM (Shuttle Radar Topography Mission) digital elevation model (DEM) data (30 m resolution), we filtered the interferograms using an adaptive filter with strength 0.4, and unwrapped using a minimum cost flow (MCF) algorithm with pixels of coherence less than 0.35 masked out.
CSK spotlight data (ascending and descending) were processed using ISCE (InSAR Scientific Computing Environment) (Rosen et al., 2012). Topographic phase contributions were removed using 30 m SRTM DEM data, we filtered the interferograms using a powerspectral filter with strength 0.4, and unwrapped using a MCF algorithm.
Interferograms were corrected for atmospheric effects using TRAIN (Toolbox for Reducing Atmospheric InSAR Noise) (Bekaert et al., 2015a). We use a phase-based linear tropospheric correction due to the unavailability of reliable atmospheric data that would enable more complex atmospheric corrections (e.g., Bekaert Table 1 Summary of west flank deformation observations at Tungurahua volcano using InSAR. asl = above sea level; LOS = line-of-sight.

Study
Observation period Uplift information Source model Interpretation Biggs et al., 2010Dec 2007-Mar 2008 Figure 1), a common issue for InSAR processing in this region (Ebmeier et al., 2013;Morales-Rivera et al., 2016). Inspection of the interferograms highlights a period of deformation in November 2015. Interferograms preceding and following this period show no deformation. Owing to the short duration of deformation, it is only captured in 8 individual interferograms (one CSK descending, one CSK ascending, two S1 ascending, and four S1 descending). Accounting for overlap in the acquisition dates across the two sensors and both viewing angles, the deformation period can be constrained to between the 3rd and 27th of November 2015 and independently viewed in 4 interferograms (one for each sensor and viewing geometry, Fig. 2). We are able to provide significantly improved temporal constraints on the deformation period compared to earlier studies due to the use of new satellites, indicating that the deformation is much faster than previously considered (e.g., Table 1). The data show a region of uplift on the west flank with maximum amplitudes in the satellite line-of-sight (LOS) of ∼3-4 cm, depending on viewing geometry and sensor.
Owing to the near-contemporaneous acquisition of both ascending and descending data for both CSK and S1 sensors, the LOS deformation can be decomposed into east-west and up-down (vertical) components (Fig. 2), using the assumption the north-south displacement is negligible due to the poor satellite sensitivity in this direction (e.g., Lloyd et al., 2018;Wright et al., 2004). Decomposition indicates that the majority of the observed uplift can be explained by vertical movement, with negligible contributions in the east-west direction. A similar observation was found for the deformation period between December 2007 and March 2008 (Biggs et al., 2010). For the current period of study (November, 2015), the deformation footprint is elongated in the east-west di-rection to a maximum footprint of ∼5 km, and the corresponding maximum LOS uplift rate is equivalent to ∼48 cm/yr. Coincident with the observed deformation, the volcano was in an eruptive phase. Throughout November 2015 there were a number of explosions, with ash columns rising between 1.5 and 4.5 km above the vent, and several secondary lahars. During this period there was also a pronounced increase in VT seismicity (Supplementary Figure 2), but these events have not been spatially located. Activity since 2015 has declined with the last recorded eruptions in 2016, and only minor seismic and fumarolic behaviour recorded in 2017. No new deformation events on the west flank have been reported.

Numerical modelling
To model and analyse the observed surface deformation we use Finite Element Analysis (FEA) with COMSOL Multiphysics v5.3. We focus on two possible mechanisms to explain the asymmetric west-flank deformation: heterogeneous subsurface material properties with a central deformation (pressure) source, or heterogeneous subsurface material properties with a deformation source directly beneath the west flank (Fig. 3). The former is based on the evidence for weak unconsolidated material infilling the 3 ka flank collapse scar (Molina et al., 2005), which may serve to amplify deformation on that flank from a deeper central deformation source (Fig. 3a), as weaker materials are known to enhance surface deformation (e.g., Hickey et al., 2013;Masterlark, 2007). The latter is based on previous modelling results for similar deformation patterns at Tungurahua (Table 1), but here we also account for known subsurface heterogeneity and surface topography which have previously been neglected (Fig. 3b). Alternative mechanisms have also been suggested to produce asymmetric volcanic surface deformation from observations of large seaward displacements in  basaltic ocean island settings, based on slow slip along a shallow detachment fault (Tridon et al., 2016) or sheared sheet intrusions (Cayol et al., 2014). However, we can discount these mechanisms at Tungurahua due to the lack of significant observed horizontal displacement (Fig. 2).

Model setup
Model boundary and initial conditions are adapted to a threedimensional (3D) geometry from the benchmarked models of Hickey and Gottsmann (2014). The deformation source is represented as a tri-axial spheroidal cavity, the surface is generated from a SRTM 30 m DEM, and the west flank is given a separate model domain (Fig. 4). The 3D geometry of the west flank domain is constrained from the surface expression of the 3 ka collapse (Hall et al., 1999) and extended to a depth of 2 km asl, as . Finite Element Analysis model geometry and boundary conditions. A pressurised cavity is used to represent the deformation source and a bordering infinite element domain prevents boundary effects. Roller conditions are assigned to all lateral boundaries, and a fixed constraint is applied to the base of the model. The top surface is free to deform. The west flank has a separate model domain in order to assign different material properties. The geometry is meshed with ∼210,000 tetrahedral elements, with greater mesh density around the source and surface.
identified in seismic P-wave imaging (Molina et al., 2005). Heterogeneous elastic material properties (varying with depth) are calculated from seismic P-wave velocities (Molina et al., 2005) by converting to a dynamic Young's Modulus, and then scaling to a static Young's Modulus, using a conservative scaling factor of 2 (e.g., Gudmundsson, 2011;Hickey et al., 2017;. The resultant values are within the expected range for volcanic regions (Gudmundsson, 2011), with the main model domain increasing from ∼5 GPa at the surface to ∼42 GPa at 20 km bsl, and the  Fig. 1) from the heterogeneous and homogeneous models that cuts through the west flank, and normalise the modelled vertical uplift to the maximum value for each profile. We plot the difference between the normalised vertical deformation profiles, which shows the change in uplift shape between the heterogeneous and homogeneous models.  Figure 3). In all models, a homogeneous value of 0.25 is assumed for the Poisson's Ratio. A number of volcano deformation studies also include the effect of viscoelastic material parameters to account for elevated deep source temperatures and/or time-dependent displacement (e.g., Hickey et al., 2016), but we omit this complexity due to the lack of temporal resolution in our data and the high probability of a shallow source above the brittle-ductile transition (e.g., Table 1).

Modelling approach
To test if material properties have a strong influence on the observed west flank deformation (Fig. 3a), we run a series of forward models comparing heterogeneous and homogeneous subsurface conditions, using a roughly centrally-located buried pressure source. The source has three equal axes of 1500 m (i.e., spherical), a pressure of 1 MPa, varies in horizontal location (longitude and latitude) up to ±3 km from the vent, and has a central depth between 0 and 5 km bsl, which produces a 3D grid of 294 source locations. These parameters are chosen to be representative of any possible deep, central, deformation source. The predicted surface displacements in the heterogeneous and homogeneous models are compared for the 294 different source locations. The heterogeneous model accounts for depth-dependent variations in Young's Modulus in both the main model and west flank domains (Supplementary Figure 3). The homogeneous model has the same constant value of Young's modulus in both domains.
To investigate a possible deformation source beneath the westflank (Fig. 3b) we run a series of FEA inversions solving for the optimum source parameters in a heterogeneous model, adapting the approach of Hickey et al. (2016) to identify the cumulative global minimum misfit to all four interferograms simultaneously. The optimised source parameters are the three axes lengths of the tri-axial spheroid, the horizontal and vertical location, dip angle, and pressure. Strike is kept constant to reflect the dominant eastwest deformation footprint, and as informed by previous studies (Table 1) and sequentially reduced the size of the parameter search grid from one optimisation to the next using the solution of one inversion as the starting point of the next inversion. This process was stopped when the model residual could no longer be minimised and ensures a robust solution (Hickey et al., 2016). Alternative starting points were also tested, as well as a number of Monte Carlo simulations, but they did not produce better solutions.

Effect of west-flank heterogeneity
To assess the influence of west flank heterogeneity on producing asymmetric surface inflation at Tungurahua, we compare deformation models that have homogeneous or heterogeneous subsurface structures. We extract a north-south cross-section from the models that cuts through the west flank, and normalise the modelled vertical uplift to the maximum value for each profile (Supplementary Figure 4). By calculating the difference between the heterogeneous and homogeneous normalised profiles for each of the different source locations, it is clear that the greatest variability in the shape of the modelled uplift patterns occurs over the west flank (Fig. 5). This indicates that the weak material (lower Young's Modulus) infilling the west flank's scar can change the shape of the modelled deformation. However, the amplitude of the change is small; it can only account for a ∼10% variation in uplift shape (Fig. 5), and this is not sufficient to explain the observed asymmetric deformation in this study (Fig. 2), or those previous (Table 1).
We test this result further by simulating two possible endmember models. In these scenarios, the west flank has a lower limit Young's Modulus of 0.1 GPa or an upper limit Young's Modulus of 10 GPa. The difference between the modelled vertical uplift is ∼0.13 cm (over a maximum value of ∼1.4 cm). These results indicate that the weak material infilling the west-flank collapse scar plays a very minor role in producing asymmetric deformation at Tungurahua. This is due to the small volume of material within the west flank compared to the total volume of material that the stress from the deformation source passes through; therefore, it only produces a minor amplification effect.
The footprint of the altered deformation pattern induced by the weak flank material is also larger than the footprint of the observed deformation. This could provide further evidence that the west flank material is not the primary cause of asymmetric deformation at Tungurahua. However, there is also the possibility  of smaller-scale heterogeneity within the west-flank that we do not capture in our model (where the west-flank is described by constant material parameters). High-resolution geophysical imaging may be able to distinguish such variations (e.g., Finn et al., 2001). Enhanced geophysical imaging may also better delineate the spatial extents of the west-flank infill and could indicate a deeper base than currently considered; this would increase the volume of the weakened material and the deformation amplification effect. Regardless of these possibilities, in the current context, we have provided the first results to show that the heterogeneity of the west-flank cannot explain the observed asymmetric deformation at Tungurahua.

Optimum deformation source model and shallow magma storage
To assess the possibility of a deformation source directly beneath the west flank as the primary cause for asymmetric surface inflation we invert our data for an optimum source model. The final best-fit source has a cigar-like shape and is oriented east-west with a shallow dip, centred at a depth of ∼2 km asl, and located beneath the west flank (Table 2 and Fig. 6). The modelled overpressure of ∼1.2 MPa is equivalent to a volume change of ∼2×10 5 m 3 , and a volume change rate of ∼3×10 6 m 3 /yr. Source parameter uncertainties were calculated using plus-or-minus one standard deviation from the optimum values, after removing outliers. Outliers were determined using a cut-off residual misfit value of two standard deviations greater than the minimum (optimum) residual.
The optimum deformation source model produces a broadly good fit to the InSAR data ( Fig. 6 (Molina et al., 2005). Given the uncertainties in defining the depth of the collapse surface, in reality the deformation source might be exploiting the weakness of the collapse surface rather than sitting directly beneath it.
al., 2016), suggesting that the same source has been repeatedly active, however, we are able to better constrain the timeframe upon which source processes may occur due to our significantly shorter observation period (∼3 weeks, Table 1). In the model, the top of the deformation source sits just beneath the bottom of the westflank domain (Fig. 7). Given the uncertainty in the seismic imaging data used to constrain the modelled geometry of the west-flank (Molina et al., 2005), in reality, the deformation source might actually be exploiting the basal collapse surface and therefore providing a new constraint on the depth and dip of that surface. The modelled source is also at the rough depth of the interface between the Paleozoic-Cretaceous regional metamorphic basement and the Tungurahua volcanics (Hall et al., 1999;Molina et al., 2005).
The nature of the deformation source can be constrained using simultaneous observations. Given the volcano was in an eruptive phase during this deformation event, it is probable that the deformation source represents a site of shallow, short-term, preeruptive magma storage (Muller et al., 2018). The cigar like shape is possibly the result of a transition from a dyke, transporting magma from deeper within the crust, to a conduit, transporting magma to the surface where it can erupt. Short-term pressurisation (on the order of weeks), and thus surface deformation, reflects a temporary volumetric imbalance between magma inflow and outflow, and also likely causes the associated increase in VT seismicity noted during the same period (Supplementary Figure 2). The rapid pressurisation and deformation will have further implications for material failure, as the stress will accumulate over a short period, producing high stress and strain rates, and reducing the opportunity for stress to dissipate over time.
We do not observe deflation following the uplift and eruption, and deflation also has not been recorded after previous inflation episodes (Table 1). The lack of observed subsidence may indicate that the deformation of the shallow crust is non-elastic (e.g., elasto-plastic), and/or that a proportion of the intruded magma is still in-situ and contributing to continued growth of the edifice (e.g., Biggs et al., 2010).
The potential for shallow magma storage has previously been suggested at Tungurahua from melt inclusion studies, indicating depths between ∼4.2 and ∼1.9 km beneath the summit, or ∼0.8 -∼3.1 km asl (Myers et al., 2014), consistent with the ∼2 km asl result from our models (Table 2). To assess the short-term nature of the shallow magma storage, a diffusion chronometry study on eruptive products following a future co-or pre-eruptive west-flank inflation might indicate possible residence times (e.g., Druitt et al., 2012) and shed further light on storage conditions. As explored later, there is an emerging picture of vertically-extensive magma storage at Tungurahua volcano.
To compare our modelled volume change (∼2×10 5 m 3 ) to magmatic volumes we must account for magma compressibility (McCormick Kilbride et al., 2016). Accordingly, the ratio of erupted volume (V e ) to modelled volume change ( V) is: where β c is chamber compressibility and β m is magma compressibility (McCormick Kilbride et al., 2016). For a first-order estimate, our chamber compressibility is calculated as 1.7 × 10 −7 Pa −1 using results from Anderson and Segall (2011) for a vertical pipe (as an analogue for our cigar-shaped reservoir, which can be thought of as a rotated vertical pipe), and magma compressibility is between 1.5 × 10 −10 and 2.5 × 10 −10 Pa −1 (Muller et al., 2018). This produces a ratio, V e / V, of 1.9 to 2.5, so our modelled volume change (∼2 × 10 5 m 3 ) predicts an erupted volume of 3.8 × 10 5 -5.0 × 10 5 m 3 . Erupted volume estimates for just the period of study, November 2015, are not available. From September to December 2015 there was 1.1 × 10 6 m 3 of dense rock equivalent (DRE) ash-fall volume recorded, but this does not account for any possible ballistics or pyroclastic density current (PDC) deposits, and we do not see other deformation episodes in this window. The lack of further constraints makes comparing the longer-term erupted volume (1.1 × 10 6 m 3 ) to our predicted volumes (3.8 × 10 5 -5.0 × 10 5 m 3 ) difficult. However, if we take the most basic assumption of equal erupted volumes for each of the months from September to December, there would be ∼2.75 × 10 5 m 3 DRE ash-fall volume for November 2015. Assuming the PDC and ballistic erupted volumes over the same period are considerably smaller than the erupted volume of ash would suggest that some intruded magma remains in the shallow crust. This would be consistent with the lack of observed subsidence following the uplift and eruption.

West-flank stability
The extent to which shallow magma storage and rapid source pressurisation might influence flank stability warrants further investigation given the previous flank collapses, observed deformation and seismicity, assumed weak infill material, ongoing eruptive activity, structural discontinuity, steep slopes, and the hazard posed by a possible future collapse event. For a first analysis we apply our optimum deformation source parameters in a series of forward models, where the source overpressure is increased from 1 MPa to 100 MPa, simulating the intrusion (or storage) of increasingly large volumes of magma (Fig. 8). Magma pressures up to 100 MPa are investigated to gain a full understanding of model behaviour, and to include the possibility of large excess pressures prior to failure (e.g., Grosfils, 2007). Results show how the shear stress along the 3 ka collapse surface linearly increases with pressure, up to a maximum of ∼26 MPa. Interestingly, the single deformation source produces two lobes of enhanced shear stress on the collapse surface, with the maximum shear stress at the centre of the lobes, and the lobes located adjacent to the point directly above the source. Increases in shear stress are known to develop shear zones and promote flank instability (e.g., Donnadieu et al., 2001;Schaefer et al., 2013), but a range of values that may trigger collapse are unknown due to the dependence on a variety of other intrinsic and extrinsic factors (e.g., Donnadieu et al., 2001;Siebert, 1984;van Wyk de Vries and Delcamp, 2015). An additional factor may arise due to the observed rapid pressurisation of the magma source (section 4.2), which reduces the opportunity for accumulated stress to be dissipated and produces higher strain rates; both are more likely to promote material failure.
The linear relationship between source pressure and shear stress (Fig. 8) reflects the use of a linear elastic rheology, however, if suitable estimates for mechanical material parameters were available an elasto-plastic rheology might be more appropriate (e.g., Schaefer et al., 2013) and would likely cause shear stress to increase non-linearly with source pressure after a yield stress is exceeded. Furthermore, the repeated deformation of the west flank (Table 1) may contribute to cyclic stressing of the flank materials which can lead to a reduction in Young's Modulus (Heap et al., 2009), and might promote elasto-plastic behaviour at lower yield stress values. Similar reductions in Young's Modulus and promotion of non-elastic deformation can be caused by hydrothermal alteration; a recent study of tilt variations suggests the deformation modulus (strongly related to Young's Modulus) for Tungurahua's edifice is on the order of 10's of MPa due to the combined effects of alteration, fracturing and porosity (Marsden et al., 2019). These values are two orders of magnitude lower than those used in our current study, and if accurate, would imply the observed InSAR deformation could be reproduced with a lower deformation source overpressure than used here, and as such, the shear stress generated on the 3 ka collapse surface would also be reduced. However, such low values of Young's Modulus have not been demonstrated in laboratory or seismic experiments.
It has previously been suggested that the 3 ka flank collapse at Tungurahua might have been triggered by a cryptodome intrusion (Hall et al., 1999), similar to Mt. St. Helens (MSH) in 1980. The MSH collapse was preceded by a large bulge on the side of the volcano and significant flank over-steepening. Intuitively, one might expect shallow magma intrusion or storage to also cause a steepening of a volcano's flank. Our models indicate this effect is minor (maximum angle change < 1 • ) for the scenarios we test using the optimum deformation source parameters and a range of overpressures (1 -100 MPa). However, if the source region was shallower (i.e., closer to the flank surface and within the weak, clastic infill material) or oriented more vertically (e.g., related to a conduit) the steepening effect would likely be much greater. This scenario might be accompanied by an increase in rock-fall events due to the unconsolidated eruptive deposits making up the west flank, allowing for additional monitoring via seismic and infrasound techniques.
The extent to which the modelled magmatic source from this study could develop into a cryptodome promoting flank failure (similar to the 3 ka event) is unknown, but continued flank deformation monitoring (GPS, InSAR, tilt and EDM) will provide additional clues. InSAR in particular is a vital tool in this regard, given its ability to remotely measure deformation at a large number of the world's volcanoes over a range of spatial footprints. We have demonstrated the link between flank deformation and potential instability at Tungurahua, but our approach is easily applicable to other volcanoes displaying asymmetric or offset deformation patterns (e.g., Arenal, Villarica, Cerro Aquihuato, Puyehue-Cordon Caulle; Ebmeier et al., 2018). The increasing detection of such signals demonstrates that the underlying processes could be widespread and effort should be continued to understand the driving forces.
In addition to magmatic scenarios, there are a number of other triggers and priming factors for flank collapse that are applicable for Tungurahua. These include large tectonic earthquakes, hydrological conditions, hydrothermal alteration, gravitational-driven creep, structural geology, and edifice stratigraphy (Donnadieu et al., 2001;Finn et al., 2001;Schaefer et al., 2013;Siebert, 1984;van Wyk de Vries and Delcamp, 2015). A full geotechnical analysis of the interplay between these factors is recommended to fully classify the west flank slope stability and possible hazard at Tungurahua; especially given the implications of a future collapse could be severe, with several communities now situated in the runout from the 3 ka event.

Vertically-extensive magma storage
Our results highlight the presence and implications of shallow (∼2 km asl) magma storage at Tungurahua, and are backed up by petrological melt inclusion studies (Myers et al., 2014). Previous large-scale deformation studies have highlighted magmatic deformation sources at depths of 11.5 km bsl (Champenois et al., 2014) and 5 km bsl (Muller et al., 2018), broadly consistent with other petrological analyses (Andújar et al., 2017;Samaniego et al., 2011). This range of observations as a whole indicates vertically-extensive magma storage , where dynamic mass transfer between different levels of magmatic storage can cause unrest episodes and trigger eruptions . At the shallowest level where we focus our study (∼2 km asl), we identify a possible site of shallow, transient, magma storage at the transition between an eruption-feeding conduit and a dyke (maybe from the 5 km bsl source), preferentially exploiting the collapse surface as a zone of weakness. However, it has been proposed that in a trans-crustal or vertically-extensive magmatic system volcanic surface deformation may be primarily caused by volatile exsolution and expansion . In this scenario, the shallow deformation source we infer might represent a region of volatile accumulation from reorganization of the deeper plumbing system as it feeds the coincident eruptions. We consider these two interpretations as possible end-members with a number of hybrid models existing between them, but currently favour a predominantly magmatic source due to the ongoing eruptive activity. Continued geophysical imaging and monitoring will be required to fully understand the dynamics beneath Tungurahua's west flank, and their implications for slope instability.

Conclusions
We observe rapid and localised surface inflation on the west flank of Tungurahua volcano, coincident with an increase in eruptive activity and VT seismicity, in November 2015. A previous (3 ka) flank collapse and recent periods of inflation in this locality indicate dynamic subsurface processes are contributing to frequent unrest and possible instability. We model our observed deformation using Finite Element Analysis to test different mechanisms causing the asymmetric surface deformation, and relate them to shallow magma storage and slope stability.
Our results indicate that the weak clastic material infilling the west-flank collapse scar plays a very minor role in producing asymmetric deformation at Tungurahua. Using these new results, we conclude that the heterogeneity associated with this sector cannot explain the current or past episodes of west-flank inflation. To reproduce the observed deformation, we invert our InSAR data and derive a best-fit cigar-shaped source model, located beneath the west flank at a depth of ∼2 km asl with a shallow dip and oriented approximately east-west, indicating shallow pressurisation. This source is inferred to represent a site of shallow, short-term, pre-eruptive magma storage, where a temporary volumetric imbalance between magma inflow and outflow causes rapid pressurisation, deformation (equivalent to ∼48 cm/yr) and seismicity. The storage region is likely exploiting the 3 ka collapse surface as a zone of preferential weakness, but the volcano is also likely underlain by a vertically-extensive magmatic system.
Using our shallow magmatic source model in a series of tests shows that the shear stress along the 3 ka collapse surface increases linearly with pressure (or volume) change, and that these stresses have little chance to dissipate owing to the rapid source processes taking place. These new findings highlight a factor that can contribute to slope instability during future volcanic unrest episodes, where greater volumes of magma could be stored or transported in this shallow region. The extent to which the modelled magmatic source from this study could develop into a cryptodome causing flank failure (perhaps similar to the 3 ka event) is unknown, and must be considered alongside other dynamic priming and triggering factors.
Shallow magma transport and storage are key dynamic processes at the heart of hazard assessment and event forecasting for eruptive activity and flank collapse. The results of this study highlight the importance of monitoring volcanoes with a history of flank collapse and how this can be incorporated into routine geophysical and geodetic volcanic surveillance, and demonstrates the efficacy of satellite InSAR for that task.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.