Abstract
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.
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 (). 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. 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 (; ). Calderas usually show long repose periods between large eruptions, but without quiescence (e.g., resurgent volcanism, subsidence, multiple pulses of uplift) (). 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) () 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 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 (). 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 (). All these circumstances influence on the eruptive history of Three Sisters.
FIGURE 1
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 (). North Sister is a monotonously mafic edifice created 120 ka ago, formed by long-lived effusive volcanism (); 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 (). However, most of the volcanic activity is identified with mafic shields and cones around the major composite volcanoes (). Geochemical anomalies suggest that episodes of intrusion may be more frequent in the Three Sisters area than the age of eruptive vents would involve (). The most recent eruptions were rhyolitic close to South Sisters, 2000 years ago (). 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 ().
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 (). 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 (; ). The uplift at Three Sisters has been aseismic except for a swarm of ∼300 small earthquakes (Mmax = 1.9) in March 2004 (; ). Previous studies (; ; ; ) 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 (). 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 (). 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 (; ; 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
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 () and unwrapped using SNAPHU (), with perpendicular baselines up to 500 m, as explained by .
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 (), 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 (). 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 (), 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 (; ). 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 (). 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 (). 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
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 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 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) (), which estimates source parameters through a Markov chain Monte Carlo method and uses, among others, analytical forward models from dMODELS software package (). 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 (), prolate spheroid () and sill-like () 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 (, ; ), 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 (). InSAR data are subsampled with an adaptative quadtree method () to reduce the computational cost of the Bayesian inversion, particularly for sill-like and prolate spheroid sources. The inversion computes 2×106 iterations for spherical point source, 5×106 iterations for sill-like source and 8×106 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, ) 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, , records 1) the incremental volume change between two acquisition dates, ti and tj from an interferogram or 2) the cumulative increment of volume since the first observation, i.e., tj = t0, being t0 the starting date of CGPS.
A linear inversion scheme using Weighted Least Squares (WLS, ) 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 () 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, , given the set of interferograms and 5,008 cumulative GPS displacements.
The least square estimator of each inversion, , is given by:with the cofactor matrix . We considered a diagonal variance-covariance matrix, Cd, 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 (; ).
4.2.2 Second Step: Volume Changes Time-series
We want to solve for the temporal evolution of volume change for each observed epoch tk from 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, ). This prevents the presence of large discontinuities in the final solution.
Let be the data vector of volume changes over the corresponding time intervals (N × 1), and the unknown vector of volume change rates (M × 1) between adjacent epochs, tj where the overdot means differentiation over time. Then,The usual method of converting the observations 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 tj 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 (tslavei) and master () epochs of each i time interval, for i = 1, ..., N. Therefore, the (i, j) component of the design matrix is Bij = (tj+1 − tj) for , and zero elsewhere. In the case of cumulative (i.e., continuous GPS), B presents lower triangular matrix blocks. For example, if volume changes are obtained over different time intervals, i.e., if tAB, tBC, and tAC (from InSAR data), and tCD, tCE, tCF, (from CGPS data) are , , , , and , the design matrix is given by:To illustrate the simple example in Eq. 4, let all dates, tA, tB, tC, tD, tE, and tF be equally spaced at time intervals of 1 year, i.e., tAB = 1 year and so on, and ×106 m3. In this case, standard least squares can be applied, given . The cumulative volume times series is then , meaning a linear inflation rate of in time interval tAB and a posterior linear inflation rate of .
However, the set of 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, , is not possible since Eq. 4 constitutes an ill-posed unstable model, with one or more eigenvalues of the normal matrix close to zero. This fact is responsible for large uncertainty on the estimated volume change rates, . Alternatives to the least square method can be proposed for an improved estimate of : Tikhonov regularization (), Bayesian and stochastic inferences (), Truncated Singular Value Decomposition (TSVD) or Principal Components (; ; ). Here, we consider TSVD as proposed by .
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 (). Furthermore, we included some estimations of the error of data () to establish some uncertainty in the estimator. To estimate the uncertainties, the Weighted Generalized Inverse method () permitted the use of the “a priori” information from the data (and optional model, CV = ATPA) 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 CV = I. is obtained by error data propagation through 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, , and Bnew = DBS−1 give the transformation of data, model parameters and forward operator in the new system of coordinates. TSVD is applied to Bnew with a specific regularization method to find the Principal Components of the observation set (). 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:
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) (; ; ). Methods like GCV sometimes fail to find the appropriate regularization parameter (flat local minima), whereas the L-curve gives a robust estimation () 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 (; ; ):where is the regularization error, and being the exact and truncated SVD solutions; P is the regularization parameter value, si are the singular values, and are called Fourier coefficients ( are data and ui 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 si.
The L-curve method is applied to resulting in turn from applying (Eq. 7) through a log-log plot of the norm of a regularized solution vs. the norm of the corresponding residual norm . As recommended by 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 (7–13)×106 m3, 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 5
TABLE 1
| Lower | Upper | Step | Start | |
|---|---|---|---|---|
| Mogi | ||||
| Xcenter (m) | −1.00 × 103 | 1.50 × 104 | 25 | 2.00 × 103 |
| Ycenter (m) | −1.00 × 103 | 1.50 × 104 | 25 | 2.00 × 103 |
| Depth (m) | 5 × 102 | 2.00 × 104 | 50 | 1.00 × 103 |
| 0.1 | 1.00 × 104 | 1.00 × 10−3 | 0.1 | |
| Sill | ||||
| Xcenter (m) | −1.00 × 103 | 1.50 × 104 | 25 | 2.00 × 103 |
| Ycenter (m) | −1.00 × 103 | 1.50 × 104 | 25 | 2.00 × 103 |
| Depth (m) | 5 × 102 | 2.00 × 104 | 50 | 1.00 × 103 |
| Radius (m) | 100 | 4,000 | 50 | 100 |
| 1 × 10−5 | 10 | 1 × 10−6 | 1 × 10−2 | |
| Spheroid | ||||
| Xcenter (m) | −1.00 × 103 | 1.50 × 104 | 25 | 2.00 × 103 |
| Ycenter (m) | −1.00 × 103 | 1.50 × 104 | 25 | 2.00 × 103 |
| Depth (m) | 5 × 102 | 2.00 × 104 | 50 | 1.00 × 103 |
| Major semi axis (m) | 100 | 4,000 | 50 | 100 |
| Aspect ratio | 0.2 | 0.99 | 0.01 | 0.5 |
| 1 × 10−5 | 10 | 1 × 10−6 | 1 × 10−2 | |
| Strike (°) | 1 | 359 | 1 | 270 |
| Plunge (°) | 0 | 90 | 1 | 45 |
Prior information for the Elastic half-space Spherical Point Pressure (), Penny-shaped sill-like () and Spheroid () deformation sources.
TABLE 2
| median | 5% | 95% | |
|---|---|---|---|
| Mogi | |||
| Lon (°) | −121.8382 | −121.8418 | −121.8350 |
| Lat (°) | 44.1055 | 44.1030 | 44.1082 |
| Depth (m) | 5,000 | 4,500 | 6,000 |
| 8.99 | 6.98 | 12.66 | |
| Sill | |||
| Lon (°) | −121.8368 | −121.8398 | −121.8339 |
| Lat (°) | 44.1018 | 44.0997 | 44.1041 |
| Depth (m) | 6,300 | 5,600 | 7,200 |
| Radius (m) | 250 | 220 | 400 |
| 0.23 | 0.05 | 0.30 | |
| Spheroid | |||
| Lon (°) | −121.8707 | −121.8791 | −121.8616 |
| Lat (°) | 44.1187 | 44.1130 | 44.1228 |
| Depth (m) | 6,100 | 5,300 | 7,400 |
| Major semi-axis (m) | 400 | 240 | 720 |
| Aspect ratio | 0.28 | 0.22 | 0.37 |
| 1.85 | 0.38 | 8.61 | |
| Strike (°) | 68 | 49 | 102 |
| Plunge (°) | 82 | 78 | 85 |
Bayesian inversion results, with the median a posterior probability solution and the 95% credible intervals, for the Elastic half-space Spherical Point Pressure (), Penny-shaped sill-like () and Spheroid () deformation sources.
Geo-reference point is [−121.9, 44.03]°.
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: dmedian = 5000 m, d5% = 4500 m and d95% = 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, , 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 (dmedian = 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, , are more uniform for the CGPS data sets (Figure 7), but more variable for the individual SAR data sets (Figure 6).
FIGURE 6
FIGURE 7
Finally, we applied the Picard Plot condition suitable to understand the conditions of the ill-posed problem (Eq. 3). Figure 8A shows how si only decays faster than the Fourier coefficients () 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
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 during 1999–2001 and a subsequent rate as much as 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 10
TABLE 3
| Inversion | Depth (m) | ||
|---|---|---|---|
| 1996–2000, InSARa | 6500 ± 400 | 23.00 ± 3.00 | ∼5.75 ± 0.75 |
| 2001–2008, Geodetic ground base datab | 5,800 | 22.20 | ∼3.14 |
| 1993–2008, InSARc | 5200 ± 100 | 57.00 [+1.95, −3.60] | ∼3.80 [+0.13, −0.24] |
| 1997–2000, InSAR [this study]d | 5000 [+1000, −500] | 8.99 [+3.67, −2.01] | ∼3.00 [+1.22, −0.67] |
Source location comparison for Three Sisters (assuming a Mogi source) from Previous Studies and the bayesian inversion carried out in this Study.
ERS Descending (08/1996–10/2000) ().
Campaign GPS, CGPS and Leveling (05/2001-late 2008) ().
ERS Track 385 Descending (08/1993-08/2008) ().
ERS Track 385 Descending (08/1997-09/2000).
TABLE 4
| Value | c (yr−1) | t0 (yr) | R–square | RMSE | |
|---|---|---|---|---|---|
| Optimal | 0.1055 | 9.48 | 1999.09 | 0.989 | 0.60 × 106 |
| Lower bound | 0.1042 | 9.60 | 1998.75 | 0.998 | 0.18 × 106 |
| Upper bound | 0.1074 | 9.31 | 1999.59 | 0.988 | 0.87 × 106 |
e-folding parameters obtained by NonLinear Least Squares, using the Least Absolute Residual method (LAR).
Upper and lower bounds for the curve fitting parameters, encompassing the results for the volume changes series related to the 5%, median and 95% percentile values of source depth estimate.
6 Discussion
6.1 Source Characterization
Studies at Three Sisters using InSAR interferometric pairs and stacks (; ), GPS (campaign and continuous) and leveling () 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 and , 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. 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 (). 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). processed three interferograms, obtaining an increment of volume change of 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 , to fit the same surface deformation. Although the time acquisition of the best interferogram of spans only 1 year more than our InSAR data interferometric pair, the important difference between our and their is mainly due to their depth estimation. Our InSAR interferogram matches one of the other two by . 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 (, Eq. 1) and the corresponding simulated observations ( 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 associated with the InSAR data display a wide dispersion, and the associated with GPS data acquire discrete values, i.e., the same 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. 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, is the characteristic relaxation time constant, here after named e-folding parameter, and t0 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 t0 parameters.
Our time-series spans 26 years and presents a characteristic time constant of years. This new result updates the value of 5.3 years obtained by from an 8-years time-series of GPS displacements from 2001 to late 2008. 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×106 m3. As late as January 2020, our inflation time-series indicated that the cumulative volume of a spherical point source with dmedian = 5000 m is 29.1×106 m3. For d5% = 4500 m and d5% = 6000 m, the values are 22.7×106 and 39.9×106 m3, respectively. This range of values is consistent with those predicted by for a prolate spheroid model, 44.9 to 51.6×106 m3 (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 ().
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 (), Long Valley (), Uturuncu (), Laguna del Maule (), Yellowstone (), Campi Flegrei (), Santorini (), Alutu (), Agung (). 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 t0 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
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 (). 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 (). 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 (). Santorini presented a source inflation process with an e-folding of 0.28 years for the period October 2011- August 2012 (). Then, its subsidence rates increased in the post-unrest period 2012–2017, suggesting the superimposition of various deformation sources (). 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 (), the general trend can be modeled as an exponential decay. After this inter-eruptive episode, a phreatomagmatic eruption occurred in July–August 2008 (). 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 (). 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 . 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′):where y0 is the displacement or volume change at t0, is the value after total relaxation (i.e., at ), and t0 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., ; ; ). 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 (; ; ; ). 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 (). Recent rhyolitic eruptions close to South Sisters, 2000 years ago (), and the evidence of a long-lived magmatic source at depth in Three Sisters, from spring geochemistry studies (), 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, 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 ∼1018 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., ; ; ; ).
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 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 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.
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.
Statements
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.
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).
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.
References
1
AcocellaV.Di LorenzoR.NewhallC.ScandoneR. (2015). An overview of recent (1988 to 2014) caldera unrest: knowledge and perspectives. Rev. Geophys.53, 896–955. 10.1002/2015RG000492
2
BackusG. E. (1988). Bayesian inference in geomagnetism. Geophys. J. Int.92, 125–142. 10.1111/j.1365-246X.1988.tb01127.x
3
BagnardiM.HooperA. (2018). Inversion of surface deformation data for rapid estimates of source parameters and uncertainties: a bayesian approach. G-cubed19, 2194–2211. 10.1029/2018GC007585
4
BattagliaM.CervelliP.MurrayJ. (2013). DMODELS: a MATLAB software package for modeling crustal deformation near active faults and volcanic centers. J. Volcanol. Geoth. Res.254, 1–4. 10.1016/j.jvolgeores.2012.12.018
5
BekaertD. P. S.HamlingtonB. D.BuzzangaB.JonesC. E. (2017). Spaceborne synthetic aperture radar survey of subsidence in Hampton roads, Virginia (United States). Sci. Rep.7, 14752. 10.1038/s41598-017-15309-5
6
BerardinoP.FornaroG.LanariR.SansostiE. (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. 10.1109/IGARSS.2002.1025900
7
BiggsJ.PritchardM. (2017). Global volcano monitoring: what does it mean when volcanoes deform?Elements.13, 17–22. 10.2113/gselements.13.1.17
8
BiggsJ.LuZ.FournierT.FreymuellerJ. 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 Earth115, 1–11. 10.1029/2010JB007577
9
BonafedeM.FerrariC. (2009). Analytical models of deformation and residual gravity changes due to a Mogi source in viscoelastic medium. Tectonophysics471, 4–13. 10.1016/j.tecto.2008.10.006
10
ChenC. W.ZebkerH. A. (2002). Phase unwrapping for large SAR interferograms: statistical segmentation and generalized network models. IEEE Trans. Geosci. Rem. Sens.40, 1709–1719. 10.1109/TGRS.2002.802453
11
DecriemJ.ÁrnadóttirT.HooperA.GeirssonH.SigmundssonF.KeidingM.et al (2010). The 2008 May 29 earthquake doublet in SW Iceland. Geophys. J. Internat.—Geophys. J. Int.181, 1128–1146. 10.1111/j.1365-246X.2010.04565.x
12
Del NegroC.CurrentiG.ScanduraD. (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. 10.1016/j.pepi.2008.10.019
13
DieterichJ. H.DeckerR. W. (1975). Finite element modeling of surface deformation associated with volcanism. J. Geophys. Res. (1896–1977).80, 4094–4102. 10.1029/JB080i029p04094
14
DzurisinD.LisowskiM.WicksC. 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. 10.1007/s00445-009-0296-4
15
DzurisinD.LisowskiM.WicksC. W.PolandM. P.EndoE. 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. 10.1016/j.jvolgeores.2005.07.011
16
DzurisinD.LisowskiM.WicksC. W. (2017). Semipermanent GPS (SPGPS) as a volcano monitoring tool: rationale, method, and applications. J. Volcanol. Geoth. Res.344, 40–51. 10.1016/j.jvolgeores.2017.03.007
17
EvansW.van SoestM.MarinerR.HurwitzS.IngebritsenS.WicksC.et al (2004). Magmatic intrusion west of Three Sisters, central Oregon, United States: the perspective from spring geochemistry. Geology32, 69–72. 10.1130/G19974.1
18
FarrT. G.RosenP. A.CaroE.CrippenR.DurenR.HensleyS.et al (2007). The shuttle radar topography mission. Rev. Geophys.45, RG2004. 10.1029/2005RG000183
19
FialkoY.KhazanY.SimonsM. (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. 10.1046/j.1365-246X.2001.00452.x
20
FiersteinJ.HildrethW.CalvertA. T. (2011). Eruptive history of South sister, Oregon cascades. J. Volcanol. Geoth. Res.207, 145–179. 10.1016/j.jvolgeores.2011.06.003
21
GalettoF.AcocellaV.CaricchiL. (2017). Caldera resurgence driven by magma viscosity contrasts. Nat. Commun.8, 1750. 10.1038/s41467-017-01632-y
22
GonzálezP. J.SamsonovS. V.PepeS.TiampoK. F.TizzaniP.CasuF.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 Earth118, 4361–4377. 10.1002/jgrb.50289
23
HansenP. C. (1992). Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev.34, 561–580. 10.1137/1034115
24
HansenP. C. (2007). Regularization tools version 4.0 for matlab 7.3. Numer. Algorithm.46, 189–194. 10.1007/s11075-007-9136-9
25
HansenP. C. (1990). The discrete picard condition for discrete ill-posed problems. Bit.30, 658–672. 10.1007/BF01933214
26
HansenP. C.O’LearyD. P. (1993). The use of the L-curve in the regularization of discrete ill-posed problems. SIAM J. Sci. Comput.14, 1487–1503. 10.1137/0914086
27
HildrethW. (2007). Quaternary magmatism in the cascades-geologic perspectives. USGS. Professional PaperVol. 1744, 125.
28
HildrethW.FiersteinJ.CalvertA. T. (2012). Geologic map of three sisters volcanic cluster, cascade range, Oregon. US Geological Survey Scientific Investigations Map, Vol. 3186, pamphlet107. p., 2 sheets, scale 1:24,000. Available at https://pubs.usgs.gov/sim/3186/.
29
HillD.Montgomery-BrownE.ShellyD.FlindersA.PrejeanS. (2020). Post-1978 tumescence at Long Valley caldera, California: a geophysical perspective. J. Volcanol. Geoth. Res.400, 106900. 10.1016/j.jvolgeores.2020.106900
30
HooperA.BekaertD.SpaansK.ArikanM. (2012). Recent advances in SAR interferometry time series analysis for measuring crustal deformation. Tectonophysics514–517, 1–13. 10.1016/j.tecto.2011.10.013
31
HossainaliM.BeckerM.GrotenE. (2010). Comprhensive approach to the analysis of the 3D kinematics deformation with appliction to the kenai Peninsula. J. Geodetic Sci.1, 59–73. 10.2478/v10156-010-0008-1
32
HutchisonW.BiggsJ.MatherT. A.PyleD. M.LewiE.YirguG.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. 10.1002/2016GC006395
33
LarsenJ.NealC.SchaeferJ.KaufmanA.LuZ. (2015). The 2008 phreatomagmatIc eruption of Okmok volcano, AleutIan Islands, Alaska: chronology, deposits, and landform changes. Bull. Volcanol.77, 69. 10.1007/s00445-015-0952-9
34
LauN.TymofyeyevaE.FialkoY. (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. 10.1016/j.epsl.2018.03.026
35
LawlessJ.WangP. (1976). A simulation study of ridge and other regression estimators. Commun. Statist. Theor. Meth, A5, 307–323. 10.1080/03610927608827353
36
Le MévelH.FeiglK. L.CórdovaL.DeMetsC.LundgrenP. (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. 10.1002/2015GL064665
37
LenglinéO.MarsanD.GotJ.-L.PinelV.FerrazziniV.OkuboP. G. (2008). Seismicity and deformation induced by magma accumulation at three basaltic volcanoes. J. Geophys. Res.: Solid Earth113. B12305. 10.1029/2008JB005937
38
LohmanR. B.SimonsM. (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. 10.1029/2004GC000841
39
MccaffreyR.QamarA.KingR.WellsR.KhazaradzeG.WilliamsC.et al (2007). Fault locking, block rotation and crustal deformation in the Pacific Northwest. Geophys. J. Internat. – Geophys. J. Iint.169, 1315–1340. 10.1111/j.1365-246X.2007.03371.x
40
MenkeW. (1989). Geophysiocal data analysis: discrete inverse theory. San Diego, Calif: Academic.
41
MogiK. (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.
42
MoranS. C. (2004). Seismic monitoring at cascade volcanic centers, 2004- status and recommendationsc. US Geological Survey Scientific Investigations Report, 2004-5211, 28.
43
NewmanA.DixonT.OfoegbuG.DixonJ. (2001). Geodetic and seismic constraints on recent activity at Long Valley Caldera, California: evidence for viscoelastic rheology. J. Volcanol. Geoth. Res.105, 183–206. 10.1016/S0377-0273(00)00255-9
44
PapageorgiouE.FoumelisM.TrasattiE.GuidoV.RaucoulesD.MouratidisA. (2019). Multi-sensor SAR geodetic imaging and modelling of Santorini volcano post-unrest response. Rem. Sens.11, 259. 10.3390/rs11030259
45
ParksM.MooreJ.PapanikolaouX.BiggsJ.MatherT.PyleD.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. 10.1002/2014JB011540
46
PolandM. P.LisowskiM.DzurisinD.KramerR.McLayM.PaukB. (2017). Volcano geodesy in the cascade arc, United States. Bull. Volcanol.79. 10.1007/s00445-017-1140-x
47
ReversoT.VandemeulebrouckJ.JouanneF.PinelV.VilleminT.SturkellE.et al (2014). A two-magma chamber model as a source of deformation at Grímsvötn Volcano, Iceland. J. Geophys. Res.: Solid Earth119, 4666–4683. 10.1002/2013JB010569
48
RiddickS. N.SchmidtD. A. (2011). Time-dependent changes in volcanic inflation rate near Three Sisters, Oregon, revealed by InSAR. G-cubed.12, Q12005. 10.1029/2011GC003826
49
RosenP. A.GurrolaE.SaccoG. F.ZebkerH. (2012). “The InSAR scientific computing environment,” in EUSAR 2012; 9th European conference on synthetic aperture radar, 23-26 April 2012, Nuremberg, Germany730–733.
50
RosenP.HenleyS.PeltzerG.SimonsM. (2004). Updated repeat orbit interferometry package released. Eos, Trans. American Geophys. Union.85, 47. 10.1029/2004EO050004
51
SchmidtM.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. 10.1093/petrology/egq094
52
ScottW. E.IversonR.SchillingS.FisherB. (2001). Volcano hazards in the three sisters region, Oregon. US Geological Survey Open-File Report, Vol. 18, 99–437. 10.3133/ofr99437
53
SegallP. (2016). Repressurization following eruption from a magma chamber with a viscoelastic aureole. J. Geophys. Res.: Solid Earth121, 8501–8522. 10.1002/2016JB013597
54
SherrodD. R.TaylorE. M.FernsM. L.ScottW. E.ConreyR. M.SmithG. A. (2004). Geologic map of the Bend 30- × 60-minute quadrangle, central Oregon. US Geological Survey Geological Investigations Series Map I-2683, 48.
55
SyahbanaD.KasbaniK.SuantikaG.PrambadaO.AndreasA.SaingU.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. 10.1038/s41598-019-45295-9
56
TikhonovA.ArseninV. (1977). Solutions of ill-posed problem. New York: John Wiley & Sons.
57
TizzaniP.CastaldoR.PepeA.ZeniG.LanariR.BattagliaM. (2015). Magma and fluid migration at Yellowstone Caldera in the last three decades inferred from InSAR, leveling and gravity measurements. J. Geophys. Res.: Solid Earth120, 2627–2647. 10.1002/2014JB011502
58
TroiseC.De NataleG.PingueF.ObrizzoF.De MartinoP.TammaroU.et al (2007). Renewed ground uplift at Campi Flegrei caldera (Italy): new insight on magmatic processes and forecast. Geophys. Res. Lett.34, L03301. 10.1029/2006GL028545
59
WalwerD.GhilM.CalaisE. (2019). Oscillatory nature of the Okmok volcano’s deformation. Earth Planet Sci. Lett.506, 76–86. 10.1016/j.epsl.2018.10.033
60
WeissJ.WaltersR.MorishitaY.WrightT.LazeckyM.WangH.et al (2020). High-resolution surface velocities and strain for Anatolia from Sentinel-1 InSAR and GNSS data. Geophys. Res. Lett.47, e2020GL087376. 10.1029/2020GL087376
61
WicksC. 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. 10.1029/2001GL014205
62
YangX.-M.DavisP. M.DieterichJ. 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 Earth93, 4249–4257. 10.1029/JB093iB05p04249
63
ZurekJ.Williams-JonesG.JohnsonD.EggersA. (2012). Constraining volcanic inflation at Three Sisters volcanic field in Oregon, United States, through microgravity and deformation modeling. G-cubed13, Q10013. 10.1029/2012GC004341
Summary
Keywords
inter-eruptive deformation, characteristic relaxation time, continuous GPS, interferometric synthetic aperture radar, geodetic time series, Three Sisters volcano
Citation
Rodríguez-Molina S, González PJ, Charco M, Negredo AM and Schmidt DA (2021) Time-Scales of Inter-Eruptive Volcano Uplift Signals: Three Sisters Volcanic Center, Oregon (United States). Front. Earth Sci. 8:577588. doi: 10.3389/feart.2020.577588
Received
29 June 2020
Accepted
23 November 2020
Published
21 January 2021
Volume
8 - 2020
Edited by
Flavio Cannavo’, National Institute of Geophysics and Volcanology, Section of Catania, Italy
Reviewed by
Maurizio Battaglia, United States Geological Survey (USGS), United States
Gilda Maria Currenti, Istituto Nazionale di Geofisica e Vulcanologia (INGV), Italy
Updates
Copyright
© 2021 Rodríguez-Molina , González, Charco, Negredo and Schmidt.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Pablo J. González, pjgonzal@liverpool.ac.uk
This article was submitted to Volcanology, a section of the journal Frontiers in Earth Science
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.