Reevaluating Volcanic Deformation Using Atmospheric Corrections: Implications for the Magmatic System of Agung Volcano, Indonesia

A major challenge in using satellite interferometry (interferometric synthetic aperture radar) for volcanic monitoring in the tropics is distinguishing volcano deformation from atmospheric noise. We reanalyzed interferometric synthetic aperture radar time series from 2007–2009 from Agung volcano, Indonesia, which had previously been used as evidence for unrest. Using uncorrected data, we find an apparent velocity of 5.0±2.7 cm/yr consistent with previous reports, but we show that this signal is consistent with predictions of atmospheric contributions derived from weather models (European Centre for Medium‐Range Weather Forecasts‐High‐Resolution). Following the correction, we find a velocity of 1.4±4.2 cm/yr and conclude that there was no significant deformation related to the inflation of a shallow magma source from 2007–2009. We discuss the implications for the inferred magma storage system at Agung and consider which other reported signals might have been wrongly attributed to deformation and should be reanalyzed.


Introduction
Deformation has been reported at over 200 volcanoes globally (Biggs & Pritchard, 2017) and the majority of reported observations are based on satellite interferometric synthetic aperture radar (InSAR). However, atmospheric artifacts can lead to the misinterpretation of volcanic deformation (i.e., Beauducel et al., 2000;Poland & Lu, 2008;Remy et al., 2015), particularly at tropical volcanoes where the water vapor content is high and rapidly varying (Ebmeier et al., 2013). Thus, a fundamental question remains, are all these volcanoes actually deforming? And if they are, what kind of deformation can we resolve with the available satellite data and ancillary products?
The Indonesian volcanic arc is home to 13% of the world's active volcanoes, including Agung, Bali, which produced one of the largest eruptions of the twentieth century in 1963 (Self & Rampino, 2012). Uplift was reported at Agung between July 2007 and February 2009 based on InSAR time series from the Advanced Land Observation Satellite (ALOS) (Chaussard & Amelung, 2012). The rate of uplift was reported as 5.7 cm/yr in Chaussard and Amelung (2012) and 7.8 cm/yr in Chaussard et al. (2013) based on different numbers of InSAR images. The deformation was attributed to magma injection into a closed-vent volcano at a depth of 1.4 km below sea level (Figure 1b).
In September 2017, the Indonesian Center for Volcanology and Geological Hazard Mitigation reported a seismic swarm and fumarole activity (Center for Volcanology and Geological Hazard Mitigation Press Release, https://magma.vsi.esdm.go.id), which was followed by eruption on 21 November 2017 . Between August and November 2017, Sentinel interferograms showed uplift between Agung and its RESEARCH LETTER 10   neighboring volcano, Batur. However, atmospheric corrections were required to properly characterize the deformation due to persistent atmospheric artifacts . The 2017 deformation was interpreted as a dike intrusion at a depth of 7-13 km b.s.l.   (Figure 1c), different to the 2007-2009 shallow magma source reported by Chaussard and Amelung (2012). This paper aims to reevaluate InSAR data at Agung for the period 2007-2009 using newly available atmospheric correction techniques and in the light of new InSAR results from Albino et al. (2019). We hypothesize that the apparent deformation signals reported by Lingyun et al. (2013), Chaussard and Amelung (2012) and Chaussard et al. (2013) are caused by topographically correlated atmospheric artifacts, suggesting that Agung was not deforming in [2007][2008][2009]. Consequently, the conceptual model of Agung's magmatic system and recharge between the 1963 and 2017 eruptions requires reanalysis.

Methods
To reevaluate the reported deformation of Agung, we processed the same interferogram pairs used by Chaussard and Amelung (2012) and Chaussard et al. (2013). The synthetic aperture radar images are acquired by the Japan Aerospace Exploration Agency ALOS Phased-Array L-band Synthetic Aperture Radar instrument. We used the Repeat Orbit Interferometry Package developed by Jet Propulsion Laboratory and the California Institute of Technology (Rosen et al., 2004) to produce 31 interferograms from ALOS ascending Track 7010, Scene 422 spanning 21 February 2007 to 26 February 2009 (Table S1 and Figure S1 in the supporting information). Three pairs used by Chaussard and Amelung (2012) had processing issues and were excluded (Table S1 and Figure S1). While Chaussard and Amelung (2012) used a reference date of 21 February 2007, we chose 9 July 2007 as the reference date for our study because it marks the beginning of the period of uplift from the uncorrected interferograms. The topographic phase was removed from each interferogram using a 30-m Shuttle Radar Topography Mission 1 Digital Elevation Model (Farr et al., 2007). All interferograms were filtered using a power spectrum filter (Goldstein & Werner, 1998) and unwrapped using the branch cut algorithm (Goldstein et al., 1988). The radar phase change was converted to millimeter range  Figure S2), which is later interpolated to plot this map. The black boxes around Agung and Batur denote LOS velocity values not included in the frequency distribution. (c) Frequency distribution of LOS velocity at Agung using data before interpolation and excluding values obtained within the black square. Before applying GACOS atmospheric correction, the highest frequency of LOS velocity at Agung is 5.0 ± 0.9 cm/yr. This is similar to the LOS velocity observed by Chaussard and Amelung (2012) at 5.7 cm/yr indicated by the red dash line. After applying GACOS correction (d), the mean LOS velocity at Agung is 1.4 ± 0.8 cm/yr. change whereby positive measurements correspond to displacement toward the satellite. Long wavelength ramps due to ionosphere or orbital error were removed using linear ramps before atmospheric correction.
We performed a least squares inversion using a Small Baseline Subset approach (i.e., Berardino et al., 2002;Schmidt & Bürgmann, 2003) to generate a time series and a line-of-sight (LOS) velocity map for Agung and Batur. We calculate time series in an 11 × 11 window and plot the mean and standard deviation (2 ) of these values. We applied a linear regression to each time series to measure the LOS velocity with errors bars stated at a 95% confidence interval (2 ).
InSAR measures relative displacement, meaning the displacement of a point of interest is dependent on the choice of reference point. It is therefore important to select a stable reference point, but this can be especially challenging where there is strong atmospheric noise. One approach is to fit a plane through the nondeforming areas to establish a reference level (i.e., Delgado et al., 2017;Henderson & Pritchard, 2013), but this is challenging for small volcanic islands, where the nondeforming area is small and poorly defined. To quantify the uncertainty associated with the choice of reference pixel at Agung, we repeat the LOS velocity calculation using different reference points on a grid of 93 m × 90 m, excluding the area around Agung and Batur. We plot the result in a map view whereby each point represents the LOS velocity of Agung measured using that point as the chosen reference pixel ( Figure S2). To help visualization, the coarse grid is interpolated using the gap filling method of Garcia (2010) and Wang et al. (2012) and the probability distribution of velocity values plotted (Figures 2a and 2b).
To estimate the atmospheric contribution to our interferograms, we use the Generic Atmospheric Correction Online Service for InSAR (GACOS). GACOS is based on the European Centre for Medium-Range Weather Forecasts High-Resolution weather model (0.125 • and 6-hr resolutions) and uses an Iterative Tropospheric Decomposition method to separate stratified tropospheric signals and turbulent signals and interpolates  (Table S3). (f) Comparison between ground velocity and topography that shows a strong correlation between InSAR and GACOS. them using the 90-m Shuttle Radar Topography Mission Digital Elevation Model (Yu et al., 2017;Yu, Li, Penna, & Crippa, 2018). We then convert the phase delays of GACOS model from zenith to the LOS of the interferograms.

Results
We apply a linear regression to the time series of uncorrected interferograms and find a LOS velocity of 5.0 ± 2.7 cm/yr (Figure 3a), which is slightly lower than that found by Chaussard and Amelung (2012) and Chaussard et al. (2013). By varying the choice of the reference pixel, we find an uncertainty of the central value of ±0.9 cm/yr (Figure 2b), consistent with the spread of previously reported values. This highlights the importance of the choice of reference pixel when reporting InSAR measurements in noisy atmospheric conditions. For this study, we select a reference pixel coordinate that corresponds to the peak of the histogram (Figure 2c), which is [115.375 • E, 8.451 • S] (Table S2).
The time series produced using the GACOS model show a similar pattern to the uncorrected time series (Figures 3d and 3e) with rates of 3.6 ± 2.0 and 5.0 ± 2.7 cm/yr, respectively. After applying the GACOS correction, the residual velocity is 1.4±4.2 cm/yr (Figure 3f). The correlation between topography and signal reduces from R 2 = 0.39 before correction to R 2 = 0.02 after correction, suggesting that the signal in the uncorrected interferograms are largely attributed to topographically correlated atmospheric artifacts rather than real deformation.
Both Agung and Batur have residual velocities on the order of 1.4 cm/yr. We investigate this further by extracting profiles across the summits of Agung and Batur (Figure 4). At both volcanoes, the uncorrected signal shows a clear relationship with topography, with a larger magnitude at the higher-relief Agung. The GACOS model shows a similar pattern, but with a lower magnitude and greater spatial extent, likely as a result of the interpolation of the coarser weather model. Thus, the correction underestimates the peak of the signal and overestimates in the surrounding area producing the characteristic donut shaped spatial pattern seen in the corrected interferograms (Figure 3c).

Importance of Atmospheric Correction in Volcano InSAR
We show that the deformation signals previously reported at Agung correlate with atmospheric artifacts estimated using weather models. This observation is consistent with previous studies, which shows that interferograms at low latitudes (Ebmeier et al., 2013) and close to coastline are prone to the effects of atmospheric noise (i.e., Albino et al., 2019;Bastin et al., 2007;Huang & Wen, 2014). Despite the similarity between the observed signal and the GACOS model, the correction is not perfect, leaving a residual signal of 1.4 cm/yr and increasing the standard deviation from from 2.7 to 4.2 cm/yr. The coarse spatial resolution of European Centre for Medium-Range Weather Forecasts High-Resolution weather model (0.125 • resolution) limits the removal of atmospheric artifacts smaller than ≈13.9 km (Zebker et al., 1997), which may leave residuals in the interferogram. This reflects the importance of using the highest available resolution weather models when correcting InSAR images (Parker et al., 2015).
Here we have demonstrated that reports of deformation at Agung in 2007-2009 were due to a misinterpretation of atmospheric signals, but how many of the 220 other volcanoes with reported deformation are also misclassified? Atmospheric artifacts are dependent on the topography and geographic location of the volcano (Parker et al., 2015) and tropical volcanoes are likely to experience the largest magnitude of artifacts (Ebmeier et al., 2013). Chaussard and Amelung (2012) and Chaussard et al. (2013) used pairwise logic (Massonnet et al., 1995) to argue that the signal at Agung volcano was real uplift. However, while this method is may still be useful in temperate regions, we find it is less useful in tropical regions where water vapor variations are continuous, rather than associated with discrete weather events. Fortunately, many deformation episodes initially detected by satellite have been confirmed by ground-based measurements (i.e., Lloyd et al., 2019;Parks et al., 2015). However, many have not, including the seven other Indonesian volcanoes where deformation was reported by Chaussard and Amelung (2012) and Chaussard et al. (2013). Of these, Lamongan and Anak Krakatau are much lower relief than Agung, and the signals in those volcanoes are unlikely to be attributed to atmospheric artifacts. However, Slamet, Lawu, Merapi, and Sinabung have summit elevations similar or greater than Agung and reanalysis using modern atmospheric correction techniques is warranted. Developing countries in Southeast Asia, such as Indonesia and Philippines, could benefit enormously from the use of satellite InSAR for monitoring volcanoes considering the large number of eruptions that occur and the lack of ground-based networks (Fernández et al., 2017).
Understanding the resemblance of atmospheric anomalies in volcanic deformation is also important for future studies, especially those conducted during volcanic crises where misidentified signals can have more serious consequences. This is particularly challenging as data become more widely available and nonspecialists share results on social media (Loughlin et al., 2015). One possibility is the use of machine learning algorithms, which have the potential to automatically distinguish between atmospheric signals and deformation with minimal expert input (Anantrasirichai et al., 2018(Anantrasirichai et al., , 2019. The output of the machine learning algorithm is the probability that the signal is caused by deformation as opposed to atmospheric artifacts. Thus, the nonexpert user is no longer required to discriminate between atmospheric artifacts and deformation. Unfortunately, it is likely that the restricted acquisition plan of ALOS-2 will continue to limit the study of volcano deformation in this part of the world, until the launch of the open data L-band NISAR mission in 2022. We also demonstrate the challenges associated with selecting a stable reference point in tropical environments where turbulent atmospheric anomalies are widespread. The differences in the rate of uplift between this study, Chaussard and Amelung (2012) and Chaussard et al. (2013) can be explained by the difference in choice of reference point. Where possible, the best approach is to integrate InSAR data with a ground-based geodetic network to provide an absolute reference frame (i.e., Delouis et al., 2002;Muller et al., 2014). However, in many cases that is not possible, and we propose a method of choosing a reference point by searching through all possible references, looking at the distribution of velocities and choosing a representative point. This is particularly useful for small islands where the nondeforming area is small.

Magma Plumbing Systems of Volcanoes
We conclude that (1) there is no evidence of significant (>2.7 cm/yr) volcano deformation at Agung from 2007-2009 and (2) there is no evidence of magma recharge into a shallow magma reservoir (<7.8 km). This conclusion requires a reassessment of our conceptual models of the plumbing system and eruptive cycle at Agung.
The lack of evidence for deformation does not disprove the existence of a shallow reservoir. The 1963 (and possibly 2017) eruption of Agung was triggered by magma mingling between a deep basaltic input and an existing basaltic-andesitic to andesitic reservoir (Self & King, 1996), suggesting such reservoirs are at least transiently present (Sparks et al., 2019). Chaussard and Amelung (2012) linked their now discredited shallow deformation source (1.4 km b.s.l.) to a transpressional tectonic regime on the basis that a compressional stress field would prevent shallow magma propagation, (i.e., Cembrano & Lara, 2009;Mériaux & Lister, 2002;Zellmer et al., 2014). In contrast, the 2017 dike intrusion between Agung and Batur indicates an extensional stress field imposed by edifice loading , suggesting that regional tectonic stresses are less important than the local stress field in guiding magma pathways consistent with observations of large edifice and high flux volcanoes elsewhere (Wadge et al., 2016). The 2017 dike intrusion also demonstrated a previously unrecognized hydraulic connection between the volcanic systems at Agung and Batur, which is consistent with observations of a small eruption at Batur following the 1963 eruption of Agung Wheller & Varne, 1986).
Similarly, the lack of measured deformation in 2007-2009 does not discount the possibility that deformation occurred during other time periods. Assuming a linear rate of magma recharge into the same reservoir that erupted in 1963 between 1963 and 2017, to recover the 400 ×10 6 m 3 of dense rock equivalent (Self & Rampino, 2012) would require a constant recharge rate of 7.4 ×10 6 m 3 /yr, with the assumption that the eruption was triggered by mass injection or deformation occurs in transient episodic pulses of deformation. For an InSAR detection threshold of 2.7 cm/yr and the assumption of Mogi point source model (Mogi, 1958), this would only be detectable for reservoirs shallower than 7.8 km ( Figure S5). Alternatively, if we consider the rate required to accumulate the 23 ×10 6 m 3 associated with the 2017 events Syahbana et al., 2019), anything deeper than 1.9 km would not be detected.
However, linear recharge of an elastic reservoir by an incompressible magma is not the only plausible model, and not even the most probable. Given the sparse geodetic record, discrete deformation episodes, exponentially decaying recharge, or incomplete emptying of the magma reservoir by the 1963 eruption or deformation below detection threshold are all plausible. Each of these could lead to an eruption in 2017 without any measurable surface deformation during the 2007-2009 period of L-band InSAR observations. Alternatively, magma input does not necessarily cause significant deformation as viscous deformation and compressible, gas-rich magmas can dissipate the pressure increase associated with magma input or have viscoelastic response, reducing surface deformation (Gottsmann et al., 2006;Kilbride et al., 2016). A lack of detectable deformation can be caused by a low rate of magma recharge, a high detection threshold, or gas-rich compressible magmas (Ebmeier et al., 2013;Kilbride et al., 2016).

Conclusion
We use newly available atmospheric correction techniques to reevaluate the volcanic deformation detected using InSAR at Agung volcano between 2007 and 2009 (Chaussard & Amelung, 2012). After atmospheric correction, we find a LOS velocity of 1.4 ± 4.2 and 1.4 ± 3.1 cm/yr at Agung and Batur, respectively. Our analysis suggests that the previously reported signals were caused by atmospheric artifacts and we conclude that there is no detectable volcano deformation related to a shallow magma source. However, the sparse geodetic data do not preclude other magma recharge scenarios, including steady recharge of a deep magma system (>7.8 km), rapid recharge in the 1960s, or pulses of recharge outside the 2007-2009 observation window. This study shows that monitoring tropical volcanoes using InSAR remains a challenge, and due to the high elevation and its proximity to the ocean, Agung is particularly prone to atmospheric artifacts (i.e., Ebmeier et al., 2013;Parker et al., 2015). However, with the right techniques and using atmospheric corrections where appropriate, spectacular examples of ground deformation at tropical volcanoes can be observed Hamling et al., 2019).