Preliminary result for crustal properties derivation related to tectonics for hazard mitigation in Eastern Indonesia using Teleseismic P Coda

Eastern Indonesia is tectonically complex, formed by different plates and microplates interactions from different origins. This complexity gives geoscientists a challenge to solve the 'jigsaw' of the complex interactions. The understanding of tectonic processes can lead to a breakthrough in both resource exploration and disaster risk reduction. We utilize teleseismic P wave coda for random coda from scattering and deterministic coda originated from the crust-mantle boundary (Moho) to derive the crustal properties, including thickness, Vp/Vs, and qualitative scattering characteristics. For the scattering properties, we apply Iterative Cross-Correlation and Stacking (ICCS) to align the waveform. At the same time, for the crust characteristic, we employ the Receiver Functions (RF) method alongside H-k stacking. The crustal thickness recovered from the RF and H-k stacking has a good correlation with the crustal origin, where the thickness in older and stable crust originated from Sundaland and Gondwana is thicker than a younger plate of the crust arc and subduction origin. The Vp/Vs is high in a region that is interpreted to be dominated by mafic lower crust originated from oceanic-oceanic subduction during Eocene, anisotropy, or by a magmatic anomaly. The P coda also correlated well with the subsurface magmatic anomaly by providing a unique pattern.


Introduction
Eastern Indonesia exhibits very complex geological and tectonic settings. This region is formed by the complex interaction between Sundaland as the eastern margin of Eurasian Plate, north moving convergence of Australian Plate, and west moving convergence of Pacific Plate [1]. This complexity is reflected by the forming of K-shaped Sulawesi and Halmahera Islands, drifting, and colliding of Sula, Buton, and other micro-continents westward, trapping and forming of the deep Banda Arc and generating major strike-slip structure of Sorong and Palu Koro [2,3,1] (Fig. 1A). The northern part features double subduction generating active volcano arc of Sangihe in the west and Halmahera in the east [1,4].
The complexities bring a challenge to understand both the crustal physical characteristics, its origin, and its geodynamics. The understanding of the properties and origin can lead to a breakthrough in resource exploration (hydrocarbon exploration, mineral, geothermal) that are dependent highly on the tectonic framework; hazard can also be minimized by knowing geodynamics and possible process of magmatism beneath the crust. Several studies with different methods have been comprehensively described the tectonic evolution of this region [3,5]. Complex plate reconstruction using computer modelling can explain * Corresponding author: ws@ugm.ac.id and model the interactions between Gondwana's microplates fragment and Sundaland from Jurassic-Cenozoic and correlate very well with the recent tectonic condition [3]. Seismic studies also provide a detailed view both in the spatial and temporal interaction of plates by using seismic tomography [5]. Seismic tomography being the bridge between mantle dynamics and crustal motions by evaluating each depth section and reconstruct by back projecting the plate motions and mantle structure variation. This method can also describe the velocity variation that is very related to the intrinsic physical characteristic of tectonic plates. Despite its powerful result, this method is very resourceextensive in the computational domain.
A simple seismic approach such as receiver function and teleseismic coda wave analysis can be the fastest and cheapest way in describing the crustal properties [6][7][8][9]. By using the phase conversion phenomenon in the velocity contrast between crust and upper mantle (Moho), we can derive and estimate the crustal thickness [8], which is important information in describing crustal and tectonic processes, including thickening, thinning, and its mantle implication. Teleseismic receiver function using vertical waveform to estimate the source waveform and remove it using deconvolution process [6,8,9]. Using the H-k stacking method, we further can estimate the bulk Vp/Vs of the crust and describe the intrinsic properties.
The coda wave also holds information on the scattering process in the subsurface; scattering is related to the distribution of small-scale heterogeneities such as magma packets as isolated strong velocity contrast [10].
Several studies have analyzed the Eastern Indonesia tectonic and crustal properties using the Receiver Function Method [11,12], yet those previous methods have not provided any analysis about the random seismic coda and the heterogeneities relationship. The relationship of the derived crustal properties with the origin of the Plate or tectonic block also has not been discussed. This study will relate the derived crustal properties and coda pattern both for deterministic and random coda with the crustal origin and its heterogeneities characteristic.

Methodology
The P coda processing will be separated into two independent processes: scattering teleseismic P coda analysis and teleseismic P Coda Receiver Function. For the first process, we collect 7 Mw > 7.0 events from WebDC3 GFZ [13] (Fig. 1B) with epicentre to station distance between 50-80 degrees to obtain nearly vertical incident of P wave propagating into the lithosphere, to ensure that the coda does not contain any other deterministic body wave phases 25 s after P arrival, the depth is constrained to be greater than 110 [14]. The sampling rate is uniform in all of the data (40 Hz), but the recording instrument varies.
Before the alignment process, the instrument response is removed to get the true ground motion using Obspy [15]. The corrected data then normalized to each maximum value, and 0.1-8.0 Hz bandpass filtered [14]. We apply the Iterative Cross-correlation and Stacking (ICCS) algorithm [16] to remove the time delay between the station recording. This algorithm ensures consistency: multitrace cross-correlation by stacking all the traces and iteratively adjusting the time picks. The algorithm is started by predicting initial time picks (arrival time) of P according to the existing velocity model IASP91 [17]. The traces are then aligned and stacked using the initial time picks. The second step is cross-correlating and shifting each trace with the stacked trace from the previous step. The process continues until the iteration does not change the stack.    The P coda for receiver function was collected for the same station with different events criterion: distance between 30-90 degrees and magnitude greater than Mw 6.0; the data is collected from 2010-01-01 until 2021-05-28 [13] (Fig. 1B). Receiver function and H-k stack calculated using RfPy [9]. The seismogram first rotated to vertical, radial, and transversal (ZRT). The receiver function utilizes the vertical component seismogram to approximate the source waveform and will be removed from the horizontal waveform using deconvolution [7][8][9]. The general deconvolution process will include the seasonal ambient noise spectrum of microseism [20,8]; this noise is reduced using the Wiener filtering approach [8]. The recovered radial and transversal waveforms stacked by back azimuth with spacing 10 degrees (Fig.  2) [6,9] (Fig. 2A). The recovered radial waveform from the deconvolution process is convolved again with the vertical waveform and cross-correlated with the observed radial waveform; the correlation coefficient quantifies the process quality [9]. For the data quality, the 0.05-1 Hz bandpass filtered vertical waveform signal-to-noise ratio is calculated in the 30s before and after the IASP91 P arrival [17,9].
Crustal thickness can be estimated directly by measuring the time difference between the first arrival P wave and its phase converted Ps. However, this simple method is very sensitive with Vp/Vs in the crust where error 0.1 s/km in Vp/Vs can lead to 4 km uncertainty in the estimated crustal thickness [19]. One way to reduce this dependency is by stacking multiple phases: Ps, PpPs (or Pps), and PsPs (or Pss) in the H and k domain by weighted sum [19]. [8] modified the stacking process by introducing the weighted product method in the amplitude corrected phases, where the negative amplitude is replaced by zero value. The maximum value obtained from stacking is the estimated Vp/Vs and crustal thickness (H).
The average Vp used is 6.0 km/s and stacked using 20 bins of slowness. The weighting combination is 0.5, 2.0, -1.0 for Ps, Pps, and Pps, respectively. The median value for every single-phase weighted stack is used to mitigate outlier that may distort the stacking process [9] (Fig. 2B). The k and H range is manually adjusted to get a reasonable value of k in crust [21] and H compared to the crustal thickness according to CRUST1.0: Global Model of Earth's Crust (http://igppweb.ucsd. edu/~gabi/crust1.html#reference).

Result and discussion
The crustal thickness of this study is in the range of 17.5-54 km, with BNDI (Banda Neira) being the area with the thinnest crust while BKB (Balikpapan) being the thickest (Fig. 4A) (Table 1). BKB is located northern of Meratus Suture zone that is resulted from the collision of East Java-West Sulawesi block to the South West Borneo block that has been amalgamated with Sundaland extended to northern in 60 Ma [18,19]. The collision might result in crustal thickening, causing the area to have the thickest crust among the study area. This station shows strong scattering in the P Coda wave that we interpret to be related with the complex melange of the ophiolite (Fig. 3A and 3B), high Vp/Vs (1.8) (Fig.  4B and Table 1) is likely to be attributed to ~10 km sediments of Kutai Basin [22] or might be the contribution of small mafic part of the concealed ophiolite complex.
TOLI2 (Toli-Toli) located in the bending of the northern arm of Sulawesi Island; this arm is interpreted to be a volcanic arc, ophiolite, and accreted formed in plate margin [3] contain mostly volcanic, volcaniclastics, and granitic stocks [1] and was an Early Eocene subduction complex [3]. This young formerly ophiolite and the volcanic arc have 29.5 km crustal thickness, far from BKB, whose block originated from much older Gondwana's block [23,3] and thickened by the accretion process. Though having a different value of crustal thickness, the Vp/Vs for this area is the same as BKB; we interpret that this high value of Vp/Vs is of mafic lower crustal contribution which forming island arc or accreted during the Early Eocene oceanic-oceanic subduction [3]. The P coda pattern is similar to BKB but showing higher fluctuation in the central part of the coda (window 585-595 s) in Fig. 3A and 3B. This difference might be attributed to more complex small-scale heterogeneities in TOLI2, either the granitic stocks or scattered magma pockets by the recent subduction zone. The surrounding LUWI also presents a past subduction zone and is the time marker of the early contact between the Australian Plate and Eurasian Plate [3] in 25 Ma. Sula Spur starts the contact as the fragment of the Australian Plate that collides with the volcanic arc of the northern Sulawesi arm. This contact formed subduction that accreted ophiolite and melange complex and formed the eastern arm [3,1]. The crustal thickness is higher than its past generation; we suggest that this slight difference is mainly drifted by the difference of subducting plate origin and rate of subduction, wherein the Sula Spur acts this case. The Vp/Vs is also higher than TOLI2, which likely is caused by the anisotropy [24] in this massively deformed region. SANI (Sanana) is part of Sula Spur that is composed mainly of granite and exposed as continental basement type [1], which might be the main contributor to its very low Vp/Vs (1.60) [21]. The crustal thickness is 37.0 km and is related to its Gondwana origin [3] and should be higher than TOLI2 and LUWI, which are much younger and are volcanic arc and ophiolite origin. The P coda shows no prominent and stable characteristic related to the observations in the other stations.
TNTI (Ternate), positioned in the Island of Maluku, has a very similar geological condition with Sulawesi, having volcanic in the western origin and ophiolite in the eastern origin [1]. TNTI lies near the currently active Gamalama Volcano; observation in this station obtains a very strange value of Vp/Vs 2.02; this anomaly can be interpreted as the existence of magma, based on its location. This high Vp/Vs can also be related to its origin that formed by oceanic-oceanic subduction and by the existence of the Molucca Sea Plate that subducts beneath this Plate. However, the latter could not be as the recovered crustal thickness is only 27 km. From the P coda pattern, magmatic activity likely to contribute to the prominent feature in the late time window (600-605 s) that is not observed in the non-volcanic related seismic stations (Fig. 3A). The late P coda prominent pattern is also observed in BNDI (Banda Neira), which is also located near a volcanically active region. Though its origin is still tried to be determined by a quantitative approach, the phenomenon is unique and spatially specific.
FAKI (Fak-Fak) being the typical characteristic of continental crust with 1.7 Vp/Vs and 35 km of crustal thickness. Tectonically, this area is part of Papua that is originated from Gondwana [3]. The coda shows a similar pattern with the other stations, but at certain time windows exhibit a weaker fluctuation.

Conclusion
The crustal thickness variation in the study area is very related to the origin of the block or plate where the observations are made. The crustal thickness for the older block originated from Gondwana (FAKI and SANI) is relatively thick, the BKB which should give the same result but thicken by the contact and amalgamation of Gondwana's microplates to the Sundaland, LUWI, TOLI2, and TNTI has thinner crust because of its young age and volcanic arc and ophiolite origin, while BNDI, the youngest arc having the thinnest crust. While the pattern Vp/Vs is harder to predict, we suggest that it is correlated with the crustal origin, anisotropy, or magmatic/fluid existence. The scattering properties observed by the P Coda can give qualitative insight into the subsurface heterogeneities; from the observation, the volcanic area relates to a unique pattern of the coda.
The crustal thickness and Vp/Vs value derived from receiver function analysis of deterministic the part of P Coda correlates well with the crustal origin and estimated subsurface condition. The random part of the coda related to the subsurface's scattering properties shows prominent and unique features that are likely to be caused by magmatic origin. Crustal characteristics will provide important information on the past and ongoing geodynamic processes, the main actor that trigger earthquake and volcanic hazard.
The plot for P Coda analysis is generated using Matplotlib; the coastline base map is obtained from Data Basin, the volcano distribution is from LAPAN, bathimetry and topography map from GEBCO plotted using QGIS, graphic annotation and layouts done using Inkscape. We thank Wakhidatun Nisa for providing digitized tectonic features in the area. This research is funded by Penelitian Dasar Unggulan Perguruan Tinggi (PDUPT) 2021 granted to Wiwit Suryanto.