# Imaging of a magma system beneath the Merapi Volcano complex, Indonesia, using ambient seismic noise tomography

T. Yudistira, J. -P. Metaxian, M. Putriastuti, S. Widiyantoro, N. Rawlinson, F. Beauducel, Z. Zulfakriza, A. D. Nugraha, A. Laurin, A. A. Fahmi, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

T. Yudistira, J. -P. Metaxian, M. Putriastuti, S. Widiyantoro, N. Rawlinson, et al. Imaging of a magma system beneath the Merapi Volcano complex, Indonesia, using ambient seismic noise tomography. Geophysical Journal International, 2021, 226, pp.511-523. 10.1093/gji/gab104. insu-03590038

HAL Id: insu-03590038

https://insu.hal.science/insu-03590038

Submitted on 26 Feb 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Imaging of a magma system beneath the Merapi Volcano complex, Indonesia, using ambient seismic noise tomography

T. Yudistira,<fnref n="1" /> J.-P. Metaxian,<fnref n="2" />,<fnref n="3" /> M. Putriastuti,<fnref n="4" /> S. Widiyantoro,<fnref n="1" />,<fnref n="5" /> N. Rawlinson,<fnref n="6" /> F. Beauducel,<fnref n="2" />,<fnref n="3" /> Z. Zulfakriza,<fnref n="1" /> A. D. Nugraha,<fnref n="1" /> A. Laurin,<fnref n="2" /> A. A. Fahmi<fnref n="2" /> and A. Budi-Santoso<fnref n="7" />

Accepted 2021 March 12. Received 2021 March 4; in original form 2020 August 10

## SUMMARY

Mt Merapi, which lies just north of the city of Yogyakarta in Java, Indonesia, is one of the most active and dangerous volcanoes in the world. Thanks to its subduction zone setting, Mt Merapi is a stratovolcano, and rises to an elevation of 2968 m above sea level. It stands at the intersection of two volcanic lineaments, Ungaran–Telomoyo–Merbabu–Merapi (UTMM) and Lawu–Merapi–Sumbing–Sindoro–Slamet, which are oriented north–south and west–east, respectively. Although it has been the subject of many geophysical studies, Mt Merapi’s underlying magmatic plumbing system is still not well understood. Here, we present the results of an ambient seismic noise tomography study, which comprise of a series of Rayleigh wave group velocity maps and a 3-D shear wave velocity model of the Merapi–Merbabu complex. A total of 10 months of continuous data (October 2013–July 2014) recorded by a network of 46 broad-band seismometers were used. We computed and stacked daily cross-correlations from every pair of simultaneously recording stations to obtain the corresponding inter-station empirical Green’s functions. Surface wave dispersion information was extracted from the cross-correlations using the multiple filtering technique, which provided us with an estimate of Rayleigh wave group velocity as a function of period. The group velocity maps for periods 3–12 s were then inverted to obtain shear wave velocity structure using the neighbourhood algorithm. From these results, we observe a dominant high velocity anomaly underlying Mt Merapi and Mt Merbabu with a strike of 152°N, which we suggest is evidence of old lava dating from the UTMM double-chain volcanic arc which formed Merbabu and Old Merapi. We also identify a low velocity anomaly on the southwest flank of Merapi which we interpret to be an active magmatic intrusion.

**Key words:** Asia; Tomography; Volcanic structure; Seismic noise.

## 1 INTRODUCTION

Several sections of the Quaternary volcanic arc in Java are characterized by the presence of two (or more) volcanoes perpendicular to the arc (Tatsumi & Eggins 1995). This feature is often referred to as a double-chain volcanic arc, which can be divided into the trench-side and backarc-side volcanoes. The Ungaran–Telomoyo–Merbabu–Merapi (UTMM) volcano group is an example of a double-chain volcanic arc. Kohno et al. (2005) examined lava that had erupted

from each volcano in the UTMM volcano group to determine its petrology and rock composition. They determined that the volcanic activity tends to migrate southward, with Ungaran acting as the back-arc side chain due to its deeper magma source compared to the Telomoyo, Merbabu and Merapi Volcanoes, which act as the trench-side chain. In this study, we focus on the adjoining Merapi and Merbabu volcanoes; Merapi is the youngest volcano in the UTMM volcanic arc and likely has strong seismic anomalies that reflect an active magma plumbing system (e.g. Ramdhan et al. 2019).

<footnote n="1">¹Global Geophysics Research Group, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung, Bandung 40132, Indonesia. E-mail: ilikwidi@gmail.com</footnote>

<footnote n="2">²Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTERRE, 38058 Grenoble, France</footnote>

<footnote n="3">³Université de Paris, Institut de physique du globe de Paris, CNRS, 75005 Paris, France</footnote>

<footnote n="4">⁴Geophysical Engineering Graduate Program, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung, 40132 Bandung, Indonesia</footnote>

<footnote n="5">⁵Faculty of Engineering, Maranatha Christian University, Bandung 40164, Indonesia</footnote>

<footnote n="6">⁶Department of Earth Sciences - Bullard Labs, University of Cambridge, Cambridge, CB30EZ, United Kingdom</footnote>

<footnote n="7">⁷Center for Volcanology and Geological Hazards Mitigation, Jl. Cendana No. 15, Yogyakarta 55166, Indonesia</footnote>
There is also evidence of changes in magma composition between Old Merapi and New Merapi (Kohno et al. 2005), which may reveal itself as a change in shear wave velocity. The adjacent Merbabu Volcano is considered to be a dormant stratovolcano (last erupted in 1797), which rises to 3142 m above sea level. Neuman van Padang (1951) showed that the Merbabu volcanic edifice consists of basalt (composed of olivine-augite minerals), andesite with augite minerals, and also andesite with hornblende hipersten-augite minerals, which is consistent with its tectonic setting.

Merapi Volcano (7°32′26″S, 110°26′48″E), one of the largest volcanoes in the Sunda Arc, is a ~3000-m-high stratocone located in Central Java. It was formed in the late Pleistocene-Early Holocene Era (Kohno et al. 2005) as a result of subduction of the Indo-Australian Plate below the Eurasian Plate. It is bounded by the intersection of two volcanic lineaments: Ungaran-Telomoyo-Merbabu-Merapi (N164°E) and Lawu-Merapi-Sumbing-Sindoro-Slamet, and also by two major faults: the Merapi-Merbabu (north-south) and the Baribis-Semarang-Kendeng, as shown in Fig. 1. In the early stages of its growth, Merapi experienced an effusive eruption of basaltic magma; this, followed by subsequent eruptions, has allowed the K₂O concentration of the lavas to be tracked over time. Old Merapi lavas have a lower K₂O concentration compared to either its younger counterpart, or the Telomoyo and Merbabu volcanoes (Kohno et al. 2006). The composition of the magma has changed in time from containing medium levels of K₂O to having high levels of K₂O, which results in a more silicic and viscous lava.

A detailed study of the evolution of Merapi Volcano has been carried out by Gertisser (2012), who provides a synopsis of all geological research from the eighteenth century. It is thought that Merapi Volcano was preceded by two Proto-Merapi volcanoes: Mt Bibi, a basaltic andesite volcano on the northeastern flank and Mt Turgo and Mt Plawangan, two almost identical basaltic edifices 5 km southwest of the New Merapi summit. Mt Batulawang and Mt Kendil also form part of the Somma-Merapi sequence.

This study describes an application of the ambient noise tomography method to produce Rayleigh wave group velocity maps and a 3-D Vs model of the Merapi-Merbabu Complex, with the goal of trying to better understand its volcanic plumbing system, and help determine whether a sizable melt source or magma chamber is present in the upper crust.

Ambient noise tomography has been carried out across a range of scales, from global (e.g. Sager et al. 2020) to regional (e.g. Yang et al. 2007; Behm et al. 2016) to local (e.g. Young et al. 2011; Hannemann et al. 2014), but it has also proven itself to be a useful tool for investigating the seismic structure beneath active volcanoes (e.g. Stankiewicz et al. 2010; Benediktsdóttir et al. 2017). For example, Huang et al. (2018) used ambient noise tomography to image beneath the Aso Caldera in Japan. By using periods as short as 1 s, they were able to recover kilometre-scale structure in the uppermost crust to a depth of 5–6 km. Low velocities in the top 1–2.5 km were interpreted to be evidence of hydrothermal reservoirs, whereas high velocities below 3 km depth were thought to be evidence of consolidated igneous rock. A return to low velocities at 5–6 km depth was interpreted to represent an accumulation of magma at the base of the upper-crustal volcanic plumbing system. Most ambient noise tomography studies exploit the surface wave dispersion characteristics of Rayleigh and/or Love waves. Typically, this dispersion information is inverted to extract either group or phase velocity maps which are then inverted for 3-D shear wave speed structure (e.g. Young et al. 2013a; Crowder et al. 2019). However, the approach of directly inverting dispersion data for 3-D

shear wave velocity structure in a single step is gradually becoming more popular (e.g. Zhang et al. 2019). In the following study, we apply the more traditional two-step method outlined above to image the mid-upper crustal structure beneath the Merapi-Merbabu volcanic complex, by exploiting ambient noise data from a temporary seismic deployment that spanned the two volcanoes.

## 2 DATA

The ambient noise tomography analysis was carried out using ten months of continuous data (October 2013–July 2014) recorded by the DOMERAPI array, which consisted of 46 three-component broad-band seismic stations distributed around Mt Merapi and Mt Merbabu (Fig. 1). Average station spacing is ~4 km, resulting in a dense array which spans approximately 50 km in latitude and 40 km in longitude. Each station included a Taurus digitizer and a Guralp 40T seismometer belonging to the French national pool of portable seismic instruments (SISMOB/RESIF). The geometry of the DOMERAPI array is well suited to the use of ambient noise tomography to image the connection between the two volcanic edifices, with the ambient noise largely generated by oceanic microseisms from the surrounding seas. The relatively narrow aperture of the array limits the maximum period of surface waves that are extracted to ~10–12 s, as determined below, which limits the imaging to mid-upper crustal depths (down to about 10 km).

## 3 METHODOLOGY

Ambient seismic noise interferometry relies on the extraction of empirical Green's functions (EGF) from simultaneously recording pairs of stations. This is typically done through cross-correlation of long-term recordings of 'background noise', which can produce interstation EGFs between all pairs of operating stations in a network.

In practice, ambient seismic noise imaging typically involves four main stages (e.g. Yudistira et al. 2017; Bensen et al. 2007; Crowder et al. 2019, among others): (1) Green's function calculation, (2) interstation group or phase velocity determination, (3) group or phase velocity tomography and (4) inversion for shear wave velocity (Vs) as a function of depth on a dense grid of points that span the geographic limits of the model in order to construct a pseudo 3-D Vs model. In this study, we choose to use group velocity, since we found its dispersion characteristics easier to extract from the cross-correlated waveforms compared to phase velocity.

### 3.1 Green's function calculation

In order to compute the EGF for each station pair we follow the procedure proposed by Bensen et al. (2007). Prior to EGF calculation, several standard processing steps were applied to the recorded data. Here, we briefly explain these steps and the EGF calculation method, although a more detailed description of the approach we adopt can be found in Bensen et al. (2007).

The first step of the data preparation stage consists of partitioning the data into 24-hr-long (or daily) segments, and then removing the mean and trend of each segment; this is a standard technique to permit efficient data handling. In the second step, in order to reduce the amplitude dominance of signals from events such as large earthquakes, one-bit time domain normalization was carried out. This step also involves spectral whitening in order to balance the amplitude in the frequency domain.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Top panel: the station distribution of the DOMERAPI network (black inverted triangles) used in this study. The mobile installation consists of three component seismometers deployed at 51 different locations surrounding Mt Merapi, Mt Merbabu and Mt Telomoyo. Bottom panel: regional location of the study area. Red lines indicate known active faults.</figcaption></figure>

The next step is cross-correlation and stacking of all data for each station pair. Cross-correlations were performed by dividing the daily data streams into 1-hr segments in order to compute cross-correlation functions (CCFs). The hourly CCFs for each station

pair are then stacked, to obtain the interstation EGF. Stacking enhances the coherent part of the EGF and suppresses incoherent signals. Due to differences in the recording periods of the seismometers in the DOMERAPI array, stations pairs generally do not have
completely overlapping data, and thus it is necessary to exclude periods in which mutual recording does not occur. After applying a manual selection process, which is based on the visual clarity of dispersive wave train signals in the cross-correlations as a function of interstation distance (e.g. Fig. 2b), we obtain about 1158 inter-station CCFs from the vertical component data, which corresponds to Rayleigh wave particle motion. Fig. 2 shows the filtered (period 3–15 s) cross-correlation gather of station ME41 with respect to all other stations in the network. The surface wave trains are clearly apparent on the positive time lag, but appear to be largely absent on the negative time lag. This corresponds to a dominant noise source from north of the DOMERAPI network, which may be related to the dynamics of the Java Sea. The same conclusion on the directivity of the dominant sources of noise is also reached by Yudistira & Widiyantoro (2016) using MERAMEX network data. Such directivity could yield artefacts in the measurement of interstation group velocities in certain period bands (e.g. Pedersen & Kruger 2007), although we see no evidence of this in our study.

## 3.2 Dispersion curve measurement

The dispersion characteristics of the EGFs are extracted by applying frequency-time analysis using a multiple filtering technique (MFT) developed by Dziewonski et al. (1969). The method involves the application of Gaussian bandpass filters to the dispersive signals across a range of central frequencies. The group arrival time of the specific center frequency is estimated by the maxima of the filtered signal envelope. A knowledge of the arrival time and the interstation distance can be used to estimate path-averaged group velocity for each station pair. Fig. 3 shows an example of a frequency-time representation of an EGF for station pair ME09-ME36; the white squares correspond to the automatic pick of the maximum amplitude of the envelope at each central frequency. Each dispersion plot is visually inspected to ensure that the picked curve is continuous from long to short periods, and represents the fundamental rather than higher modes (which are often not visible at longer periods and are faster than the fundamental mode at the same period).

Bensen et al. (2007) found that a reliable dispersion measurement for a given period requires an interstation spacing of at least three wavelengths. The maximum interstation distance in the DOMER-API array is about 50 km, and the estimated group speed in the region is ~1.5 km s⁻¹; thus, the maximum period we can safely exploit is ~12 s. We can clearly observe the degradation of the dispersion measurement at periods greater than this maximum cut-off value, at least for group velocities. Therefore, we use a 3–12 s period range for this study, or about 0.1–0.3 Hz in frequency, which still includes the typical frequency range of oceanic microseisms (0.03–0.3 Hz). One important factor in this analysis is the choice of the Gaussian constant, which controls the resolution of the time and frequency domain in MFT. It is impossible to simultaneously obtain the best frequency and time resolution: we can only manage an optimum trade-off value between both domains. There is no ideal strategy for obtaining the best Gaussian constant value, so trial and error was used. Next, we carefully selected the inter-station dispersion curve by visually inspecting the causal, acausal and symmetric (average of the acausal and causal signal) dispersion characteristics for each station pair for the period range 3 to 12 seconds, and choosing the one that is most clear. The selection process is based on a visual assessment of the signal-to-noise ratio of the CCF (see Fig. 3) and the continuity of the resulting dispersion curve. During the selection process, we discarded some of the station pairs

which did not have a clear pattern. This process led to a different number of dispersion curves for each period since every station pair will produce different results. As a consequence, the final data set consists of 206 causal, 433 symmetric and 481 acausal dispersion curves. This means that from all possible station pairs, more than 50 per cent had poor quality dispersion curves that were rejected.

## 3.3 Surface wave tomography

The selected interstation group velocity dispersion data are inverted to obtain a set of group velocity maps in the period range 3–12 s, with an interval of 0.5 s. The inverse problem was solved using a regularized least squares method. The group velocity model was parametrized by 2-D rectangular blocks, and both damping and smoothing were applied in the inversion to obtain parsimonious velocity models. Similar to many previous studies (e.g. Ritzwoller & Levshin 1998; Barmin et al. 2001; Yoshizawa & Kennett 2002; Yudistira et al. 2015) we assume that rays follow a great circle path such that the inverse problem is linear.

The solution of the inverse problem involves minimizing a least-squares objective function of the form

$$ \|W_d^{1/2} (d_{\text{obs}} - d_{\text{cal}})\|_2^2 + \|W_m^{1/2} m\|_2^2, \quad (1) $$

where the first term of the right-hand side is the data residual norm and the second term is the model norm. $W_d = \text{diag}(\sigma_{d1}^{-1} \sigma_{d2}^{-1} \dots \sigma_{dN}^{-1})$ and $W_m = \text{diag}(\sigma_{m1}^{-1} \sigma_{m2}^{-1} \dots \sigma_{mN}^{-1})$ are diagonal matrices where $\sigma_{d1}$ and $\sigma_{m1}$ are the data and model covariances respectively, and $W_m$ regulates the characteristics of the solution model. $d_{\text{obs}}$ is the data observation vector which contains the inter-station traveltimes and $d_{\text{cal}}$ is vector of predicted traveltimes through the current model for a specific frequency. For a purely linear inverse problem, the data perturbations are related to the model perturbations $m$ as $Gm$ where $G$ is the forward operator matrix.

The Thikonov solution for minimizing equation (1) can be expressed as $m = (G^T W_d G + W_m)^{-1} G^T W_d d_{\text{obs}}$. $W_m$ consists of a model damping and a smoothing matrix, where $W_m^{1/2} = \lambda_1 L_1 + \lambda_2 L_2$ with $L_1$ a diagonal matrix of ray hit counts in each block such that more damping is applied to cells with fewer rays passing through them, and $L_2$ a matrix of the discrete Laplacian operator (which imposes smoothing). To reduce the number of tuning parameters, we use the same value of $\lambda_1$ and $\lambda_2$ so that the task is to find a single regularization value that balances the model variance and the data variance by manual investigation of trade-off curves. The final solution can then be written as:

$$ \mathbf{m} = (G^T \mathbf{W}_d \mathbf{G} + \lambda^2 [\mathbf{L}_1^T \mathbf{L}_1 + \mathbf{L}_2^T \mathbf{L}_2])^{-1} G^T \mathbf{W}_d \mathbf{d}_{\text{obs}} \quad (2) $$

which in matrix form is equivalent to the least-squares solution of

$$ \begin{bmatrix} W_d^{1/2} G \\ L_1 \\ L_2 \end{bmatrix} m = \begin{bmatrix} W_d^{1/2} d_{\text{obs}} \\ 0 \\ 0 \end{bmatrix}. \quad (3) $$

The tomographic solution model is dependent on a number of initial parameters, including choice of the starting velocity model, and the level of damping and smoothing that is applied. In order to investigate the influence of these parameters, trial and error testing was performed with different sets of damping, smoothing, and number of blocks to represent structure.

To illustrate the potential of the data to recover structure, a simple synthetic checkerboard test (SCT) was carried out. The checkerboard model consists of a pattern of low and high velocity perturbations and a synthetic data set is generated by solving the forward
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The bandpass filtered (3–15 s period) cross-correlation functions as a function of inter-station distance from station ME41 to all other stations. The dispersive wave trains correspond to fundamental mode Rayleigh waves and can be extracted to determine seismic velocity.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Frequency–time analysis of the cross-correlation for station pair ME09-ME36. The picked high amplitude envelope represents the fundamental mode group velocity dispersion. The left-hand panel corresponds to the normalized filtered (period 3–15 s) cross-correlation function (CCF) and the top panel is its normalized period-dependent signal-to-noise ratio (SNR). The CCF and SNR were amplitude-normalized between -1 to 1 and 0 to 1 respectively, relative to their peak values.</figcaption></figure>

problem in the presence of this model using the same set of receiver pairs as the observational data set. This synthetic data set is then inverted using the same input parameters as the real data inversion to provide insight into which regions of the model are robustly recovered. In our SCT, the checkerboard pattern is represented by 2×2 blocks of the base parametrization. In order to determine an optimal number of cells, we performed the SCT by varying the number of blocks in each dimension. Fig. 4 shows the results for 16 × 16, 24 × 24 and 32 × 32 blocks for period 3 and 7 s. Based

on these results we choose to use 24 × 24 blocks to describe the model, since it is the finest discretization that is still able to recover a clear checkerboard pattern. The SCT results for all periods using this parametrization choice are shown in Fig. 5. As expected, the central area, where the ray path coverage is dense, produces a good recovery of the input checkerboard. Although the span of the resolved area varies slightly between periods due to differences in ray path coverage, we use the same region of each model in order to produce the 3-D shear wave velocity model.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The checkerboard test using different cell sizes for periods 3 and 7 s. Based on these results, we choose to represent structure in our group velocity maps using 24 × 24 cells. The interstation paths for the corresponding periods are shown in the third row of the figure.</figcaption></figure>

The regularization parameter controls the trade-off between data fit, model complexity, and perturbations from the starting model. The most common approach for selecting the value of such parameters is via *ad hoc* trade-off tests (Barmin *et al.* 2001). We found that trade-off curves for each period do vary, but not by a large amount (Fig. 6) and we choose regularization values ranging between 4–7. The resultant group velocity maps are shown in Fig. 7 and will be discussed together with the shear wave velocity model.

## 3.4 Shear velocity inversion

The last stage of the modelling procedure is the inversion of the group velocity maps to obtain the 3-D distribution of shear wave velocity ($V_s$). This requires discretely sampled pseudo dispersion curves to be extracted from a regular grid of points which span the group velocity maps. For each group velocity dispersion curve

produced by this process, a 1-D inversion for shear velocity as a function of depth can be carried out (e.g. Young *et al.* 2013b). Here, this was achieved by using the daver toolbox from the Geopsy package by Wathelet *et al.* (2004), which implements the nonlinear Neighbourhood Algorithm (NA) inversion scheme of Sambridge (1999). NA iteratively samples the model space by partitioning it into Voronoi cells based on a random distribution of points. The distribution of misfit values in these cells is then used to guide the search process for data fitting models. The best-fitting 1-D model is chosen at each gridpoint, and a pseudo 3-D model is obtained by combining them through simple linear interpolation. Based on the 3–12 s period range of dispersion curves, the maximum depth of good recovery will be around 10 km. The 1-D $V_s$ models are parametrized in terms of constant velocity layers of variable thickness. In this case, after running a series of experiments using different numbers of layers, we chose to define the crust by five layers, which results in a total of nine parameters (five velocity and four
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Checkerboard resolution tests for selected periods (3–12 s) using 24 × 24 cells. Based on this synthetic test the central region, which is the focus of our interpretation, is considered well resolved.</figcaption></figure>

thickness). In general, adding additional layers to the model did not result in an improved fit to the data. $V_p$ and density are also permitted to vary in the inversion, but due to their low sensitivity to the observables, they are coupled to $V_s$ via a linear relationship, and hence not considered in the interpretation. Furthermore, since

we have no prior information about $V_s$ structure in our study area, we allow the NA to perform a relatively broad search, from 100 to 2500 m s⁻¹, across the full depth range of the model.

The two primary NA ‘tuning’ parameters are the number of model samples (ns) and the number of subsequent resamples (nr) at each iteration. Following detailed experimentation using a range of values, we finally settled on using $nr = ns = 100$ for 1000 iterations, so that in the end we have 100 000 sampled models. The best-fitting model of the complete ensemble at each point was then chosen for input into the pseudo 3-D model. Fig. 8(a) shows a typical example of a 1-D shear wave velocity model ensemble and the associated group dispersion (both observed and predicted). The best-fitting model—shown as a black dashed line—is used as input in the pseudo 3-D model. Here, the region of good fit is generally well defined, and is narrowest in the upper 1–2 km of the model. The largest velocity jump (~0.3 km s⁻¹) occurs at a depth of ~4 km; based on the spread of models, this seems to be a robust feature, although we note that the parametrization requires velocity jumps rather than smooth gradients. Since surface waves cannot discriminate between either scenario, we acknowledge that this is an end-member model, and a smoothed version of this model that does not contain discontinuities is equally plausible as far as the data are concerned. In the pseudo 3-D model, this choice of model parametrization is evident through the presence of sharp changes in velocity with depth, but overall, we feel that the pattern of anomalies is robust, even if the nature of the transition between them (smooth or sharp) is not. With the period range of dispersion curves spanning 3s–12 s, the model should be well constrained between ~1 and 10 km depth; in order to confirm this we calculate the sensitivity kernels of a representative model shown in Fig. 8(b).

## 4 RESULTS AND DISCUSSION

The group velocity maps for periods of 3–12 s are shown in Fig. 7, and vertical slices through the 3-D shear wave velocity model are shown in Fig. 9. A 3-D rendering of the final $V_s$ model is also shown in Fig. 10. The two sets of results reveal strong group- and/or shear velocity contrasts and several coherent structures beneath the Merapi–Merbabu–Telemoyo Volcano lineament. One of the primary features of interest is a clear high shear velocity anomaly (body H1) with a NW–SE orientation between 5 and 10 km depth capped by a more superficial anomaly (body H2) with the same orientation (Figs 9a–c). This high shear velocity anomaly (also evident in the group velocity maps in Fig. 7) lies beneath the Merapi–Merbabu Complex and is surrounded by much lower velocities to the east and west. This anomaly appears to connect the two volcanoes at depth, and has a higher value beneath Merbabu, which is dormant. This result is consistent with a single igneous body underlying both volcanoes that is cooler (and hence faster) beneath Merbabu (cf. Widiyantoro et al. 2018). Another significant feature is the low shear velocity anomalies (mirrored to some extent in the group velocity maps) on the west and south flank of Merapi Volcano (bodies L1 and L2). In Figs 9(a), (b) and (d), body L1 appears to continue down to 10 km depth and could be an extension of the main magma chamber identified by Widiyantoro et al. (2018), which is also positioned south of Merapi, between 10 and 20 km depth. In Fig. 9(b), body L1 appears to be connected to the less prominent body L2. The position of this superficial low velocity body is consistent with the recent eruptive activity, which is located in the south region of the volcano. Analyzing seismic signals produced by phreatic explosive activity that occurred in 2014, Métaxian et al. (2020)
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 6.** Trade-off curves used to determine the optimal regularization parameter. Based on these curves we use values in the range 4–7 for the group velocity maps in the period range 3–12 s.</figcaption></figure>

was able to track the migration of microseismicity from 900 to 2700 m above sea level, which corresponds to the elevation of the crater floor. The microseismicity was interpreted to be the result of fluid influx from the shallow magma system into the hydrothermal

system; this caused a sudden decompression due to the transition of the volcanic fluid to a gaseous phase, and produced high-frequency microseismicity and the migration of a low-frequency tremor source towards the surface.
<figure><img src="image_19.png" /><figcaption>**Figure 7.** Rayleigh wave group velocity perturbation maps for selected periods. The velocity at all periods is relative to a group velocity of 1.11 km s⁻¹. The dashed line denotes the region of good resolution according to the checkerboard tests.</figcaption></figure>

Although there are differences in resolution, we compare our new results with those of Koulakov *et al.* (2007), Zulfakriza *et al.* (2014) and Koulakov *et al.* (2016) who exploit MERAMEX data. We also compare our results with geological analyses from Van Bemmelen (1949), geochronology and petrogenesis by Kohno *et al.* (2005), and other additional geophysical studies such as gravity by Tiede *et al.* (2005) and magnetotellurics by Müller & Haak (2004), Koulakov *et al.* (2007) performed body wave tomography using data from all 134 MERAMEX stations. The stations operated from May to October 2014 and covered an area of approximately 150 × 200 km in Central Java, Indonesia, with only about six stations in the Merapi–Merbabu Complex. Although their coverage of the Merapi–Merbabu Complex is limited, they noted that even though Merapi lies atop a large area of low velocities which correspond to the Merapi Lawu Anomaly (MLA), their shear wave velocity model reveals a deeper high velocity anomaly directly beneath Merapi and Merbabu volcanoes. The top of the anomaly appears at a depth of about 5 km and becomes weaker in the depth range 10–25 km. This appears to be approximately consistent (apart from the onset of high velocities at 5 km depth) with our new shear wave velocity model, which shows the presence of a large high velocity structure (bodies H1 and H2 in Fig. 9 and Fig. 10) beneath Merapi and Merbabu volcanoes which continues down to at least 10 km depth.

Zulfakriza *et al.* (2014) and Koulakov *et al.* (2016) both performed ANT in Central Java using MERAMEX data. Their results show a low-velocity anomaly located beneath Merapi and Lawu Volcanoes, which is also present in the body wave tomography results of Koulakov *et al.* (2007). Koulakov *et al.* (2016) then used their results to interpret the magma plumbing system of Central Java and the shallow structures present beneath the Merapi Complex; in particular, they image a clear low velocity anomaly beneath the summit of Merapi that extends to a depth of 8 km. In our case, the low velocity anomalies L1 and L2 (Fig. 9) occur further to the south and west. These differences are perhaps unsurprising, since resolution tests undertaken by Koulakov *et al.* (2016) indicate a horizontal resolution of the order of 30 km, whereas our peak horizontal resolution is closer to 5 km due to the DOMERAPI data, which have much denser path coverage.

Geochronology and petrogenesis research in the UTMM Volcanic Arc by Kohno *et al.* (2005) shows that there is a relationship between the Telomoyo, Merbabu, and Old Merapi Volcanoes in terms of their lava composition. Kohno *et al.* (2005) suggested that these three volcanoes are part of a volcanic arc that act as the trench side chain, while the Ungaran Volcano acts as the backarc chain. Thus, there is a distinct possibility that the Merbabu and Old Merapi Complex actually came from the same primary magma. The high shear velocity anomalies (bodies H1 and H2 in Fig. 9 and Fig. 10) discovered in this study also appears to reinforce this argument.

Tiede *et al.* (2005) produce a 3-D density model for the Merapi–Telomoyo region, based on the inversion of 443 gravity measurement points, which shows a strong correlation with our results. The density model reveals a NW–SE high-density structure ranging from 2 km above and 7 km below mean sea level. This high-density structure decreases strength to the southeast. This anomaly may originate from the older basaltic lava underlying the newer deposits, with the intrusive complex starting from the large, positive density anomaly body on the NW border and extending SE across the volcano. Furthermore, Tiede *et al.* (2005) also found a low-density body at depths of about 1–8 km to the southwest of Merapi Volcano which coincides with a region of high conductivity discovered by Müller & Haak (2004).

By integrating the results of this study with previous studies (Kohno *et al.* 2005; Tiede *et al.* 2005; Koulakov *et al.* 2007) we
<figure><img src="image_20.png" /><figcaption>**Figure 8.** (a) An example of a 1-D *V*<sub>S</sub> model ensemble and associated dispersion curves from the NA inversion for one gridpoint. The red colours indicate the best fit, while yellow, green and blue denote increasingly poor fits to the data (see panel on right). The black squares in the right-hand panel correspond to the group velocity map values at the selected gridpoint and the dashed black line in the left-hand panel is the selected best-fitting model. (b) The Rayleigh wave fundamental model group velocity sensitivity kernels for a representative model of the study area.</figcaption></figure>

can create a unified interpretation to explain the existence of the
high velocity anomaly (bodies H1 and H2 in Fig. 9 and Fig. 10).
The magma, which was formed by partial melting from subduction,
migrated through pre-existing transverse faults and produced the
Merbabu and Old Merapi Complex. Eventually, all of this magma

became inactive due to cooling and solidification, and the Merbabu
and Old Merapi Complex became an intrusive complex. This the-
ory is supported by the analysis of Merapi's evolution using strati-
graphic field data, radiocarbon ages, and whole-rock geochemical
data at locations such as Mt Bibi as a Proto-Merapi, and also Mt
<figure><img src="image_21.png" /><figcaption>**Figure 9.** Slices taken through the 3-D shear velocity model. Centre: horizontal slice taken at a depth of 6 km: (a) 110.41°E red dashed line across the western flanks of Merapi and Merbabu volcanoes; (b) 110.43°E green dashed line across Merapi and Merbabu volcanoes, (c) 7.45°S green line across Merbabu volcano; (d) 7.55°S red line across Merapi volcano. L1 and H1 indicate the low and high velocity anomalies ranging from ~4 km down to >10 km depth beneath Merapi and Merbabu, respectively. L2 is a near surface low velocity anomaly to the south of Merapi and H2 high is a high velocity anomaly, which underlies both Merapi and Merbabu at depth. These four anomalies are also seen in the 3-D representation shown in Fig. 10.</figcaption></figure>

Kendil and Mt Lawang (Gertisser 2012) as a part of Somma-Merapi, which lies above body H2, northeast of New Merapi. There is also Mt Turgo and Mt Plawangan on Merapi's southern flank that can be interpreted as either basaltic vents of Old Merapi or erosional remnants of a Proto-Merapi (Newhall et al. 2000). The new tomography results reveal a deep high velocity anomaly beneath Mt Turgo and Mt Plawangan as shown in Fig. 10.

While the active magma of New Merapi has a very different composition compared to that of Old Merapi, as discussed by Kohno et al. (2005), the location of the active magma may be constrained

in this study. The low velocity anomalies (bodies L1 and L2 in Fig. 9) correspond to both the low-density anomaly (Tiede et al. 2005) and the high conductivity body (1 Ωm) found by magnetotellurics (Müller & Haak 2004). Several alternative explanations for this feature are discussed by Müller & Haak (2004). The first suggestion is that the SW-conductor corresponds to pyroclastic and lahar deposits produced by New Merapi. However, the estimated production rate of about 10⁶ m³ yr⁻¹ over 1350 yr are too small to explain the volume of this anomaly (Müller & Haak 2004). The second alternative explanation is clay derived from accumulations of
<figure><img src="image_22.png" /><figcaption>**Figure 10.** 3-D image of the shear wave velocity perturbation model represented by isosurfaces taken at -22 and 16 per cent relative to an average shear velocity of 1.44 km s⁻¹.</figcaption></figure>

volcanic ash. However, a resistivity of 1 Ωm is too low for clay which usually has a resistivity of about 5–100 Ωm. Such a high conductivity body would require one million years at ambient temperatures to form (Risk *et al.* 1998). Therefore, we add another suggestion: that bodies L1 (4–10 km) and L2 (<3 km) are the possible location of magmatic intrusive bodies, which is consistent with the low shear wave velocities, low gravity anomaly, high conductivity and a highly active Merapi Volcano (Tiede *et al.* 2005; Widiyantoro *et al.* 2018).

## 5 CONCLUSIONS

In this study, we apply ambient noise tomography to data collected by the DOMERAPI passive seismic experiment, which consisted of 46 broad-band stations deployed around Merapi Volcano in Central Java. Rayleigh wave group velocity dispersion was extracted from the long-term cross-correlation of interstation waveforms, and inverted for period dependent group velocity maps between 3 and 12 s period. These maps were then inverted using the Neighbourhood Algorithm to produce a pseudo 3-D shear wave speed model, which is well resolved between the surface and 10 km depth (~5 km lateral resolution and 1–2 km vertical resolution). The new shear wave velocity model reveals a number of coherent anomalies that can be interpreted. This includes a large high velocity anomaly which underlies both Merapi and Merbabu at depth; we interpret this to be evidence of a massive igneous intrusion which used to be the upper crustal magma source for both Old Merapi and Merbabu Volcanoes. Furthermore, a distinct and pervasive low velocity anomaly located beneath the southwest flank of Merapi Volcano may represent a newly discovered magma source that supplies magma to New Merapi.

## ACKNOWLEDGEMENTS

The authors are grateful to the French National Research Agency (ANR) for their funding as well as to the Institut de Recherche pour le Développement for their support of the DOMERAPI project. Gratitude is also offered to the French Pool of Portable Seismological Instruments (SISMOB/RESIF) for providing the seismic stations and all of the DOMERAPI project partner organizations. We thank the two reviewers for their constructive comments on the original version of this manuscript. SW would like to thank RISTEK-BRIN for WCR funding, contract no. 127/SP2H/AMD/LT/DRPM/2020. This research is also partially funded by the Indonesian Ministry of Research and Technology/National Agency for Research and Innovation, and the Indonesian Ministry of Education and Culture under the World Class University (WCU) Program managed by Institut Teknologi Bandung.

## DATA AVAILABILITY

The DOMERAPI data used in this study are available at http://www.fdsn.org/networks/detail/YR_2013/with citation information provided by https://doi.org/10.15778/RESIF.YR2013.

## REFERENCES

* Barmin, M.P., Ritzwoller, M.H. & Levshin, A.L., 2001. A fast and reliable method for surface wave tomography, *Pure appl. Geophys.*, **158**, 1351–1375.

* Behm, M., Nakata, N. & Bokelmann, G.T.J., 2016. Regional ambient noise tomography in the eastern Alps of Europe, *Pure appl. Geophys.*, **173**, 2813–2840.
Benediktsdóttir, Á., Gudmundsson, Ó., Brandsdóttir, B. & Tryggvason, A., 2017. Ambient noise tomography of Eyjafjallajökull volcano, Iceland, *J. Volc. Geotherm. Res.*, **347**, 250–263.

Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *Geophys J. Int.*, **169**, 1239–1260.

Crowder, E., Rawlinson, N., Pilia, S., Cornwell, D.G. & Reading, A.M., 2019. Transdimensional ambient noise tomography of Bass Strait, south-east Australia, reveals the sedimentary basin and deep crustal structure beneath a failed continental rift, *Geophys. J. Int.*, **217**(2), 970–987.

Dziewonski, A., Bloch, S. & Landisman, M., 1969. A Technique for the analysis of transient seismic signals, *Bull. seism. Soc. Am.*, **59**(1), 427–444.

Gertisser, R., 2012. The geology evolution of Merapi Volcano, Central Java, Indonesia. *Bull. Volcanol.*, **74**, 1213–1233.

Hannemann, K., Papazachos, C., Ohrnberger, M., Savvaidis, A., An-thymidis, M. & Lontsi, A.M., 2014. Three-dimensional shallow struc-ture from high-frequency ambient noise tomography: new results for the Mygdonia basin-Euroseistest area, northern Greece, *J. geophys. Res.*, **119**, 4979–4999.

Huang, Y., Ohkura, T., Kagiyama, T., Yoshikawa, S. & Inoue, H., 2018. Shallow volcanic reservoirs and pathways beneath Aso caldera revealed using ambient seismic noise tomography, *Earth Planets Space*, **70**, 169, doi:10.1186/s40623-018-0941-2.

Kohno, Y., Setijadji, L.D., Pecskay, Z., Harijoko, A., Utami, P., Imai, A. & Watanabe, K., 2006. Geochronology and petrogenetic aspects of qua-ternary across arc magmatism on Merapi-Merbabu- Telomoyo-Ungaran Volcanoes, Central Java, Indonesia, in *Proceedings of the 3rd Int. Symp. Earth Resources and Geological Engineering Education*, Gadjah Mada University, Yogyakarta, Indonesia.

Kohno, Y., Setijadji1, L.D., Utami, P., Harijoko, A., Pecskay, Z., Imai, A. & Watanabe, K., 2005. Geochronology and Petrogenetic Aspects of Merapi- Merbabu-Telomoyo-Ungaran Volcanoes, Central Java, Indonesia, in *Proceedings of the Joint Conv. HAGI-IAGI-PERHAPI Annual Conference and Exhibition*, Surabaya, Indonesia.

Koulakov, I., Bohm, M., Asch, G. & Luhr, B.G., 2007. P and S velocity structure of the crust and the upper mantle beneath Central Java from local tomography inversion, *J. geophys. Res.*, **112**, B08310.

Koulakov, I., Maksotova, G., Jaxybulatov, K., Kasatkina, E., Shapiro, N.M., Luehr, B.-G., El Khrepy, S. & Al-Arifi, N., 2016. Structure of magma reservoirs beneath Merapi and surrounding volcanic centers of Cen- tral Java modeled from ambient noise tomography, *Geochem. Geophys. Geosyst.*, **17**, 4195–4211.

Métaxian, J.-P., et al., 2020. Migration of seismic activity associated with phreatic eruption at Merapi volcano, Indonesia, *J. Volcanol. Geotherm. Res.*, 396:106795, DOI: 10.1016/j.jvolgeores.2020.106795, *J. Volcanol. Geotherm. Res.*, **396**, 106795, doi: 10.1016/j.jvolgeores.2020.106795.

Müller, A. & Haak, V., 2004. 3-D Modelling of a conductivity structure: integrating magnetotellurics, induction vectors and the effect of steep topography—a case study from Merapi Volcano (Central Java), *J. Volc. Geotherm. Res.*, **138**(3–4), 205–222.

Neuman van Padang, M., 1951, *Catalogue of the Active Volcanoes of Indonesia*, Catalogue of the Active Volcanoes of the World Series, International Volcanological Association, 271 pp.

Newhall, C.G., Bronto, S. & Alloway, B., 2000. 10,000 years of explosive eruptions of Merapi Volcano, Central Java: archaeological and modern implications, *J. Volc. Geotherm. Res.*, **100**, 9–50.

Pedersen, H. A. & Kruger, F., 2007. Influence of the seismic noise charac-teristics on noise correlations in the Baltic shield. *Geophys. J. Int.*, **168**, 197–210.

Ramdhan, M. *et al.*, 2019. Detailed seismic imaging of Merapi volcano, Indonesia, from local earthquake travel-time tomography, *J. Asian Earth Sci.*, **177**, 134–145.

Risk, G.F., Bibby, H.M. & Caldwell, T.G., 1998. Resistivity structure of the central Taupo volcanic zone, New Zealand, *J. Volc. Geotherm. Res.*, **90**, 163–181.

Ritzwoller, M.H. & Levshin, A.L., 1998. Eurasian surface wave tomography: group velocitites, *J. geophys. Res.*, **103**, 4839–4878.

Sager, K., Boehm, C., Ermert, L., Krischer, L. & Fichtner, A., 2020. Global-scale full-waveform ambient noise inversion, *J. geophys. Res.*, **125**, e2019JB018644, doi:10.1029/2019JB018644

Cambridge, M., 1999. Geophysical inversion with a Neighbourhood Algorithm - I. Searching a parameter space, *Geophys. J. Int.*, **138**, 479–494.

Stankiewicz, J., Ryberg, T., Haberland, C., Fauzi, F. & Natawidjaja, D.H., 2010. Lake Toba volcano magma chamber imaged by ambient seismic noise tomography, *Geophys. Res. Lett.*, **37**, L17306, doi:10.1029/2010GL044211.

Tatsumi, Y. & Eggins, S., 1995. Subduction zone magmatism, in *Frontiers in Earth Sciences*, 1st edn, Blackwell Science, 211pp.

Tiede, C., Camacho, A.G. & Gerstenecker, C., 2005. Modelling the density at Merapi volcano area, Indonesia, via the inverse gravimetric problem, *Geochem. Geophys. Geosyst.*, **6**(9), 679, 1–13.

Van Bemmelen, R.W., 1949. *The Geology of Indonesia*, Government Printing Office.

Wathelet, M., Jongmans, D. & Ohrnberger, M. (2004). Surface-wave inversion using a direct search algorithm and its application to ambient vibration measurements. *Near Surf. Geophys.*, 2. doi:10.3997/1873-0604.2004018.

Widiyoranto, S. *et al.*, 2018. Seismic imaging and petrology explain highly explosive eruptions of Merapi Volcano, Indonesia, *Sci. Rep.*, **8**(1), 13656.

Yang, Y., Ritzwoller, M.H., Levshin, A.L. & Shapiro, N.M., 2007. Ambient noise Rayleigh wave tomography across Europe, *Geophys. J. Int.*, **168**, 259–274.

Yoshizawa, K. & Kennett, B.L.N., 2002. Non-linear waveform inversion for surface waves with a neighbourhood algorithm—application to multimode dispersion measurements, *Geophys. J. Int.*, **149**(1), 118–133.

Young, M.K., Rawlinson, N., Arroucau, P., Reading, A. & Tkalčić, H., 2011. High-frequency ambient noise tomography of southeast Australia: new constraints on Tasmania’s tectonic past, *Geophys. Res. Lett.*, **38**, L13313, doi:10.1029/2011GL047971.

Young, M.K., Cayley, R.A., McLean, M.A., Rawlinson, N., Arroucau, P. & Salmon, M., 2013a. Crustal structure of the east Gondwana margin in southeast Australia revealed by transdimensional ambient seismic noise tomography, *Geophys. Res. Lett.*, **40**, 4266–4271.

Young, M.K., Rawlinson, N. & Bodin, T., 2013b. Transdimensional inversion of ambient seismic noise for 3D shear velocity structure of the Tasmanian crust, *Geophys.*, **78**, WB49-WB62.

Yudistira, T., Paulssen, H. & Trampert, J., 2017. The crustal structure beneath the Netherlands derived from ambient seismic noise, *Tectonophysics*, **721**, 361–371.

Yudistira, T. & Widiyoranto, S., 2016. Characteristics of seismic noise in Central Java, Indonesia, *AIP Conf. Proc.*, **1730**, 020013, doi:10.1063/1.4947381.

Zhang, X., Hansteen, F., Curtis, A. & de Ridder, S., 2019. 1-D, 2-D, and 3-D Monte Carlo ambient noise tomography using a dense passive seismic array installed on the North Sea Seabed, *J. geophys. Res.*, **125**(2), doi:10.1029/2019JB018552.

Zulfakriza, Z., Saygin, E., Cummins, P.R., Widiyoranto, S. & Nugraha, A.D., 2014. Upper crustal structure of Central Java, Indonesia, from transdimensional seismic ambient noise tomography, *Geophys. J. Int.*, **197**(1), 630–635.

SUPPORTING INFORMATION

Supplementary data are available at GJI online.

VIDEO-2020-12-15-10-54-23 (1).mp4

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.