From quiescence to unrest: 20 years of satellite geodetic measurements at Santorini volcano, Greece

Periods of unrest at caldera‐forming volcanic systems characterized by increased rates of seismicity and deformation are well documented. Some can be linked to eventual eruptive activity, while others are followed by a return to quiescence. Here we use a 20 year record of interferometric synthetic aperture radar (InSAR) and GPS measurements from Santorini volcano to further our understanding of geodetic signals at a caldera‐forming volcano during the periods of both quiescence and unrest, with measurements spanning a phase of quiescence and slow subsidence (1993–2010), followed by a phase of unrest (January 2011 to April 2012) with caldera‐wide inflation and seismicity. Mean InSAR velocity maps from 1993–2010 indicate an average subsidence rate of ~6 mm/yr over the southern half of the intracaldera island Nea Kameni. This subsidence can be accounted for by a combination of thermal contraction of the 1866–1870 lava flows and load‐induced relaxation of the substrate. For the period of unrest, we use a joint inversion technique to convert InSAR measurements from three separate satellite tracks and GPS observations from 10 continuous sites into a time series of subsurface volume change. The optimal location of the inflating source is consistent with previous studies, situated north of Nea Kameni at a depth of ~4 km. However, the time series reveals two distinct pressure pulses. The first pulse corresponds to a volume change (ΔV) within the shallow magma chamber of (11.56 ± 0.14) × 106 m3, and the second pulse has a ΔV of (9.73 ± 0.10) × 106 m3. The relationship between the timing of these pulses and microseismicity observations suggests that these pulses may be driven by two separate batches of magma supplied to a shallow reservoir. We find no evidence suggesting a change in source location between the two pulses. The decline in the rates of volume change at the end of both pulses and the observed lag of the deformation signal behind cumulative seismicity, suggest a viscoelastic response. We use a simple model to show that two separate pulses of magma intruding into a shallow magma chamber surrounded by a viscoelastic shell can account for the observed temporal variation in cumulative volume change and seismicity throughout the period of unrest. Given the similarities between the geodetic signals observed here and at other systems, this viscoelastic model has potential use for understanding behavior at other caldera systems.


Introduction
Although caldera-forming volcanic systems spend most of their lifetime in a state of quiescence, they may also undergo extended periods of restlessness, which may or may not be accompanied by eruptive activity [e.g., Newhall and Dzurisin, 1988;Orsi et al., 1996;Wicks et al., 2006]. The past two decades have seen notable examples of seismic unrest and deformation without eruption at calderas including Campi Flegrei, Italy [e.g., Del Gaudio et al., 2010;Aiuppa et al., 2013], Yellowstone [Wicks et al., 2006], and Long Valley [Hill and Prejean, 2005]. Rabaul caldera (Papua New Guinea) showed caldera-wide uplift and increasing seismicity from 1971 to 1984 but did not erupt until 10 years later, in 1994later, in [McKee et al., 1985. The deformation histories of restless calderas are often temporally complex and may or may not correlate with earthquake seismicity (e.g., Santorini [Konstantinou et al., 2013], Campi Flegrei [D'Auria et al., 2011], and Rabaul [Saunders, 2001]). The underlying causes of these episodes of unrest have been attributed to both hydrothermal and magmatic processes, or combinations of the two (e.g., Nisyros [Gottsmann et al., 2005], Yellowstone [Waite and Smith, 2002], and Campi Flegrei [Chiodini et al., 2012] • Readme • Table S1 • Figure S2 • Figure S3 • Text S4 • Figure S5 Correspondence to: M. M. Parks, michelle@hi.is Since episodes of unrest may or may not lead to an eruption, learning how to best interpret the geodetic, seismic, and geochemical data in near real time during a volcanic crisis is a key aim for volcanologists seeking to understand such systems. To this end, here we analyze an extended data set of geodetic measurements at Santorini volcano, a caldera-forming volcano in Greece, spanning a 20 year period, to obtain an improved understanding of the long-term deformation history. We also use a recent phase of unrest (2011-2012) as a test case to determine whether or not it may be possible to infer additional information on source parameters and melt supply during this episode, which may be applicable at other volcanoes during similar phases of unrest. We employ a joint inversion technique to convert geodetic measurements into a time series of subsurface volume change and interpret this time series in terms of multiple phases of increasing pressure to a shallow magma chamber surrounded by a viscoelastic shell.
Between 1950 and 2010, the volcano remained quiet, with no documented reports of seismic activity. Although reports of minor inflation were documented in the northern part of the caldera between 1994 and 2000 [Stiros and Chasapis, 2003;Stiros et al., 2010;Saltogianni and Stiros, 2012], other studies [e.g., Papageorgiou et al., 2010;Foumelis et al., 2013] have argued that the observations are not consistent with the inflation of a shallow magma chamber. However, multiple measurements of surface deformation between 1992 and 2010 have confirmed a slow subsidence signal in the southwest region of Nea Kameni [Lagios et al., 2005;Papageorgiou et al., 2010Papageorgiou et al., , 2012Foumelis et al., 2013 andLagios et al., 2013]. The location and velocity of this signal remained relatively stable throughout this period, with an observed average subsidence rate on the order of 5 to 6 mm/yr [Foumelis et al., 2013].
In January 2011, the volcano entered a period of unrest, characterized by the onset of detectable seismicity and caldera-wide uplift [Newman et al., 2012;Parks et al., 2012]. Figure 1 shows the location of volcano-tectonic (VT) earthquakes with local magnitude (M L ) greater than 2.0, which occurred within the caldera between October 2010 and September 2012 [Aristotle University of Thessaloniki, 2005]. The earthquakes are clustered close to the Kameni line-a NE-SW trending normal fault/fracture zone-which is thought to control vent locations for historic dome-forming eruptions on Nea Kameni [Pyle and Elliott, 2006].
The unrest that began in late December 2010/January 2011 marked the first significant activity detected at Santorini since the end of the last eruption in 1950. The results of several geodetic studies covering the period of unrest have already been published (Table 1).
Prior studies have used both GPS and interferometric synthetic aperture radar (InSAR) measurements to model the location of the source of the inflation observed at Santorini between January 2011 and April 2012 (Table 1); however, only one previous study [Foumelis et al., 2013] has attempted to model the earlier subsidence signal observed at Nea Kameni by assuming a Mogi source [Mogi, 1958] and modeling the subsidence with a deflating pressure source at a depth of 800 ± 200 m beneath the southwest corner of Nea Kameni, deflating at (À14 ± 3) × 10 3 m 3 yr À1 .
Modeled locations of the later source of inflation typically agree within~500 m (Table 1), yielding a median longitude, latitude, and depth to source of 25.389°E, 36.426°N, and 4.2 km, respectively, and a median rate of inflation of 10.4 × 10 6 m 3 yr À1 . The one exception to the consensus on source location is the work of Papoutsis et al. [2013], and we suggest that this may derive from their use of data from a single satellite track (descending Envisat track 93). Most recently, Saltogianni et al. [2014] explored whether multiple pressure sources were required to explain the GPS times series and patterns of seismicity throughout the period of unrest.  [Parks et al., 2012] 25.389 36.430 4.4 11.8 GPS [Papoutsis et al., 2013] 25.384 36.429 3.5 12.4 InSAR [Papoutsis et al., 2013] 25.403 36.426 6.3 24.2 InSAR [Foumelis et al., 2013] 25.389 36.429 4.3 9.4 GPS [Foumelis et al., 2013] 25.389 36.424 4.0 10.3 Joint InSAR-GPS [Foumelis et al., 2013] 25.389 36.425 3.8 9.2 InSAR [Lagios et al., 2013] 25.395 36.426 4.5 9.2 GPS [Lagios et al., 2013] 25 While these studies all attempted to determine the source locations and volume changes, there has been no attempt to bring these interpretations together to develop a broader picture of the evolution of the shallow magmatic system throughout the period of unrest. Here we analyze deformation measurements derived from C band (5.65 cm wavelength) and X band (3.11 cm wavelength) satellite synthetic aperture radar (SAR) data acquired between 1993 and 2012 (Table S1 in the supporting information) and continuous global positioning system (cGPS) data from 2010 to 2012, collected from a network of 10 sites, installed on the caldera complex ( Figure 1). We use modeling to explore the causes of the slow subsidence on Nea Kameni and a joint inversion technique to interpret both interferometric synthetic aperture radar (InSAR) and cGPS measurements during the period of unrest as a time series of subsurface volume change.
These previous studies lay the foundations of our current work and allow us to use data from this well-studied system to gain a detailed understanding of the geodetic signals seen, specifically by the following: 1. Using an extended data set and combining geodetic observations from multiple sources (GPS and InSAR) covering two distinct episodes of deformation (localized subsidence from 1993 to 2010 and caldera-wide inflation from 2011 to 2012). 2. Applying a joint inversion technique to simultaneously analyze independent measurements of the surface displacement and provide a time series of volume change. This facilitates the detection of any changes in source characteristics with time. 3. Exploring alternative composite models to explain the slow subsidence signal observed at Nea Kameni (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). 4. Proposing a new viscoelastic model for the shallow magma chamber that explains both the variation in deformation and seismicity during the period of unrest (2011)(2012).
This work has immediate application to volcano monitoring for the improved understanding of magma supply rates and the temporal evolution of deformation source parameters during unrest at other caldera systems.

Interferometric Synthetic Aperture Radar Processing
InSAR is an established technique routinely used to measure centimeter-scale ground movements at volcanoes between repeat passes of a satellite [e.g., Massonnet and Sigmundsson, 2000;Pritchard and Simons, 2004;Hooper et al., 2012]. In this study, we used SAR data from ERS, Envisat, and TerraSAR-X satellites acquired between 1993 and 2012 (Table 2). Interferograms were processed using the Stanford Method for Persistent Scatterers software (StaMPS) small baseline (SB) technique [Hooper et al., 2007;Hooper, 2008]. This method provided improved coherence on the island of Thera and minimized small unwrapping errors observed in the previous two-pass processing [see Parks et al., 2012, Figure 2]. We used a 15 m digital elevation model (DEM) (Figure 1), provided by P. Moore of IntegralGIS, to remove topographic fringes.

Continuous Global Positioning System Processing
Prior to the 2011-2012 unrest, four permanent cGPS stations were operating on Santorini (SNTR, NOMI, KERA, and PKMN). After July 2011, a further six cGPS sites were installed across the Santorini archipelago (MOZI, SANT, MKMN, DSLN, WNRY, and RIBA). Data from each of these 10 cGPS sites were processed with Bernese GPS Software V5.0 [Dach et al., 2007], using the double differenced carrier-phase observables approach and  Table S1 in the supporting information.
Journal of Geophysical Research: Solid Earth

10.1002/2014JB011540
final International Global Navigation Satellite Systems Service products [Dow et al., 2009]. We solved for daily site coordinates and hourly tropospheric parameters using Niell mapping functions [Niell, 1996], with the wet part of the tropospheric delay estimated as a piecewise linear function every 2 h. Horizontal tropospheric gradient parameters were estimated for each station in 24 h intervals, using the "tilting" model [Meindl et al., 2004].
The time series data for each of the cGPS sites were not detrended prior to incorporation into the volumetric inversion (sections 3.3 and 4.2.1). The majority of the sites were only in operation since July 2011 (or later), and we consider it ambiguous to estimate annual and semiannual signals over such short time periods. Instead, we accounted for these residual trends in the volumetric inversion by solving for several nuisance parameters, including a constant offset and velocity for the north, east, and up (NEU) components at each cGPS site.

Joint Inversion of Geodetic Observations for Determining the Temporal Evolution of Source Parameters
We used a joint inversion technique outlined by Biggs et al. [2010] and adapted from Berardino et al. [2002] to convert cGPS and InSAR observations to a subsurface volume change by treating the displacements as arising from a point source pressure variation at depth within an elastic crust. Unlike standard joint inversion techniques, this method provides a time series of the change in source volume through time, which is more useful for interpreting changes in the behavior of a volcano over a selected time period. We used a simple Mogi source [Mogi, 1958] to convert our displacement time series to volume change. We estimate the incremental volume change of the Mogi source V, using equations (7) and (8)   The primary strength of this method is that it enables us to jointly invert InSAR measurements from multiple tracks with different look angles and three component GPS measurements from multiple sites. The incremental volumes are then integrated to provide a time series of volume change [Biggs et al., 2010]. We solved for several different sources of noise or nuisance parameters: those that were temporally and spatially correlated and those that were uncorrelated. The correlated nuisance parameters included a constant offset and velocity for the NEU components at each cGPS site and a constant phase shift in the interferograms. The random noise component was estimated using the relative errors on each observation and a diagonal variance-covariance matrix to weight the inversion [Menke, 1989;Sudhaus and Jónsson, 2009]. The cGPS uncertainties were estimated during the initial Bernese processing, and the error on each interferogram was estimated by calculating the standard deviation beyond a radius of 4 km from the source location [e.g., Biggs et al., 2010]. More detailed information regarding both the inversion technique and error estimation may be found in Biggs et al. [2010].

Period of Subsidence (1993-2010)
The mean line of sight (LOS) velocity maps from 1993 to 2010 (generated using ERS and Envisat SB interferograms) display a slow subsidence signal on the central volcanic island of Nea Kameni (the blue colors

Slow Subsidence (Nea Kameni)
Scenarios that might explain the steady slow subsidence observed on Nea Kameni from 1993 to 2010 include the cooling and contraction of a shallow subsurface magma body [e.g., Foumelis et al., 2013, although the authors suggest that this may result from variations within the hydrothermal system], contraction of historic lava flows, and loading by flows inducing viscoelastic relaxation of the substrate. Since there is no independent evidence for a shallow magma body at this location beneath Nea Kameni, whereas this location coincides with the most voluminous of the historic lava flows (1866-1870) [Nomikou et al., 2014], we limit the modeling and discussion to investigate whether thermal contraction and continued loading from historic lava flows can together plausibly account for the observed signal.

Thermal Contraction of Historic Lava Flows
Postemplacement contraction of lava flows may result in surface deformation by two processes: (1) repacking of clasts and compression of voids within underlying flows or strata and (2) thermal contraction. The loss of void space in underlying lava flows/strata is thought to occur rapidly following emplacement [Stevens et al., 2001], so we do not consider (1) a likely ongoing phenomenon for flows emplaced over 65 years ago like those at Santorini. Thus, we address only (2) below.
The slow subsidence signal observed at Nea Kameni corresponds spatially to the locations of the 1866-1870 and 1939-1941 lava flows (Figure 6a). We first explore whether the observed subsidence on Nea Kameni between 1993 and 2010 may be related to the cooling and contraction of this lava pile by modeling each flow as a slab of thickness d, with temperature T (relative to ambient) and coefficient of thermal expansion α. If the individual flows are emplaced at time t before present, on a semi-infinite medium of thermal diffusivity κ [Carslaw and Jaeger, 1959], then the rate of contraction, ω, is Numerical models have been applied at other volcanoes to determine the rates of lava flow cooling, resulting from multiple thermodynamic processes [e.g., Peck et al., 1977;Shaw et al., 1977;Patrick et al., 2004]. However, models incorporating transient, short-term, posteruptive processes are only relevant for the first 2 years of cooling, after which time thermal conductivity and porosity become the dominant factors governing heat loss [Patrick et al., 2004]. Given the extended period since lava extrusion, we do not see the need to use a more complex model other than the thermal contraction model presented here.
We computed the expected contraction rate as a function of flow thickness (Figure 6b) for the 1866-1870, 1925-1928, and 1939-1941 lava flows, using equation (1) and the input variables: α = 3 × 10 À5 K À1 , κ = 10 À6 m 2 s À1 [Brown and Marco, 1958;Henderson and Henderson, 2009] and T = 1223 K [Pyle and Elliott, 2006]. The time since emplacement (t) is~140 years for the 1866-1870 flow,~80 years for the 1925-1928 flow, and~65 years for the 1939-1941 flow. The difference between the three curves (displayed in Figure 6b) results from the different emplacement times. The average thickness of the 1886-1870 flow is approximately 90 m, which corresponds to~3 mm/yr of contraction. The average thickness of the 1939-1941 flow is 20 m [Nomikou et al., 2014], which results in <1 mm/yr contraction, while the 40 m thick 1925-1928 flow [Nomikou et al., 2014] would correspond to <2 mm/yr of contraction ( Figure 6b). The 1866-1870 flow has the largest continuing contraction rate based on its substantial thickness and could account for~50% of the observed subsidence rate for the period of 1993-2010. Although the combined rates from the three lava flows modeled (1866-1870, 1925-1928, and 1939-1941) would add up to the 6 mm/yr of observed average subsidence, the 2 mm/yr from the 1925-1928 flow cannot contribute to the subsidence observed in the southern portion of the island, as this flow was only extruded in the northeast of Nea Kameni (Figure 6a). The 1939-1941 flows were partially extruded on the northwestern section of the 1866-1870 flows, so may contribute to the observed subsidence in this area but at a much lower rate (<1 mm/yr).
Since the observed subsidence exceeds that due to thermal contraction alone, we now investigate other processes that may contribute to the signal. The most common co-contributing process proposed in prior work to account for subsidence on active volcanoes is loading from historic flows and induced viscoelastic relaxation of the substrate [Ebmeier et al., 2012;Ofeigsson et al., 2011;Fournier et al., 2010].

Loading From Historic Lava Flows and Viscoelastic Relaxation of the Substrate
Since the 1866-1870 flow was the most voluminous of the historic lava flows extruded on Nea Kameni [Pyle and Elliott, 2006;Nomikou et al., 2014], we consider that loading from this flow and viscoelastic relaxation of the substrate may account for at least part of the observed subsidence signal. To test this, we modeled a cylindrical load placed on an elastic lid of thickness H, on a viscoelastic half-space [Nakiboglu and Lambeck,  England et al., 2013]. The load has a radius of 1 km, a height of 100 m, and a density of 2500 kg/m 3 . The elastic lid has a density of 2800 kg/m 3 and a Poisson's ratio of 0.25. We use a shear modulus of 5 GPa, and it is assumed that loading commenced~140 years ago. Our model (Figure 7) suggests that >60% of the observed rate of subsidence at Nea Kameni may be induced by the loading of a 100 m thick lava flow, assuming a elastic lid thickness between 2 and 3 km. The processes of cooling and contraction and loading are not mutually exclusive. Given that the location of the signal agrees closely with the extent of the 1866-1870 lava flow and significant seafloor subsidence occurred in the area of the extruded flows following emplacement, we suggest that a combination of thermal contraction and load-induced relaxation is sufficient to account for the observed subsidence, and we see no need to infer a shallow subsurface deformation source for which we have no additional evidence. This is the first study to provide an explanation for the long-term subsidence observed at Nea Kameni, whereby the geodetic observations are supported by a combination of modeling results and bathymetric data. Similar mechanisms have been proposed by Briole et al. [1997] and Stevens et al. [2001] to explain posteruptive deformation at Mount Etna (Italy) and by Ofeigsson et al. [2011] at Hekla volcano (Iceland).

Caldera Wide Uplift 4.2.1. Joint Inversion
We have used the Joint Inversion of Geodetic Observations technique (section 3.3) to convert deformation measurements derived from the 10 cGPS sites (Figure 1) and 199 SB interferograms (Table S1 in the supporting information) covering the period of unrest, into subsurface volume change. Prior to running the inversion, the SB interferograms were resampled to a grid with a spacing of 200 m up to an inner radius of 4 km from the source, with a coarser spacing of 400 m outside this area. This increased the number of sample points close to the source and reduced the number of sample points in the inversion that are more likely to be affected by atmospheric delay on the eastern edge of Thera. For each data set (InSAR-only, GPS-only, and combined InSAR-GPS), we searched a three-dimensional grid of 200 m spacing in x, y, and z to determine the optimal source location that minimized the misfit between the data and the model. The error on the total volume estimates for each of the inversions is the cumulative sum of the weighted root-mean-square (WRMS) residual for individual observations (computed during the inversion). The optimum depth was located using the minimum WRMS derived during the inversion. The bounds on the depth estimate were derived using the volume change error estimates and the trade-off plot between the depth and the volume. This is the first application of a simultaneous inversion technique at Santorini volcano, capable of generating a time series of subsurface volume change. The time series of cumulative volume change for the InSAR-only, GPS-only, and joint InSAR-GPS data sets is shown in Figure 8 and the best fit source locations in Figure 9. The combined InSAR, GPS, and joint inversion results suggest that the best fitting point source has an effective  (Table 1 and Figures 8, 9, and 1). The depth to the source is in good agreement with those reported in previous studies (Table 1), as is the geographic location (longitude: 25.388-25.390, latitude: 36.428) and the rate of volume change within the shallow chamber~11-18 × 10 6 m 3 /yr. The estimated total volume change from the InSAR-only inversion (from January 2011 to April 2012) is (15 ± 1) × 10 6 m 3 . The volume change between January and March 2011 (data gap in InSAR measurements) is estimated using the average rate displayed in parentheses in Table 1. The total volumes computed from the GPS-only and joint inversions over the same time period (January 2011 to April 2012) are slightly higher (19 ± 4) × 10 6 m 3 . These differences may result from the varying sensitivity of InSAR and GPS measurements to vertical and horizontal movements or may be suggestive of a higher horizontal to vertical expansion of the source, as suggested by Foumelis et al. [2013]. The misfit between the input observations and model is displayed in Figures S2 and S3 in the supporting information.
Figures 9a-9c display the results of a grid search for source location using the RMS misfit-a measure of goodness of fit rather than the formal error discussed in section 3.3. The red filled circles represent the geographic location of the minimum WRMS error, representing the optimum source location (e.g., Figures 9d-9f). Surprisingly, the InSAR data (Figure 9a) alone places a weaker constraint on the source location than the GPS data ( Figure 9b). The key advantages of InSAR are usually the greater spatial resolution and ability to measure deformation without ground-based equipment. For the 2010-2012 period of unrest at Santorini, the joint inversion of InSAR and GPS data produces a result that is identical to the GPS-only inversion. In this case, the 10 cGPS sites provide adequate spatial coverage given the simplicity of the spatial pattern of deformation, and the InSAR coverage is restricted by the configuration of the islands. As a result, the three-dimensional daily solutions for the 10-station cGPS network provide better constraints on the time series of displacement than the InSAR.
In contrast, the subsidence from 1993 to 2010 occurs during a time period with no suitable cGPS coverage, a situation common to the majority of the world's volcanoes, even now. Furthermore, the spatial pattern of deformation is limited to Nea Kameni and would not have been covered, even if the current 10-station GPS network had been deployed since 1993. In several past eruptions, ground-based stations have been damaged or rendered inoperable. For example, in the Eyjafjallajökull 2010 eruption, a semipermanent GPS installation was carried away in a jökulhlaup, and several campaign benchmarks were buried by flood deposits (S. Hreinsdóttir, personal communication, 2014). It is also common for unrest to evolve into more complex spatial patterns, for example, recent activity at Bárðarbunga volcano Iceland, involved simultaneous dyke inflation, and caldera subsidence [Sigmundsson et al., 2014]. Therefore, in the long term, a combination of both GPS and InSAR remains the optimal strategy for measuring ground deformation at restless volcanoes. It is also useful to compare the information that each technique can provide in isolation.
An interesting temporal variation in volume change of the shallow chamber, throughout the period of unrest is observed in Figure 8. To determine whether this variation may be related to instability in the source location, we ran a series of tests subdividing the period of unrest into a series of time intervals and rerunning the inversion for the different time periods. Our results show no evidence for any significant source migration (in depth or geographic location) throughout the period of unrest. This agrees with the single source modeling results obtained from Saltogianni et al. [2014]. However, these authors preferred a combination of seismic and aseismic, multiple double sources (deep and shallow) at different locations beneath the caldera to explain the variations in deformation and seismicity observations. Interestingly, they found no rapid change in direction or center of origin of the displacement vectors, and their Coulomb stress models appear to be in disagreement with an earlier study undertaken by Feuillet [2013].
Although in reality, the magmatic plumbing system beneath Santorini volcano is likely to be significantly more complex than either a single or double Mogi source, and a complete analysis of stress distribution and deformation at caldera volcanoes might be amenable to finite element modeling (such as that described by Geyer and Gottsmann [2010] and Hickey et al. [2013]); as a simple approximation, we model the variation in volume change within a single shallow source, believed to be responsible for the bulk of the deformation at Santorini observed during the period of unrest. Assuming a single shallow source of deformation, the time series plot (Figure 8) suggests that the shallow (~4 km) chamber may have undergone two episodes of increasing pressure during the period of unrest. A steady rate of volume change was observed between December 2010 and August 2011 (first pulse). We interpret this first pulse as an increase in pressure within the shallow magma chamber, associated with an initial magmatic intrusion. After August 2011, there is a slight reduction in the rate of cumulative volume change, followed by an additional increase in volume change rate from October 2011 to April 2012. Each pulse corresponds to a volume change of~1 × 10 7 m 3 (Figure 8). We suggest that these pulses may relate to separate batches of melt being supplied to the shallow chamber.

Viscoelastic Modeling
To investigate the likely cause of these signals, we looked first at the variable rate of melt supply to a shallow magma reservoir. This was achieved by modeling linear pressure changes within a sphere, surrounded by an elastic medium. However, this model could not replicate the distinctive rollover observed in the volume change curve at the end of each pulse (e.g., around August 2011 and April 2012) identified on Figure 8, without invoking a more complex melt supply history. This could be explained by melt supply that ramps up at the start of each pulse and then tapers at the end, but both pulses must ramp up and taper with the same decay rate. This type of exponential decay often indicates stress relaxation through viscoelastic processes. In the event of a magmatic intrusion, heat is transferred to adjacent wall rocks, which may result in the formation of a ductile aureole (also referred to as a viscoelastic shell). This shell promotes a region of stress relaxation in the crust immediately adjacent to the magma reservoir, which effectively acts as a filter, modifying the stress transfer between the walls of the magma chamber, and the surrounding crust. Field evidence and thermal models demonstrate that thermal aureoles may develop in the host rock surrounding magma intrusions [Hodge, 1974;Gasparini et al., 1981]. There is strong evidence for the existence of a thermal or metamorphic aureole beneath the Kameni islands, both from basement rock inclusions in Kameni lavas [Nicholls, 1971] and from the carbon-isotopic compositions of Kameni CO 2 emissions [e.g., Tassi et al., 2013].
Although this model is more complex than a Mogi source in an elastic half-space, the pressure history that it requires to fit our data is simpler, and it explains the symmetry of the signal. Below we derive a new model based on this viscoelastic concept (following the approach of Segall [2010]) to explore whether this scenario might cause the observed variation in the volume change curve displayed in Figure 8. We modeled the observed temporal variation in volume change (derived from the joint inversion), associated with both pulses, in terms of the ductile deformation of a concentric viscoelastic shell surrounding a shallow magma chamber, driven by two separate increases in pressure. Our viscoelastic model consists of a magma chamber located at a depth, d, beneath the surface, with a radius R 1 . The magma chamber is surrounded by a viscoelastic shell of radius of R 2 -R 1 and rigidity μ, embedded in an elastic half-space.
Several studies at other active volcanoes have used viscoelastic modeling to better constrain source parameters responsible for observed deformation [e.g., Bonafede et al., 1986;Dragoni and Magnanensi, 1989;Newman et al., 2001;Hickey et al., 2013]. Bonafede et al. [1986] were the first to explore the effect of the anelastic properties of the Earth's crust in relation to volcanic deformation. These authors presented analytical expressions to model the deformation produced by a pressure source embedded in a viscoelastic half-space. Dragoni and Magnanensi [1989] expanded on this study by considering a magma chamber surrounded by a viscoelastic shell embedded in an infinite medium (full space). The analytical solutions of Dragoni and Magnanensi [1989] are very useful in determining the variation in the shape of the deformation pattern over time; however, the output may not be directly compared with GPS or InSAR measurements, as the elastic full-space solution would underestimate the observed deformation (by a factor of 3 if using a Poisson's ratio of 0.25) due to the exclusion of the free surface [Newman et al., 2001]. More recent studies [Hickey et al., 2013;Masterlark et al., 2010;Newman et al., 2006] have used a finite element approach to model viscoelastic deformation within an elastic half-space.
Although the finite element method enables modeling of magma chambers with increased complexity, these models are often time consuming, computationally intensive, and difficult to validate. Segall [2010] presented two analytical solutions to model vertical displacements due to a pressurized sphere in an elastic half-space, based on a modification of Dragoni and Magnanensi's [1989] original full-space solution. The first of Segall's [2010] solutions models the pressure change within the magma chamber as an instantaneous step increase (Heavyside function), whereas the second models the pressure change as a gradual increase followed by an exponential decay [Segall, 2010, pp. 243-248]. Although we applied both solutions in an attempt to model the variation in volume change/deformation observed at Santorini between December 2010 and October 2012, the solutions did not provide a good fit to the observations. This suggests that the pressure increase within the shallow magma chamber at Santorini was not instantaneous nor did it gradually decay with time. In light of this, we derived a new analytical solution to model the variation in volume change within a shallow magma chamber surrounded by a viscoelastic shell, considering an alternate pressure history.
We suggest that the observed variations in seismicity and the rate of volume change of the shallow magma chamber ( Figure 10) are intrinsically connected and result from two separate phases of increasing pressure within the shallow reservoir. To test this idea, we followed the same equations presented in Segall [2010, section 7.6] but modify equation (7.96) (p. 243) to use an alternate pressure function. We then transform the final equation so the output is given in terms of volume change rather than deformation, using an equivalence relation to the time dependent Mogi model.
Our pressure function corresponds to two consecutive ramp functions-this represents a linearly increasing pressure to the shallow chamber between times t 0 and t 1 , a fixed pressure from t 1 to t 2 , linearly increasing pressure between t 2 and t 3 and a constant pressure from t 3 onward ( Figure 10). Several studies support a linear ramp-like pressure increase, and a very similar deformation history has been observed at multiple volcanoes including Campi Flegrei [Chiodini et al., 2012], Mount Etna [Heap et al., 2011], Mount Pinatubo [Smith and Kilburn, 2010], and Rabaul [Saunders, 2001]. The pressure change within the shallow magma chamber may be represented in terms of volume change [McTigue, 1987]. The time-dependent volume change is given by equation (2):  Figure 4). It should be noted that these CO 2 values are derived from a small area on the summit of Nea Kameni and can only be seen as suggestive of trends in the broader flux of CO 2 through the Santorini volcanic system. Temporal variation in the gas chemistry from fumaroles, reported by Tassi et al. [2013], is also consistent with a period of magma intrusion to the shallow system from early 2011 to early 2012. The green lines represent the normalized pressure model, displaying a linearly increasing pressure from t 0 to t 1 and from t 2 to t 3 . In between t 1 and t 2 and after t 3 , the pressure within the shallow magma chamber is held constant.
Journal of Geophysical Research: Solid Earth 10.1002/2014JB011540 and The pressure function is given by equation (5), and the full derivation may be found in Text S4 in the supporting information.
The predicted volume change is fitted to the data, weighted by the errors between consecutive observations, using a least squares minimization routine in Mathematica [Wolfram, 2013]. We solve for the ratio of the outer to inner radius (R 2 /R 1 ), the characteristic relaxation time (t R ), the shut-off times (t 1 and t 3 ), the onset time for the second pressure pulse (t 2 ) and a scaling coefficient, which is a combination of the pressure increase within the magma chamber and the chamber radius (R 1 ). The best fit parameters minimizing the misfit between the volume change derived from the joint inversion and the viscoelastic model are displayed in Table 3. The effective viscosity of the shell was calculated using the relaxation time scale (derived during the modeling) and equation (4), specifying both the shear modulus (5 GPa) and Poisson's ratio (0.25).
The best fit model provides a good match to the joint inversion result, and the derived shell viscosity (4.7 × 10 15 Pa s) and R 2 /R 1 ratio (1.73) are within the expected range, considering values reported from other studies [e.g., Newman et al., 2001, and references therein]. The pressure and radius terms in the Mogi approximation are coupled [Mogi, 1958], and it is not possible to individually determine either parameter. However, a trade-off curve between pressure change and radius may be used to determine a range of physically viable parameters. The computed trade-off curve is displayed in Figure S5 in the supporting information (for example, if the overpressure in the shallow chamber was increased to 10 MPa, this would dictate a chamber radius (R 1 ) of~1200 m and an outer radius (R 2 ) of~2000 m), given the R 2 /R 1 ratio of 1.73). The derived shut-off times (t 1 and t 3 ) for the pressure pulses also provide a good fit to the observed cumulative seismicity (Figure 10), both coinciding with a lull in earthquake activity, following a spike in cumulative seismic moment. Between t 1 and t 2 and after t 3 , the pressure within the magma chamber is held constant (this may be achieved via a number of processes-one of which could involve vesiculation and concurrent degassing), but the cumulative volume change continues to increase at an exponentially decaying rate.
This viscoelastic response explains the lag between the cessation in earthquake activity and the flattening of the cumulative volume change curve shown in Figures 8 and 10, i.e., the VT earthquakes cease following the removal of the source of increasing pressure; however, the viscoelastic shell surrounding the magma chamber acts as a buffer, resulting in a delayed decline in deformation driven by the pressure shut off.
The timing of the second pressure shut off (t 3 ) corresponds to a significant earthquake swarm at the end of January 2012, which was followed by an almost complete cessation in microseismic activity in the center of the caldera. The swarm also corresponded to a change in CO 2 emissions on the summit of Nea Kameni . Although CO 2 soil gas measurements were undertaken at regular intervals throughout the period of unrest, the peak of soil-CO 2 outgassing was not observed until a year after the onset of unrest, following the January 2012 swarm   (Figure 10). This observation may also be consistent with brittle creep laboratory experiments, whereby axial strain and rock damage (cracking) are seen to accelerate prior to brittle rock failure (formation of fractures) [Heap et al., 2011]. The combined deformation, seismicity, and gas observations at Santorini may in fact be indicative of the final acceleration stage in these experiments, during which crack propagation and coalescence produce fractures providing new gas 11.56 ± 0.14 × 10 6 9.73 ± 0.10 × 10 6 a Parameter η is the viscosity of the viscoelastic shell surrounding the magma chamber. R 2 /R 1 is the ratio of the outer-inner radius. Parameters t 0 and t 2 are the onset times of the increasing pressure pulses. The shut-off times are t 1 and t 3 . Parameter t R is the relaxation time scale.
Journal of Geophysical Research: Solid Earth 10.1002/2014JB011540 migration pathways to the surface [Kilburn, 2003;Wright et al., 2009]. Konstantinou et al. [2013] postulated that gas migration may be responsible for the observed seismicity, highlighting the fact that the seismicity and source of deformation were not co-located. Konstantinou et al. [2013] suggested that the intrusion north of Nea Kameni resulted in an increase/redistribution of stress within the crust, triggering earthquakes along preexisting zones of weakness-e.g., Kameni line and possible ring faults. In this scenario, fracture propagation may be triggered by a combination of processes including increasing pore pressure [e.g., Martens and White, 2013] and stress corrosion [Kilburn and Voight, 1998, and references therein]. Temporal variations in V p /V s ratios inferred by Konstantinou et al. [2013] correspond broadly to our two distinct pressure pulses displayed in Figures 8 and 10. During February-August 2011, the authors observed lower V p /V s values (<1.8), indicative of gas rather than melt-filled cracks, followed by an increase between September 2011 and December 2011 and subsequent decrease in February 2012.
The distinctive rollover observed in the volume change curve at the end of each pulse has been observed in the deformation time series at other caldera volcanoes (e.g., Long Valley [Newman et al., 2006]). Therefore, the viscoelastic model presented here has immediate potential application to other caldera-forming systems during episodes of unrest to help better understand the source mechanism and temporal relationship between the deformation and seismicity observations.  Table 4). Several similarities can be drawn between the recent unrest at Santorini and the episodes of minor uplift and seismicity that occurred at Campi Flegrei between 1989 and 2012 [Chiodini et al., 2012] (Table 4). These include (i) the inferred existence of multiple reservoirs (shallow and deep) in both settings [Cottrell et al., 1999;Gottsmann et al., 2006], (ii) the observation that seismicity and deformation are not spatially coincident [D'Auria et al., 2011;Konstantinou et al., 2013], and (iii) the observed delay between the onset of unrest and an increase in gas emissions [Chiodini et al., 2012;Aiuppa et al., 2013]. Several authors [e.g., D' Auria et al., 2011;Chiodini et al., 2012] have suggested that the recent phases of minor uplift and seismicity at Campi Flegrei are related to the repeated transfer of magmatic fluids to a shallow hydrothermal reservoir. A similar cause of unrest has been inferred for the pulses of uplift, accompanied by earthquake swarms and increased CO 2 emissions that occurred at Long Valley between 1978 and 1998 ( Table 4). The role of (hydrous) fluid migration in the recent unrest at Santorini remains unknown. Although both geodetic and gas-chemical data sets can be interpreted in terms of the supply of new magma to a shallow reservoir beneath Santorini [Parks et al., 2012Tassi et al., 2013], and this model is consistent with viscoelastic modeling (section 4.2.2), a detailed gravity survey would be required to constrain the mass changes associated with this phase of unrest.

Conclusions
Prior to the start of 2011, a slow subsidence signal persisted at Nea Kameni. Our modeling suggests that this subsidence could be accounted for by the combination of thermal contraction and the viscoelastic relaxation of the substrate associated with recent lava flows. There is reasonable correspondence between the signal and the location of the 1866-1870 lava flow, the most voluminous historical flow, and our modeling suggests that the observed subsidence rate between 1993 and 2010 may be accounted for by a combination of both processes. This is the first study to explore alternative mechanisms for the measured subsidence at Nea Kameni and provide an explanation, supported by both modeling and independent observations.
A new phase of unrest commenced at Santorini volcano in January 2011, characterized by the onset of microseismic activity focused along the Kameni line, as well as caldera wide uplift. We have used a joint inversion of InSAR and GPS data to determine the temporal evolution of source parameters by converting both InSAR and GPS observations into a time series of subsurface volume change. The best fitting parameters derived from the combined InSAR-GPS inversions are in good agreement with most other studies and suggest that the point pressure source is centered north of Nea Kameni at a depth of~4 km beneath the surface and has an effective volume of~14-23 × 10 6 m 3 . The time series of subsurface volume change has identified two separate pressure pulses during the period of unrest. A lag response is also identified between the lulls in microseismicity within the caldera and the exponential decay observed in the rate of subsurface Journal of Geophysical Research: Solid Earth 10.1002/2014JB011540 volume change (several months later). We propose that the two pulses observed are associated with episodes of magma supply to a shallow reservoir and that the resulting pressure increase within this reservoir is responsible for the VT earthquake swarms focused along the preexisting zone of weakness-the Kameni line. We have modified the model of the shallow magma chamber beneath Santorini volcano by considering the magma chamber as a sphere inside a viscoelastic shell (representing a ductile aureole) in order to account for the lag between the seismicity and the shape of the volume change time series. This characteristic pattern has also been observed at other caldera volcanoes (e.g., Long Valley); thus, the viscoelastic model presented here has immediate application to other caldera-forming systems.
It is currently not possible to determine when or if a new pulse of magma will be supplied to the existing shallow chamber, or indeed if the next pulse will bypass this chamber completely and be intruded at an alternate location. However, this study does provide additional insight into the subsurface rheology and the temporal evolution of magma supply during periods of unrest at Santorini. The next phase of unrest will provide important constraints on the time period between these rapid pulses. Continuous deformation, seismicity, and gas monitoring are essential to identify the onset of renewed unrest. IntegralGIS for providing the 15 m DEM used in the InSAR processing. This work was supported by the National Environmental Research Council (NERC) through an urgency grant NE/J011436/1, a NERC small grant NE/I015760/1, and through COMET. M. P. was supported through an NCEO studentship and the Leverhulme trust. TSX data were acquired through DLR proposal GEO1228. SAR data used in this study may be obtained through the relevant satellite agencies. Continuous GPS data are made available through UNAVCO.