Unrest at the Laguna del Maule volcanic field 2005–2020: renewed acceleration of deformation

The Laguna del Maule volcanic field in Chile has been exhibiting unrest since 2005. New GPS and InSAR data reveal a second episode of accelerated deformation beginning in late 2016 and continuing through May 2020, with an uplift rate > 290 mm/year between 2019 and 2020. To explain the spatial and temporal pattern of deformation, we apply a dynamic model of viscous magma flowing through a conduit into a fluid-filled reservoir surrounded by a heterogeneous, viscoelastic crust. A Monte Carlo procedure optimizes the ellipsoid reservoir geometry and the inlet pressure history. The two episodes of accelerating uplift are each modeled with a pressure increase rate of ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\sim $\end{document}9 MPa/year. Since 2016, 0.10 km3 of magma was injected into the system for a total of 0.37 km3 since 2005.


Introduction
The volcanic field at Laguna del Maule (LdM) in the southern Andes exhibits recent silicic volcanism and ongoing, rapid deformation, but no historical eruption Andersen et al. 2019). It has been uplifting faster than 200 mm/year since 2007 Le Mével et al. 2015;Le Mével et al. 2016). Episodes of unrest characterized by ground uplift are commonly observed at large silicic caldera systems (Sandri et al. 2017), as illustrated by the recent unrest at Domuyo volcano in Argentina (Astort et al. 2019;Lundgren et al. 2020), Long Valley caldera (Montgomery-Brown et al. 2015), and Yellowstone caldera (Chang et al. 2010). Episodes of uplift without subsequent eruptions have also been observed at mafic calderas such as Alcedo volcano in the Galapagos  (Galetto et al. 2019). Surface uplift signals pressurization of the underlying magma reservoir(s) which can result from emplacement of new magma and/or perturbations to pore-fluid pressure in the associated hydrothermal system (e.g., Mogi 1958;Gottsmann et al. 2006;Dzurisin 2007;Troiano et al. 2011;Coco et al. 2016). Surface uplift can also presage eruptive activity (Sigmundsson et al. 2010;Peltier et al. 2008). One of the outstanding challenges in interpreting measurements of volcanic deformation thus lies in determining when a pre-eruptive stage is reached. Since the long, multi-sensor time series of geodetic measurements at LdM includes episodes of quiescence, acceleration, and deceleration, this data set provides a unique opportunity to study magma emplacement at long-lived magmatic systems during an inter-eruptive time interval.
The LdM volcanic field sits in the Southern Volcanic Zone of Chile (Fig. 1). Previous InSAR and GPS studies have revealed three phases of ground deformation for the ongoing geodetic unrest at LdM: (1) No displacement was observed between 2003 and 2004 (Fournier et al. 2010); (2) Ground uplift accelerated from some time before 2007 through 2010 with rates exceeding 200 mm/year Le Mével et al. 2015), and (3) Uplift continued with slowly decelerating rates from 2010 to 2016 (Le Mével et al. 2016). The deformation affects an area of ∼500 km 2 centered on the western shore of the lake basin. This area includes ∼40 km 3 of rhyolite erupted since the last glaciation (ca. 20 ka) and as late as ca. 2000 ka from at least Antuco -SAR frame outlines (color-coded rectangles with track number), and MAUL continuous GPS station used as a reference for the GPS time series (white circle). Lower-right inset shows location of LdM (red dot) within South America 24 vents distributed around the lake basin Andersen et al. 2019). The source of uplift has been modeled as an opening sill at a depth of 5 km . Le Mével et al. (2016) introduce analytic solutions for a dynamic model of magma injection into a spheroid chamber to explain nonlinear evolution of uplift observed in the displacement time series.
Previous studies have imaged the magmatic system at LdM, including: a low-density body at a depth of 2 to 4 km sitting directly above the modeled deformation source (Miller et al. 2017), a conductor extending laterally northward and downward to a depth of 15 km (Cordell et al. 2018), and a low seismic velocity zone of ∼500 km 3 (Wespestad et al. 2019;Bai et al. 2020). A network of active normal faults striking NE-SW was also revealed onshore and in the lake sediments Garibaldi et al. 2020). Petrologic and geochemical findings at LdM support the physical model of volatile-rich magma injection into an existing reservoir (Klug et al. 2020;Andersen et al. 2018;Andersen et al. 2019).
In this study, we analyze geodetic data at LdM to calculate time series of displacement from 2007 to 2020 and model them in terms of magma input to the reservoir.

GPS
Five continuously recording GPS stations have been installed since 2012 (Fig. 2). Precise point-positioning (Zumberge et al. 1997) using GIPSY-OASIS (release 6.4) provides daily estimates of station coordinates referred to the stable South America reference frame with uncertainty better than 1 mm. To remove correlated regional noise and leftover regional postseismic deformation from the 2008 Mw 8.8 Maule earthquake, we refer each site's relative position to the continuous station MAUL, located ∼30 km NW of the study area (Fig. 1). For each GPS station, we calculate the mean velocity between February 2019 and February 2020. Figure 2 shows uplift and radial motion outward from a point on the western shore of the lake basin. Before 2016, the uplift rate was slowly decreasing (Le Mével et al. 2016). Since then, however, the new GPS results show an increase in velocity that is still ongoing as of June 2020 ( Fig. 2 and S1). The GPS station MAU2 is moving upward at 291 ± 2 mm/year between 2019 and 2020, the fastest uplift rate measured since the installation of the GPS network in 2012.

InSAR
All available SAR images covering LdM from 2007 to 2020 ( Fig. 1) acquired by ALOS-1 PALSAR (tracks 112 and 113), TerraSAR-X/TanDEM-X (tracks 28, 111, and 119), and SENTINEL-1 Copernicus (paths 18 and 83) are analyzed using ISCE (v2.0; Rosen et al. (2012). The resulting interferograms map the displacement along the line of sight (LOS) between target pixels on the ground and the radar sensor aboard the satellite. These maps are multilooked, unwrapped, and coregistered to a common SAR acquisition date (Table S1). Interferograms with average temporal coherence lower than 0.7 are excluded. The phase delay due to a stratified troposphere is modeled and removed using the atmospheric model ERA5 from ECMWF and the PyAPS software (Jolivet et al. 2011;Jolivet et al. 2014). A bi-linear (planar) phase ramp is estimated and removed. The resulting seven stacks are then inverted to calculate a time series of displacement using the weighted, network-inversion approach of Yunjun et al. (2019) implemented in MintPy. We estimate the uncertainty for the displacement at each epoch as the sample standard deviation over a 3-by-3-pixel area around the pixel of interest. The maximum uncertainty for each track varies between 1.3 and 13.6 mm (0.21% and 1.98% of the displacement value, respectively) which is smaller than the symbols in Fig. 3. The InSAR time series agrees with the GPS (projected on each LOS vector for comparison) confirming the acceleration in late 2016 (Fig. 3).

Methods
In a first step, we fit the mean GPS velocities observed between 2019 and 2020 ( Fig. 2). In our model, a constant change in pressure is applied at the boundary of a 3dimensional ellipsoidal cavity embedded in a heterogeneous viscoelastic crust with temperature-dependent viscosity (Del Negro et al. 2009;Gregg et al. 2012). The resulting ground displacement is calculated numerically using the finite element method in COMSOL Multiphysics (v.5.4). This type of numerical solution has been described in detail and compared with analytic solutions by Le Mével et al. (2016). We run 60,000 Monte Carlo simulations, randomly sampling the 9 model parameters: location (easting, northing, depth), dip, strike, principal semi-axes (a, b, c) of the ellipsoid and the applied overpressure P (Fig. 4). We calculate and minimize the misfit to the GPS data as a χ 2 statistic: χ 2 = n i=1 ((u obs − u mod )/σ obs ) 2 i , where n = 15 is the number of GPS measurements of displacement (3 components at 5 stations), u obs is the annual displacement between 2019 and 2020 measured by GPS with uncertainty σ obs , and u mod is the modeled annual displacement. A statistical F-test with 6 degrees of freedom (Bevington et al. 1993) yields the critical value of χ 2 that delimits the 68% confidence intervals (Fig. S2). The geometry of the best-fitting model and the modeled annual displacement at each GPS site are shown in Fig. 4.
Second, we describe the temporal evolution of the surface displacement using a dynamic model of magma injection into the reservoir. As described in Le Mével et al. (2016) and sketched in Fig. 4c, we assume viscous magma flowing upward through a vertical conduit into a fluid-filled magma reservoir. We assume the magma to be an incompressible Newtonian fluid. A time-dependent overpressure P i (t) is applied at the conduit's inlet located at the base of the crust. The governing equations account for laminar flow in the fluid domain, stress in a linear viscoelastic (Maxwell) solid, and the two-way coupling at the interface between them. The numerical solution calculated using COMSOL is fully 3-dimensional and time-dependent. This model includes seven parameters: the starting time t 0 of the first magma injection, the time t 1 of the first inflection point, the maximum pressure P 1 reached at t 1 , the starting time t 2 of the second increase in inlet pressure, the maximum inlet pressure P 3 reached at the end of the simulation, the conduit radius a c , and the magma viscosity η m (Fig. 5).
To optimize the 7 model parameters, we consider the LOS displacement time series from ALOS T112 (Fig. S3) and the GPS time series of vertical displacement at MAU2, decimated to 42 data points (one measurement every 60 days) (Fig. 2). After running 6000 Monte Carlo trial simulations, we retain the set of values that minimizes the misfit to the data using the χ 2 statistic.  (Figs. 4 and S2). The best-fit ellipsoid is centered offshore between the western peninsula and the rhyolite flow of Las Nieblas (star on Fig. 4). The source is oriented NE-SW (azimuth φ = N25 • E, φ ∈ [10, 42] • E), slightly dipping to the NW (plunge θ = −3 • , θ ∈ [−21, 10] • ). This orientation resembles that of the Troncoso normal fault, one of the major structural features in the basin, and of the faults mapped under the lake Peterson et al. 2020). The location and spatial extent of the deformation pattern do not appear to have changed through time, as evidenced by the similarity of the new model geometry to the best-fitting dislocation model estimated previously by Feigl et al. (2014) and Le Mével et al. (2015). This ellipsoid explains the GPS velocities well (Fig. 4). The χ 2 misfit decreases from ∼10 5 to 10 3 (Fig. S2). The residual velocities are less than 14 mm/year and 33 mm/year for the horizontal and vertical components, respectively, at all 5 GPS stations (i.e. 11% of the displacement rate). While most of the volcanic inflation signal is explained by the model of a pressurized ellipsoidal source, the residuals to the InSAR velocity maps (Fig. S4) reveal additional localized deformation, including subsidence on the rhyolite flow SW of the lake basin and possibly motion along the Troncoso fault. Modeling these additional deformation sources is beyond the scope of this paper. The geometry of our bestfitting model for the deformation source is in accordance  (Table S1) spanning 2019 to 2020, and GPS site MAU2 (black triangle). (Right column) LOS displacement time series for each SAR track calculated for pixel at MAU2 (blue circles) and GPS observations projected onto each LOS (beige dots)

Dynamic modeling results
We assume that the two inflection points identified on the uplift time series indicate two periods of magma injection and model them by an increase in pressure at the inlet of the conduit (Fig. 5, time intervals I and III). During time interval II, the inlet pressure is constant and the magma continues to flow to the reservoir but with decreasing flow rates, leading to decreasing surface uplift rates. The best-fitting model shows: an increase in conduit inlet pressure from P 0 = 0 to P 1 = 48.  [2010.18, 2010.97]), and a second increase from P 2 = P 1 to P 3 = 79.9 MPa (P 3 ∈ [77.3, 84.7] MPa) between t 2 = 2016.73 (t 2 ∈ [2016.54, 2017.03]) and t 3 = 2020.11, with magma viscosity η m = 7.1 × 10 9 Pa s (η m ∈ [5.8 × 10 9 , 7.5 × 10 9 ] Pa s), and conduit radius a c = 82 m (a c ∈ [78, 84] m). The two episodes of magma injection are thus described with very similar rates of pressure increase, 9.0 and 9.2 MPa/year, respectively. Integrating the magma flow rate rising through the conduit to the reservoir Q(t) over time gives the cumulative magma injection volume V (t) as a function of time (Fig. 5). In total, 369 × 10 6 m 3 of magma was added to the ellipsoidal reservoir. The fit to the time series is good with residual displacements less than 52 mm (3%) and 46 mm (5%) for the GPS and InSAR-derived displacements, respectively, at each of the 51 data points included in the inversion (Fig. 5). The χ 2 misfit decreases from 0.024 to 0.009 (Fig. S5).

Discussion
The seismic monitoring network at LdM consists of 6 broadband seismometers installed between 2011 and 2013 (Cardona et al. 2018). The seismicity recorded between 2011 and 2020 occurs as repetitive swarms of volcanotectonic (VT) earthquakes of small magnitude (M L 3.3), located mostly to the southwest of the volcanic field near the Troncoso fault at depths between 2 and 5 km (Cardona et al. 2018) (Fig. 4b). The southwest cluster of VT earthquakes coincides with the southwestern edge of the ellipsoidal reservoir modeled in this study. In the same area, recent measurements of diffuse CO 2 degassing through soil have revealed an anomaly with a flux as high as 478.4 g/m 2 /day in March 2020 (Fig. 4b, OVDAS volcanic activity report of March 20, 2020). In addition, gravimetric changes measured between 2013 and 2016 align with the western edge of the modeled ellipsoid but indicates shallower mass injection at 1.5 to 2 km depth (Miller et al. 2017). The changes in deformation rates (inflection points) do not coincide with the occurrence of a particular earthquake swarm (Fig. 5). Nonetheless, the temporal evolution of the cumulative number of earthquakes seems to follow the longterm trend of the deformation rates. The cumulative number of earthquakes has decreased between 2014 and 2017 and conversely, have recently increased during interval III. We also note an increase in earthquake magnitude with more M L 2-3 earthquakes over the last two years. The most energetic seismic episode ever recorded at LdM happened in June 2020, with ∼2500 seismic events recorded during three earthquake swarms on June 11, 13, and 16. All of these observations are consistent with a pressurizing magma reservoir at depth affecting the stress distribution in the surrounding crust and hence facilitating fluid motion and/or the reactivation of local faults, producing the observed  InSAR data  from ALOS T112 and T113  (black circles), and GPS (gray dots) observations at MAU2. During time intervals I and III the measured uplift rates increase while they decrease during time interval II. For plotting purposes, each scalar LOS measurement from InSAR has been decomposed into the three components of displacement using the method outlined in Text S1; this calculation infers the vertical component of displacement from the observed LOS displacement by using the orientation from the modeled displacement vector. Seismic events recorded at LdM (black bars, maximum daily number annotated on right axis) and cumulative number (blue line). (Middle) Inlet pressure applied at base of conduit (brown) and resulting fluid pressure in the reservoir (pink). (Bottom) Temporal changes in volume flow rate feeding the reservoir (purple) and cumulative volume of magma added to the reservoir (orange) spatial and temporal patterns in seismicity and deformation Zhan et al. 2019).
The results from modeling the deformation field measured by GPS and InSAR indicate two distinct episodes of increased pressure feeding an ellipsoidal reservoir under the center of the lake basin. The temporal changes in uplift rates indicate that magma supply to the LdM magmatic system varies over short time scales of the order of several years. The total volume of magma added to the system between 2005 and 2020 is 0.37 km 3 , that is 0.45% of the total volume of the modeled reservoir (81.7 km 3 ). The modeled rates of inlet pressure change are similar in each of the two intervals of accelerated deformation (2005-2010 and 2016-2020), indicating a similar underlying mechanism. This identified pattern might repeat regularly. The magma flow rate varies between 0 and 0.0377 km 3 /year, yielding an average of 0.0243 km 3 /year. The high magma recharge rates associated with volcanic unrest observed on the decadal timescale supports the hypothesis of episodic growth of magma bodies. Some ∼35 episodes of non-eruptive intrusions of magma like the current one would be required to make up the 13 km 3 estimated by Singer et al. (2018) to explain the 62-m uplift of the paleoshoreline at LdM during the Holocene.
The high recharge rate at LdM is comparable to the volume change rate modeled during the ongoing unrest at Domuyo volcano (Lundgren et al. 2020), a silicic volcanic center ∼50 km to the SE of LdM. In their statistical study of caldera unrest, Sandri et al. (2017) show that the majority of pre-eruptive unrest lasts less than 18 months and display high rates of seismicity and degassing. The long duration of unrest at LdM (∼15 years) and relatively low seismic activity is typical of non-eruptive caldera unrest episodes, as also observed at Long Valley, Campi Flegrei or Izu-Oshima (Sandri et al. 2017). Nevertheless, the stress models of Zhan et al. (2019) show that the trajectory of the overpressure building up in the reservoir controls the occurrence and timing of potential failure at LdM. Our results showing renewed acceleration at LdM therefore have important implications for the future evolution of unrest and the hazard assessment at this volcano.

Conclusions
In this paper, we have analyzed and modeled new GPS and InSAR data revealing a renewed acceleration in uplift at LdM starting in late 2016 and reaching more than 290 mm/year between February 2019 and February 2020. The 2019-2020 GPS velocities are modeled using a 7.1 by 1.9 by 1.3 km pressurized ellipsoid at a depth of 4.8 km in a viscoelastic crust. We then use a dynamic model of magma injection to describe the uplift time series from 2005 to 2020 and find that two episodes of magma injection have accumulated a total volume of 0.37 km 3 . Both intervals of accelerating uplift are explained with a similar rate of inlet pressure increase (∼9 MPa/year). If the new episode follows the previous pattern observed between 2005 and 2010, we could expect the uplift rates to keep increasing for two additional years. The rapid deformation, together with increasing seismicity and CO 2 degassing, demonstrates the need for continued monitoring at the Laguna del Maule volcanic field.
Acknowledgements The authors thank the editors Drs. Joël Ruch and Andrew Harris, and two anonymous reviewers for providing helpful comments that helped to improve the manuscript. The authors also thank Brad Singer and Nathan Andersen for helpful discussions. We thank the following organizations for facilitating access to data: DLR (research project RES1236), JAXA, ESA, UNAVCO, WinSAR, and ASF. This research was partially supported by grants from the U.S. National Science Foundation (EAR1411779 -HLM & KLF) and Jet Propulsion Laboratory (SRTD 1641138 to KLF). KLF gratefully acknowledges sabbatical support from a Merle Tuve Fellowship at Carnegie, the Center for Academic Partnership at JPL, and the College of Letters & Science at the University of Wisconsin-Madison. GNSS and seismic data can be requested to OVDAS (SERNAGEOMIN). The COMSOL MPH files are available upon request. InSAR time series products are available at https://zenodo.org/record/4662908.

Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.