Transient deformation associated with explosive eruption measured at Masaya volcano (Nicaragua) using Interferometric Synthetic Aperture Radar. Journal Volcanology

Deformation caused by processes within a volcanic conduit are localised, transient, and therefore challeng- ing to measure. However, observations of such deformation are important because they provide insight into conditions preceding explosive activity, and are important for hazard assessment. Here, we present measurements of low magnitude, transient deformation covering an area of ∼ 4 km 2 at Masaya volcano spanning a period of explosive eruptions (30th April–17th May 2012). Radial uplift of dura- tion 24 days and peak displacements of a few millimeters occurred in the month before the eruption, but switched to subsidence ∼ 27 days before the onset of the explosive eruption on 30th of April. Uplift resumed during, and continued for ∼ 16 days after the end of the explosive eruption period. We use a ﬁnite element modelling approach to investigate a range of possible source geometries for this deformation, and ﬁnd that the changes in pressurisation of a conduit 450 m below the surface vent (radius 160 m and length 700 m), surrounded by a halo of brecciated material with a Young’s modulus of 15 GPa, gave a good ﬁt to the InSAR displacements. We propose that the pre-eruptive deformation sequence at Masaya is likely to have been caused by the movement of magma through a constriction within the shallow conduit system. Although measuring displacements associated with conduit processes remains challenging, new high resolution InSAR datasets will increasingly allow the measurement of transient and lower magnitude defor- mation signals, improving the method’s applicability for observing transitions between volcanic activity characterised by an open and a closed conduit system. © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Volcanoes that transition between effusive and explosive behaviour pose a threat to nearby communities, but the dynamics behind this transition are not well-defined (Dingwell, 1996;Gonnermann and Manga, 2003). Characteristic eruptive behaviour at a volcano depends on the magma ascent rate, which is in turn determined by subsurface geometries, physical properties of the magma and overpressures within the system (Woods and Koyaguchi, 1994;Dingwell, 1996;Gonnermann and Manga, 2007). Eruptive activity can also sometimes be related to changes in the conduit structure, and particularly whether the magma plumbing system is open or closed to the atmosphere. Open conduit systems have an unrestricted connection between the magma reservoir and the surface, allowing the lava lake height to act as a pressure gauge (Patrick et al., 2015, e.g. Kılauea summit lava lake). Volcanoes that have a conduit geometry that prevents the free flow of gas and magma between the magma reservoir and the surface are known as closed conduit systems (Worster et al., 1993;Tazieff, 1994). Vulcanian eruptions, which are violent and short-lived, are thought to be due to the temporary closure of the conduit by a plug of cold, solidified magma in the conduit that allows for a build-up of pressure (Francis and Oppenheimer, 2004). These plugs are most commonly observed at silicic volcanoes, which have high viscosity magmas (Albino et al., 2011, e.g. Colima), but have also been found to form within basaltic volcanic systems (Lyons and Waite, 2011, e.g. Fuego). The presence of an ephemeral lava lake indicates that a volcano transitions between an open and closed conduit system (Tazieff, 1994), with the lake appearing or disappearing in a cyclic (periodic) or acyclic manner (Barker et al., 2003;Witham and Llewellin, 2006;Hirn et al., 2008).

ARTICLE IN PRESS
Masaya volcano in Nicaragua is an excellent example of a volcano that switches between effusive and explosive behaviour. It has a persistent high gas flux and contains an ephemeral lava lake, which suggests that the conduit system can transition between open and closed conditions. In this study, we used high temporal (1-14 days) resolution satellite data to capture transient conduit processes. COSMO-SkyMed (CSK) Synthetic Aperture Radar (SAR) acquisitions (ASI) were obtained over Masaya volcano during a period of explosive activity and no lava lake (30 April-17 May 2012). We present measurements of deformation at Masaya over 7 months spanning the explosive eruption. We use Finite Element Modelling (FEM) to investigate potential deformation source geometries and discuss potential deformation mechanisms.
Models developed to explain ground deformation resulting from conduit processes involve changes in flow dynamics within the volcanic conduit, which generate shear and normal stresses (e.g. Chadwick et al., 1988;Bonaccorso and Davis, 1999;Beauducel et al., 2000;Green et al., 2006), and/or changes in surface loads due to dome growth or collapse (e.g. Beauducel et al., 2000). Early analytical models developed to describe inflation of a vertical conduit include a closed pressurised pipe, a dislocating pipe, and inflationdeflation patterns of a vertical conduit, which were shown to be due to the removal of a blockage within the pipe (Bonaccorso and Davis, 1999). Numerical experiments performed by Albino et al. (2011) to examine the deformation field due to plug evolution within a conduit showed that the contribution of magma conduit flow to the radial deformation field is only significant for distances of less than 0.2 of the conduit length (e.g. 200 m away from the vent for a 1 km length conduit), with almost negligible contribution to the vertical displacement field at distances greater than a few hundred meters from the vent.
In the case of plug emplacement within a conduit, there are 2 potential mechanisms for plug removal (Albino et al., 2011): (1) stick-slip transition resulting from changes in overpressure beneath the plug (Lensky et al., 2008) or brittle failure conditions of the magma (Collier and Neuberg, 2006), and (2) changes in conduit wall permeability that prevents lateral degassing and removes the plug either by reducing its length or its viscosity contrast (Edmonds et al., 2003;Iguchi et al., 2008).

Previous studies of pre-eruptive displacement
Deformation related to conduit processes has primarily been observed days to seconds before the onset of explosions (e.g. Dzurisin et al., 1983;Chadwick et al., 1988;Rowe et al., 1998;Anderson et al., 2010;Voight et al., 2010;Albino et al., 2011;Lyons and Waite, 2011;Iguchi et al., 2008;Salzer et al., 2014;Mothes et al., 2015). This has included inflation, marking a transition between an open and closed conduit, for example ∼7 h before the start of eruptions at Colima (Salzer et al., 2014), but only a few seconds before the eruption at Erebus (Rowe et al., 1998). Another common feature at several volcanoes is contraction within the conduit prior to or at the onset of explosive activity (e.g. Voight et al., 2010;Iguchi et al., 2008;Lyons and Waite, 2011;Mothes et al., 2015). At Sakurajima, Semeru and Suwanesojima volcanoes, inflation occurred minutes before explosive activity, followed by contraction a few minutes to immediately before the explosions took place (Iguchi et al., 2008). Similarly, cycles of pressurisation-depressurisation prior to explosive activity have been observed at Soufrière Hills Volcano, at Mount St. Helens and at Fuego (Anderson et al., 2010;Voight et al., 2010;Lyons and Waite, 2011). The general interpretation of these pressurisation cycles is of plug emplacement and subsequent removal when the overpressure exceeded the plug strength. Plug ascent within the conduit has also been recorded by tiltmeters at elevation as an uplift-subsidence pattern (Mothes et al., 2015, e.g. Tungurahua). Other studies have found that shear stresses provide a better explanation for within the conduit radial displacements at Mount St. Helens (observation period 1981-1982), Mount Merapi and Soufrière Hills Volcano (Dzurisin et al., 1983;Chadwick et al., 1988;Beauducel et al., 2000;Green et al., 2006).

Geological setting of Masaya volcano
Masaya volcano (11.984N, 86.161W, 635 m) is a basaltic caldera located in south-western Nicaragua, and contains Masaya Lake, and several nested crater pits (namely San Pedro, Nindirí, Santiago, Masaya) and cones (namely Nindirí, Masaya) (Rymer et al., 1998;Delmelle et al., 1999;Williams-Jones, 2001;Duffell et al., 2003 see Fig. 1). Ring faults seen in an exposed section of the Nindirí crater pit dip outwards from the crater at an angle of ∼80 • (Rymer et al., 1998;Harris, 2009). The red dashed lines in Fig. 1 indicate the inferred location of a larger ring fault within the caldera, along which the most recent vents at Masaya are also located (Crenshaw et al., 1982;Rymer et al., 1998;Mauri et al., 2012;Caravantes Gonzalez, 2013;Zurek, 2016). The exact location and orientation of this larger ring fault, however, has not been well-constrained. Several profiles constructed using geophysical field methods (magnetics and Very Low Frequency (VLF) electromagnetics) suggest that the ring fault is inclined inwards into the caldera (Caravantes Gonzalez, 2013).
Masaya has rapidly transitioned between sudden explosive activity and quiescent degassing in the past (e.g. 23 April 2001 VEI 1 explosive eruption), and thus poses a hazard to tourists visiting the volcano (Global Volcanism Program, 2001). This increases the importance of using observational methods that can provide information on transient processes occurring within the ephemeral system.

Previous observations & structural interpretation
Previous studies of Masaya have proposed several subsurface geometries based on deformation, seismicity, gas geochemistry observations and geochemical methods, and geophysical techniques (self-potential, electromagnetic and potential field methods) ( Table 1). The models agree that there is a shallow magma reservoir present below Masaya, with estimated depths ranging between 0.4-5 km below the active vent (Métaxian et al., 1997;Rymer et al., 1998;Williams-Jones et al., 2003). Early models from seismicity suggest an open system at Masaya, in which the source locations from permanent tremor indicate a shallow source at 400 m below  Crenshaw et al. (1982, pp. 346) and Caravantes Gonzalez (2013)). Blue stars indicate the locations of the closest towns. White box indicates region plotted in (b). b) Masaya crater pit area. Red dashes indicate the inferred extent of the larger ring fault from (a). Thick white dashes indicate extent of the Masaya and Nindirí cones, while thin white dashes indicate extent of the crater pits. Green line indicates location of the Nindirí ring fault, and green dashes indicate the inferred extent of the Nindirí ring fault. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the active vent between 1992 and 1993 (Métaxian et al., 1997). Multidisciplinary observations between 1993-2001 (microgravity, SO 2 flux and ground deformation) suggest that the volcanic system extends under the San Pedro, Nindirí and Santiago crater pits, with a shallow reservoir at ∼1 km below sea level (which corresponds to 1-2 km below the active vent) (Rymer et al., 1998;Williams-Jones et al., 2003). This shallow subsurface region is thought to consist of magma pathways within a brecciated region extending to ∼800 m depth, with an accumulation of bubble-rich magma near the surface (Rymer et al., 1998;Williams-Jones et al., 2003). Another model developed using these microgravity and SO2 2 flux datasets (Rymer et al., 1998;Williams-Jones et al., 2003) simply considers a shallow reservoir with a foam layer at the top connected to a deeper reservoir by a narrow vertical conduit (2-6 m radius) (Stix, 2007). Self-potential, soil CO 2 and ground temperature observations suggest that the hydrothermal system at Masaya is more complex than previously envisioned, and that hidden structures, such as faults and dykes, influence the movement of groundwater flow within the caldera (Mauri et al., 2012). A SW-NE profile across the active Santiago crater suggests the presence of a shallow magma plumbing system beneath Santiago and a shallow groundwater flow less than 150 m below the surface.
Very few geodetic observations have been made at Masaya due to lack of ground instrumentation, and spatial and temporal coverage of satellite datasets. While GPS data at Masaya is scarce and poorly sampled, long-term subsidence (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) has been observed from re-measuring eGPS sites within the caldera (LaFemina, P. & Geirsson, H., pers. comm., 2016). Previous InSAR observations at Masaya have been performed using ALOS-1 SAR acquisitions between 2007 and 2010 (Ebmeier, S.K., pers. comm., 2017) (see Supplementary Table  3 and Supplementary Fig. 1). A steady subsidence on the order of ∼20 mm/year was observed within the inferred large ring fault. This subsidence occurred during a period of lava lake activity, as well as various ash, steam and phreatomagmatic eruptions from the Santiago crater pit (Global Volcanism Program, 2009;Global Volcanism Program, 2011;Ebmeier et al., 2013b). Any transient deformation that may have been associated with lava lake activity and conduit processes was not captured due to the poor temporal resolution of the ALOS-1 data (46 days). Lower magnitude deformation may also have been undetected as the average minimum detection threshold for this dataset was estimated to be 13 mm/year over Masaya (Ebmeier et al., 2013b).

Masaya activity 14th January-1st August 2012
During our period of InSAR observations, there was no lava lake present, and no magma was erupted during the explosive period (30 April-17 May 2012) (INETER, 2012a,e). Observations made by Instituto Nicaragüense de Estudios Territoriales (INETER) showed that the average flux of SO 2 for the SW-trending plume had been increasing since December 2011, reaching a peak of ∼ 1000 tonnes per day (t/d) on 20 March 2012. Incandescence from the crater was reported on the evening of the 27 April, and at 0500 the following morning  (Global Volcanism Program, 2012a). At 0800, a strong VEI 1 explosion of ash and gas was seen from the Santiago crater, reaching up to 1 km above the summit. Pink and yellow sandsized ash, lapilli and blocks up to 10 cm in diameter were ejected from the vent. Sporadic explosions of gas and ash were recorded throughout the explosive period (30 April-17 May), with a decreasing SO 2 flux from 760 t/d to 530 t/d observed through April and May (Global Volcanism Program, 2012b). In some areas, the ash fall thickness was 2 mm, and a total ejecta volume of 736 cm 3 was calculated by INETER. Masaya returned to normal high degassing volumes ( >1000 t/d SO 2 ) after the end of the explosive period, with no explosive activity or incandescence observed at the Santiago crater pit, and no notable seismic activity (INETER, 2012c,d). The SO 2 flux again increased to ∼ 1970 t/d by the end of July, dropping to 500 t/d in mid-December (INETER, 2012b,c).

Interferogram processing
40 ascending CSK SAR acquisitions were obtained over the Masaya region through the Committee on Earth Observation Satellites (CEOS) Volcano Pilot Project for the period 14 January-1 August 2012 (X-band wavelength, 3.1 cm). Interferograms were processed using the JPL/Caltech InSAR Scientific Computing Environment (ISCE) software (Rosen et al., 2011). Topographic phase contributions were removed using a 30 m Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM). A total of 154 interferograms were selected using the criteria that the perpendicular baseline separation between the satellite positions at the time of acquisition was 230 m or less. The interferograms were unwrapped using the snaphu network-flow algorithm (Chen and Zebker, 2001). The unwrapped interferograms were downsampled further to 110 m resolution in order to decrease computation time, and a coherence threshold of 0.1 was applied to remove decorrelated areas.

Correcting for atmospheric artefacts, orbital ramps and unwrapping errors
Topographically-correlated atmospheric artefacts have the potential to mask true deformation at volcanic centers (e.g. Ebmeier et al., 2013a), and thus pairwise logic was used to visually inspect for atmospheric artefacts (e.g. Massonnet and Feigl, 1995). It was found that approximately half of the interferograms were affected by atmospheric contributions, either at large-or small-scales. North American Regional Reanalysis (NARR) atmospheric model data (specifically temperature, specific humidity and geopotential datasets) were used to create atmospheric phase screens (APS) for the dry delay and wet delay for each CSK acquisition date (e.g. Parker et al., 2015). These phase screens were downsampled and applied to the downsampled unwrapped interferograms. Interferograms were then re-referenced by taking an average of the cumulative displacement within the caldera, with the crater pits, where deformation was localised, masked out. Since the spatial resolution of the weather model is low (32 km), only large-scale atmospheric artefacts were removed from the interferograms. Small-scale turbulent atmospheric artefacts that were identified using pairwise logic were manually removed. For example, a small turbulent artefact was observed in interferogram pairs associated with date 15 June 2012 to the south of the crater pits.
Linear and quadratic ramps were applied to remove orbital errors from the final dataset. Obvious unwrapping errors were either removed or manually corrected according to their spatial extent (i.e. a single anomalous pixel was removed, while a region of 4 or more pixels was corrected for). Manual correction of unwrapping errors was performed by selecting the region of interest, and either increasing or decreasing the displacement by a factor of k/(4p), where k is the wavelength of the satellite radar.

Time-series analysis
Small BAseline Subset (SBAS) time-series analysis was performed using a linear least squares inversion for a network of well-connected interferograms over the whole caldera in order to examine and characterise the relative ground displacement around the crater pits before, during and after the explosive period (e.g. Lundgren et al., 2001;Berardino et al., 2002). The displacement of the first acquisition date is explicitly assumed to have no deformation (e.g. Lundgren et al., 2001). Since the magnitude of the deformation is low, we referenced the interferograms to a small annulus region centered over the crater pits (inner radius ∼0.49 km, outer radius ∼1.5 km) to reduce the contribution of any remaining spatially correlated atmospheric artefacts to the time-series over Santiago. This means that the displacements presented here are all relative to displacements a distance of less than 1.5 km away. For estimation of time-series errors, the standard deviation of the average of the caldera, with a broad elliptical region over the crater pits masked out (radius 1.33 km), was computed using Monte Carlo inversion for 200 runs (as described by Ebmeier et al., 2010).
Bootstrapping was performed on an initial subset of 64 interferograms (selected according to network connectivity) in two parts: (1) bootstrapping of individual interferograms, and (2) bootstrapping of individual satellite acquisition dates. This was performed to examine the reliability of interferograms in the time-series network, as well as to assess the presence of any remaining atmospheric artefacts in each interferogram. 54 interferograms (constructed from 22 CSK acquisition dates, see Supplementary Table 4 and Supplementary Fig. 2) were chosen after the bootstrapping analysis according to coherence within the Masaya caldera and around the crater pits, connectivity of the interferogram network, and the reliability of the interferograms and acquisition dates.

Preliminary assessment of corrected interferograms
The coherence across all the interferograms was found to depend on the time of the year, the interferogram timespan, lava flows, presence of vegetation and water bodies. Within the caldera, there is good coherence along the 1670 and 1772 lava flows and the pre-Holocene flows on the caldera floor. The contribution of the topography was examined for several data points in the area of interest around the crater pits by comparing perpendicular baseline to phase (see Supplementary Fig. 3). Linear regression values for a plot of perpendicular baseline against phase were less than 0.5 and thus it was concluded that there are probably no residual topographic contributions and that the DEM errors do not dominate the phase in these locations. Fig. 2 demonstrates the different corrections that were applied to the deformation field for a single unwrapped interferogram (11 March-3 April). These results demonstrate the importance of careful application of corrections to avoid introducing additional artefacts and over-interpretation of the data due to the sub-centimeter magnitudes observed (Saepuloh et al., 2013). The application of NARR APS to the interferograms appears to have removed large-scale atmospheric artefacts (3-10 km scale), and also greatly reduced apparent topography-related tropospheric contributions around the crater pits in several interferograms (e.g. Fig. 2b). The removal of orbital ramps for several interferograms was shown to be useful in removing broad deformation patterns that persisted across an entire interferogram (e.g. Fig. 2c). Since our measurements of displacements near the pit-craters were made relative to a reference area less than 1.5 km away, we do not expect any residual errors in estimation of orbits to affect our results.

Observed displacements at Masaya
Throughout the bootstrap analysis of our InSAR-derived timeseries, there are three phases of deformation that persist. We find the turning points in the interpolated time-series of displacement within 2 km from Santiago crater, and use the dates of peak uplift (3 April) and peak subsidence (28 April) to define separate phases of deformation for our analysis and modelling. From approximately 11 March-3 April, radial uplift is centered over Santiago crater. After 3 April, subsidence centered on Santiago continued until 28 April when uplift resumes. The observation that the second turning point in our time-series occurred within 48 h of explosive eruption onset, and on the day that volcanic tremor was first recorded, gives us confidence that our measurements are dominated by deformation rather than atmospheric artefacts.
The time-series plots shown in Fig. 3 are for 54 selected interferograms (22 CSK dates) over two regions as shown in Fig. 3a: (1) the Santiago crater pit (radius 0.55 km, dark red circle) (Fig. 3b), and (2) the northern flank of Santiago crater pit (radius 0.44 km, bright pink circle) (Fig. 3c). The active Santiago crater pit region was chosen in order to examine potential contributions of conduit processes related to lava lake activity (or lack thereof) to the deformation field, while the northern flank was chosen to demonstrate the region of maximum deformation observed throughout the deformation phases. The magnitude of cumulative displacement for both time-series plots is notably higher for the maximum deformation region, however we still observe the three persistent deformation phases. As the theme of this paper is on conduit processes related to lava lake activity, the discussion will focus mainly on the time-series results over the active Santiago crater pit region. The first uplift (11 March-3 April) occurred over a period of ∼3 weeks, reaching peak uplift of ∼2 mm roughly 4 weeks prior to the explosive period (Fig. 4a). This peak uplift extends radially ∼1.5 km from the active vent. Subsidence of ∼5 mm then occurred over ∼27 days leading up to the explosive period (3 April-28 April) (Fig. 4c). The subsidence pattern is observed in the same spatial region as the first uplift period. The second period of uplift (28 April-16 July) occurred gradually over 2 months during and after the explosive period, appearing to stabilise ∼16 days after the end of the explosive period (Fig. 4e). The deformation observed for this period has a smaller observed areal extent of less than 1 km radially from the active vent. Estimated errors for the time-series data points vary between ∼2.5 mm and 4 mm. Assuming constant deformation, we calculated the estimated minimum detection threshold for this study to be 2.2 mm/year. This was performed by taking the expected magnitude of noise of the cumulative displacement for a point within the caldera divided by the total timespan of the interferograms.

Profiles for modelling
We used stacking (e.g. Schmidt and Bürgmann, 2003) to improve the signal-to-noise ratio for each deformation period (Fig. 4a, c, e). Profiles (length 2.2 km) for each time period were obtained by taking an average of the cumulative displacement at a radial distance interval of 110 m within a 105-degree arc from the central vent location in Santiago (Fig. 4b, d, f), starting at a distance of about 500 m to the north of the Santiago crater pit (Fig. 4a). We chose this region as we have more confidence in apparent deformation observed to the north of the crater pits (i.e. less turbulent atmospheric artefacts), and previous studies have suggested that the majority of the magmatic system is located to the north of the crater pits (e.g. Métaxian et al., 1997;Williams-Jones, 2001;Williams-Jones et al., 2003). The error for each averaged radial distance along the profile was calculated by taking the standard deviation for each radial bin. The standard deviation for each radial bin was found to exceed the variance used to estimate the time-series errors (as described in Section 3.3).
As a preliminary test to examine the sources for each deformation period, we computed correlation plots of the subsidence profile against the first uplift episode profile (11 March-3 April) and against the second uplift episode profile (28 April-16 July) (see Supplementary Fig. 4). For both plots, the linear correlation was high, with a linear regression of 0.945 for the first uplift correlation plot (11 March-3 April), and 0.989 for the second uplift correlation plot (28 April-16 July). These high correlations suggest that the source for all the periods of deformation is the same.
We focus our efforts on modelling the pre-eruptive subsidence (∼5 mm subsidence between 3 April and 28 April) (Fig. 4d).

Modelling subsidence (3 April-28 April)
We used COMSOL Multiphysics 5.1, a Finite Element Analysis software, to numerically model the observed InSAR profiles at Masaya. The decompression of a magma reservoir in homogeneous media was simulated to replicate the observed InSAR deformation by imposing a pressure change on a buried cavity. We tested for various geometries and depths of reservoirs since the shape and amplitude of the deformation is dependent on the size, shape and depth of the deformation source, the rock mechanics and the imposed pressure change. We tested two different numerical models: a homogeneous (single domain) model and a heterogeneous (two domain) model. The model is a 2D axisymmetric domain, 30 km in radius and in height, with boundary conditions as set out by Hickey and Gottsmann (2014). The size of the model and the mesh were tested for convergence and found to be satisfactory. We used a simple ellipse shaped cavity which becomes an ellipsoid by rotation during computation, and started with the simplest model setting of a homogenous elastic halfspace. It is important to state, however, that the assumption of symmetry in the deformation field would not necessarily hold for all origins of deformation. Ignoring topography and heterogeneities in near surface deformation modelling has been shown to have a detrimental impact on the results (Young and Gottsmann, 2015), however the topography at Masaya is complicated and the simplifications necessary to render it in an axisymmetric setting could themselves introduce errors. As discussed in Section 5.1, the profile used for modelling is an average of a 105-degree swath, and does not take into account variations in topography that are not axisymmetric. Additionally, when tested, the effect of topography on predicted displacements is within the expected InSAR error. We therefore neglect the topography for this model setting.
There is little available information about Masaya's crustal properties (e.g. seismic tomography), and we therefore use constant material properties for the model (i.e. the properties do not change with depth). The chosen material properties of the model are based upon the natural rock strength of basalts. Estimates of Young's modulus range from 22 GPa for Etnean basalts (Heap et al., 2009(Heap et al., , 2007 to 36 GPa for fresh Icelandic basalt (Heap et al., 2007), and a Poisson's ratio for volcanoes in Nicaragua of 0.25 (Cailleau et al., 2007). However, experimental pressurisation and depressurisation cycles on Etnean basalts have the effect of lowering the Young's modulus by 30% (from 33 to 22 GPa) and raising the Poisson's ratio by a factor of 3 (from 0.2 to 0.5) (Heap et al., 2009). Mechanical testing of basalts from Pacaya volcano shows that rock strength decreases with porosity, increases with strain rate and temperature, and has a variable response to thermal stressing (Schaefer et al., 2015). We therefore selected 20 GPa as an appropriate lower limit of Young's modulus for the whole model domain (E1), and 0.25 for Poisson's ratio.
The pressure change imposed is constrained by the tensile strength of the rocks. Subsidence of the volcano implies a depressurisation event, and therefore tensile deformation of the surrounding rock. Natural rocks have a maximum tensile deformation of up to −10 MPa (Gudmundsson, 2011). As the subsidence was not accompanied by an eruption, and consequently rock failure, we assumed that the tensile strength limit was not reached and selected −10 MPa as the maximum pressure.
For the homogeneous model, we performed parametric sweeps of the reservoir length, reservoir radius and reservoir top depth (meters below the surface) in COMSOL. Incremental sweeps of 200 m for reservoir length, 100 m for reservoir top depth and 50 m for reservoir radius were performed, with the parameter limits expanded until the upper and lower limits of each parameter was found (see Supplementary Table 5).
Resultant models were only selected if they fell within the expected InSAR error at every data point using a weight sum of squared errors as per Eq. (1): where O is the observed (COMSOL) data and E is the expected (InSAR) data. w 2 = 0 means a perfect fit to the InSAR data, and at each data point, if w 2 i is <1, the observed data fits within the InSAR error. Therefore, if w 2 is <n and each w 2 i is <1, the model fits the InSAR data, and the closer w 2 is to 0, the better the fit. We then explored modelling a heterogeneous halfspace as the contribution of the shallow subsurface required further examination according to previous Masaya models (Rymer et al., 1998;Williams-Jones et al., 2003). The heterogeneous model was constructed to include the shallow vertical reservoir surrounded by a cylindrical domain with a lower Young's modulus to represent fractured rocks from repeated pit collapse, which has been previously suggested by Rymer et al. (1998) and Williams-Jones et al. (2003). The effect of a halo of weakened rock around a shallow conduit style reservoir was tested using the concepts given by Young and Gottsmann (2015). The halo radius is fixed at 500 m, approximating the extent of recent pit craters. The halo extends from the model surface to 200 m below the base of the conduit to prevent boundary errors. Incremental sweeps of 100 m for reservoir length, 10 m for reservoir depth, and 20 m for reservoir radius were performed, again expanding the parameter limits until the upper and lower limits of each parameter was found. However, for this model, the upper radius limit was set to 250 m, as this is the Santiago crater radius. The halo Young's modulus was also tested (E2) at 5, 10 and 15 GPa (see Supplementary Table 5). Again, the resultant model was only selected if it fell within the expected InSAR error.

Source modelling results
The results for the homogeneous and heterogeneous domain models for the subsidence period are shown in Fig. 5, with the final model parameters given in Table 2. For the homogeneous models, the best fit radius and length (according to the w 2 criteria) for both the shallowest and deepest possible models are shown, as well as the best fit model with best fit radius, length and depth (Fig. 5a, c and Table 2a). Since the model is a homogeneous halfspace, the depths are with reference to the reference surface at 0 m. The horizontal and vertical reservoirs vary the shape and amplitude of the resultant deformation signal, as does the depth of the modelled reservoir. All of these models fall within the expected InSAR error, with the bestfit model being a sill at ∼ 1500 m depth with a w 2 of 1.6. However, because we expect magma to be present in the shallow subsurface (e.g. Métaxian et al., 1997, Rymer et al., 1998, Williams-Jones et al., 2003, we examined the possibility of a shallower deformation source further. The homogeneous shallow vertical reservoir fit is greatly improved by the addition of the halo domain (Fig. 5b, c and Table 2b), where a w 2 of 1.9 approaches the w 2 obtained for the best-fit homogeneous sill model (Fig. 5c). The conduit radius of 160 m is wellconfined to the width of the Santiago crater pit (radius 250 m), and the depth to the top of the vertical reservoir, 450 m, is deeper than the 400 m depth suggested by Métaxian et al. (1997) and shallower than ∼800 m modelled by Williams-Jones et al. (2003). The length of the reservoir, 700 m, however corresponds well with the suggested length of 600 m for the "conduit region" connected to a shallow reservoir from the same model by Williams-Jones et al. (2003).
All 4 models presented above fit within the expected InSAR error, and highlight the non-uniqueness of modelling solutions and the importance of modelling using a priori knowledge from previous studies to find the most plausible best-fit. While we recognise that the deformation field observed at a volcano can be the superposition of multiple sources at different depths and locations (Bonaccorso and Davis, 1999), we do not consider that the CSK interferograms contain information sufficient to constrain any deeper source. We consider the heterogenous shallow source model as the best model because it both accounts for the majority of the displacement signal, and agrees well with previous models of the shallow subsurface at Masaya (e.g. Métaxian et al., 1997;Rymer et al., 1998;Williams-Jones et al., 2003).

Potential origins of deformation
Our estimation of the turning point in our InSAR time-series of deformation places the transition from uplift to subsidence on 3 April (Fig. 3b). The first signs of activity above the background persistent degassing were observed on 21 April, when jetting sounds could be heard from deep within the crater, and there was evidence of rockfall collapse on the western wall of Santiago crater pit (Global Volcanism Program, 2012a). The lack of explosive activity during the subsequent period of subsidence suggests that the depressurisation of the system must have resulted from a non-explosive mechanism that caused the system to destabilise. Depressurisation of the system at Masaya may have occurred in several ways (either individually or (1) The start of (or increase in) flow through a constriction in the conduit (magma and/or gas), (2) the removal of a plug or obstacle within the system, (3) a decrease in shear stresses between the magma and reservoir walls (e.g. Anderson et al., 2010;Chadwick et al., 1988;Iguchi et al., 2008;Salzer et al., 2014), or (4) slip on the inferred larger ring fault (e.g. Fig. 1). The presence of an ephemeral lava lake at Masaya suggests regular changes in the subsurface magmatic plumbing. The absence of a lava lake during our period of observation implies that there was either a relative lack of magma being supplied into the system or there was a change in the conduit geometry preventing the flow of magma (Tazieff, 1994;Worster et al., 1993). As Masaya is persistently degassing, we assume that there was no lack of gas-rich magma within the system, and that rather there was an obstruction (such as a constriction or a blockage) preventing the free flow of magma to the surface. We know from previous models that there is a shallow brecciated region (500-800 m depth) that probably contains magma pathways to the surface and fault structures due to repeated pit collapse that spans the width of the San Pedro, Nindirí and Santiago pit craters (Rymer et al., 1998;Williams-Jones et al., 2003). This brecciated region then appears to narrow with depth and is connected to a magma reservoir 1-2 km below the surface vent (Williams-Jones et al., 2003). Our preferred FEM model places

ARTICLE IN PRESS
the deformation source at depths of between 450 m and 1150 m. We therefore expect any obstruction to flow to be located at the top of this complicated brecciated region. Stix (2007) proposed that the shallow reservoir is connected to the deeper reservoir by a narrow (2-6 m radius) conduit (depth not specified). Although our preferred shallow conduit source would lie above this constriction, it may be more consistent with our best-fit homogeneous halfspace sill model (depth 1.5 km).
An alternative explanation is that the pre-eruptive subsidence was caused by a decrease in shear stress. Shear stresses are due to pressure gradients that develop along the conduit walls, and an abrupt decrease in shear stresses has been shown to result in subsidence at the surface (e.g. Green et al., 2006). Subsidence has been attributed to drops in shear stress at Mount St. Helens (Chadwick et al., 1988;Dzurisin et al., 1983), Mount Merapi (Beauducel et al., 2000), and Soufrière Hills Volcano (Green et al., 2006). However, shear stresses have only been shown to cause deformation at high viscosity systems (e.g. Beauducel et al., 2000;Chadwick et al., 1988;Dzurisin et al., 1983;Green et al., 2006), and most examples of conduit pressurisation attributed to abrupt decreases in shear stress have occurred rapidly (days to seconds before explosive activity), while subsidence at Masaya was observed over ∼27 days. It is also unlikely that the deformation we measure records the gradual ascent of a discrete body of magma (e.g. Mothes et al., 2015), given that the spatial pattern of deformation is similar in the months before the eruption.
The timing of deformation switching from subsidence to uplift at the onset of eruptive activity suggests that the primary cause of deformation is related to a change in magma movement in the conduit. However, ring faults have been shown to accommodate uplift and subsidence at other volcanic centers during eruptive activity (e.g. Bathke et al., 2015Bathke et al., , 2013Lipman, 1997;Walker, 1984) (see Supplementary Table 2). Previous studies from ALOS-1 interferograms suggests subsidence of ∼20 mm/year over a three year period in the southern sector of the fault during a period of lava lake activity and phreatomagmatic explosions. In our study, we do not observe features in the deformation at the nested crater pit ring faults (Fig. 1b, green lines), or along the inferred larger ring fault (Fig. 1, red dashed  lines). It is possible that conduit pressurisation had an impact on the ring fault slip, perhaps through static stress changes, or by causing changes in pore-fluid pressure, and that may have contributed to the observed deformation field. The location and geometry of the larger ring fault is poorly constrained, and the poor coherence within the crater pits region, through which the larger ring fault is inferred to pass, present challenges for assessing any ring fault contribution to the deformation field.
Our observations are in some respects analogous to plug emplacement (associated with Strombolian and Vulcanian-style eruptions) at other basaltic volcanoes (Lyons and Waite, 2011, e.g. Fuego), although the processes we observe took place on much longer timescales (pre-eruptive uplift over ∼24 days, and pre-eruptive subsidence over ∼27 days) and the deformation source is at a greater depth (500-800 m) (Rymer et al., 1998;Williams-Jones et al., 2003). The general concept of pressurisation due to plug emplacement suggests that once a plug is emplaced, the overpressure increases beneath the plug and the edifice inflates (Albino et al., 2011). Once the overpressure exceeds the strength of the plug, it is released by the escape of gas, and the edifice begins to deflate (e.g. Albino et al., 2011;Anderson et al., 2010;Bonaccorso and Davis, 1999;Iguchi et al., 2008). Previous studies found that minutes to seconds after this blockage failure, eruptive activity was observed at the surface and the edifice continued to deflate (e.g. Iguchi et al., 2008). At Masaya, however, the gradual depressurisation before the explosive activity suggests an increase in magma (or gas) flow rate over time as opposed to the abrupt opening of a previously blocked pathway. This suggests that the magma ascent is more likely to have gradually increased through a constriction than reached an overpressure that suddenly exceeded the strength of the blockage.
If the deformation pattern across the crater pits were asymmetric, we would require a more complicated model than the simple axially symmetric conduit pressurisation used here, but poor coherence and prevalence of atmospheric artefacts to the south of the crater pits make this difficult to assess. Additionally, the presence of an active hydrothermal system at Masaya (e.g. Mauri et al., 2012) suggests that the movement and/or trap of fluids near the active vent may have contributed to the deformation field. Testing for this contribution would require measurements of ground water movement over the same period of observational time, geological observations on the ash products released during the 2012 explosive period, as well as a modelling approach incorporating fluid flow mechanics.

Implications
As the resolution of satellite data improves, more studies have been able to explore conduit processes and geometries at volcanoes worldwide (e.g. Ebmeier et al., 2014;Hamling et al., 2016;Pinel et al., 2011;Salzer et al., 2014;Wang et al., 2015). We have shown that it is possible to detect low magnitude, transient deformation associated with shallow conduit processes at a volcano with an ephemeral lava lake. Our results also suggest that conduit processes may contribute to the radial deformation field at distances greater than 0.2 of the conduit length (i.e. up to 0.8 of the conduit length).
Distinguishing between the different depressurisation mechanisms occurring at a volcanic system is important for prediction and hazard assessment. Previous studies at Mount St. Helens have shown that different depressurisation mechanisms can be operating at different times during the volcano's lifetime (Anderson et al., 2010;Chadwick et al., 1988). The timescales over which depressurisation occurs also provides important information about the pre-eruptive conditions that lead to an eruption.

Conclusions
We analysed InSAR data spanning a period of explosive activity at Masaya volcano, and our results demonstrate that the April-May 2012 explosive sequence was associated with local, low magnitude displacements before and during the explosive period. Uplift over the active Santiago crater pit was observed over a period of ∼24 days reaching a peak uplift of ∼2 mm, followed by subsidence of ∼5 mm over ∼27 days leading up to the explosive period. Uplift was observed again over ∼2 months from the onset of the explosive activity, and appeared to stabilise ∼16 days after the end of the explosive period.
We performed finite element modelling, and our results showed that a heterogeneous model of a shallow vertical reservoir (∼450 m depth) surrounded by a weaker brecciated region gave a good fit to the observed InSAR displacements. The depressurisation leading up to the explosive period followed by the repressurisation of the system at the onset of explosive activity is evidence of transitions between open and closed conduit behaviour. Using these results and models from previous studies (e.g. Métaxian et al., 1997;Rymer et al., 1998;Stix, 2007;Williams-Jones et al., 2003), we proposed that the pre-eruptive uplift-subsidence may have resulted from magma flowing through a constriction within the shallow conduit system. The constriction creates an increase in pressure as magma flows through it, and the subsequent gradual decrease in pressure occurs as the magma flow rate increases to overcome the constriction. Future models will require more detailed observations and measurements of the exact location and geometry of the larger ring fault in order to assess the potential contribution of motion on the ring fault to the deformation field.

ARTICLE IN PRESS
The results presented in this study highlight the potential of high temporal and spatial resolution satellite data for detecting shallow, transient conduit processes. Observations of the timescales over which shallower source, pre-eruptive deformation takes place are important for identifying the changes in magmatic plumbing that can lead to an eruption.

Author contributions
KJS processed and analysed InSAR data with assistance from SKE. NKY performed all COMSOL FEM modelling. This work originates in KJS's MSc thesis at the University of Bristol, supervised by SKE and JB. All authors read and contributed to improving the manuscript.