Tomographic Imaging of the Agung-Batur Volcano Complex, Bali, Indonesia, From the Ambient Seismic Noise Field

The Agung-Batur Volcanic Complex (ABVC), part of the Sunda volcanic arc, is the source of some of the most hazardous volcanic activity in Indonesia. The ABVC has undergone many small (VEI 1-2) eruptions since historical records began in the early 19th century, but Mt. Agung has experienced much larger (VEI 5) eruptions, both in the modern (1963) and historical (1843) eras, as well as several times during the past 2000–3000 years. The 1963 eruption caused more than 1000 deaths, and a more recent eruption in 2017 caused the evacuation of 140,000 people. Delineating the magma structure beneath ABVC is an important first step in understanding the physics of these eruptions. This paper presents the first local-scale study of Rayleigh wave group velocity structure and the seismic velocity structure beneath the ABVC using ambient seismic noise tomography. Seismic data were collected using 25 seismometers deployed across the ABVC during early January to March 2019. The local seismic network provides good resolution beneath both Mt. Agung and Mt. Batur. We obtained 158 Rayleigh Green’s Functions, extracted by cross correlating noise simultaneously recorded at available station pairs. We used sub-space inversion to calculate group velocity at different periods and to estimate the lateral variations in group velocity for given periods. 2-D tomographic maps obtained from the inversion of the group velocity of Rayleigh waves clearly showed some pronounced velocity anomalies beneath the ABVC. We applied the Neighbourhood Algorithm (NA) technique to invert the Rayleigh wave dispersion curves to obtain shear wave velocity (Vs) vs. depth profiles. These profiles indicate a low Vs of about 1 km/s underlying the volcanic complex between Mt. Agung and Mt. Batur at depths up to 2 km, which we suggest is due to a combination of low-Vs volcanic deposits as well as a shallow hydrothermal fluids system associated with magma fluids and/or gases produced by magma intrusion at depths >7 km.


INTRODUCTION
Knowledge of the structure and dimensions of a volcano's plumbing system is important for understanding its eruptive behavior (Pinel and Jaupart, 2003). Local-scale seismic arrival time tomography is an effective way to obtain such knowledge, as shown by Monteiller et al. (2005) and Widiyantoro et al. (2018) for Kilauea, Hawaii and Merapi, Indonesia, respectively. Seismic ambient noise tomography (ANT) (see, e.g., Shapiro and Campillo, 2004;Shapiro et al., 2005) has also been successfully applied to obtain the 3D upper crustal structure in many volcanic regions around the world, such as Toba in Sumatra, Indonesia (Stankiewicz et al., 2010), Asama in Japan (Nagaoka et al., 2012) and Colima in Mexico (Escudero and Bandy, 2017). This study aims to use ANT to delineate upper crustal structure beneath Agung Batur Volcanic Complex (ABVC), Bali, Indonesia using a local seismic network.
Understanding the eruptive behavior of the ABVC is important because of the threat it poses to populations in Bali and beyond. The ABVC has experienced at least 29 historic eruptions since the early 19th century (Global Volcanism Program, 2013). Almost all of the eruptions have Volcanic Explosivity Index (VEI) 1-2. The VEI 5 eruption of Mt. Agung in 1963, however, is regarded as one of the most significant eruptions in the 20th century, having a global climate impact and causing over 1000 fatalities, making it one of Indonesia's most deadly eruptions (Self and Rampino, 2012). Fontijn et al. (2015) used tephrostratigraphic analysis to show that the 1843 eruption of Mt. Agung, as well as at least 4 other events 2000-3000 BP, were likely as big as the 1963 eruption. The composition of eruptive products ranges is from basaltic andesite to andesite. The latter characterized the large 1963 and 1843 eruptions.
Since December 2018, a temporary seismographic network has been deployed around the ABVC that we refer to as the Agung Seismic Experiment. The Agung Seismic Experiment aims to monitor the heightened seismicity around the ABVC after the eruption in November 2017, and to use ANT to determine seismic velocity structure that can improve our understanding of the ABVC's magmatic plumbing and hydrothermal systems.

GEOLOGICAL SETTING AND RECENT ERUPTION RECORD
Bali, an island in the Indonesian archipelago, is part of the Sunda arc ( Figure 1A), which stretches from Sumatra to Nusa Tenggara through Bali. The Agung Batur Volcanic Complex lies above the subduction zone where the Indo-Australian plate dives below the Sunda block, resulting in the production of melt that drives volcanic activity. Chaussard and Amelung (2012)  Mt. Anak Krakatau. The evidence of inflation can be seen in this study, which showed that these volcanoes have shallow magma reservoirs at depths 1 -3 km, and significant vertical deformation rates of 3-8 cm/yr were detected, indicating in particular inflation of Mt Agung during 2006-2009. The Agung Batur Volcanic Complex is a stratovolcano formed by volcanic deposits left by past eruptions that produced andesite lava, volcanic breccia, volcanic ash, and pyroclastic debris. Previous studies have inferred the geological structure of the volcanic layers and depths of magma chambers beneath the ABVC. Geiger et al. (2018) described a multi-level magma plumbing based on mineral-melt equilibrium thermobarometry of lavas produced by the 1963 eruption of Agung and eruptions of Batur in 1963 and 1974. They established that, for both Agung and Batur, magmas were sourced from both an upper crustal reservoir at 3-5 km depth and a lower crustal reservoir near the Moho at about 20 km depth (just above the Moho at 18 km depth for Batur, and right at the Moho for Agung). Using satellite InSAR estimates of surface deformation and 3D numerical modeling, Albino et al. (2019) suggested that a sub-vertical magma intrusion midway between Agung and Batur, at 7-13 km depth, might play an important role in connecting these magmatic systems. Zen and Hadikusumo (1964) showed that the 1963 Mt. Agung eruption had a major impact on residential and agricultural areas in Bali. Monitoring of volcanic activity in Mt. Agung at that time did not use the type of sophisticated instrumentation that would be used today, so no one knows in detail the changes in seismic activity, inflation and heat flow. The increased volcanic activity at that time was known only through the tremors that were felt by people living within 6 km from the peak of Mt. Agung. A gust of dust billowed and formed a column of hot clouds from March 17 to May 16 at an altitude of >20 km above the crater peak (Self and Rampino, 2012). About 7 days after the earthquake tremors commenced, a thick lava flow emerged from the crater and flowed northward. Increased activity continued until mid-March 1963 (Surjo, 1981). Self and Rampino (2012) explain the chronology of the eruption of Mt. Agung in 1963-1964 as a continuous eruption. The lava flow produced in the eruption of Mt. Agung has some unusual attributes. The total volume of basaltic andesite material erupted is estimated to be <0.5 km 3 (0.1 km 3 lava; <0.4 km 3 pyroclastic fall deposits; 0.05 km 3 pyroclastic flow deposits -all dense rock equivalent volumes).
The Mt. Agung crisis in 2017 began with an increase in earthquakes swarm in mid-May 2017, this time recorded by the PVMBG seismometer network . On September 22, 2017 there was an increase in seismicity around Mt. Agung with the number of detected earthquakes reaching more than 800. On the same date, PVMBG gave a warning about the danger of eruption of Mt. Agung and raised the alert to its highest level, 4 ("Awas"). The announced danger zone was within 9-12 km from the top of the crater and about 140,000 people evacuated, including about 70,000 who self-evacuated from outside the exclusion zone . On November 21 2017, a phreatic eruption produced a hot cloud that rose up to 700 m above the crater. This condition resulted in major disruption of flight schedules in Bali. Syahbana et al. (2019) chronology of the 2017 eruption includes analysis of seismic activity and geodetic measurements that indicate a cluster of earthquakes as well as inflation centered in an area between Mt. Agung and Batur Cladera. From this evidence Syahbana et al. (2019) inferred that magma intruded into a mid-crustal dike between Mount Agung and Batur Caldera, and that this initiated the earthquake swarm in late September.

Seismic Data Acquisition
The temporary local seismograph network at ABVC was installed in mid-December 2018 and has been recording up to August, 2019. We deployed 25 seismographs, covering the area around the ABVC (Figure 1) with minimum and maximum interstation distances of 5 km and 7 km, respectively. There are 300 potential station pairs for cross correlation. We deployed 3 Geobit, 2 Guralp, and 18 Nanometrics Compact Trillium broadband seismometers, and 2 Lennartz shortperiod seismometers (Supplementary Table S1). The 3 Geobit sensors are in boreholes and the remaining sensors are surface deployments. Geobit and Guralp instruments are integrated sensor-recorder packages, while data from the remaining sensors are recorded by low-power dataloggers manufactured by the Australian National University. All waveform data are recorded with a 100 Hz sampling rate. We dug pits to keep the seismometer under the ground surface for thermal stability and use solar panels to recharge the battery voltage so that the device can record continuously. A photograph of one of these seismic recorders is shown in Figures 1C-G.
The Geobit and Nanometrics Trillium seismometers are broadband sensors, both with response flat to velocity in the period range 0.1-20 s. The response of the Lennartz short-period sensor, on the other hand, is flat only over the range 0.03-1 s. The manufacturer of the Lennartz claims that the instrument response is sufficiently precise to be reliably deconvolved to periods as long as 20 s (Lennartz Electronic, 2017), longer than the maximum period used here. An example of a cross correlogram between short-period seismogram and broad-band seismogram is shown in Supplementary Figure S1. If some correlogams using data from the Lennartz seismometers are contaminated with instrument noise for the longer periods used in our study (up to 12 s), we expect they will be rejected when we screen the correlograms to keep only those with high signal-to-noise ratio (see Section "Ambient Seismic Noise Cross Correlation").
Although each station records 3-components of motionvertical, North and East -for the analysis discussed here we utilized only the vertical component data. Because several different instrument types were used in the Agung Seismic Experiment, we first corrected the instrument response before calculating cross correlograms between each station pair to be used in our ANT analysis. We used a pole-and-zero representation of the frequency response of each different instrument to obtain the corresponding displacement.

Ambient Seismic Noise Cross Correlation
Ambient noise cross correlation was introduced by Shapiro et al. (2005) as a way of extracting from seismic "background" noise signals representing the elastic impulse response of the subsurface. This technique is now widely used to image subsurface structure using group velocity inversion. ANT used to delineate the subsurface structure in Australia Kennett, 2010, 2012), Netherlands (Yudistira et al., 2017). Furthermore, in Indonesian region the ANT used to understand upper crustal structure beneath central Java (Zulfakriza et al., 2014), Jakarta the capital city of Indonesia (Saygin et al., 2016), eastern Java (Martha et al., 2017), western Java (Rosalia et al., 2020), and Bandung Basin (Pranata et al., 2019). On the other hand, ANT also used to delineate structures within volcano bodies, such as Lake Toba, Indonesia (Stankiewicz et al., 2010), Mt. Asama in Japan (Nagaoka et al., 2012) and Colima volcano complex, Mexico (Escudero and Bandy, 2017). In this study, we applied the ANT technique to delineate the Rayleigh wave group velocity structure of the ABVC.
Cross-correlating the ambient seismic noise recorded simultaneously at a pair of stations over many days (∼60 days in our study) makes use of the coherent noise produced by the interaction of the ground surface with the atmosphere and ocean wave activity, to produce a signal that approximates the impulse response of the medium, or Green's Function, between the two stations (Bensen et al., 2007). This estimate of the Green's Function obtained from cross correlation of interstation noise is referred to here as the Time Domain Empirical Green's Function, or TDEGF (Snieder, 2004). When the vertical component of motion at the two sites is cross correlated, the resulting TDEGF can be used to measure Rayleigh wave group velocities at different periods (Shapiro and Campillo, 2004;Sabra et al., 2005;Shapiro et al., 2005). The procedure used to estimate Rayleigh-wave TDEGFs follows the steps carried out by Saygin and Kennett (2010). The cross correlation between two stations A and B can be written as follows, Where V(X A , t) and V(X B , t) are signals recorded continuously by two different stations at positions X A and X B , respectively. We used 60 days of observed data and divided each into several segments. Each segment has a duration of 3,600 s with overlaps of 1,200 s for each segment. After that, we stacked the result to get an estimate of the TDEGF for each day of observation, so that 1-day TDEGFs were finally obtained for all observation days, and these were finally stacked to get a single TDEGF for each interstation pair.
Although our dataset was of only 60 days duration, we found that a collection of TDEGFs associated with a particular station, when sorted according to inter-station distance, displays a clearly identifiable arrival of energy with an average moveout of about 2 km/s, which we assume to be the fundamentalmode Rayleigh wave (Figure 2A). Following Saygin and Kennett (2010) and Zulfakriza et al. (2014), for each correlogram we selected a window that encompassed this energy, while rejecting correlograms that were too noisy to unambiguously identify the onset and duration of the Rayleigh wave. We note that more than half (158) of the 300 available correlograms were rejected as too noisy, including some with data at long period from the short-period Lennartz seismometers. A narrow-band Gaussian is applied to each these TDEGFs, and group velocity arrivals for each period are manually picked from time-frequency plots over periods of 0.5-15 s, as indicated with white dots in Figures 2B,C.

Sub-Space Inversion of Group Velocity
The TDEGF for each station pair can provide measurements of group velocity for different periods by applying the Multiple Filter Technique (MFT), popularized by Dziewonski et al. (1969). This technique has been widely applied in previous studies such as Saygin and Kennett (2010), Zulfakriza et al. (2014), and Martha et al. (2017). In this way group velocity dispersion curves are obtained for each interstation pair in the network of seismographs, giving rise to a web of intersecting inter-station paths along which dispersion is measured. Tomographic techniques can then be applied to provide an image of group velocity for select periods. In many cases these maps of group velocity anomalies can be directly interpreted in terms of lowvelocity sedimentary basins or high-velocity metamorphic rock formations.
In this study, we used the surface wave tomography approach developed by Rawlinson (2005), which uses a subspace inversion method to obtain the optimal group velocity map along with the Fast Marching Method (FMM) to trace surface wave raypaths in heterogeneous media. Subspace inversion is an iterative method in which an objective functional is minimized by using successive quadratic approximations of the function in an n-dimensional subspace (Kennett et al., 1988). Once the minimum of the quadratic approximation is found in the model subspace, a new quadratic approximation is made, either in the same or a different model subspace, and the process repeated. At each step the rays are re-traced using FMM, so the non-linear relationship between velocity and travel time is accounted for. The method has been successfully applied to obtain the group velocity structure in various regional studies (Saygin and Kennett, 2012;Martha et al., 2017). According to Rawlinson et al. (2010) the incorporation of FMM and subspace inversion can provide stable inversion results even for strongly heterogeneous media.
The inversion problem is solved by optimization the objective function S (m), which includes misfit, damping and smoothing terms as indicated in equation (2) below.
where m is a vector of model parameters (i.e., the group velocity map), g(m) are the predicted group traveltimes from the model, d are the observed group traveltimes, C e is the data error covariance matrix, m 0 is the reference model, C m is the model parameter covariance matrix, D is a flatness/smoothness matrix, ε is the damping parameter, and η is the smoothing parameter.

Shear Wave Velocity Inversion
There are several ways to estimate velocity-depth profiles from surface wave dispersion curves. Shapiro et al. (2000) estimated the model uncertainty using a Monte Carlo scheme, which performs non-linear inversion of dispersion curves using a random search of the model space. Herrmann (2013) used a linearized method to find the stack of horizontally homogeneous layers that best fit dispersion curve. Throughout the inversion, the thickness and Vp/Vs for each layer remained fixed, and the density was estimated from the Vp. We applied the scheme of neighborhood algorithm introduced by Sambridge (1999a,b) to address the inversion problem. This algorithm uses a random search that exploits information accumulated over the history of the model space search to adapt the sampling, concentrating it in parts of the model space surrounding the lower misfit models. We used the AK135 model developed by Kennett et al. (1995) as an initial upper crustal mode; over depths 500 -15000 m, with maximum Vp 5.8 km/s, and maximum Vs 3.46 km/s.

RESULTS
Here we consider the application of tomographic inversion to the fundamental-mode Rayleigh wave TDEGFs estimated from inter-station cross correlograms of ambient seismic noise on the ABVC, as well as Neighbourhood Algorithm (NA) inversion of the resulting dispersion curves. We first consider an analysis of the lateral resolution achievable with our experimental configuration, followed by a discussion of the images of Rayleigh wave group velocity obtained by tomographic inversion of the observed travel times at each period. Finally, we apply NA inversion to the dispersion curves obtained from the tomographic images, to obtain 1D Vs to depth profiles at each point of a regular grid covering the study area, and discuss the implications for 3D Vs structure of the ABVC.

Checkerboard Resolution Test
We conducted tests using "checkerboard" group velocity patterns of alternating low-and high-velocity cells to generate synthetic Rayleigh wave travel times for the inter-station paths of our dataset, which we then invert to assess how well the original velocity pattern is recovered. From such tests we determine the optimum parameterization and achievable resolution of the inversion. In this paper, the checkerboard resolution tests were used with noiseless data.

Tomographic Group Velocity Model
Following the inversion procedure described above, Rayleigh wave group velocity tomograms at different periods were obtained. The damping and smoothing hyperparamters -and η, respectively, in Eq. 2 -were each assigned a value of 50, which we found provided acceptable trade-offs between the final model's deviation from its initial value and its smoothness, and the data misfit. Figure 4 shows 2-D maps of tomographic group velocity models obtained for periods of 1, 3, 5, 8, 10, and 12 s, with insets depicting the ray-path coverage. In this study, longer periods (more than 10 s) have poor resolution due to sparse ray-path coverage.
The most prominent feature of the tomograms at different periods plotted in Figure 4 is the low velocity anomaly between Mt. Agung and Mt. Batur. This is a significant low velocity anomaly, with group velocity about 500 m/s at period 3 s (Figure 5). The position of this anomaly coincides with the cluster of earthquake epicenters between Mt. Agung and the Batur Caldera at depths less than 10 km estimated by Syahbana et al. (2019).

Shear Wave Velocity Model
We applied the NA as described by Sambridge (1999a,b) to invert group velocity dispersion curves. The group velocity curves were extracted for a set of regularly spaced grid points covering the study area, and inverted to create 1-D depth profiles of Vs (Supplementary Figure S2). There are 59 points that covered the study area at 5 km intervals. The group velocity at these spatial grid points retrieved from tomograms at periods between 0.5 and 12 s. The dispersion curve at each grid point was inverted using the program "dinver" in the Geopsy software package developed by Wathelet et al. (2008), to produce Vs profiles having a depth variation of velocity consistent with the observed dispersion curves. We allowed dinver to consider 4-layer model, with each layer of variable thickness and velocity range, and each having five "sublayers" over which the layer's top and bottom velocities are linearly interpolated. The allowed layer parameters are Vs ranges of 500-2000, 1000-3000, and 1000-3000 m/s and thickness ranges of 500-2000, 2000-5000, and 5000-8000 m, respectively, for the top three layers, and 2000-4000 m/s for the Vs of the underlying half-space.
We interpolated the minimum-misfit 1-D Vs profiles produced by the NA inversion at each grid point to create a 3-D model of the Vs structure of the ABVC. The resulting maps of Vs for several depths (0.5-8 km) are plotted in Figure 5. The spatial pattern of low and high velocity at each depth is roughly consistent with that of the Rayleigh wave group velocity at different periods shown in

DISCUSSION
The group velocity and shear velocity model as shown in Figures 4, 5 show the variation of velocity structure beneath the ABVC. A prominent low-velocity anomaly can be seen at periods u to 3 s for the group velocity maps, and depths up to 2 km for Vs, but disappears at greater periods and depth. We did synthetic test using initial model then conducted the inversion based on initial model (Supplementary Figure S3). The low velocity anomaly appears in between Mt Batur and Mt Agung. Vertical cross-sections of velocity structure in W-E and S-N directions can be seen in Figure 6. Both of these vertical cross-sections intersect the low velocity (∼0.5-2.2 km/s) anomaly between Mt. Agung and Mt. Batur that extends from the surface to about 2 km depth. In Figure 6 we overlay the Vs tomographic map with event hypocenters in the interval December 2018 to June 2019, to show the positive correlation between the low velocity anomaly and the hypocenter distribution. The hypocenters were estimated using the Non-Linloc method (Lomax et al., 2014), with manually picked P-and S-wave phases (Supplementary Figure S4) and the 1D seismic velocity model of Ak 135 (Kennett et al., 1995). We estimated the locations of 307 events that occurred during the seismometer deployment, with average location uncertainty: 1.70 km (EW), 1.75 km (NS), and 1.99 km (Z), for which the uncertainty histograms are shown in Supplementary Figure S5. Some of the hypocenters are at shallow depth near the low velocity anomaly, but much of the seismicity extends to greater depth beneath the low-velocity anomaly (Figure 6).
On the one hand, the low velocity observed at such shallow depth between Mt. Agung and Mt. Batur could simply reflect the accumulation of low-Vs volcanic deposits in this topographic "saddle" between the volcanoes. This could explain how the anomalous Vs could be so low, at 1000 m/s. However, this interpretation fails to explain the 2 km thickness of the low Vs anomaly, or the clustering of seismicity that extends from depth to the low Vs anomaly at 2 km depth (Figure 6).
Another explanation for the low velocity anomaly between Mt. Agung and Mt. Batur is that it could reflect magmatic processes. From an InSAR analysis of Sentinal-1 SAR data prior to the September 2017 eruption, Albino et al. (2019) inferred magma intrusion into a subvertical dike aligned with and located about midway along the axis connecting Mt. Batur and Mt. Agung -i.e., directly below the Vs anomaly in Figure 5, but at 7-13 km depth. Syahbana et al. (2019) found that both preeruptive seismicity and displacements at GNSS stations on the ABVC are consistent with this interpretation, and also note the potential for pressurization of groundwater above the dike. The shallow, low-Vs anomaly we observed is much shallower than the depth of the dike, so it is unlikely to be caused by magma intrusion, but it may be explained by overpressurization of groundwater. As described by Syahbana et al. (2019), meteoric water seeping down from the summit of Mt. Agung might interact with magmatic fluid and/or gases rising from the deep magmatic intrusion, leading to over-pressurization. Such overpressurization is consistent with both reduced Vs (Vanorio et al., 2005;Brenguier et al., 2008) and the occurrence of earthquakes (Coulon et al., 2017).

CONCLUSION
The evolution of our knowledge of the ABVC has involved several studies to understand volcanic processes beneath Mt. Agung and Mt. Batur. This paper presents the first localscale study of Rayleigh wave group velocity structure of the ABVC, delineating the structure using ambient seismic noise tomography, which provides good resolution beneath and between Mt. Agung and Mt. Batur.
We estimated inter-station Rayleigh wave Green's functions using ambient seismic noise cross correlation from available station pairs across the ABVC. Observed Rayleigh waves contain energy in the period range of 1-10 s, and arrival times of these surface waves were picked for available station pairs (see Figure 3).
2-D tomographic maps obtained from the inversion of the group velocity of Rayleigh waves and the Vs inverted from the resulting dispersion curves clearly delineate a strong low velocity anomaly as deep as 2 km depth between Mt. Agung and the Batur Caldera, and this coincides with a clustering of seismicity observed during the experiment. We suggest this low Vs anomaly can be explained by a combination of: (1) the accumulation of low-Vs volcanic deposits in the topographic saddle between Mt. Batur and Mt. Agung, and (b) the presence of over-pressurized hydrothermal fluids, possibly sourced from magma intrusion into a dike or system of dikes at depth >7 km inferred by Albino et al. (2019) and Syahbana et al. (2019). Of these explanations, only the presence of over-pressurized fluids provides a mechanism for generating the seismicity near the low Vs anomaly we imaged at <2 km depth.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
ZZ, AN, SW, PC, and DPS conceived the survey and seismic study on ABVC. ZZ, DPS, and SR contributed to the writing of the manuscript. All authors contributed to the preparation of the manuscript, and read and approved the final manuscript.