First Onset of Unrest Captured at Socompa: A Recent Geodetic Survey at Central Andean Volcanoes in Northern Chile

We report the first detection of unrest at Socompa volcano during our recent survey of Central Andean volcanos in Northern Chile using Interferometric Synthetic Aperture Radar measurements spanning January 2018 to October 2021. We find that Socompa volcano, whilst initially undeforming and no recorded eruptions for 7.2 kyr, shows a steady uplift (17.5 mm/yr) from November 2019, independently recorded by near‐field continuous Global Positioning System data. The deformation pattern can be fitted with pressure increase in an ellipsoidal source region stretching from 2.1 to 10.5 km, with a volume change rate of ∼6.2 × 106 m3/yr. Our observations of the onset of uplift suggest it is unlikely that a nearby Mw 6.8 deep intraslab earthquake on 3 June 2020 triggered the unrest. The deformation signal we detect indicates the initiation of unrest at Socompa, after at least two decades without measurable deformation, and many thousands of years without volcanic activity.

Magma Body (APMB, Ward et al., 2014). The deformation signal near Cerro Overo, which has no ground-based monitoring, is intriguing as it transitioned from subsidence to uplift in ∼2003-2005, which involves fluid accumulation and loss within the crust at ∼10 km depth . The deformation pattern of uplift at Lazufre (Lastarria and Azufre) has been interpreted to represent magma accumulation in the mid-upper crust with source depth <10 km (Díaz et al., 2015;Pearse & Lundgren, 2013;Remy et al., 2014;Ruch & Walter, 2010), in addition to a source ∼1 km under the Lastarria summit observed between 2003and 2005(Froger et al., 2007. Other deformation, for example, during the 2010 unrest at Lascar (which erupted in 2015-2017 without deformation, Gaete et al., 2020, andmost recently in December 2022) can be linked to crater evolution processes such as gravitational slumping or piston-like subsidence (Richter et al., 2018). Putana and Sillajhuay showed short-lived, low-magnitude uplift in 2009 and 2007-2010, respectively, related to hydrothermal or seismic activity Pritchard et al., 2014;Stebel et al., 2014).
The triggers for episodes of magmatic uplift (or subsidence) in the Central Andes are obscure but could potentially include (a) variations in flux from lower crustal bodies of melt (e.g., the APMB) or (b) changes within shallow reservoirs such as crystal mushor degassing as inferred in other settings . These processes cause pressure changes within reservoirs, thus controlling the initiation and cessation of inflation and deflation (or even possible episodic inflation and deflation as shown at Cerro Overo and Uturuncu, Walter & Motagh, 2014). The initiation of deformation may be linked to external events like earthquakes. For example, large subduction earthquakes in the Southern Andes and Japan caused stress field changes that triggered episodes of subsidence at multiple volcanoes Takada & Fukushima, 2013), and regional earthquakes are also thought to have triggered delayed uplift through surface waves (e.g., Lupi et al., 2017).
Here, we analyze ∼4 years of Sentinel-1 InSAR time series data, spanning January 2018 to October 2021, in the region of Antofagasta, Chile (Figure 1). Similarly to those previously reported, we observe uplift at Uturuncu, and Cerro Overo, and Azufre. However, we also find a previously unreported deformation signal centered on the Socompa volcano, where no deformation has previously been observed from regional InSAR studies (1992. We measure a steady linear uplift (rate of 17.5 ± 3.7 mm/yr) starting from November 2019, which continues through the rest of our InSAR observation time (until October 2021).
Socompa is a large stratovolcano (peak elevation 6,031 m) and is the site of a trainline and manned border control between Chile and Argentina. It is known for the failure of the northwestern flank that produced a 600 km 2 debris-avalanche deposit and triggered post-collapse eruptions ∼7, 200 years ago (Wadge et al., 1995). As a result of its remote location and presumed quiescence, it lacks targeted ground monitoring, although it was selected as the site of a single Global Positioning System (GPS) station (ID: SOCM) installed in 2011 as part of the NSF PLUTONS network (Pritchard et al., 2018) due to its location halfway between Lazufre and Uturuncu and previously used as a far-field non-deforming reference station. A small lake at the foot, and several warmspots near the summit of the volcano form a complex microbial ecosystem (Costello et al., 2009;Farías et al., 2013;Halloy, 1991) where both water and CO 2 degassing have been observed during field studies (but not from a satellite IR survey; Jay et al., 2013), implying the presence of active hydrothermal and therefore magmatic systems.
We determine the precise onset time of uplift at Socompa using a time-dependent parameterized model fitting the GPS time series from SOCM, and investigate the temporal relationship between the onset of Socompa uplift and nearby earthquakes to explore the potential trigger mechanisms. We reconstruct the cumulative deformation fields using an InSAR time series approach (Liu et al., 2021), and combine both InSAR and GPS data as inputs to assess several potential geodetic source models to explain the Socompa deformation. Finally, we discuss the sudden onset of uplift at Socompa in the context of the long timescales of unrest as observed with InSAR in the Central Andean volcanoes since the 1990s.

InSAR Time Series Analysis
Here we process Sentinel-1 InSAR time series using the LiCSAR processing chain (Lazeckỳ et al., 2020) and StaMPS software (Hooper et al., 2007, processing details in Text S1 in Supporting Information S1). Compared to single interferograms, InSAR time series analysis provides frequent estimates of surface displacement through time (every 6 or 12 days for Sentinel-1), reducing measurement uncertainties from noise (Osmanoğlu et al., 2016). In addition, if we assume a deformation model appropriate for displacements due to Socompa unrest, we can reconstruct the post-onset cumulative deformation field via a time-dependent parameterized model fitted to the InSAR time series. As the observed velocity change at Socompa is approximately linear, we assume that surface displacement at time t following the onset time t 0 can be decomposed as follows: where Η(*) is a Heaviside step function, V 1 is the background long-term linear deformation rate, V 2 is the linear velocity change after the onset time, and b is a constant reference offset in observations. We do not fit the seasonal signals because the already-applied Generic Atmospheric Correction Online Service (GACOS, Yu et al., 2018) should suppress the seasonality, and it is also difficult to model it accurately considering the noise level within the InSAR data in this region (particularly from ionospheric delay, Liang et al., 2019). After fitting the data, we reconstruct the cumulative pre-and post-onset deformation field in the line of sight (LOS) direction via the difference between the points at both ends of the fitting lines ( Figure 2).
As the Central Andes predominately lacks vegetation, coherence is very high, significantly lowering the impact of unwrapping errors and fading signals (Agram & Simons, 2015). The main InSAR error sources arise from atmospheric noise, including both tropospheric and ionospheric components. Although the applied GACOS corrections (Yu et al., 2018) improve the data quality (with average standard error reductions of 16.9% and 45.7% for ascending and descending interferograms, respectively, Figure S1 in Supporting Information S1), the  (Perkins et al., 2016). The M w 6.8 earthquake epicenter (3 June 2020 with a 112 km centroid depth), is marked by the focal mechanism, while all the M w > 6.0 historical earthquakes since 1976 (where a precise global seismic network was established) in this region are shown by blue circles colored by centroid depth (records from the USGS, United States Geological Survey).
ionospheric noise is very strong and could not be ignored, especially on ascending track. We therefore remove a linear ramp that spans the whole interferogram to reduce ionospheric noise, and other long wavelength signals associated with orbit errors and plate motion. Overall, the noise level in the ascending data is much higher than for descending, and noticeable atmospheric artifacts remain in high topography areas.

Onset Time Implications
Determining the onset time for Socompa uplift is important not only for reconstructing the cumulative deformation field, but also for investigating potential causes for initiating unrest at Socompa. Due to historically (a) The pre-and post-onset cumulative deformation fields using ascending track data. The focal mechanism, black dashed polygon, and green square represent the epicenter of the M w 6.8 earthquake, the approximate boundaries of the salar (salt pan) regions (Text S2 in Supporting Information S1), and the InSAR reference points, respectively. The InSAR time series plots of some peak displacement pixels near the volcanoes are shown in (c). (b) Same as (a) but for descending track data. In all figures, positive values mean movement toward the satellite. Note the displacement signal associated with Cerro Overo is ∼40 km southeast of the volcano and falls outside of the data coverage on the descending frame.
lower temporal resolution satellite imagery and the typically long duration of unrest, it has not previously been possible to determine precisely the initiation time of deformation at a Central Andean volcano. For example, while uplift at Sabancaya is known to have started in 2013, the distribution of Synthetic Aperture Radar (SAR) acquisitions means that it could have taken place at any point over several months (Macqueen et al., 2020). Here we investigate potential triggering mechanisms from earthquakes in this region by exploring all major historical events (M w > 6.0 since 1976, see Figure 1). We find the closest event in space and time is a M w 6.8 intraslab earthquake with a 112 km centroid depth, which occurred on 3 June 2020, and whose epicenter is ∼120 km northwest of Socompa ( Figure S2 in Supporting Information S1). Initial visual inspection of the deformation signal pointed to the potential for a causal relationship given the close correlation in time, but the InSAR time series are noisy and contain seasonal signals that overprint changes in long-term trends.

Onset Time Determination
To determine the exact onset time, we collect data from a previously installed GPS station (SOCM, , which is located ∼8 km southwest of the Socompa volcano and captures the deformation signal ( Figure 3c). We use a time series with average daily positions processed by the Nevada Geodetic Laboratory in a South American Plate reference frame (Blewitt et al., 2018) to do a time-dependent parameterized fitting, using the trajectory model: ( ) = 1 + ( − 0) 2 + 1 sin(2 + 1) where the unit of time t is year, A 1 , A 2 and φ 1 , φ 2 are the amplitudes and phases of annual and semi-annual terms respectively, t eq(i) and C i are historical earthquake event times and corresponding coseismic offsets that are close to the station (based on the database of the Nevada Geodetic Laboratory). We use a Markov Chain Monte Carlo approach to determine that the optimal onset time t 0 is 19 November 2019 (197 ± 12 days ahead of the earthquake event time, Text S3 in Supporting Information S1), using all three components of GPS time series data and weighting them by the noise level within each component (Figure 3a).
Such analysis highlights strong seasonal effects in the GPS time series, especially in the North direction, leaving some uncertainty about the onset time in the data fitting. To reduce the influence of seasonal signals, we further use a novel method of vector decomposition, transforming the East-North vectors into another orthogonal coordinate system, aligned along the movements parallel and perpendicular to the direction of seasonal motion (Text S4 in Supporting Information S1). The decomposition results (Figure 3b) clearly show an onset time half a year ahead of the earthquake, suggesting that unrest at Socompa was unlikely dynamically triggered by seismic waves induced by this earthquake.

Volcanic Geodetic Source Modeling
As we detect unrest at Socompa for the first time, we explore several source models to explain the observed deformation, including a point pressure source (Mogi, 1958), prolate spheroid (Yang et al., 1988), dipping dike with uniform opening (Okada, 1985), and a point Compound Dislocation Model (pCDM, Lundgren et al., 2017;Nikkhoo et al., 2017). We use the reconstructed post-onset cumulative deformation fields from the InSAR time series (November 2019 to October 2021), and cumulative GPS deformation at SOCM station that has the same time scale as InSAR data, as it improves the signal-to-noise ratio (SNR) of input data and subsequently provides more robust modeling results.
We first use a nested uniform downsampling of the reconstructed post-onset cumulative deformation fields, with a greater pixel density in the deformation area ( Figure S3 in Supporting Information S1). Then we use the Geodetic Bayesian Inversion Software (GBIS, Bagnardi & Hooper, 2018), a Bayesian approach for the inversion of multiple geodetic data sets that provides the posterior probability density functions of source model parameters, to invert the model parameters. We embed the code of pCDM  into the GBIS software so that all models run in the same environment, and use the data uncertainty within the InSAR and GPS observations to weight them during the inversion ( Figure S4 in Supporting Information S1).
To obtain the equivalent volume change of pCDM, we further use the point Ellipsoidal Cavity Model (pECM, Nikkhoo et al., 2017), a special case of pCDM i.e., constrained to represent a pressurized ellipsoidal cavity,  Table S1 in Supporting Information S1, the values of reduced chi-squared, Text S5 in Supporting Information S1 are 0.42, 0.43, 0.48, and 0.65 for pCDM, Yang, Okada, and Mogi, respectively), and these two models give similar source depth (6.3 and 7.3 km for pCDM and Yang, respectively), geometry, and volume change (∼1.1 × 10 7 m 3 for both models).

Other Volcanic Deformations
Several volcanoes have been reported to be deforming in the past few decades ( Figure 2) and we tie our InSAR observation to the GPS network in this region (Text S6 and Figure S10 in Supporting Information S1) to compare our results to these earlier studies. Starting in the north, Uturuncu previously showed a deformation rate of ∼15 mm/yr in the 1990s (Fialko & Pearse, 2012;, but this gradually slowed in the 2010s (Gottsmann et al., 2017) to a rate of 3-5 mm/yr in ∼2017 (Lau et al., 2018) and ∼3 mm/yr in ∼2020 (Eiden et al., 2023). In agreement with previous studies, we observe an uplift rate of 2.5 ± 1.8 mm/yr on Uturuncu (2018-2021, Figure S11a in Supporting Information S1). Putana displayed short-lived uplift totaling 40 mm displacement in 2009-2010 , whilst we find potential subsidence of −3.9 ± 2.1 mm/yr (November 2019 to October 2021), with an onset which seems coincident with the deformation at Socompa (Figure S11b in Supporting Information S1). The deformation signal close to the Cerro Overo, which previously changed from subsidence of −4 mm/yr (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) to uplift of 5 mm/yr through 2010 on descending track , continues to uplift at a rate of 3.8 ± 2.6 mm/yr (ascending LOS velocity, 2018-2021). Lazufre volcano shows uplift rates of 11.2 ± 1.7 mm/yr (2018-2021, Figure S11c in Supporting Information S1), consistent with the trend of surface deformation slowing at this volcano Remy et al., 2014).

Discussion
Since uplift at Socompa started months before the M w 6.8 intraslab earthquake, we instead consider a plausible explanation for the sudden uplift to be the ascent of magmatic fluids from a deeper melt source into a shallower reservoir. A magnetotelluric study (Ślęzak et al., 2021) in the Atacama region found a high conductivity zone at Socompa (∼5 km west of the volcano and spanning 2 km to over 30 km depth), although there is significant uncertainty on the existence of the high conductivity zone as it is constrained by only one measurement point at Socompa, whereas the high resistivity region suggests that Socompa is unlikely connected to APMB (Figure 4c). The crust beneath Socompa and Cerro Overo has not been subject to the same level of study as Uturuncu and the APMB (e.g., Comeau et al., 2016). Deformation at Socompa has some first order similarities to that at Lazufre (∼90 km to the South): they have similar source depth (<10 km) and rate of volume change in order of 10 6 m 3 / yr Remy et al., 2014). The shallow reservoir and hydrothermal system beneath Lazufre have been suggested to be linked to the Southern Puna Magma Body (Budach et al., 2013;Stechern et al., 2017;Ward et al., 2017), but there is no independent evidence for this at Socompa.
An interesting question here is whether the initiation of uplift at Socompa will maintain a linear rate, decrease exponentially like at Lazufre or whether it will gradually slow and eventually cease. The current geodetic observations show no trace of deceleration, which in a purely elastic system would imply a constant pressure increase. Alternatively, it may be too early to detect any decrease in magnitude of a hydraulic connection to a deeper magma supply. Note that the current pressure we obtain from the model is ∼13.2 MPa, which is far less than the overpressure required for chamber wall failure (Gerbault et al., 2018).
Deformation at Socompa is very consistent with other observations of unrest in the Central Andes (e.g., Pritchard & Simons, 2004a). The deformation rate is low (usually <30 mm/yr), and uplift started after an apparently very long period of quiescence (like Lazufre), consistent with deformation taking place on much longer timescales that other parts of the Andes (inter-eruptive and co-eruptive deformation rates are much higher in both the Northern and Southern Andes, e.g., Pritchard & Simons, 2004b;Fournier et al., 2010; Morales Rivera et al., 2016). This means that Holocene activity is not necessarily a good basis for assessing whether Central Andean volcanoes have active magmatic systems or are likely to enter a phase of unrest. InSAR measurements of deformation are therefore critical for the detection of volcanic unrest. However, volcano deformation in the Central Andes is generally not associated with eruption (except at Sabancaya and Lascar), reflecting lower rates of reservoir pressurization and therefore lower rates of magma flux that are more conducive to intrusion growth than brittle failure, dyke propagation and magma ascent (Biggs & Pritchard, 2017).

Conclusion
Our observations update the volcanic monitoring at Central Andean volcanoes in Northern Chile. We first detected unrest at Socompa volcano, contribute to a picture of low-rate, episodic deformation in this region, indicative of magmatic processes that take place on very long time scales of decades. We determine the onset time of Socompa uplift in November 2019, with a linear rate of 17.5 mm/yr up to January 2023, using InSAR and GPS observations. We test several geodetic source models, finding the best-fit for an ellipsoidal source located at a depth of ∼6.3 km and volume change of ∼6.2 × 10 6 m 3 /yr. We capture the onset of deformation at a Central Andean volcano for the first time at high temporal resolution, which suggests earthquake triggering is unlikely in this case. This provides a potentially important dataset for assessing the temporal development and therefore origin of such deformation.

Data Availability Statement
The Sentinel-1 SAR data are copyrighted by the European Space Agency and additionally distributed by the Alaska Satellite Facility Distributed Active Archive Center (https://earthdata.nasa.gov/eosdis/daacs/asf). All the GPS data we use, including SOCM (https://www.unavco.org/data/doi/10.7283/T5TT4P27), is processed by the Nevada Geodetic Laboratory (http://geodesy.unr.edu), which are based on services provided by the GAGE Facility, operated by UNAVCO, Inc., with support from the National Science Foundation, the National Aeronautics and Space Administration, and the U.S. Geological Survey under NSF Cooperative Agreement EAR-1724794. The InSAR time series fitting and geodetic modelling results are available on Zenodo (https://doi.org/10.5281/ zenodo.7688945). This work is supported by the UK Natural Environment Research Council (NERC) through the Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET, http://comet. nerc.ac.uk), and the Looking into the Continents from Space (LiCS) large Grant (NE/K010867/1). Figures were made using the Generic Mapping Tools (GMT) (Wessel et al., 2013). FL acknowledges support from the Great Britain-China Educational Trust. JE and TC acknowledge support from the Royal Society through University Research Fellowships (UF150282 and URF\R1\180088). SE is funded by a NERC Independent Research Fellowship (NE/R015546/1). AH and CNL acknowledge support from the European Research Council (ERC) through the EU Horizon 2020 project DEEPVOLC (Grant 866085). FD thanks FONDECYT Iniciacion 11220513 research grant for funding. We acknowledge Kristina Butler, Catherine Gagnon, Joaquín Castillo, Sofía Parra, Milton Quinteros, and Gabriela Herrera for collecting the GPS data of the SOCM station. We thank Scott Henderson, Julie Elliott, and especially Matthew Pritchard (for installing the SOCM station) for their help in retrieving and processing the data.