Early Results of Time Domain Receiver Function Data Processing in Mt Merapi and Mt Merbabu

The receiver function method is a method to image the earth subsurface by utilizing Ps conversion waves that are formed due to the velocity boundary. In general, the receiver function method estimates depth of structures such as the mantle-crust boundary by deconvoluting the vertical component from the horizontal component. Typical receiver function data processing is done in the frequency domain where the deconvolution process can be seen as a division between two components. In this study, we tried to reprocess the data using a deconvolution technique in time domain, popularly known as iterative time-domain deconvolution. The principle of iterative time domain deconvolution consists of iterative cross-correlation between the horizontal and vertical component. We use data from the DOMERAPI seismic station network located in the vicinity of Mt Merapi and Mt Merbabu. Mt Merapi is one of the most active volcanoes in the world with frequent eruptions and located at the ring of fire chain volcano in Indonesia. Note that the previous receiver function study in this region showed complex signals at some stations that may be related to sediment at shallow sediment and possible layers of low velocity zone that interfering main signal for a crust-mantle boundary. Our current results show iterative time domain RFs have clearer and smoother signal than the frequency domain that help interpreting the waveform signals. We estimate a range of crust thickness between 26-31 km near Mt Merapi. However, we noticed that iterative time domain calculation requires longer computation time and input signal.


Introduction
Central Java's location is in the middle part of Java island, Indonesia. This area has an active tectonic activity where the Indian-Australian plate subducts the southern part of the Sundaland. The subduction of Indo-Australian plate beneath Sundaland began after SE Asia evolution in the middle Miocene about 42 MA [1] and existed to date with a convergence rate of 67 mm per year. Central Java is characterized by many volcanic activities and arc magmatism due to the subduction setting's influence. There are two main volcanic arcs in this area: The Southern Mountain Arc (SMA) and Modern Volcanic Arc (MVA). Southern Mountain Arc (SMA) formed in the island's southern coastal region started at the middle Eocene about 45-20 Ma [2]. The subduction in this area ceased in the Cretaceous but then resumed in the late Miocene and formed Modern Volcanic Arc (MVA), which located about 50 km north of SMA [3]. On the edge of Sundaland, a large basin formed in the north of MVA, named as the Kendeng Zone. In this study, the focus area is Mt. Merapi and Mt. Merbabu which are part of MVA.
Mt. Merapi is one of the most active volcanoes in the world and has an eruption frequency between 2 to 6 years [4]. The type of Mt. Merapi is andesitic basaltic volcanoes with a crater containing a lava dome and has pyroclastic flows [5]. [4] studied seismicity discovered a non-seismic zone within the volcano cone and suggested two magma reservoirs at shallow (1-2 km) and deeper depth (5 km). This hypothesis was supported by the petrological composition study [6].
Interestingly, Mt. Merbabu, located 9 km to the North of Mt. Merapi shows lower activities compared to Mt. Merapi. In this study, we tried to reprocess data from DOMERAPI network with receiver function analysis. Receiver function analysis is a method to image the Earth's distinct layer such as Crust-Mantle boundary using 3-components waveform information from far-distant events. [7] showed complex RFs signals in this area by processing the data in a frequency domain. Alternatively, the deconvolution method in RF can also be done in time domain that could enhance the signal. The objective of this study is to reprocess the same data set using time domain deconvolution to better RFs signal. We compare the results between the frequency domain and time domain process to investigate the benefit between each method and finally help the interpretation.

Data and Method
We used recorded waveform event from DOMERAPI temporary seismic network carried out from October 2013 to April 2015. The network consisted of 53 three component seismometers with 4 km station spacing spread out near Mt. Merapi and Mt. Merbabu ( Figure 1). We selected recordings of teleseismic earthquake with magnitude higher 6 with epicentral distance between 30-90° and gather more than 130 events. In general, after removing some events that has lower to signal ratio, we collected 50 good RFs per station.

Receiver Function
Receiver function technique is a method to generate P-S converted waves from teleseismic data. Ps converted wave formed when the P wave pass through velocity boundaries or discontinuities underneath seismic station. To process the data in RFs, we rotated vertical, North-South and East-West components into vertical, radial and tangential components using back azimuth information between source to events. To improve the Ps signals, further rotation can be done into P, SV and SH component using incidence angle information estimated using ak135 velocity model. In RFs, to isolate Ps waves, the radial component is deconvolved from the vertical component or SV components is deconvolved from the P component.
The deconvolution method can be done either in time or frequency domain. In this study we used the deconvolution technique called iterative time domain deconvolution. The basic of this method is crosscorrelation process between radial and vertical component to analyze similarity from both components. Before the deconvolution, first the data will be filtered to reduce the noise. After that, the filtered vertical component is cross-correlated with the radial component to estimate the lag of the first and largest spike in receiver function [8]. Zero time is the expected P wave time arrival.
From the cross-correlation result will be obtained the first iteration spike which is a sign of similarity between the signal of the vertical component and the radial component. Then the spike will be filtered with Gaussian width factors that will control the bandwidth of the signal to be estimated ( Figure 2). Figure 2 display that typically the first iteration in deconvolution will create the first spike at zero second, and followed by later signals. The result of the filter process will be convoluted with vertical components to produce the prediction of the first iteration radial component. Figure 5 show the difference between radial component estimated from RFs and observation data. Depending of the complexity of the signals, Most RFs show good results after 50 iterations. The whole process will be repeated to estimate spike lag and other amplitudes through iterations until they get results with small errors.

Iteration 150
The  Figure 3. Radial observation compared to radial prediction. Horizontal axis is time at 0.01 s per tick.
Zero time represents the arrival of P wave phase.
In this study we also performed H-k stacking method [9]. H-k stacking approach is grid search technique from Zhu and Kanamori to estimate the thickness of crust as well as Vp /Vs ratio below a station. Theoretical delay times of Ps, PpPs and PpSs phases are calculated for each pair of H, k according to a given mean P-wave velocity above the discontinuity of interest. The best estimation of the depth to the discontinuity and the Vp/Vs ratio is expected to coincide with the maximum of the amplitude stack.

Result and Discussion
Each receiver function from event was stacked to image subsurface beneath the station. The RF stacked result makes determining the arrival times of p and s waves easier to analyse. H-k stacked was applied to final receiver function stack to observe crustal thickness and Vp/Vs ratio around study area.     H-k contour stack from two methods showed by figure 5b (time domain) and 5c (frequency domain). RF time domain method estimate crustal thickness beneath ME29 station is 28.85 km with 1.76 Vp/Vs ratio. In other hands the RF frequency domain method estimate crustal thickness at 28.56 km with 1.77 Vp/Vs ratio. Crustal thickness generated from both methods have fairly similar result.

Conclusion
As the conclusion, we proposed that the receiver function time deconvolution method can be used as alternative method for deconvolution data processing since it can generate similar results with frequency domain deconvolution and show simpler signals with less wiggles in the later time. However, still visual inspection is recommended at each RFs to enhance overall data quality.