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 ∼1018 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.

A classical inflation-eruption-deflation cycle of a volcano is useful to conceptualize the timeevolving 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 intereruptive 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 18 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.

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 intereruptive period.

BACKGROUND
The Three Sisters Volcanic Complex forms a N-S chain of stratovolcanoes located at 44.1 + 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.
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 rhyoliticintermediate 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. Microgravity data collected between 2002 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.

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.

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.

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 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 summerto-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)   (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.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 577588 5 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 timeseries. 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-tonoise 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).

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 twostep 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.

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 halfspace is assumed in all the approaches with Poisson's ratio 0.25. We assumed, as previous studies (Dzurisin et al., 2006(Dzurisin et al., , 2009Riddick 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 6 iterations for spherical point source, 5×10 6 iterations for sill-like source and 8×10 6 iterations for prolate spheroid, which stabilizes the inversion procedure.

Temporal Evolution of the Source Volume Changes
Once the magmatic source is fixed, we perform the quasidynamic 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.

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 Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 577588 intervals (increments of volume changes, Δ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, Δ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 0 , being t 0 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 timeseries) 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 (Δ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, ΔV ij , given the set of interferograms and 5,008 cumulative GPS displacements.
The least square estimator of each inversion, m, is given by: 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).

Second Step: Volume Changes Time-series
We want to solve for the temporal evolution of volume change for each observed epoch t k from Δ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 ΔV be the data vector of volume changes over the corresponding time intervals (N × 1), and _ V the unknown vector of volume change rates (M × 1) between adjacent epochs, t j where the overdot means differentiation over time. Then, The usual method of converting the observations ΔV on volume change rates is: 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 slave i ) and master (t masteri ) 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 slave i ≤ t j < t master i , and zero elsewhere. In the case of cumulative Δ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 Δv AB , Δv BC , Δv AC , Δv CD , Δv CE and Δv CF , the design matrix is given by: 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 ΔV [2, 1, 3, 1, 2, 3] ×10 6 m 3 . In this case, standard least squares can be applied, given _ V [2, 1, 1, 1, 1] × 10 6 m 3 /year. The cumulative volume times series is then V [2, 3, 4, 5, 6] × 10 6 m 3 , meaning a linear inflation rate of 2 × 10 6 m 3 /year in time interval t AB and a posterior linear inflation rate of 1 × 10 6 m 3 /year.
However, the set of Δ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, _ V (B T PB) − 1 B T ΔV, is not possible since Eq. 4 constitutes an ill-posed unstable model, with one or more eigenvalues of the normal matrix B T PB close to zero. This fact is responsible for large uncertainty on the estimated volume change rates, _ V. Alternatives to the least square method can be proposed for an improved estimate of _ 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).

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 (ΔV) to establish some uncertainty in the _ V estimator. To estimate the uncertainties, the Weighted Generalized Inverse method (Menke, 1989) permitted the use of the "a priori" information from the data C ΔV (and optional model, C _ V A T PA) covariance matrix. Such matrices can be decomposed as: 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.C ΔV is obtained by error data propagation through Δ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, ΔV new DΔV, _ V new S _ V and B new DBS −1 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 (Δ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:

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): where _ V − _ V T1 2 is the regularization error, _ V and _ V T1 being the exact and truncated SVD solutions; P is the regularization parameter value, s i are the singular values, and u T i ΔV i are called Fourier coefficients (Δ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 _ V T 1 resulting in turn from applying (Eq. 7) through a log-log plot of the norm of a regularized solution _ V T 2 2 vs. the norm of the corresponding residual norm B _ V T 2 − ΔV 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.

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 ΔV (7-13)×10 6 m 3 , 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 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 bimodality (slightly unstable inversion result). For the prolate spheroid source, the aspect ratio between semi-axes and the Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 577588 9 excess pressure PDFs display also a bi-modal shape. The major semi-axis PDF exhibits multi-modality.

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, Δ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, ΔV i , are more uniform for the CGPS data sets (Figure 7), but more variable for the individual SAR data sets ( Figure 6).
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 ( u T i ΔV i ) 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).
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 ∼ 1.60 × 10 6 m 3 /yr during 1999-2001 and a subsequent rate as much as ∼ 0.75 × 10 6 m 3 /yr for the period 2015-January 2020. Data since   (Mogi, 1958), Penny-shaped sill-like (Fialko et al., 2001) and Spheroid (Yang et al., 1988) deformation sources.

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   (Mogi, 1958), Penny-shaped sill-like (Fialko et al., 2001) and Spheroid (Yang et al., 1988)  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 ΔV 23 × 10 6 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 Δ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 ΔV rate ∼ 3.00 × 10 6 m 3 /yr and their ΔV rate ∼ 5.75 × 10 6 m 3 /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.

Time Series of Volume Changes: Regularization Using the Truncated SVD
To assess the effect of the regularization, we compared the increments of volume change (ΔV i , Eq. 1) and the corresponding simulated observations ( ΔV ireconstruction  V (i+1) inversion − V iinversion 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 ΔV i associated with the InSAR data display a wide dispersion, and the ΔV i associated with GPS data acquire discrete values, i.e., the same Δ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).

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: where a, b and d are constants, 1/c ϵ is the characteristic relaxation time constant, here after named e-folding parameter, and t 0 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    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 timeseries 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 6 m 3 . 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 6 m 3 . For d 5% 4500 m and d 5% 6000 m, the values are 22.7×10 6 and 39.9×10 6 m 3 , 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 6 m 3 (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 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).
Frontiers in Earth Science | www.frontiersin.org 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 0 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).
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 noneruptive 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′): where y 0 is the displacement or volume change at t 0 , y t∞ is the value after total relaxation (i.e., at t t ∞ ), and t 0 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][2003][2004][2005][2006][2007][2008][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 18 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).

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 threecomponents 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 [ ± 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.

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 cowrote the manuscript. AN co-supervised SR-M and contributed to the manuscript. DS contributed to the 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).