Time-Scales of Inter-Eruptive Volcano Uplift Signals: Three Sisters Volcanic Center, Oregon (United States)

A classical inflation-eruption-deflation cycle of a volcano is useful to conceptualize the time-evolving deformation of volcanic systems. Such a model predicts accelerated uplift during pre-eruptive periods, followed by subsidence during the co-eruptive stage. Some volcanoes show puzzling persistent uplift signals with minor or no other geophysical or geochemical variations, which are difficult to interpret. Such temporal behaviors are usually observed in large calderas (e.g., Yellowstone, Long Valley, Campi Flegrei, Rabaul), but less commonly for stratovolcanoes. Volcano deformation needs to be better understood during inter-eruptive stages, to assess its value as a tool for forecasting eruptions and to understand the processes governing the unrest behavior. Here, we analyze inter-eruptive uplift signals at Three Sisters, a complex stratovolcano in Oregon (United States), which in recent decades shows persistent inter-eruptive uplift signals without associated eruptive activity. Using a Bayesian inversion method, we re-assessed the source characteristics (magmatic system geometry and location) and its uncertainties. Furthermore, we evaluate the most recent evolution of the surface deformation signals combining both GPS and InSAR data through a new non-subjective linear regularization inversion procedure to estimate the 26 years-long time-series. Our results constrain the onset of the Three Sisters volcano inflation to be between October 1998 and August 1999. In the absence of new magmatic inputs, we estimate a continuous uplift signal, at diminishing but detectable rates, to last for few decades. The observed uplift decay observed at Three Sisters is consistent with a viscoelastic response of the crust, with viscosity of ∼10¹⁸ Pa s around a magmatic source with a pressure change which increases in finite time to a constant value. Finally, we compare Three Sisters volcano time series with historical uplift at different volcanic systems. Proper modeling of scaled inflation time series indicates a unique and well-defined exponential decay in temporal behavior. Such evidence supports that this common temporal evolution of uplift rates could be a potential indicator of a rather reduced set of physical processes behind inter-eruptive uplift signals.

1 Introduction

Many volcanoes follow a common deformation pattern consisting of uplift during inter-eruptive periods and subsidence in co-eruptive stages, occasionally interrupted by periods of quiescence or subsidence. Some other volcanoes do not however exhibit this simple behavior (Biggs and Pritchard, 2017). Part of them show puzzling non-steady persistent uplift signals that can last from days to years with minor or no other geophysical or geochemical variations, which are difficult to interpret. Therefore, uplift during inter-eruptive episodes cannot be only interpreted as a pre-eruptive precursory indicator. Such temporal behavior is usually observed in large calderas (e.g., Yellowstone, Long Valley, Campi Flegrei, Rabaul), but less commonly for stratovolcanoes.

An important goal in volcano eruption forecasting is to find how the deformation time-series can distinguish among physical processes, especially during inter-eruptive periods leading to a pre-eruptive scenario. The latter are characterized by new injections of magma/increment of volatiles, viscoelastic relaxation of the media, or a mixing of different coeval processes. Therefore, we must constrain what controls usually long-lived or persistent uplift at volcanic centers. Le Mével et al. (2015) show that the temporal evolution of deformation surprisingly follows the same pattern for different volcanic systems at specific analyzed periods (Yellowstone, Long Valley, Laguna de Maule and Three Sisters). This is consistent with the hypothesis that similar processes may be at work in similar volcanoes. After this stage, these volcanoes presented an eventual pause and/or change to subsidence (related to seismic events and/or hydrothermal changes), but did not produce an eruption.

Stratovolcano behavior contrasts with better documented restless calderas (Acocella et al., 2015; Galetto et al., 2017). Calderas usually show long repose periods between large eruptions, but without quiescence (e.g., resurgent volcanism, subsidence, multiple pulses of uplift) (Biggs and Pritchard, 2017). The Three Sisters complex volcano is a good example of a system with long lasting monotonous inter-eruptive uplift without associated eruptive activity or significant seismicity and might correspond with different magmatic system behaviors.

In this work, we studied the decadal deformation time-series of the Three Sisters volcano in order to understand the processes governing its unrest behavior and to find out whether it is still inflating. For this purpose, we needed to perform a consistent analysis of the volume change time-series underlying the Three Sisters uplift signal, which is a challenging task. We analyze available continuous GPS data since 2001 and multiple satellite interferometric data spanning the 1993–2020 period. We proposed combining multiple geodetic data-sets using an improved linear regularization method based on Truncated Singular Value Decomposition (TSVD) (González et al., 2013) to find an optimal regularization criteria for Three Sister data-set combination. We obtained a seamless, continuous time series of volume change (and its uncertainties) with which to rigorously assess changes over the 26-years studied period. Finally, we compared the Three Sisters temporal behavior to other well-known examples of uplifting volcanoes to understand 1) whether a variety or not of physical mechanism are at work behind deformation and, if so, 2) if uplift time-scales are informative of whether a certain volcano is on a late or early stage of the inter-eruptive period.

2 Background

The Three Sisters Volcanic Complex forms a N-S chain of stratovolcanoes located at ${44.1}^{\circ }$N in the central Oregon Cascades (Figure 1), an active volcanic arc produced by the subduction of the Juan de Fuca and Gorda plates beneath the North American plate. In addition to this convergent motion, there is an oblique relative plate motion and northward push of the Sierra Nevada-Great Basin microplate, favoring a N-S stress orientation of the vents within the Oregon Cascades (Mccaffrey et al., 2007). South Sisters is near the propagating tip of a crustal melting anomaly westward across Oregon, progressing since the mid-Miocene, going through the Cascades in the Quaternary (Fierstein et al., 2011). All these circumstances influence on the eruptive history of Three Sisters.

FIGURE 1

FIGURE 1. (A) Shaded relief map for the Cascade range, with locations of representative Cascade volcanoes. Three sisters volcano is highlighted with black squares. (B) Regional map of Three Sisters volcano complex. The continuous GPS station HUSB, located at ∼5 km west of South Sister volcano, is shown as a purple circle with black outline.

The Three Sister area includes shields, composite volcanoes and cinder cones, with basaltic to rhyolitic volcanism. The three eponymous volcanoes are progressively younger from north to south and exhibit little family resemblance (Hildreth et al., 2012). North Sister is a monotonously mafic edifice created 120 ka ago, formed by long-lived effusive volcanism (Schmidt and Grunder, 2011); Middle Sister is an andesite-basalt-dacite cone constructed between 48 ka and 14 ka and South Sister is a bimodal rhyolitic-intermediate edifice built between 50 ka and 2 ka, both with histories of explosive volcanism (Scott et al., 2001). However, most of the volcanic activity is identified with mafic shields and cones around the major composite volcanoes (Hildreth, 2007). Geochemical anomalies suggest that episodes of intrusion may be more frequent in the Three Sisters area than the age of eruptive vents would involve (Evans et al., 2004). The most recent eruptions were rhyolitic close to South Sisters, 2000 years ago (Hildreth et al., 2012). There is a potential volcanic hazards threat if future eruptions are similar to South Sister’s recent past. Tephra fallout might accumulate to 1–2 cm thick in the Bend area, and small-volume lahars and pyroclastic flows could pose a hazard to nearby areas (Sherrod et al., 2004).

Before 2001, Three Sisters was considered a dormant volcano. Nonetheless, ERS-1/2 satellite interferometric synthetic aperture data (InSAR) analysis from 1992 to 2000 revealed active uplifting located 6 km west of South Sister (Wicks, 2002). Geodetic information for Three Sisters volcanic center has been accumulating since this discovery. Nowadays, deformation is continuously monitored through Leveling surveys since 2002, Continuous GPS (CGPS), GPS campaign since 2001 and Semi Permanent GPS (SPGPS) since 2009 (Dzurisin et al., 2009; Dzurisin et al., 2017). The uplift at Three Sisters has been aseismic except for a swarm of ∼300 small earthquakes (M_max = 1.9) in March 2004 (Moran, 2004; Dzurisin et al., 2009). Previous studies (Wicks, 2002; Dzurisin et al., 2006; Dzurisin et al., 2009; Riddick and Schmidt, 2011) show evidence that observed uplift can be described by a spherical point source within an homogeneous isotropic elastic half-space. Nevertheless, deformation source geometry is non-unique and sources as horizontal crack, vertical prolate spheroid, and sill-like have been proposed at Three Sisters to fit geodetic data. Interpretation of the temporal evolution of InSAR, leveling and GPS data suggests the beginning of deformation in late 1997 or 1998, with a maximum uplifting rate of 3‐5 cm/year during 1998–2001. Microgravity data collected between 2002 and 2009 show no significant change in the mass flux across the deforming area (Zurek et al., 2012). No studies have been published about the uplift evolution over the last decade. The uplift process was still on-going in January 2020, when this manuscript was prepared.

3 Geodetic Data Processing

We aim to extend the detailed uplift history at Three Sisters already mentioned above to the present (2020) by using the available CGPS and InSAR data-sets.

3.1 Continuous GPS

In May 2001, the U.S Geological Survey (USGS), in collaboration with the U.S. Forest Service, installed a continuous GPS station (HUSB) near the actively deforming area. It was strategically installed at a location approximately ∼2 km northwest of the detected uplift center. HUSB is part of the USGS Pacific Northwest Network, so it is automatically processed to obtain daily coordinates. No other regional and local continuous GPS station falls within the deformation area. Hence, the HUSB time-series is particularly important to understand the surface deformation time-scales at Three Sisters.

Daily GPS data (coordinates and their uncertainties) are analyzed by the USGS using GIPSY/OASIS II software. Common-mode daily biases are estimated and removed using QOCA (Dzurisin et al., 2009). Three Sisters is located near the actively deforming Cascadia margin, so any geodetic data and coordinates must consider the wider regional deformation patterns. The motion of a background steady rigid-body with a rotation pole situated near the eastern limit of Oregon must therefore be removed from the time-series, and a correction for predicted horizontal tectonic motion should be applied. Here, we remove a linear trend of 4.29 mm/yr for the North component and 1.50 mm/yr for the East component. This model prediction represents an update and improved version of the horizontal displacements at HUSB (Dzurisin et al., 2006; Dzurisin et al., 2009; M. Lisowsky personal communication; Cascades Volcano Observatory, 2017).

Figure 2 shows the resulting GPS displacements between July 2001 and January 2020. CGPS data reveals several gaps that occurred due to snowfall in the winter seasons. Furthermore, CGPS shows a gap and a posterior data offset during August 2017-August 2018. USGS data site reports some readjustment of HUSB permanent station during this period and these could explain some of the gaps and offsets in the time series.

FIGURE 2

FIGURE 2. North (U_North), east (U_East) and vertical (U_up) components of the continuous GPS displacements at HUSB.

3.2 Interferometric Synthetic Aperture Radar (InSAR)

Our InSAR data set includes 85 interferometric pairs, with temporal baselines from 35 to 2,894 days, from four satellite missions (ERS, ENVISAT, ALOS-1 and Sentinel-1). ERS and ENVISAT SAR images were acquired during summer and fall between 1993 and 2010 (descending orbits, tracks 113 and 385; ERS and ENVISAT look angles 20.2^° and 19.8^°, respectively). We used 51 interferograms processed with the ROI PAC software (Rosen et al., 2004) and unwrapped using SNAPHU (Chen and Zebker, 2002), with perpendicular baselines up to 500 m, as explained by Riddick and Schmidt (2011).

To improve the temporal coverage of InSAR observations, we also analyzed data from the ALOS-1 and Sentinel-1 SAR data missions. The mean line-of-sight velocity of ALOS-1 data (path 219, ascending orbit, look angle 38.7^°) was obtained during January 2007- March 2011. Most individual interferograms in the Cascades range show poor coherence because of vegetation and seasonal snow coverage, hence we also processed 4 Sentinel-1 summer-to-summer and summer-to-late spring interferometric pairs, between September 2015 and May 2018, for descending (path 115, look angle 39.8^°) and ascending (path 137, look angle 38.8^°) orbits. To provide deformation data during the GPS gap mentioned above, two Sentinel-1 interferometric pairs cover that period. We used JPL InSAR Scientific Computing Environment (ISCE) software (Rosen et al., 2012), processing Level-0 raw ERS-1 and ALOS-1, and SLC-level Sentinel-1 radar data. All interferograms were corrected for orbital and topographic contributions using precise orbit information and the SRTM digital elevation model (Farr et al., 2007). We also reprocessed a highly coherent interferometric pair for ERS-1 track 365 (descending orbit, corresponding to August 1997 - September 2000). This interferogram was essential to further re-evaluate the magmatic source location and constrain its uncertainties.

L-band data (wavelength ∼24 cm) from ALOS-1 were very useful to avoid decorrelation owing to the vegetation encompassing the Three Sisters area. Although the LOS deformation rate from 2006 to 2010 is small (about 6–8 mm/yr) making it difficult for a single L-band interferogram to detect the deformation signal (Riddick and Schmidt, 2011), a cumulative LOS deformation time-series can detect such changes in rate. The corresponding time-series was processed by StaMPS Version 3.3.b1 to study the surface deformation, applying the Small Baselines method for 30 interferometric pairs (Hooper et al., 2012; Bekaert et al., 2017). Small Baselines method minimizes decorrelation in natural terrains. So, it is an appropriate method for the Three Sisters area, which lacks man-made structures and hence offers few dominant persistent scatterers.

Due to the small deformation rate (5.8 mm/year for the period June 2015 - January 2020) and low signal-noise ratio in the Cascades, geodetic data must be analyzed carefully (Poland et al., 2017). Following this recommendation, we consider a 1 cm standard deviation for neighboring pixels in all interferograms. Only four good quality interferometric pairs were used for the Sentinel-1 observation period. Adding more interferograms did not significantly improve the analysis of the volume change time-series. Moreover, the analysis of a 6-years long Sentinel-1 dataset in Turkey indicates that surface displacement rate uncertainties are mostly dominated by length of observation, rather than larger numbers of available interferograms (Weiss et al., 2020). Hence, we consider Sentinel-1 summer-to-summer and summer-to-late spring InSAR data to avoid decorrelation due to snow coverage, and to fill a noticeable GPS time-series gap. However, the deformation rates could be reexamined in future, using longer Sentinel-1 datasets.

Figure 3 represents the mean LOS velocity (mm/year) for the ascending path 219 ALOS-1 from January 2007 to March 2011. Due to the high signal-to-noise ratio, the StaMPS LOS velocity results were noisy and we post-processed them to reduce undesirable oscillations of non-volcanic origin. We applied a band pass filter to retain spatial deformation signals between 10 and 0.8 km, using a median filter (GMT blockmedian). Although close to the signal-to-noise ratio value, results indicate a 6 km circular uplift pattern west of South Sister with a mean LOS velocity of ∼5–10 mm/year. This mean velocity is consistent with a value obtained for the HUSB CGPS station during the same period (5.2 mm/year).

FIGURE 3

FIGURE 3. The linear surface deformation LOS velocity (mm/year) obtained for the ascending Path 219 ALOS-1, using the small baselines method implemented in StaMPS (Hooper et al., 2012) for the period January 2007 to March 2011. Positive LOS velocity values corresponds to displacements toward the satellite, i. e., uplifting. Black triangles and star represent the Three Sister complex volcano system and the approximate center of the uplifting area. StaMPS LOS velocity results were noisy and we post-processed to reduce undesirable oscillations of non-volcanic origin. We applied a band pass filter to retain spatial deformation signals between 10 and 0.8 km using a median filter (GMT blockmedian). Results indicate a 6 km circular uplift pattern west of South Sister with a mean LOS velocity of approximately 5–10 mm/year, consistent with a value obtained for the Husband CGPS station during the same period (5.21 mm/year), shown as circle with a black outline.

4 Methods

Here, we introduce a mathematically rigorous strategy for the joint inversion of time-dependent InSAR (different look angles and sensors, high spatial resolution) and continuous GPS (daily sampling) data to achieve a complete timeline of volcanic activity and quantify a single time series of volume flux rates. The strategy captures the benefits of each system avoiding the time evolution determination on a point-by-point basis. It is based on the two-step approach proposed by González et al. (2013) that produces time series of volume from sets of different look angles and satellite sensors once an active source is determined for an inflation period. In this section, we provide a description of the González et al. (2013) algorithm and extent it 1) to include continuous GPS data and therefore to combine different components and/or look vectors using a unique source model and 2) to afford a defined method of truncation of the TSVD technique used to regularize the inverse problem, with the goal of finding the time evolution of volume and therefore to improve the accuracy on the estimation of volume time series. First, we show how the active source location (horizontal position and depth) and geometry is determined using a Bayesian inversion approach. Subsequently, we solve for temporal evolution of volume.

4.1 Source Characterization

First of all, we infer the active magmatic source through a Bayesian inversion, using InSAR data spanning the period of maximum displacement. The horizontal location, depth and geometry of the inflation source at Three Sisters were computed using the MATLAB^®-based software package GBIS (Geodetic Bayesian Inversion Software) (Bagnardi and Hooper, 2018), which estimates source parameters through a Markov chain Monte Carlo method and uses, among others, analytical forward models from dMODELS software package (Battaglia et al., 2013). It obtains the posterior probability distributions (PDFs) for all model parameters by taking into account uncertainties in the data (e.g., data errors). To achieve this, considering the pattern of surface deformation, we employ simple elastic models such as point source (Mogi, 1958), prolate spheroid (Yang et al., 1988) and sill-like (Fialko et al., 2001) models. An elastic, homogeneous and isotropic half-space is assumed in all the approaches with Poisson’s ratio 0.25. We assumed, as previous studies (Dzurisin et al., 2006, 2009; Riddick and Schmidt, 2011), a stationary source and we used an interferogram spanning August 1997 - September 2000 to look for source parameters. This interferogram fulfills two important conditions to determine the best-fitting static displacements: 1) it spans the shortest time during the period of maximum deformation; 2) it shows acceptable signal-to-noise ratio. InSAR spatially correlated error (caused mainly by the “wet” tropospheric delay) is estimated by modeling experimental semivariograms in deformation-free regions (Bagnardi and Hooper, 2018). InSAR data are subsampled with an adaptative quadtree method (Decriem et al., 2010) to reduce the computational cost of the Bayesian inversion, particularly for sill-like and prolate spheroid sources. The inversion computes 2×10⁶ iterations for spherical point source, 5×10⁶ iterations for sill-like source and 8×10⁶ iterations for prolate spheroid, which stabilizes the inversion procedure.

4.2 Temporal Evolution of the Source Volume Changes

Once the magmatic source is fixed, we perform the quasi-dynamic time-dependent model using a linear inversion scheme to look for the volume changes at each interferogram’s period and the cumulative volume changes since the first GPS observation. Both volume changes at each interferogram and cumulative volume changes from the GPS data are used to solve for the time evolution of volume using TSVD.

4.2.1 First Step: Piecewise Volume Changes Over Temporal Data Periods

Once the location and geometry of the inflation source are fixed, we determine the volume changes over the corresponding time intervals (increments of volume changes, $\text{Δ}V$) for both, InSAR and CGPS data sets, which are assumed to be uncorrelated. In this way, observations from several interferograms and GPS sites can be combined to estimate increments of volume changes assuming a unique source model. Each volume change, $\text{Δ}{V}_{ij}$, records 1) the incremental volume change between two acquisition dates, t_i and t_j from an interferogram or 2) the cumulative increment of volume since the first observation, i.e., t_j = t₀, being t₀ the starting date of CGPS.

A linear inversion scheme using Weighted Least Squares (WLS, Menke, 1989) is applied. The inversion is constrained by 55 interferometric pairs (ERS, ENVISAT, Sentinel-1), one ALOS-1 interferogram (cumulative LOS deformation time-series) and a 3-component GPS time series. The forward problem is defined by d = Gm, where d is the data vector (InSAR or GPS), m is the model parameter ($\mathbf{\Delta }\mathbf{V}$) vector, and G is the Green’s function matrix representing the impulse response for the specific elastic source, projected into the three components of GPS or the satellite line-of-sight. Therefore, a total of 5,064 independent linear inversions were performed to find the increments of volume changes, $\text{Δ}{V}_{ij}$, given the set of interferograms and 5,008 cumulative GPS displacements.

The least square estimator of each inversion, $\widehat{\mathbf{m}}$, is given by:

$\widehat{\mathbf{m}}={\left[{\mathbf{G}}^T{\mathbf{C}}_{\mathbf{d}}^{-1}\mathbf{G}\right]}^{-1}{\mathbf{G}}^T{\mathbf{C}}_{\mathbf{d}}^{-1}\mathbf{d},(1)$

with the cofactor matrix ${\mathbf{C}}_{\widehat{m}}={\left[{\mathbf{G}}^T{\mathbf{C}}_{\mathbf{d}}^{-1}\mathbf{G}\right]}^{-1}$. We considered a diagonal variance-covariance matrix, C_d, assuming that all data are independent, which significantly reduces the computation time of the inversions. Hence, we ignore the possible spatial and temporal correlation noise in InSAR data (e.g., pixel correlation due to atmospheric artifacts, topography structures, repeated acquisitions) and between GPS components (Lohman and Simons, 2005; Biggs et al., 2010).

4.2.2 Second Step: Volume Changes Time-series

We want to solve for the temporal evolution of volume change for each observed epoch t_k from $\text{Δ}{V}_{ij}$ obtained on first step considering both InSAR and CGPS data. Instead of volume change itself, the rate of volume changes is inverted as a function of time by applying the Short Baseline Subset Approach (SBAS, Berardino et al., 2002). This prevents the presence of large discontinuities in the final solution.

Let $\mathbf{\Delta }\mathbf{V}$ be the data vector of volume changes over the corresponding time intervals (N × 1), and $\dot{\mathbf{V}}$ the unknown vector of volume change rates (M × 1) between adjacent epochs, t_j where the overdot means differentiation over time. Then,

$\dot{\mathbf{V}}=\left[\begin{array}{c}\dot{v_1}=\frac{v_1}{t_1-t_0}\\ ⋮\\ \dot{v_M}=\frac{v_M-v_{M-1}}{t_M-t_{M-1}}\end{array}\right].(2)$

The usual method of converting the observations $\mathbf{\Delta }\mathbf{V}$ on volume change rates is:

$\mathbf{B}\dot{\mathbf{V}}=\mathbf{\Delta }\mathbf{V},(3)$

where B is the design matrix (N × M). To determine the components of B, we define a (M × 1) vector E, containing the single epochs t_j present in all time intervals and sorted in chronological order, for j = 1, ..., M; and a (N × 2) vector F, whose columns are the slave (t_slavei) and master ($t_{{\text{master}}_i}$) epochs of each i time interval, for i = 1, ..., N. Therefore, the (i, j) component of the design matrix is B_ij = (t_j+1 − t_j) for $t_{{\text{slave}}_i}\le t_j<t_{{\text{master}}_i}$ , and zero elsewhere. In the case of cumulative $\mathbf{\Delta }\mathbf{V}$ (i.e., continuous GPS), B presents lower triangular matrix blocks. For example, if volume changes are obtained over different time intervals, i.e., if t_AB, t_BC, and t_AC (from InSAR data), and t_CD, t_CE, t_CF, (from CGPS data) are $\text{Δ}{v}_{AB}$, $\text{Δ}{v}_{BC}$, $\text{Δ}{v}_{AC}$, $\text{Δ}{v}_{CD}$, $\text{Δ}{v}_{CE}$ and $\text{Δ}{v}_{CF}$, the design matrix is given by:

$\begin{array}{c}\left(\begin{array}{ccccc}\left(t_B-t_A\right)& 0& 0& 0& 0\\ 0& \left(t_C-t_B\right)& 0& 0& 0\\ \left(t_B-t_A\right)& \left(t_C-t_B\right)& 0& 0& 0\\ 0& 0& \left(t_D-t_C\right)& 0& 0\\ 0& 0& \left(t_D-t_C\right)& \left(t_E-t_D\right)& 0\\ 0& 0& \left(t_D-t_C\right)& \left(t_E-t_D\right)& \left(t_F-t_E\right)\end{array}\right)\left(\begin{array}{c}{\dot{v}}_{AB}\\ {\dot{v}}_{BC}\\ {\dot{v}}_{CD}\\ {\dot{v}}_{DE}\\ {\dot{v}}_{EF}\end{array}\right)=\left(\begin{array}{c}\text{Δ}{v}_{AB}\\ \text{Δ}{v}_{BC}\\ \text{Δ}{v}_{AC}\\ \text{Δ}{v}_{CD}\\ \text{Δ}{v}_{CE}\\ \text{Δ}{v}_{CF}\end{array}\right)\end{array}.(4)$

To illustrate the simple example in Eq. 4, let all dates, t_A, t_B, t_C, t_D, t_E, and t_F be equally spaced at time intervals of 1 year, i.e., t_AB = 1 year and so on, and $\mathbf{\Delta }\mathbf{V}=\left[\mathrm{2,1,3,1,2,3}\right]$ ×10⁶ m³. In this case, standard least squares can be applied, given $\dot{\mathbf{V}}=\left[\mathrm{2,1,1,1,1}\right]\times {10}^6{\text{m}}^3/\text{year}$. The cumulative volume times series is then $\mathbf{V}=\left[\mathrm{2,3,4,5,6}\right]\times {10}^6{\text{m}}^3$, meaning a linear inflation rate of $2\times {10}^6{\text{m}}^3/\text{year}$ in time interval t_AB and a posterior linear inflation rate of $1\times {10}^6{\text{m}}^3/\text{year}$.

However, the set of $\text{Δ}{V}_{ij}$ forms, in general, an unconnected set of observations with at least one time step not directly related to data, making Eq. 4 an ill-posed problem without solution even in the least square sense. Thus, given P a definite positive matrix, a least-square solution, $\widehat{\dot{\mathbf{V}}}={\left({\mathbf{B}}^T\mathbf{P}\mathbf{B}\right)}^{-1}{\mathbf{B}}^T\mathbf{\Delta }\mathbf{V}$, is not possible since Eq. 4 constitutes an ill-posed unstable model, with one or more eigenvalues of the normal matrix ${\mathbf{B}}^T\mathbf{P}\mathbf{B}$ close to zero. This fact is responsible for large uncertainty on the estimated volume change rates, $\widehat{\dot{\mathbf{V}}}$. Alternatives to the least square method can be proposed for an improved estimate of $\widehat{\dot{\mathbf{V}}}$: Tikhonov regularization (Tikhonov and Arsenin, 1977), Bayesian and stochastic inferences (Backus, 1988), Truncated Singular Value Decomposition (TSVD) or Principal Components (Lawless and Wang, 1976; Hansen, 1990; Hansen, 1992). Here, we consider TSVD as proposed by González et al. (2013).

4.2.3 Regularized Linear Joint Inversion

A key difficulty in applying the TSVD method is how to set up proper criteria to truncate eigenvalues due to the lack of a theoretically solid foundation to discard small nonzero eigenvalues. We developed a strategy to circumvent this difficulty based on Picard condition and L-curve methodology. In such way, we are assured a good balance, filtering out enough noise without losing too much information in the computed solution (Hansen and O’Leary, 1993). Furthermore, we included some estimations of the error of data ($\mathbf{\Delta }\mathbf{V}$) to establish some uncertainty in the $\widehat{\dot{\mathbf{V}}}$ estimator. To estimate the uncertainties, the Weighted Generalized Inverse method (Menke, 1989) permitted the use of the “a priori” information from the data ${\mathbf{C}}_{\mathbf{\Delta }\mathbf{V}}$ (and optional model, C_V = A^TPA${\mathbf{C}}_{\dot{\mathbf{V}}}={\mathbf{A}}^T\mathbf{P}\mathbf{A}$) covariance matrix. Such matrices can be decomposed as:

$\begin{array}{c}{\mathbf{C}}_{\mathbf{\Delta }\mathbf{V}}={\left[\mathbf{D}{\mathbf{D}}^T\right]}^{-1}\\ {\mathbf{C}}_{\dot{\mathbf{V}}}={\left[\mathbf{S}{\mathbf{S}}^T\right]}^{-1}\end{array}(5)$

where the D (N × N) and S (M × M) matrices are determined from the eigenvalue problem of each covariance matrix. In our case, no model covariance information is used, so C_V = I. ${\mathbf{C}}_{\dot{\mathbf{V}}}=\mathbf{I}.{\mathbf{C}}_{\mathbf{\Delta }\mathbf{V}}$ is obtained by error data propagation through $\widehat{\mathbf{\Delta }\mathbf{V}}$ estimator. The utility of D and S is the introduction of a transformed coordinate system where data (and optional model) parameters each have uncorrelated errors and unit variance. Therefore, $\mathbf{\Delta }{\mathbf{V}}_{\mathbf{n}\mathbf{e}\mathbf{w}}=\mathbf{D}\mathbf{\Delta }\mathbf{V}$, ${\dot{\mathbf{V}}}_{new}=\mathbf{S}\dot{\mathbf{V}}$ and B_new = DBS⁻¹ give the transformation of data, model parameters and forward operator in the new system of coordinates. TSVD is applied to B_new with a specific regularization method to find the Principal Components of the observation set ($\mathbf{\Delta }\mathbf{V}$). Then, the problem is back to the original coordinates to achieve the solution and finally, the volume change time-series is obtained by integrating the volume change rate in time:

$V\left(t\right)={{\displaystyle \sum }}^{\text{}}\dot{V}\delta t(6)$

4.2.4 Regularization: Techniques Used for Truncation of Small Eigenvalues

Some workable criteria for truncation in interdisciplinary problems include L-curve, Discrete Picard condition and Generalized cross validation (GCV) (Hansen, 1992; Hansen, 2007; Hansen and O’Leary, 1993). Methods like GCV sometimes fail to find the appropriate regularization parameter (flat local minima), whereas the L-curve gives a robust estimation (Hansen and O’Leary, 1993) and the appropriate smoothing solution, which is very attractive from a mathematical point of view. We thus designed a strategy based on a L-curve to set up proper criteria to truncate the eigenvalues.

First, we considered the Discrete Picard condition to explain the instability of the transformed linear inverse problem (Eq. 3) and disregarded the smallest singular values (Hossainali et al., 2010; Hansen, 1990; Hansen, 2007):

${\Vert \dot{\mathbf{V}}-{\dot{\mathbf{V}}}_{T_1}\Vert }_2\le P^{\frac{1}{2}}\underset{1\le i\le M-P}{\max }\left\{\frac{\left|u_i^T\text{Δ}{V}_i\right|}{s_i}\right\}(7)$

where ${\Vert \dot{\mathbf{V}}-{\dot{\mathbf{V}}}_{T_1}\Vert }_2$ is the regularization error, $\dot{\mathbf{V}}$ and ${\dot{\mathbf{V}}}_{T_1}$ being the exact and truncated SVD solutions; P is the regularization parameter value, s_i are the singular values, and $\left|u_i^T\text{Δ}{V}_i\right|$ are called Fourier coefficients ($\text{Δ}{V}_i$ are data and u_i the corresponding eigenvectors of the data space). The Discrete Piccard condition is satisfied if, for all singular values larger than P, the corresponding Fourier coefficients decay faster on average than s_i.

The L-curve method is applied to ${\dot{\mathbf{V}}}_{T_1}$ resulting in turn from applying (Eq. 7) through a log-log plot of the norm of a regularized solution $\Vert {\dot{\mathbf{V}}}_{T_2}{\Vert }_2$ vs. the norm of the corresponding residual norm ${\Vert \mathbf{B}{\dot{\mathbf{V}}}_{T_2}-\mathbf{\Delta }\mathbf{V}\Vert }_2$. As recommended by Hansen (1992) we fit the log-log plot of discrete points with some curves, choosing a 2D spline curve and then search for the truncation parameter by computing the L-corner (maximum curvature point). This corner of the spline curve is approximated to the closest discrete point. The L-corner is located exactly where the solution changes in nature from being dominated by regularization errors to being dominated by the residual size. This regularization filters out the contribution of small singular values and noisy data.

5 Results

5.1 Re-Evaluation of Source Location and Geometry

We performed the Bayesian inversion for point source, prolate spheroid and sill-like sources, with similar results. The results are represented in Figures 4, 5 and Tables 1, 2 report the prior information and the PDFs with the 95% credible intervals for all model parameters, respectively. The inversion reveals that the surface displacements can be explained by a spherical point source with depth (4500–6000)m and $\text{Δ}V$ (7–13)×10⁶ m³, by a sill-like source with depth (5600–7200)m, radius (220–400)m and dimensionless excess pressure 0.05–0.30 and by a prolate spheroid source with depth (5300–7400)m, major semi-axis (240–720)m, aspect ratio (0.22–0.37)m, dimensionless excess pressure 0.38–8.61, strike angle (49–102)^° and plunge (78–85)^°. The descending ERS Track 385 wrapped interferogram reveals a near axisymmetric deformation pattern, with a maximum LOS surface displacement of ∼5 cm recorded at the center of the uplifting pattern (Figures 4A, D and G). Figures 4B, E and H present the predicted spherical point, sill-like and prolate spheroid forward models, using the median value of the PDF of the model parameters. As expected from the deformation pattern, spherical point and sill-like models are very similar, suggesting that the geometry of the source is far from unique. The extra modeling parameters of the prolate spheroid do not improve the misfit. Therefore, we favor the simplest spherical point source model over a sill-like and prolate spheroid source, to fit deformation pattern displayed in Figures 4A, D and G. Blue stars represent the horizontal location of spherical point, sill-like and prolate spheroid estimated sources (Table 2). Figures 4C, F and I show the residuals between observed LOS displacement, and spherical point, sill-like and prolate spheroid model predicted displacements. The residual is larger close to the Three Sisters complex volcano (green triangles), due to orbital and topographic contributions, and also in the western half of the uplift pattern, where data are less dense.

FIGURE 4

FIGURE 4. Wrapped InSAR data and model results for Track 385 of the ERS satellite. (A, D, G) Line-of-sight (LOS) deformation observed in a nearly 3 years period from August 24, 1997 and September 17, 2000, considering only pixels with coherence >0.2. Green triangles represent the Three Sister complex volcano system. (B, E, H) Bayesian model and horizontal location for the median a posteriori probability solution (blue star) for a predicted Mogi, sill-like and prolate spheroid source, respectively. (C, F, I) Residual maps for Mogi, sill-like and prolate spheroid source. The model parameters results are presented in Figure 5 and Tables 1, 2.

FIGURE 5

FIGURE 5. Posterior probability distributions for the Mogi, sill-like and prolate spheroid source models. Black solid lines show the optimal value for the corresponding model parameter (Tables 1, 2).

TABLE 1

TABLE 1. Prior information for the Elastic half-space Spherical Point Pressure (Mogi, 1958), Penny-shaped sill-like (Fialko et al., 2001) and Spheroid (Yang et al., 1988) deformation sources.

TABLE 2

TABLE 2. Bayesian inversion results, with the median a posterior probability solution and the 95% credible intervals, for the Elastic half-space Spherical Point Pressure (Mogi, 1958), Penny-shaped sill-like (Fialko et al., 2001) and Spheroid (Yang et al., 1988) deformation sources.

Figure 5 displays the histograms of marginal PDF for the four spherical point, five sill-like and eight prolate spheroid source parameters. Black solid lines show the optimal values for the corresponding model parameters. For the sill-like source, the radius and dimensionless excess pressure PDFs exhibit bi-modality (slightly unstable inversion result). For the prolate spheroid source, the aspect ratio between semi-axes and the excess pressure PDFs display also a bi-modal shape. The major semi-axis PDF exhibits multi-modality.

5.2 Source Inflation Time-Series

We performed a CGPS and InSAR joint inversion to obtain the time-line of volume changes, considering the best fitting source location for the point source geometry to better characterize the time-dependent inflation of the magma source at Three Sisters.

To apply our two-step-approach (section 3.2), we use the median, and 5% and 95% percentile values of the PDF of depth estimated by the Bayesian inversion. The corresponding values are: d_median = 5000 m, d_5% = 4500 m and d_95% = 6000 m. The volume change time-series is determined using InSAR data from four satellite missions (ERS, ENVISAT, ALOS-1 and Sentinel-1), on five different tracks and look angles, and CGPS data from HUSB.

First, we obtain the increments of each volume change, $\text{Δ}{V}_i$, relating the Green’s functions (representing the source impulse for a buried point source) to the LOS deformation data observed along each satellite track and the three component CGPS data. Figures 6, 7 respectively show results regarding estimation of the median value of source depth (d_median = 5000 m) for InSAR and CGPS data. The cumulative increments of volume changes detected at HUSB show gaps due to ice and snow accumulation during winter. By means of the daily GPS measurements, the corresponding increments of volume change, $\text{Δ}{V}_i$, are more uniform for the CGPS data sets (Figure 7), but more variable for the individual SAR data sets (Figure 6).

FIGURE 6

FIGURE 6. Increments of volume change obtained for all interferograms (ERS, ENVISAT, ALOS-1 and SENTINEL-1), according with the median value for the source depth.

FIGURE 7

FIGURE 7. Cumulative increments of volume changes for the CGPS station, Husband. The figure shows the results according with the median value for the source depth.

Finally, we applied the Picard Plot condition suitable to understand the conditions of the ill-posed problem (Eq. 3). Figure 8A shows how s_i only decays faster than the Fourier coefficients ($\left|u_i^T\text{Δ}{V}_i\right|$) for the smallest nonzero singular values. Hence, the problem can be considered stable, discarding the last 10% of the singular values. Due to the stability of the problem, the Picard Plot provides no clues about the appropriate level of truncation (Eq. 7). Therefore, we use L-curve to determine the truncation level. L-curve criterion is fulfilled when L-corner = 1198, i.e., when only the first 24.1% nonzero singular values are used in the inversion (Figure 8B).

FIGURE 8

FIGURE 8. (A) Discrete Picard Plot condition, suited for the analysis of ill-posed problems. The solution is stable when the Fourier coefficients, $\left|u_i^T\text{Δ}{V}_i\right|$, on average decay to zero faster than the reciprocal singular values, s_i. In this case, the problem can be consider stable, discarding the last 10% singular values. (B) L-curve showing the trade-off between minimizing the residual norm ($\Vert \mathbf{B}{\dot{\mathbf{V}}}_{T_2}-\mathbf{\Delta }\mathbf{V}{\Vert }_2$) and minimizing the regularized solution size (${\left|\left|{\dot{\mathbf{V}}}_{T_2}\right|\right|}_2$). The L-corner (represented in red) is located exactly where the solution changes in nature from oversmoothing (i. e, dominated by regularization errors) to being dominated by residual size.

The analysis and results of the final inflation time-series are shown in Figures 9 and 10, and Tables 3 and 4. The inflation time-series associated with the median estimation of the source depth suggests a maximum volume change rate of $\sim 1.60\times {10}^6{\text{m}}^3/\text{yr}$ during 1999–2001 and a subsequent rate as much as $\sim 0.75\times {10}^6{\text{m}}^3/\text{yr}$ for the period 2015- January 2020. Data since 2018 show a subtle, but significant change in the trend, instead of following asymptotic behavior.

FIGURE 9

FIGURE 9. $\text{Δ}{V}_i$ data (increment of volume change obtained from InSAR and CGPS deformation data, Eq. 1), vs. Simulated observations ($\text{Δ}{V}_{i_{\text{reconstruction}}}=V_{\left(i+1\right)}-V_i$, Eq. 6). (A, B) Inversion solutions with regularization, using only the first 0.2 and 2% nonzero singular values, respectively (oversmoothing solution). (C) Regularized solution using the L-curve criteria corresponding with 24.1% nonzero singular values (appropriate smoothing solution). (D) No regularization case (undersmoothing solution).

FIGURE 10

FIGURE 10. Volume change time series for a Mogi source (constrained by both the GPS and InSAR data), at different depth estimations (Table 1): median value (light green), 95 and 5% percentile values (yellow and dark green). Cian solid lines represent predicted curves through Non Linear Least Squares (NLS), Least Absolute Residuals (LAR) method, following the e-folding characteristic shape function of Eq. 8. Stars show the location of the beginning of the exponential function: for data associated with median depth value estimation (light green), 95 and 5% percentile values (yellow and dark green).

TABLE 3

TABLE 3. Source location comparison for Three Sisters (assuming a Mogi source) from Previous Studies and the bayesian inversion carried out in this Study.

TABLE 4

TABLE 4. e-folding parameters obtained by NonLinear Least Squares, using the Least Absolute Residual method (LAR).

6 Discussion

6.1 Source Characterization

Studies at Three Sisters using InSAR interferometric pairs and stacks (Wicks, 2002; Riddick and Schmidt, 2011), GPS (campaign and continuous) and leveling (Dzurisin et al., 2009) assessed various source geometries such as spherical point, sill-like or crack and ellipsoidal. These different sources can all fit the data in a satisfactory way. Our results are consistent with previous findings (Table 3). However, volume change rates and depths vary slightly, possibly due to the fact that: 1) there may be a poorly resolved deeper magma source, 2) inversions were limited to purely kinematic models, 3) the source not being a stationary, pressurized cavity in an isotropic elastic half-space, thus producing bias due to spatial or temporal considerations, 4) diverse inversion techniques and related possible mathematical artifacts, 5) different types of data sets and 6) ambiguity of source geometries. We assumed a simple, stable, purely kinematic model as a valid approach, following the results of Dzurisin et al. (2009) and Riddick and Schmidt (2011), for estimating volume time series. Now, we focus on discussing the implications of 3), 4), 5), and 6).

We revisited the assumption of location stationarity of the inflation source, with focus on the most recent periods. Riddick and Schmidt (2011) already showed that the temporal evolution of the uplift signal can be represented by a stationary volcanic source geometry and location with a decreasing inflation rate at least from 1992 to 2010. The extent of deformation pattern remains constant over 1992–2010 time period providing a source depth range compatible with the uncertainties of inversion models, as it is expected given that the extent of deformation pattern mainly depend on the source depth and strength. The inversion for source parameters using Sentinel-1 interferograms (2014–2018) reveals very large uncertainties on the parameters. Such results could be due to the interferograms’ low signal-to-noise ratio caused by slower uplift rates during this period. Nevertheless, the best fitting spherical sources are not able to predict the observed displacements in the HUSB CGPS displacement time series. Moreover, the lack of substantial changes in the trends of each component of the HUSB CGPS time series indicates that the source might not have changed position. Therefore, we assumed the source has not changed significantly either in shape nor in location since the onset of deformation.

A range of common techniques to estimate source location has been used, like forward modeling, grid search by iteratively fixing of one parameter, arithmetic mean to obtain range values, or grid search. However, the Bayesian approach presents important advantages: 1) robust inversion for a single or more InSAR interferograms with an acceptable signal-to noise-ratio and/or GPS data, 2) rapid simultaneous inversion of all model parameters; 3) use of data uncertainty and prior model information; and 4) efficient sampling of posterior PDFs to estimate optimal model parameters and the associated range of error. Bearing in mind such advantages, to obtain a robust estimate of source geometry and location we only need geodetic data with high spatial coverage, spanning the most appropriate period (shortest time, high deformation). For this, we use the period undergoing maximum deformation, displaying as much as ∼5 cm of line-of-sight deformation (Figure 4). Selection of the ERS-1 track 365 (descending orbit) interferometric pair spanning August 1997- September 2000 satisfies both criteria. To reduce the signal-to-noise ratio, the interferogram is filtered for a pixel coherence threshold of 0.2. Other ERS-1 InSAR data were also processed for similar periods of time, but not used due to the low signal-to-noise ratio. No GPS data were available until 2001 and cannot be used to study the maximum uplift rate period.

We favor the simplest source, a spherical point source, to infer volume time series at Three Sisters. Our inversion results for spherical, prolate spheroid and a sill-like sources showed quantitatively similar results, in terms of data misfit and surface displacement pattern. We noted that there is slightly elongated pattern of the InSAR data in the North-South direction (Figure 4). This pattern cannot be perfectly reproduced with an axisymmetric source geometry. The elongation could be also be caused by the topography of Three Sisters area at the east side of the deforming area or tropospheric delays in the interferometric data. To robustly distinguish between different source geometries, we must have full three dimensional surface displacements (Dieterich and Decker, 1975). Therefore, in our case, the almost symmetric shape of the 2D deformation pattern implies that source geometry remained far from being uniquely resolved. Nevertheless, to accept more complex geometries than spherical, we should have obtained significantly lower data misfit values. Furthermore, the spheroid modeling residuals are not fully consistent with the observed pattern in the western area of deformation (Figures 4I). In this case, we also consider that the topography could have a minor effect because the highest topographic relief is concentrated on the far field area of deformation signal. Therefore, we assume that the noted asymmetry in the InSAR data should not affect the overall interpretation of the time series of volume changes.

The differences between inversion methods and data selection might explain that our optimal inversion results suggest a slightly shallower source with a corresponding smaller increment of volume change. Despite that, considering that the models fit the data well and yield similar misfit values, we conclude that it is reasonable to assume the spherical point source as the simplest kinematic model that explains the signal. Furthermore, the values for depth and increments of volume change lie within our 95% credible intervals (Table 3). Wicks (2002) processed three interferograms, obtaining an increment of volume change of $\text{Δ}V=23\times {10}^6{\text{m}}^3$ and depth of 6500 m for the one with the largest apparent signal-to-noise ratio. Ultimately, a deeper source will trade-off with a greater $\text{Δ}V$, to fit the same surface deformation. Although the time acquisition of the best interferogram of Wicks (2002) spans only 1 year more than our InSAR data interferometric pair, the important difference between our $\text{Δ}{V}_{rate}=\sim 3.00\times {10}^6{\text{m}}^3/\text{yr}$ and their $\text{Δ}{V}_{rate}=\sim 5.75\times {10}^6{\text{m}}^3/\text{yr}$ is mainly due to their depth estimation. Our InSAR interferogram matches one of the other two by Wicks (2002). For that interferogram, their model gives depths ∼1 km shallower, being closer to our depth estimation.

6.2 Time Series of Volume Changes: Regularization Using the Truncated SVD

To assess the effect of the regularization, we compared the increments of volume change ($\text{Δ}{V}_i$, Eq. 1) and the corresponding simulated observations ( $\text{Δ}{V}_{i_{\text{reconstruction}}}=V_{{\left(i+1\right)}_{\text{inversion}}}-V_{i_{\text{inversion}}}$ from the results of Eq. 6) for different levels of regularization. Figures 9A, B show the solution for an extreme regularization, using only the first 0.2% and 2% nonzero singular values (oversmoothing solutions). For the 0.2% case, the values of $\text{Δ}{V}_i$ associated with the InSAR data display a wide dispersion, and the $\text{Δ}{V}_i$ associated with GPS data acquire discrete values, i.e., the same $\text{Δ}{V}_i$ value is obtained for many different time intervals. Figure 9D represents the case of non-regularization (less filtering, maximum solution size and minimum misfit). The residual in Figure 9D is minimized because the solution reproduces even the seasonal perturbations of the GPS data (undersmoothing solution). There is a seasonal deformation of the crust associated with the surface load of the snow cover. It is possible that the magma reservoir’s internal pressure also fluctuates seasonally in response to this effect. However, from CGPS data alone we cannot resolve the cause of these fluctuations. Therefore, a smoother solution is preferred to depict the long-term deformation time series. This is given by the combination of Picard condition and L-curve criteria; it corresponds to the appropriate smoothing solution, i.e, 24.1% nonzero singular values (Figure 9C).

6.3 Time-Scales of Inter-Eruptive Uplift Signals: Three Sisters and Other Volcanoes

The continuous and extended regularized time-series of volume change allows us to study the inflation process in detail. Riddick and Schmidt (2011) proposed a piece-wise linear parametrization with two changes in rate provides a good fit to uplift rates till 2009 explaining two different inflation processes beginning at 1998 and 2004 at Three Sisters. This model was supported by the detection of a seismic swarm in 2004. Our denser, longer time-series clearly shows a smooth and continuous function, which we interpret as a fast inflation followed by relaxation of the crust (Figure 10). We are specifically interested on the interval of decaying rates. Consequently, the time-series is divided into two main intervals separating increasing and decaying behavior of volume rates. An exponential function can reasonably reproduce the relaxation process. Therefore, we propose a piece-wise parametric of the form:

$V\left(t\right)=\{\begin{array}{cc}d,& t<t_0\\ a-b\text{exp}-\left(c\cdot \left(t-t_0\right)\right),& t\ge t_0\end{array}(8)$

where a, b and d are constants, $1/c=\mathit{ϵ}$ is the characteristic relaxation time constant, here after named e-folding parameter, and t₀ is the start of the exponential trend. We solved the parameters of this model using a non-linear least-square fit. To minimize the influence of outliers, we used regression method: the Least Absolute Residuals (LAR) and Bisquare weights methods, considering also the data uncertainties (weighted). Four methods (Bisquare, Bisquare-Weighted, LAR, LAR-Weighted) show very similar fit, LAR performing the best (Figure 10; Table 4). Time-series from the median, and 5% and 95% percentiles of the PDF distribution for depth, along with the variance of the curve fitting, permit constraining the lower and upper limits of e-folding and t₀ parameters.

Our time-series spans 26 years and presents a characteristic time constant of $9.48[+0.12,-0.17]$ years. This new result updates the value of 5.3 years obtained by Dzurisin et al. (2009) from an 8-years time-series of GPS displacements from 2001 to late 2008. Riddick and Schmidt (2011) hyphotesized that one injection of magma started between June 1996 and July 1997, given StaMPS results for T385 ERS. The updated volume time-series presented in this study shows a clear exponential decay trend. We estimate a start date for the exponential trend between October 1998 and August 1999 (Table 4). These results suggest that the continuous uplift signal will be detectable for a few decades, considering volume change rates as low as 0.1×10⁶ m³. As late as January 2020, our inflation time-series indicated that the cumulative volume of a spherical point source with d_median = 5000 m is 29.1×10⁶ m³. For d_5% = 4500 m and d_5% = 6000 m, the values are 22.7×10⁶ and 39.9×10⁶ m³, respectively. This range of values is consistent with those predicted by Dzurisin et al. (2009) for a prolate spheroid model, 44.9 to 51.6×10⁶ m³ (where the uncertainty is 10–20% of those values). This estimation of the cumulative volume is 10–20 times less than the volume erupted from Mount St. Helens in May 1980 (Wicks, 2002).

To compare the characteristic relaxation time for other volcanoes with recent and well-studied unrest episodes, we compiled and modeled the available geodetic and volume time-series of the following volcanic systems: Okmok (Biggs et al., 2010), Long Valley (Hill et al., 2020), Uturuncu (Lau et al., 2018), Laguna del Maule (Le Mével et al., 2015), Yellowstone (Tizzani et al., 2015), Campi Flegrei (Troise et al., 2007), Santorini (Parks et al., 2015), Alutu (Hutchison et al., 2016), Agung (Syahbana et al., 2019). We have selected unrest periods showing deformation that deviates from the background trend and are characterised by an exponential decay. These unrest periods span from the beginning of the deformation series or from a new increment of uplift, until it occurs a new change in the trend, indicated by a departure from the flattened part of the exponential behavior. Therefore, the analyzed periods are based solely on the exponential shape of data regardless of the unrest outcome (non-eruptive unrest or pre-eruptive unrest) and of the duration of the entire inter-eruptive period. Only a posteriori can it be known whether a period of uplift with exponential decay trend can be taken as indicative or related to volcanic eruptions. Considering this, and to avoid possible bias, we analyze a heterogeneous group of volcanoes (different types, diverse geographical locations, distinct deformation/volume change time-series and different stages within each inter-eruptive period), in order to discern whether the e-folding parameter may be or not a useful parameter for understanding posterior volcanic activity. We estimated the e-folding parameter and t₀ following Eq. 8 in order to compare different volcanic systems. Greater e-folding indicates longer characteristic relaxation times, as shown in Figure 11A. The estimated e-folding parameters vary between 0.033 and 10 years, for uplift and inflation episodes lasting between 60 days and 26 years. Our selected volcanoes, particularly those of the North and South American volcanic arcs, present the longest e-folding times (Figure 11B).

FIGURE 11

FIGURE 11. (A) Estimated e-folding parameter of other available geodetic and volume time series from volcanoes with recent and well-studied unrest episodes: Okmok (Biggs et al., 2010), Long Valley (Hill et al., 2020), Uturuncu (Lau et al., 2018), Laguna del Maule (Le Mével et al., 2015), Yellowstone (Tizzani et al., 2015), Campi Flegrei (Troise et al., 2007), Santorini (Parks et al., 2015)), Alutu (Hutchison et al., 2016), Agung (Syahbana et al., 2019). (B) Geographic location and classification according to the type of Volcano. (C, E, F) Normalized uplift or volume change, y′ (Eq. 9) as function of normalized time, t′ (Eq. 10).

From the small selection of volcanoes, neither the type nor composition of the volcanoes seems to be decisive for the characteristic relaxation time associated with their inter-eruptive periods. However, those with a shorter e-folding display significant changes in their behavior. For instance, Campi Flegrei exhibited an increment in displacement from 1968–1972, with an e-folding of 1.38 years. Next, the subsidence rate was small until 1982, and the deformation trend changed due to a new uplift episode during 1982–1984, followed by subsidence during 1985–2004. Those features were related to an overpressure in the magmatic source and fracturing of the rock volume between the magmatic fluids and the aquifer (Troise et al., 2007). Alutu underwent two inflation pulses, the latest showing a short e-folding of 0.033 years, during the period October-December 2008, then a slow deflation took place. These short time-scale suggest the migration of hydrothermal or magmatic fluids or volatiles (Hutchison et al., 2016). Agung went through an uplift from August-October 2017, characterized by an e-folding of 0.038 years, then in late November, a phreatomagmatic eruption and stronger explosions took place (Syahbana et al., 2019). Santorini presented a source inflation process with an e-folding of 0.28 years for the period October 2011- August 2012 (Parks et al., 2015). Then, its subsidence rates increased in the post-unrest period 2012–2017, suggesting the superimposition of various deformation sources (Papageorgiou et al., 2019). The Okmok inflation episode during the period May 2002 - September 2007 had an e-folding of 1.24 years. Although there was a small amount of deflation (Biggs et al., 2010), the general trend can be modeled as an exponential decay. After this inter-eruptive episode, a phreatomagmatic eruption occurred in July–August 2008 (Larsen et al., 2015). Long Valley deformation series presents an e-folding of 5.31 years for the period 1978–1988. No significant seismicity was detected during the interval ∼1986–1988. After the exponential decaying trend, there was a renewed unrest, characterized by recurring earthquake swarms and tumescence of the resurgent dome (Hill et al., 2020). Yellowstone went through an uplift during 2005–2010, with an estimated e-folding of 2.1 years. Since 2015, subsidence of Yellowstone caldera has occurred at an average rate of 2–3 cm per year, as reported by Yellowstone Volcano Observatory (USGS). Laguna del Maule is the only volcano that yields a high relaxation time value for a short inter-eruptive period (2010–2014), according to Le Mével et al. (2015). However, the fit of the data for this period could also be due to a linear inflation pulse superimposed on an exponentially decaying deformation rate. On the other hand, Uturuncu and Three Sisters present the longest e-folding (8.93 years and 9.48 years), without showing significant changes in their volcanic activity. Although the e-folding parameter seems to be an informative variable in the magnitude of the inter-eruptive period time scales, it does not provide any parametric information that differentiates pre-eruptive unrest (Agung, Okmok) from non-eruptive unrest (e.g. Alutu, Campi Flegrei).

We re-scaled the observed time-series to properly emphasize similarities on the exponential decay. We normalize displacement or volume change (y′), as a function of normalized time by means of the e-folding parameter (t′):

$y^{\prime }\left(t\right)=\left(y\left(t\right)-y_0\right)/\left(y_{t_{\infty }}-y_0\right)(9)$

$t^{\prime }=\left(t-t_0\right)/\mathit{ϵ}(10)$

where y₀ is the displacement or volume change at t₀, $y_{t_{\infty }}$ is the value after total relaxation (i.e., at $t=t_{\infty }$), and t₀ is the onset of the exponential function. It is worth noting that using the e-folding parameter accurately represents the characteristic relaxation times, and hence re-scales invariantly the observations.

Figures 11C,D,E show the resulting normalized time-series for each volcano, revealing a strikingly similar pattern. This behavior seems to be independent of the e-folding value or duration of the inter-eruptive episode. Accordingly, the temporal invariance could indicate that there is a limited set of physical scenarios underlying inter-eruptive inflation episodes. This evidence seems to support inter-eruptive physical processes with exponential time-dependent solution. Several physical models could explain deformation with this pattern (e.g., Lengliné et al., 2008; Reverso et al., 2014; Walwer et al., 2019). Dzurisin et al. (2009) put forward several mechanistic explanations for Three Sisters: 1) hydraulic pressurization or instantaneous response of the crust to continued intrusion at depth; 2) pressurization of the hydrothermal system; 3) viscoelastic response of the crust due to an intrusion emplaced at depth. Our analyses estimate the start of the exponential decay around 1998–1999, in agreement with previous studies (Wicks, 2002; Dzurisin et al., 2006; Dzurisin et al., 2009; Riddick and Schmidt, 2011). Although not unique on which physical mechanism is acting at Three Sisters, the observed uplift decay is consistent with a viscoelastic response of the crust.

The lack of hot springs and thermal pools rules out that the deformation could be due to an active hydrothermal system (Wicks, 2002). Recent rhyolitic eruptions close to South Sisters, 2000 years ago (Hildreth et al., 2012), and the evidence of a long-lived magmatic source at depth in Three Sisters, from spring geochemistry studies (Evans et al., 2004), are indicative of magma evolving during thousand of years before eruption. This is compatible with the formation of viscoelastic aureoles as a result of the alteration of the mineralogical and rheological properties of the surrounding rocks. Previously, Zurek et al. (2012) concluded, based on the lack of gravity change from 2002–2009, that deformation at Three Sisters reflects a viscoelastic response of the crust to an intrusion of magma. The estimated e-folding indicates viscosities of ∼10¹⁸ Pa s of a viscoelastic shell surrounding the magmatic source with a pressure change which increases in finite time from 0 to a constant value, considering theoretical estimations on the time constants and viscoelastic medium parameters (e.g., Newman et al., 2001; Bonafede and Ferrari, 2009; Del Negro et al., 2009; Segall, 2016).

7 Conclusion

The evolution of volume change time-series at active volcanoes can be studied by combining heterogeneous geodetic datasets. For Three Sisters volcano, we combined high spatial coverage from multiple InSAR satellite data and long term temporal information on the three-components of the only available continuous GPS. We improved a previous two-step approach to volume time-series reconstruction, by including a non-arbitrary truncation level. The cut-off criterion for truncation (i.e., type of filter) is necessary to obtain a solution without too much loss of resolution affecting the stability of the inversion. We proposed a method that combines the Discrete Picard Condition and the L-curve. Furthermore, our approach takes propagation errors into account in all inversion steps. The final time-series is determined considering volume change rates instead of increments of displacement, avoiding problems deriving from the amplification of uncorrelated noise between adjacent GPS data or propagation through the time-series of the uncertainty of the first acquisition.

The inflation time-series of Three Sisters since 2018 shows a noticeable change in the trend, which departs from the previous asymptotic trend toward a constant decay rate. This change can be explained by a fixed step in the position, such as that caused by a change in the instrumentation or monument stability. However, we cannot rule out a minor injection of magma or fluid pressurization beneath Three Sisters. Considering the wide range of erupted magma compositions and eruption styles, any changes in the Three Sisters background uplift behavior should be evaluated as an important indicator of future volcanic activity.

The Three Sisters volcano uplift is still on-going. The Bayesian inversion of source parameters gives 95% credible intervals, the depth for a spherical point source being between (4500–6000)m. Parametric modeling of the inflation time-series associated with the median, and 5% and 95% percentiles of source depth allows us to constrain the onset of the exponential trend to between October 1998 and August 1999 and the characteristic relaxation time to $9.48[+0.12,-0.17]$ years. Therefore, in the absence of different or new unrest signals, we estimate a continuous uplift signal, at diminishing but detectable rates, lasting for few decades (currently estimate to $2054[\pm 2$ years$]$). The observed uplift decay is consistent with a viscoelastic response of the crust. Our future efforts will be focused on elucidating whether the deformation could be a viscoelastic response to a very rapid magma emplacement or to several years of active intrusion of magma.

This analysis is a step toward understanding the time-scale of inter-eruptive processes. Inter-eruptive uplift/volume change signals of analyzed volcanoes show rather simple and time-scale invariant behavior, after a proper scaling. We interpret this observation as pointing to a rather reduced set of physical mechanisms underlying inter-eruptive inflation episodes that are consistent with exponential decay (viscoelastic response and/or hydraulic pressurization). Furthermore, we suggest that the magnitude of the characteristic relaxation time can be indicative of significant changes of the background behavior on volcanoes. Temporally persistent, long-lasting and overlapping uplift signals are potential confounding indicators for the classic inflation-eruption-deflation cycle model. We highlight the importance of high-temporal and continuous surface deformation monitoring to identify any departures from background temporal behavior (potentially very complex), as an indicator of future eruption hazard in persistent uplifting volcanoes. In regards to eruption forecasting, the uplift/inflation itself cannot be used as a pre-eruptive precursor without knowing what controls it through the combination of petrological and/or geophysical data.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: JAXA/METI ALOS PALSAR L1.0 was accessed through ASF DAAC, https://asf.alaska.edu [June 2018] doi:10.5067/J4JVCFDDPEW1. ERS-1, ESA [1997–2000] downloaded from ASF DAAC, https://asf.alaska.edu [June 2018]. Copernicus Sentinel data was retrieved from ASF DAAC, https://asf.alaska.edu [June 2018]. This work contains modified Copernicus Sentinel data [2014–2018], processed by ESA. HUSB CGPS data and errors were obtained from the Pacific Northwest network on the USGS website https://earthquake.usgs.gov/monitoring/gps/Pacific_Northwest/husb.

Author Contributions

SR-M processed radar images, performed modeling, carried out the comparison among different volcanic systems, prepared the figures and co-wrote the manuscript. PG and DS processed radar images. MC and PG conceptualized the study, interpreted results and co-wrote the manuscript. AN co-supervised SR-M and contributed to the manuscript. DS contributed to the manuscript.

Funding

Part of this work was carried out during research visits of SR-M to the University of Liverpool and the U.S. Geological Survey (USGS) Cascades Volcano Observatory (CVO). PG was supported by the United Kingdom NERC project Centre for the Observation and Modeling of Earthquakes, Volcanoes and Tectonics (COMET, GA/13/M/031, comet.nerc.ac.uk). SR-M, MC, and AN research was supported by MINECO projects CGL2014‐58821‐C2‐1‐R and PGC2018‐095154‐B‐I00. SR-M was also supported by MINECO fellowship BES-2015-074228.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Mike Lisowsky for providing the predicted horizontal tectonic motion at HUSB. We would like to thank Yu Jiang for his support in InSAR processing and Xiao Wang for his help with the University of Liverpool, School of Environmental Sciences computing resources. This manuscript was edited by Guido Jones, currently funded by the Cabildo de Tenerife, under the TFinnova Program supported by MEDI and FDCAN funds. We also acknowledge Maurizio Battaglia, Gilda Gurrenti and Valerio Acocella their thoughtful comments and suggestions. This is a contribution of the CSIC Thematic Platform Volcanism and Society (www.ptivolcan.es).

References

Acocella, V., Di Lorenzo, R., Newhall, C., and Scandone, R. (2015). An overview of recent (1988 to 2014) caldera unrest: knowledge and perspectives. Rev. Geophys. 53, 896–955. doi:10.1002/2015RG000492

CrossRef Full Text | Google Scholar

Backus, G. E. (1988). Bayesian inference in geomagnetism. Geophys. J. Int. 92, 125–142. doi:10.1111/j.1365-246X.1988.tb01127.x

CrossRef Full Text | Google Scholar

Bagnardi, M., and Hooper, A. (2018). Inversion of surface deformation data for rapid estimates of source parameters and uncertainties: a bayesian approach. G-cubed 19, 2194–2211. doi:10.1029/2018GC007585

CrossRef Full Text | Google Scholar

Battaglia, M., Cervelli, P., and Murray, J. (2013). DMODELS: a MATLAB software package for modeling crustal deformation near active faults and volcanic centers. J. Volcanol. Geoth. Res. 254, 1–4. doi:10.1016/j.jvolgeores.2012.12.018

CrossRef Full Text | Google Scholar

Bekaert, D. P. S., Hamlington, B. D., Buzzanga, B., and Jones, C. E. (2017). Spaceborne synthetic aperture radar survey of subsidence in Hampton roads, Virginia (United States). Sci. Rep. 7, 14752. doi:10.1038/s41598-017-15309-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Berardino, P., Fornaro, G., Lanari, R., and Sansosti, E. (2002). A new algorithm for monitoring localized deformation phenomena based on small baseline differential SAR interferograms. IEEE Internat. Geosci. Remote Sens. Symp. 2, 2375–2383. doi:10.1109/IGARSS.2002.1025900

CrossRef Full Text | Google Scholar

Biggs, J., and Pritchard, M. (2017). Global volcano monitoring: what does it mean when volcanoes deform? Elements. 13, 17–22. doi:10.2113/gselements.13.1.17

CrossRef Full Text | Google Scholar

Biggs, J., Lu, Z., Fournier, T., and Freymueller, J. T. (2010). Magma flux at Okmok Volcano, Alaska, from a joint inversion of continuous GPS, campaign GPS, and interferometric synthetic aperture radar. J. Geophys. Res.: Solid Earth 115, 1–11. doi:10.1029/2010JB007577

CrossRef Full Text | Google Scholar

Bonafede, M., and Ferrari, C. (2009). Analytical models of deformation and residual gravity changes due to a Mogi source in viscoelastic medium. Tectonophysics 471, 4–13. doi:10.1016/j.tecto.2008.10.006

CrossRef Full Text | Google Scholar

Chen, C. W., and Zebker, H. A. (2002). Phase unwrapping for large SAR interferograms: statistical segmentation and generalized network models. IEEE Trans. Geosci. Rem. Sens. 40, 1709–1719. doi:10.1109/TGRS.2002.802453

CrossRef Full Text | Google Scholar

Decriem, J., Árnadóttir, T., Hooper, A., Geirsson, H., Sigmundsson, F., Keiding, M., et al. (2010). The 2008 May 29 earthquake doublet in SW Iceland. Geophys. J. Internat.—Geophys. J. Int. 181, 1128–1146. doi:10.1111/j.1365-246X.2010.04565.x

CrossRef Full Text | Google Scholar

Del Negro, C., Currenti, G., and Scandura, D. (2009). Temperature-dependent viscoelastic modeling of ground deformation: application to Etna volcano during the 1993–1997 inflation period. Phys. Earth Planet. In. 172, 299–309. doi:10.1016/j.pepi.2008.10.019

CrossRef Full Text | Google Scholar

Dieterich, J. H., and Decker, R. W. (1975). Finite element modeling of surface deformation associated with volcanism. J. Geophys. Res. (1896–1977). 80, 4094–4102. doi:10.1029/JB080i029p04094

CrossRef Full Text | Google Scholar

Dzurisin, D., Lisowski, M., and Wicks, C. W. (2009). Continuing inflation at Three Sisters volcanic center, central Oregon cascade range, United States, from GPS, leveling, and InSAR observations. Bull. Volcanol. 71, 1091–1110. doi:10.1007/s00445-009-0296-4

CrossRef Full Text | Google Scholar

Dzurisin, D., Lisowski, M., Wicks, C. W., Poland, M. P., and Endo, E. T. (2006). Geodetic observations and modeling of magmatic inflation at the Three Sisters volcanic center, central Oregon cascade range, United States. J. Volcanol. Geoth. Res. 150, 35–54. doi:10.1016/j.jvolgeores.2005.07.011

CrossRef Full Text | Google Scholar

Dzurisin, D., Lisowski, M., and Wicks, C. W. (2017). Semipermanent GPS (SPGPS) as a volcano monitoring tool: rationale, method, and applications. J. Volcanol. Geoth. Res. 344, 40–51. doi:10.1016/j.jvolgeores.2017.03.007

CrossRef Full Text | Google Scholar

Evans, W., van Soest, M., Mariner, R., Hurwitz, S., Ingebritsen, S., Wicks, C., et al. (2004). Magmatic intrusion west of Three Sisters, central Oregon, United States: the perspective from spring geochemistry. Geology 32, 69–72. doi:10.1130/G19974.1

CrossRef Full Text | Google Scholar

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., et al. (2007). The shuttle radar topography mission. Rev. Geophys. 45, RG2004. doi:10.1029/2005RG000183

CrossRef Full Text | Google Scholar

Fialko, Y., Khazan, Y., and Simons, M. (2001). Deformation due to a pressurized horizontal circular crack in an elastic half-space, with applications to volcano geodesy. Geophys. J. Int. 146 (1), 181–190. doi:10.1046/j.1365-246X.2001.00452.x

CrossRef Full Text | Google Scholar

Fierstein, J., Hildreth, W., and Calvert, A. T. (2011). Eruptive history of South sister, Oregon cascades. J. Volcanol. Geoth. Res. 207, 145–179. doi:10.1016/j.jvolgeores.2011.06.003

CrossRef Full Text | Google Scholar

Galetto, F., Acocella, V., and Caricchi, L. (2017). Caldera resurgence driven by magma viscosity contrasts. Nat. Commun. 8, 1750. doi:10.1038/s41467-017-01632-y

PubMed Abstract | CrossRef Full Text | Google Scholar

González, P. J., Samsonov, S. V., Pepe, S., Tiampo, K. F., Tizzani, P., Casu, F., et al. (2013). Magma storage and migration associated with the 2011–2012 El Hierro eruption: implications for crustal magmatic systems at oceanic island volcanoes. J. Geophys. Res.: Solid Earth 118, 4361–4377. doi:10.1002/jgrb.50289

CrossRef Full Text | Google Scholar

Hansen, P. C. (1992). Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev. 34, 561–580. doi:10.1137/1034115

CrossRef Full Text | Google Scholar

Hansen, P. C. (2007). Regularization tools version 4.0 for matlab 7.3. Numer. Algorithm. 46, 189–194. doi:10.1007/s11075-007-9136-9

CrossRef Full Text | Google Scholar

Hansen, P. C. (1990). The discrete picard condition for discrete ill-posed problems. Bit. 30, 658–672. doi:10.1007/BF01933214

CrossRef Full Text | Google Scholar

Hansen, P. C., and O’Leary, D. P. (1993). The use of the L-curve in the regularization of discrete ill-posed problems. SIAM J. Sci. Comput. 14, 1487–1503. doi:10.1137/0914086

CrossRef Full Text | Google Scholar

Hildreth, W. (2007). Quaternary magmatism in the cascades-geologic perspectives. USGS. Professional Paper Vol. 1744, 125.

Google Scholar

Hildreth, W., Fierstein, J., and Calvert, A. T. (2012). Geologic map of three sisters volcanic cluster, cascade range, Oregon. US Geological Survey Scientific Investigations Map, Vol. 3186, pamphlet 107. p., 2 sheets, scale 1:24,000. Available at https://pubs.usgs.gov/sim/3186/.

Google Scholar

Hill, D., Montgomery-Brown, E., Shelly, D., Flinders, A., and Prejean, S. (2020). Post-1978 tumescence at Long Valley caldera, California: a geophysical perspective. J. Volcanol. Geoth. Res. 400, 106900. doi:10.1016/j.jvolgeores.2020.106900

CrossRef Full Text | Google Scholar

Hooper, A., Bekaert, D., Spaans, K., and Arikan, M. (2012). Recent advances in SAR interferometry time series analysis for measuring crustal deformation. Tectonophysics 514–517, 1–13. doi:10.1016/j.tecto.2011.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hossainali, M., Becker, M., and Groten, E. (2010). Comprhensive approach to the analysis of the 3D kinematics deformation with appliction to the kenai Peninsula. J. Geodetic Sci. 1, 59–73. doi:10.2478/v10156-010-0008-1

CrossRef Full Text | Google Scholar

Hutchison, W., Biggs, J., Mather, T. A., Pyle, D. M., Lewi, E., Yirgu, G., et al. (2016). Causes of unrest at silicic calderas in the East African Rift: new constraints from InSAR and soil-gas chemistry at Aluto volcano, Ethiopia. G-cubed. 17, 3008–3030. doi:10.1002/2016GC006395

CrossRef Full Text | Google Scholar

Larsen, J., Neal, C., Schaefer, J., Kaufman, A., and Lu, Z. (2015). The 2008 phreatomagmatIc eruption of Okmok volcano, AleutIan Islands, Alaska: chronology, deposits, and landform changes. Bull. Volcanol. 77, 69. doi:10.1007/s00445-015-0952-9

CrossRef Full Text | Google Scholar

Lau, N., Tymofyeyeva, E., and Fialko, Y. (2018). Variations in the long-term uplift rate due to the Altiplano–Puna magma body observed with Sentinel-1 interferometry. Earth Planet Sci. Lett. 491, 43–47. doi:10.1016/j.epsl.2018.03.026

CrossRef Full Text | Google Scholar

Lawless, J., and Wang, P. (1976). A simulation study of ridge and other regression estimators. Commun. Statist. Theor. Meth, A5, 307–323. doi:10.1080/03610927608827353

CrossRef Full Text | Google Scholar

Le Mével, H., Feigl, K. L., Córdova, L., DeMets, C., and Lundgren, P. (2015). Evolution of unrest at Laguna del Maule volcanic field (Chile) from InSAR and GPS measurements, 2003 to 2014. Geophys. Res. Lett. 42, 6590–6598. doi:10.1002/2015GL064665

CrossRef Full Text | Google Scholar

Lengliné, O., Marsan, D., Got, J.-L., Pinel, V., Ferrazzini, V., and Okubo, P. G. (2008). Seismicity and deformation induced by magma accumulation at three basaltic volcanoes. J. Geophys. Res.: Solid Earth 113. B12305. doi:10.1029/2008JB005937

CrossRef Full Text | Google Scholar

Lohman, R. B., and Simons, M. (2005). Some thoughts on the use of InSAR data to constrain models of surface deformation: noise structure and data downsampling. Geochem. Geophys. Geosyst. 6, Q01007. doi:10.1029/2004GC000841

CrossRef Full Text | Google Scholar

Mccaffrey, R., Qamar, A., King, R., Wells, R., Khazaradze, G., Williams, C., et al. (2007). Fault locking, block rotation and crustal deformation in the Pacific Northwest. Geophys. J. Internat. – Geophys. J. Iint. 169, 1315–1340. doi:10.1111/j.1365-246X.2007.03371.x

CrossRef Full Text | Google Scholar

Menke, W. (1989). Geophysiocal data analysis: discrete inverse theory. San Diego, Calif: Academic.

Google Scholar

Mogi, K. (1958). Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Bull. Earthquake Res. Inst. Univ. Tokyo. 36, 99–134.

Google Scholar

Moran, S. C. (2004). Seismic monitoring at cascade volcanic centers, 2004- status and recommendationsc. US Geological Survey Scientific Investigations Report, 2004-5211, 28.

CrossRef Full Text | Google Scholar

Newman, A., Dixon, T., Ofoegbu, G., and Dixon, J. (2001). Geodetic and seismic constraints on recent activity at Long Valley Caldera, California: evidence for viscoelastic rheology. J. Volcanol. Geoth. Res. 105, 183–206. doi:10.1016/S0377-0273(00)00255-9

CrossRef Full Text | Google Scholar

Papageorgiou, E., Foumelis, M., Trasatti, E., Guido, V., Raucoules, D., and Mouratidis, A. (2019). Multi-sensor SAR geodetic imaging and modelling of Santorini volcano post-unrest response. Rem. Sens. 11, 259. doi:10.3390/rs11030259

CrossRef Full Text | Google Scholar

Parks, M., Moore, J., Papanikolaou, X., Biggs, J., Mather, T., Pyle, D., et al. (2015). From quiescence to unrest: 20 years of satellite geodetic measurements at Santorini volcano, Greece. J. Geophys. Res. B: Solid Earth. 120, 1309–1328. doi:10.1002/2014JB011540

CrossRef Full Text | Google Scholar

Poland, M. P., Lisowski, M., Dzurisin, D., Kramer, R., McLay, M., and Pauk, B. (2017). Volcano geodesy in the cascade arc, United States. Bull. Volcanol. 79. doi:10.1007/s00445-017-1140-x

CrossRef Full Text | Google Scholar

Reverso, T., Vandemeulebrouck, J., Jouanne, F., Pinel, V., Villemin, T., Sturkell, E., et al. (2014). A two-magma chamber model as a source of deformation at Grímsvötn Volcano, Iceland. J. Geophys. Res.: Solid Earth 119, 4666–4683. doi:10.1002/2013JB010569

CrossRef Full Text | Google Scholar

Riddick, S. N., and Schmidt, D. A. (2011). Time-dependent changes in volcanic inflation rate near Three Sisters, Oregon, revealed by InSAR. G-cubed. 12, Q12005. doi:10.1029/2011GC003826

CrossRef Full Text | Google Scholar

Rosen, P. A., Gurrola, E., Sacco, G. F., and Zebker, H. (2012). “The InSAR scientific computing environment,” in EUSAR 2012; 9th European conference on synthetic aperture radar, 23-26 April 2012, Nuremberg, Germany 730–733.

Google Scholar

Rosen, P., Henley, S., Peltzer, G., and Simons, M. (2004). Updated repeat orbit interferometry package released. Eos, Trans. American Geophys. Union. 85, 47. doi:10.1029/2004EO050004

CrossRef Full Text | Google Scholar

Schmidt, M., and Grunder, (2011). Deep mafic roots to arc volcanoes: mafic recharge and differentiation of basaltic andesite at north Sister volcano, Oregon cascades. J. Petrol. 52, 603–641. doi:10.1093/petrology/egq094

CrossRef Full Text | Google Scholar

Scott, W. E., Iverson, R., Schilling, S., and Fisher, B. (2001). Volcano hazards in the three sisters region, Oregon. US Geological Survey Open-File Report, Vol. 18, 99–437. doi:10.3133/ofr99437

Google Scholar

Segall, P. (2016). Repressurization following eruption from a magma chamber with a viscoelastic aureole. J. Geophys. Res.: Solid Earth 121, 8501–8522. doi:10.1002/2016JB013597

CrossRef Full Text | Google Scholar

Sherrod, D. R., Taylor, E. M., Ferns, M. L., Scott, W. E., Conrey, R. M., and Smith, G. A. (2004). Geologic map of the Bend 30- × 60-minute quadrangle, central Oregon. US Geological Survey Geological Investigations Series Map I-2683, 48.

Google Scholar

Syahbana, D., Kasbani, K., Suantika, G., Prambada, O., Andreas, A., Saing, U., et al. (2019). The 2017–19 activity at Mount Agung in Bali (Indonesia): intense unrest, monitoring, crisis response, evacuation, and eruption. Sci. Rep. 9, 8848. doi:10.1038/s41598-019-45295-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Tikhonov, A., and Arsenin, V. (1977). Solutions of ill-posed problem. New York: John Wiley & Sons.

Google Scholar

Tizzani, P., Castaldo, R., Pepe, A., Zeni, G., Lanari, R., and Battaglia, M. (2015). Magma and fluid migration at Yellowstone Caldera in the last three decades inferred from InSAR, leveling and gravity measurements. J. Geophys. Res.: Solid Earth 120, 2627–2647. doi:10.1002/2014JB011502

CrossRef Full Text | Google Scholar

Troise, C., De Natale, G., Pingue, F., Obrizzo, F., De Martino, P., Tammaro, U., et al. (2007). Renewed ground uplift at Campi Flegrei caldera (Italy): new insight on magmatic processes and forecast. Geophys. Res. Lett. 34, L03301. doi:10.1029/2006GL028545

CrossRef Full Text | Google Scholar

Walwer, D., Ghil, M., and Calais, E. (2019). Oscillatory nature of the Okmok volcano’s deformation. Earth Planet Sci. Lett. 506, 76–86. doi:10.1016/j.epsl.2018.10.033

CrossRef Full Text | Google Scholar

Weiss, J., Walters, R., Morishita, Y., Wright, T., Lazecky, M., Wang, H., et al. (2020). High-resolution surface velocities and strain for Anatolia from Sentinel-1 InSAR and GNSS data. Geophys. Res. Lett. 47, e2020GL087376. doi:10.1029/2020GL087376

CrossRef Full Text | Google Scholar

Wicks, C. W. (2002). Magmatic activity beneath the quiescent Three Sisters volcanic center, central Oregon cascade range, United States. Geophys. Res. Lett. 29, 26-1–26-4. doi:10.1029/2001GL014205

CrossRef Full Text | Google Scholar

Yang, X.-M., Davis, P. M., and Dieterich, J. H. (1988). Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing. J. Geophys. Res.: Solid Earth 93, 4249–4257. doi:10.1029/JB093iB05p04249

CrossRef Full Text | Google Scholar

Zurek, J., Williams-Jones, G., Johnson, D., and Eggers, A. (2012). Constraining volcanic inflation at Three Sisters volcanic field in Oregon, United States, through microgravity and deformation modeling. G-cubed 13, Q10013. doi:10.1029/2012GC004341

CrossRef Full Text | Google Scholar