A preliminary results: crustal structure in celebes region based on receiver function analysis

14004


Introduction
The Indonesian archipelago is located at the convergence of three major tectonic plates that are highly active, namely the Eurasian Plate, the Pacific Plate, and the Indo-Australian Plate, as well as the micro Philippine Plate.[1].Therefore, Indonesia's territory is highly vulnerable to natural disasters such as tectonic earthquakes and tsunamis, both on land and at sea.Celebes Island is located at the convergence of three tectonic plate movements: the Indian Plate, the Australian Plate, and the Eurasian Plate.Celebes Island is geologically vulnerable, with frequent movements and collisions of tectonic plates resulting in geological phenomena such as earthquakes and tsunamis.In the northern part of Eastern Indonesia, the Pacific Plate collides with the northern part of Papua Island, the northern part of the Maluku Islands, the northern part of Celebes Island, and the Maluku Sea region.Celebes Island and its surroundings, especially North Celebes, are one of the most complex active margins in terms of geology, structure, and tectonics.This region is the convergence center of three converging crustal plates due to the interaction of three major crustal plates during the Neogene era [2].
The Celebes Region, both on land and at sea, is close to sources of tsunamis caused by earthquakes and tectonic processes.There are active faults in the Celebes region that cause to geological phenomena.The earthquake sources in the Celebes Sea usually caused by the subduction of the North Celebes Plate, the Mayu Trench Plate, and the Sangihe Plate located to the east of North Celebes.Previous studies using geodetic, geological, and seismological data have shown that the Indonesian tectonic region can be divided into several smaller plates, namely the Burma Plate, Sunda Plate, Bandameri Plate, Maluku Sea Plate, Timor Plate, Bird's Head Plate, Maoke Plate, and Woodlark Plate [3].
Based on historical earthquake data from the Incorporated Research Institutions for Seismology (IRIS) and the Meteorology, Climatology, and Geophysics Agency (BMKG) over the past 5 years (2017-2022), there were 12,391 earthquake events in the research area.From 1821 to 2017, 31 damaging earthquakes and 2 significant tsunami-triggering earthquakes were recorded (Earthquake and Tsunami Center, 2018).According to the seismicity map in (Figure 1), the research area (Celebes) experiences a high level of seismic activity Due to the high level of seismic and volcanic activity, it is essential to study the structure of the Earth's crust in Celebes region.The method used for this research is the receiver function method.In this study, the investigation of the Earth's crust will include the models of P-wave and Swave velocities, Vp/Vs ratios, and the depth of the Moho Discontinuity beneath receiver stations in the Celebes region.This research utilizes 8 permanent seismic stations owned by BMKG located in different tectonic and geological zones.It is expected that this research can provide insights into the models of P-wave and S-wave velocities, Vp/Vs ratios, and the depth of the Moho Discontinuity in the Celebes region, which can improve the accuracy of earthquake hypocenter determination.Research on the crustal structure, P and S waves, and Vp/Vs in the Celebes region has been conducted by several researchers before [4].Research on the Earth's crust in the Eastern Indonesian region, specifically at the TMSI station in North Celebes, the variation in crustal thickness is estimated to be around 19 km in the northern part of Celebes.The research employs the 1-D velocity wave structure model method in the Minahasa region, which concludes that the 1-D velocity wave structure model in the Minahasa region ranges from (2.95 -9.07) km/s at a depth of 36 km with an RMS residual of 3.7 [5].
A local tomography study was conducted in the Palu region of Central Celebes using LOTOS-12 (Lotos Tomography Software-12) to analyze the data of P-wave and S-wave velocity distribution anomalies in the Palu region of Central Celebes.The research results indicate that the Central Celebes region has several layers in the upper crust (0-20 km) with a Vp (P-wave velocity) value of 5.961 km/s and a Vs (S-wave velocity) value of 3.589 km/s.In the lower crust (20-40 km), Vp >6.423 km/s and Vs 3.816 km/s were obtained, and in the upper mantle (> 40 km), Vp values were around 6.423 km/s, and Vs values were 3.816 km/s.The minimum Vp/Vs ratio anomaly value was found to be approximately 1.6, while the maximum Vp/Vs ratio anomaly value was around 1.888 [6].A 3D wave velocity modeling study using the double-difference method and travel time tomography for the identification of the Palu-Koro fault using a 1D P-wave velocity model that yields velocity values ranging from 3.28 km/s to 9.07 km/s in the depth range of 0 -36 km [7].Receiver function studies at the crust-mantle boundary (Moho) are being conducted to determine crustal properties, including thickness and Vp/Vs ratio.At the TOLI-TOLI station, the crustal thickness is 29.5 km with a Vp/Vs ratio of 1.8, and at the LUWI station, the crust thickness is 32.5 km with a Vp/Vs ratio of 1.88 [8].In this study, receiver function analysis was performed using M>6 teleseismic earthquakes with an epicenter distance of 30°-90° recorded at 8 broadband seismic stations belonging to the Meteorology, Climatology and Geophysics Agency was carried out to observe crustal structure, P, S and Vp/Vs wave velocities at the regional level.The stations are located in 4 main geological areas in Celebes.

Receiver function
The teleseismic receiver function method is a technique used to depict the properties of the Earth's structure using teleseismic seismic data recorded at 3-component receiver stations.This method relies on the phase changes of seismic waves that occur as seismic waves pass through layer of different impedance, including the conversion of P-waves to S-waves (Ps) generated when P-waves traverse discontinuity layers beneath the receiver station.By eliminating source effects and instrument response, seismic recordings are obtained that contain only information about the medium through which the waves travel from the source to the receiver.The diagram illustrates the scheme of seismic wave phase shift (Figure 2).

Receiver function rotation
Seismic data recorded on a seismogram generally use the ZNE coordinate system, which stands for vertical (Z-axis), North-South horizontal (North-South axis), and East-West horizontal (East-West axis).To isolate the energy and make the wave phases more clearly visible, seismogram rotation is performed.This rotation utilizes the back azimuth (γ) of the arriving wave.The vertical component (Z) remains unchanged, the N-S component is rotated into the radial (R) component, and the E-W component is rotated into the tangential (T) component.This coordinate system is known as the ZRT coordinate system.The rotation process is carried out using the equations below: From the equation above, it can be noted that the vertical component (Z) remains constant, while the horizontal component (NE) is rotated into radial components (R) parallel to the direction of wave arrival, and tangential components (T) perpendicular to the direction of wave arrival.Here, (γ) represents the back-azimuth of the arrival direction of teleseismic waves, measured from the 0° direction (north) from the observing station [10].

Iterative time domain deconvolution
The seismic data received at the three-component station consists of earthquake source function, local station structure response, and instrument response [10].The function of the ZRT coordinate components on a seismogram in the time domain can be expressed using the equation: The seismogram components Z(t), R(t), T(t) represent the vertical component, radial component, and tangential component, respectively.These are the results of the convolution (*) calculation involving the instrument response I(t), seismic source function S(t), and the functions EZ(t), ER(t), ET(t), which are the structural responses in the vertical, radial, and tangential components.
Iterative Time Domain Deconvolution, or deconvolution using an iterative time-domain approach, is based on minimizing the least-squares difference between the components of observed signals and synthetic signals.This is achieved by iteratively performing a forward modeling process of convolution using a wavelet shape model for the vertical Z(t) component, as well as the radial and tangential components.Additionally, a signal model with minimal misfit is utilized as a response to the local structure, comparing the convolution results of synthetic signals to the observed signals [9,11,12].
Ri(t) and Ti(t) being seismograms of the radial and tangential components resulting from the iterative timedomain convolution model, EiR(t) and EiT(t) represent the iterative models of the receiver function for the radial and tangential components.The iteration EiR(t) is equal to the observation ER(t) if the misfit value between the observation R(t) and the synthetic R_i(t) resulting from the convolution between EiR(t) and Z(t) is minimized.To eliminate high-frequency noise from the receiver function and isolate teleseismic signals at low frequencies, a filter is used in the iterative calculations, namely a low-pass Gaussian filter with the equation: With "α" as the control parameter for the bandwidth of the Gaussian filter for high-frequency noise filtering and "ω" as the gain of the Gaussian filter, where at ω=0 there is unity (unit-area impulse).

Inversion
Inverse modeling, also known as inversion modeling, is a technique used to adjust a model so that it matches observed data.The process involves modifying specific parameters within the model until the calculated data aligns with the observed data.This approach is useful for analyzing complex systems and improving the accuracy of model predictions.However, it should be noted that inverse modeling is not always necessary or appropriate, and objective evaluation is recommended before deciding to use this technique [13].The inversion process starts with evaluating the initial model and then conducting forward modeling to determine the optimal fit between observed data and the model to obtain a velocity model.The initial velocity model used is the AK135 model [14].The formula for modeling the velocity of a receiver function is expressed as follows [15].In the equation as follows: Specifically, the above non-linear equation is inverted by linearizing it, and the iterative approach is expressed as follows:

Stacking H-k
The H-κ stacking method, developed by [16], aims to determine the crustal thickness and calculate the average value of ��/�� from the receiver function time series data.The mathematical formulation for estimating crustal thickness is as follows : With � being the number of receiver functions, ∑ �� represents the weighting factor subject to the condition ∑ � = 1, and (�) is the amplitude of the receiver function at �1, �2, and �3, where t, respectively, is the time difference of the arrival times of the Ps, PpPs, and PpSs+PsPs phases relative to the direct P-wave.

Data
The research data consists of teleseismic earthquake waveform with a magnitude greater than 6 (M>6) and an epicenter distance to the receiving station in the range of 30°-90°.These recordings, available in miniSEED format, were obtained from 8 broadband three-component permanent receiving stations owned by the Meteorology, Climatology, and Geophysics Agency (BMKG), The earthquake catalog link at https://geof.bmkg.go.id/webdc3/ , provided by BMKG, allows for downloading of waveform and instrument response data.The selection of stations used in this study was based on their corresponding zones.The Ensialic Volcanic Arc Zone was represented by TMSI and SMSI stations, the Ensimatic Volcanic Arc Zone was represented by KRSI and MKS stations, the Metamorphics Zone was represented by WKCM and SRSI stations, and the Ophiolites Zone was represented by BBSI and LUWI stations (refer to Figure 3).

Methods
The teleseismic waveforms are obtained by downloading data from the BMKG network, with the epicenter distance between 30°-90° from the receiving station and the magnitude greater than 6 (M>6), using miniSEED data format.Furthermore, the instrument response data from each receiving station is also obtained.The data is downloaded during the time window of 5 minutes before and 20 minutes after the P-wave arrival.
The downloaded data is then rotated from ZNE coordinates (vertical, north-south, east-west) to ZRT coordinates (vertical, radial, and tangential) by applying the back-azimuth based on the receiving station coordinates and the earthquake epicenter coordinates.After this, the ZRT rotated data is deconvolved to remove the earthquake source effects, which result in the Earth's structure response function.Additionally, Gaussian filtering with a width of α=1.5 is used in this process to eliminate noise with frequencies above 0.75 Hz.The receiver function calculation is performed on the signal with the smallest misfit, within a time window of 5 seconds before and 20 seconds after the P-wave.
The deconvolved results are migrated into the depth domain and stacked to enhance and strengthen the receiver function signals while removing background noise.The migration and stacking results will reveal seismic phases representing the boundaries of Earth's crustal structures.The migration is performed using the velocity model resulted from the inversion outcome and the AK-135 velocity model.The velocity model determination is accomplished through inversion based on the processed data.The determination of P-wave and S-wave velocity models are determined from the initial AK-135 velocity model and the receiver function model of the radial component.Inversion is then carried out based on these two models to obtain the local velocity profiles beneath the station.
The Earth's crust thickness (H) and the Vp/Vs ratio are estimated using the stacking H-κ method.The Splitrflab program runs this technique in Matlab program developed by [18].The radial receiver function data, which is the output of the CPS330 program, serves as input data in this program.Subsequently, the wave velocity and slowness values or wave parameters are determined.The function (H, κ) used in the study is defined for 20 to 40 kilometers with a 0.01 km interval and for κ values ranging from 1.5 to 2.7 using a 0.01 grid.The next step involves computing time differences for the arrival phases , , and + relative to the phase.The time differences for each phase are then stacked with weighting factors of 0.7, 0.2, and 0.1 for each respective time [16].The maximum value of s (H, κ) is then indicated as the location of the Moho discontinuity boundary.

Results and discussion
The chapter's results depict the waveform phase of the receiver function was plotted based on the back azimuth of the earthquake.The profile of P-wave and S-wave velocities was then determined by inverting the receiver function and finding the depth of the Moho discontinuity.The models for P-wave and S-wave velocity, as well as the Vp/Vs ratio, were determined using the H-κ stacking method.The research results are based on the classification of the characteristics of the basement rocks in the Celebes region, including the Ensialic Volcanic Arc Zone, Ensimatic Volcanic Arc Zone, Ophiolites Zone, and Metamorphics Zone [17].

Ensialic Volcanic Arc Zone
The Ensialic Volcanic Arc Zone is represented by 2 three-component broadband seismic stations from BMKG.The earthquake event data used covers the period from March 2019 to March 2023.The first station, Tondano Timur, Minahasa Regency (TMSI), is located at coordinates 1.29° S 124.93° E, and the second station, Sumalata, North Gorontalo Regency (SMSI), is located at coordinates 0.99° S 122.37° E. (Figure 5 Phase at the TMSI station appears to not align precisely at second 0, but there is a delay of approximately 1 second.According to [19], The occurrence of a delay in the arrival time of the phase is suspected to be due to the presence of sediment layers beneath the station, creating a strong impedance contrast beneath the low-velocity sediment layer.Meanwhile, the conversion of the wave into the wave or wave is observed at around ~6 seconds, and the PpPs phase is observed at around ~9 seconds.The direct arrival phase of the P wave ( ) at the SMSI station is seen at 0 seconds.Subsequently, the conversion phase of the P wave into the S wave (Ps) from the Moho layer is observed at approximately 6 seconds.Furthermore, the observed PpPs phase is relatively small at 10 seconds.In Fig. 5, a strong negative phase is also observed at 2 seconds and remains consistent in each quadrant.According to [20] Interpreting the low-velocity zone.The low-velocity zone of this negative phase is associated with the presence of partial melting at Mount Mahawu and Mount Lokon-Empung and active subduction between the Maluku Sea Plate and the Halmahera Plate beneath the Sangihe Arc [21].
The observed receiver function signals, represented by the blue line, have a match of 58.30% with the inversion results represented by the red line can be seen in (Figure 6).The profile of P-wave velocity at TMSI station can be seen in.The P-wave velocity has a minimum value of 4.1 km/s and a maximum value of 12.1 km/s.The estimated crust thickness is approximately 35-40 km, indicated by an increase in P-wave velocity from 7.2 km/s to 9.7 km/s can be seen in (Figure 7).
The S-wave velocity at TMSI station can be seen in (Figure 8).It has a minimum value of 2.3 km/s and a maximum value of 7 km/s.The crust thickness is estimated to be around 35-40 km, marked by an increase in S-wave velocity from 4.07 km/s to 5.4 km/s.The inversion results of the receiver function yield a matching value of 95.10% between the observed signal and the inverted signal at SMSI station (Figure 9).
The P-wave velocity profile at the SMSI station can be seen in (Figure 10).The P-wave velocity has a minimum value of 4.4 km/s and a maximum value of 9.9 km/s.The estimated crustal thickness is approximately 35-40 km, indicated by the increase in P-wave velocity from 5.9 km/s to 6.3 km/s.
The S-wave velocity at the SMSI station is shown in (Figure 11).It has a minimum value of 2.4 km/s and a maximum value of 4.3 km/s.The crustal thickness is estimated to be around 35-40 km, marked by the increase in S-wave velocity from 3.3 km/s to 3.5 km/s.Lastly, to obtain depth values and Vp/Vs ratios, the H-k Stacking method was employed using the Matlab application.The analysis results indicate that the depth of the Mohorovičić Discontinuity layer at TMSI Station is approximately 30.89 km, with a Vp/Vs ratio of 2.07 km/s (Figure 13).At SMSI Station, the Moho layer was found at a depth of around 34.17 km, with a Vp/Vs ratio of 1.99 km/s (Figure 14).

Ensialic Volcanic Arc Zone
The Ensimatic Volcanic Arc Zone is represented by The results of the radial component teleseismic receiver function signals from various back azimuths and the stacking results were performed to clarify the phases of the receiver function, namely the direct P-wave arrival (Pp) phase at KRSI station, observed at 0 seconds.Subsequently, the conversion phase from P-wave to S-wave (Ps) from the Moho layer is observed around 5 seconds.The direct P-wave (P) and the P-wave converted to S-wave (Ps) phases are also observed.The Pp phase at MKS station appears to be slightly delayed, not precisely at 0 seconds, but with approximately a 1-second delay.According to [19], The occurrence of delay in the arrival time of the phase is suspected to be due to the presence of sediment layers beneath the station, creating a strong impedance contrast beneath the low-velocity sediment layer.Meanwhile, the conversion of the wave into the or wave is observed at approximately 5 seconds.The inversion results of the receiver function yield a matching value of 96.46% between the observed signal and the inverted signal at station KRSI (Figure 16).
The profile of P-wave velocity at Station KRSI can be seen in (Figure 17).The P-wave velocity has a minimum value of 6.9 km/s and a maximum value of 11 km/s.The estimated crustal thickness is approximately 38-45 km, indicated by the increase in P-wave velocity from 7.8 km/s to 9.1 km/s.
The S-wave velocity at Station KRSI can be observed in (Figure 18).It has a minimum value of 3.8 km/s and a maximum value of 6.3 km/s.The crustal thickness is estimated to be around 38-45 km, marked by the increase in S-wave velocity from 4.3 km/s to 5.3 km/s.The inversion results of the receiver function yield a velocity model for P and S waves at the MKS station with a matching value of 89.39% between the observed and inverted signals (Figure 19).
The velocity profile for P waves at the MKS station can be seen in (Figure 20).The minimum P-wave velocity is approximately 3.5 km/s, and the maximum is 9 km/s.The estimated crust thickness is around 30-35 km, indicated by the increase in P-wave velocity from 5.8 km/s to 6.3 km/s.
The S-wave velocity at the MKS station can be observed in (Figure 21), with a minimum value of 3.5 km/s and a maximum of 4.5 km/s.The crust thickness is estimated to be approximately 30-35 km, as indicated by the increase in S-wave velocity from 3.2 km/s to 3.5 km/s.Finally, to obtain depth values and Vp/Vs ratios, the H-k Stacking method was used with the Matlab application.The analysis results indicate that the depth of the Mohorovičić Discontinuity layer at KRSI Station is approximately 39.21 km, with a Vp/Vs velocity ratio of 1.75 km/s (Figure 23).At MKS Station, the Moho layer was found at a depth of approximately 33.44 km, with a Vp/Vs velocity ratio of 1.87 km/s (Figure 24).

Metamorphics Zone
The crustal structure in the Metamorphic Zone is represented by 2 three-component broadband seismic stations from BMKG.The earthquake event data used spans from March 2019 to March 2023 and includes Station Langgomali in Wolo Subdistrict, Kolaka Regency, Southeast Celebes (WKCM) located at coordinates 3.83° S, 121.30°E, and an elevation of 116 meters, as well as Station Cendana Hitam in Tomoni Timur Subdistrict, East Luwu Regency, South Celebes (SRSI) situated at coordinates 2.53° S, 120.88°E, and an elevation of 36 meters.In (Figure 25(a)), 27 events of the radial component teleseismic receiver function are displayed, along with stacking results at WKCM station with a fit value ≥ 90%, and in (Figure 25(b)), 46 events of the radial component teleseismic receiver function are shown, along with stacking results at SRSI station with a fit value ≥ 90%.Both stations exhibit the same initial signal pattern.
The results of the radial component receiver function signal from various back azimuths and the stacking results are conducted to clarify the phases of the receiver function, namely the direct P-wave arrival (Pp) observed at the WKCM station at 0 seconds.Subsequently, the phase conversion of the P-wave to the S-wave (Ps) from the Moho layer is observed at around 4 seconds.The direct P-wave (P) phase and the P-wave converted to Swave (Ps) are evident.The Pp phase at the SRSI station is observed at 0 seconds, while the phase conversion of the P-wave to the S-wave (Ps) is observed at approximately 4 seconds.The inversion results of the receiver function yield a matching score of 98.96 % between the observed signal and the inverted signal at station WKCM (Figure 26).
The P-wave velocity profile at station WKCM is shown in (Figure 27).The P-wave velocity has a minimum value of 4.5 km/s and a maximum value of 10 km/s.The estimated crustal thickness is approximately 35-38 km, indicated by an increase in Pwave velocity from 7.1 km/s to 7.9 km/s.
The S-wave velocity at station WKCM can be seen in (Figure 28).It has a minimum value of 2.6 km/s and a maximum value of 5 km/s.The estimated crustal thickness is around 28-34 km, marked by an increase in S-wave velocity from 3.9 km/s to 4.4 km/s.The result of the receiver function inversion is a model of P and S wave velocities at the SRSI station, which shows a matching value of 96.34% between the observed signal and the inverted signal (Figure 29).
The P-wave velocity profile at the SRSI station can be seen in (Figure 30).The P-wave velocity has a minimum value of 5 km/s and a maximum value of 9 km/s.The estimated crustal thickness is approximately 25-30 km, indicated by an increase in P-wave velocity from 6 km/s to 6.7 km/s.
The S-wave velocity at the SRSI station can be observed in (Figure 31).It has a minimum value of 2.8 km/s and a maximum value of 4 km/s.The crustal thickness is estimated to be around 28-34 km, marked by an increase in S-wave velocity from 3.9 km/s to 4.4 km/s.Finally, to obtain depth and Vp/Vs velocity values, the H-k stacking method was employed using the Matlab application.The analysis results indicate that the depth of the Mohorovicic Discontinuity layer at WKCM Station is approximately 34.38 km, with a Vp/Vs velocity ratio of 1.75 km/s (Figure 33).At SRSI Station, the Moho layer was found at a depth of approximately 21.01 km, with a Vp/Vs velocity ratio of 2.07 km/s (Figure 34).The results of the receiver function signal from various back azimuths and the stacking process were performed to clarify the phases of the receiver function, namely the direct P-wave arrival ( ) phase at station LUWI, observed at 0 seconds.Furthermore, the conversion phase of P-wave to S-wave (Ps) from the Moho layer is observed at around 5 seconds.The phase at station BBSI is observed at 0 seconds, while the conversion of the P-wave into S-wave ( ) is observed at approximately 5 seconds.The results of the receiver function inversion yielded a model of P and S wave velocities at the LUWI station with a matching accuracy of 90.13% between the observed signals and the inverted signals (Fig. 36).
The P-wave velocity profile at the LUWI station can be seen in (Figure 37).The P-wave velocity has a minimum value of 4.2 km/s and a maximum value of 10 km/s.The estimated crustal thickness is approximately 35-38 km, indicated by an increase in Pwave velocity from 8.3 km/s to 9.3 km/s.
The S-wave velocity at the LUWI station is shown in (Figure 38).It has a minimum value of 2.5 km/s and a maximum value of 5 km/s.The crustal thickness is estimated to be around 35-38 km, as evidenced by the increase in S-wave velocity from 4.6 km/s to 5.2 km/s.The inversion results of the receiver function yield a matching value of 90.39% between the observed signal and the inverted signal at station BBSI (Figure 39).
The P-wave velocity profile at Station BBSI can be seen in (Figure 40).The P-wave velocity ranges from a minimum of 4 km/s to a maximum of 9 km/s.The estimated crustal thickness is approximately 35-38 km, as indicated by the increase in P-wave velocity from 6.7 km/s to 7.8 km/s.
The S-wave velocity at station BBSI can be observed in (Figure 41).It has a minimum value of 2 km/s and a maximum of 4.5 km/s.The crustal thickness is estimated to be around 35-38 km, marked by the increase in S-wave velocity from 3.7 km/s to 4.3 km/s.To obtain depth and Vp/Vs velocity ratios, the H-k stacking method was employed using the Matlab application.The analysis results indicate that the depth of the Mohorovičić Discontinuity layer at the LUWI station is approximately 35.67 km, with a Vp/Vs velocity ratio of 2.51 km/s (Figure 43).At the BBSI station, the Moho layer was found at a depth of approximately 22.07 km, with a Vp/Vs velocity ratio of 2.25 km/s (Figure 44).A summary of research findings on the structure of the Earth's crust, including Moho depth, P and S wave velocity, vp/vs ratio, and crustal thickness map at 8 BMKG seismic stations on Celebes Island can be displayed in Table 1 and Figure 45.We would like to thank BMKG for providing the data and UP2KM STMKG for funding the publication of this research.

Fig. 2 .
Fig. 2. Phases of the teleseismic receiver function (left), and a scheme depicting the phase changes of P and S waves when encountering different mediums [9].

Fig. 3 .
Fig. 3.The station distribution map used is based on the main geological zones in Celebes [17].
(a)) displays 41 events showing the radial component of teleseismic receiver functions from various back azimuths and the stacking results at TMSI station with a fit value ≥ 90%.(Figure 5(b)) shows 42 teleseismic receiver function radial components from various back azimuths and the stacking results at the SMSI station with a fit value ≥ 90%.

Fig. 5 .
Fig. 5.The results of the receiver function plotted based on back azimuth (left) and stacking results (right) at TMSI station (a) and SMSI station (b).

Fig. 6 .
Fig. 6.Comparison of inversion results (red) and observations (blue) of the stacking signal receiver function at TMSI Station.

Fig. 9 .
Fig. 9. Comparison of inversion results (red) and observations (blue) of the stacking signal receiver function at SMSI Station.

Fig. 11 .
Fig. 11.S-wave velocity profile at SMSI Station.To obtain depth information using migration, the AK 135 velocity model is employed to migrate the receiver function from the time domain to the depth domain.The migration results using the AK-135 velocity model indicate that the depth of the Mohorovicic Discontinuity layer beneath Station TMSI is approximately 40-45 km, with a negative polarity phase (Low Velocity Zone) at a depth of around 18 to 22 km (Figure 12(a)).Meanwhile, at Station SMSI, the Mohorovicic Discontinuity layer is located at a depth of approximately 35-42 km, with a negative polarity phase (Low Velocity Zone) at a depth of about 14 to 20 km (Figure 12(b)).

Fig. 12 .
Fig. 12.The results of the migration of receiver functions to depth at TMSI station (a) and SMSI station (b) located in the Ensialic Volcanic Arc Zone.

Fig. 13 .
Fig. 13.The results of the Moho depth and Vp/Vs velocity using H-k stacking at TMSI station.

Fig. 14 .
Fig. 14.The results of the Moho depth and Vp/Vs velocity using H-k stacking at SMSI station.
2 three-component broadband seismic stations from BMKG.The earthquake event data used spans from March 2019 to March 2023, specifically from the Lembah Hopo Station in Karossa District, Mamuju Regency, West Celebes (KRSI) located at coordinates 1.82° S, 119.42°E, with an elevation of 60.3 meters, and the Tamarunang Station in Gowa Regency, South Celebes (MKS) situated at coordinates 5.22° S, 119.47°E, with an elevation of 80.32 meters.In (Figure 15(a)), 27 teleseismic receiver function events with radial components from various back azimuths are displayed, along with the stacking results for station KRSI with fit values ≥ 90%.(Figure 15(b)) shows 39 teleseismic receiver function events with radial components from various back azimuths and the stacking results for station MKS with fit values ≥ 90%.The initial signals displayed at both stations have different patterns.

Fig. 15 .
Fig. 15.The results of the receiver function plotted based on back azimuth (left) and stacking results (right) at KRSI station (a) and MKS station (b).

Fig. 16 .
Fig. 16.Comparison of inversion results (red) and observations (blue) of the stacking signal receiver function at KRSI Station.

Fig. 19 .
Fig. 19.Comparison of inversion results (red) and observations (blue) of the stacking signal receiver function at MKS Station.

Fig. 21 .
Fig. 21.S-wave velocity profile at MKS Station To obtain depth information using Migration, the AK-135 velocity model is employed to migrate receiver functions from the time domain to the depth domain.The results of migration using the AK-135 velocity model indicate that the depth of the Mohorovičić Discontinuity layer beneath Station KRSI is approximately 39-42 km, and there is a negative polarity phase (Low Velocity Zone) at a depth of about 5-22 km (Figure 22(a)).Meanwhile, at Station MKS, the Mohorovičić Discontinuity layer is located at a depth of around 35-40 km, with a negative polarity phase (Low Velocity Zone) at a depth of approximately 20-35 km (Figure 22(b)).

Fig. 22 .
Fig. 22.The results of the migration of receiver function to the depth at KRSI station (a) and MKS station (b) located in the Ensialic Volcanic Arc Zone.

Fig. 23 .
Fig. 23.The results of the Moho depth and Vp/Vs velocity using H-k stacking at KRSI station.

Fig. 24 .
Fig. 24.The results of the Moho depth and Vp/Vs velocity using H-k stacking at MKS station.

Fig. 25 .
Fig. 25.The results of the receiver function plotted based on back azimuth (left) and stacking results (right) at WKCM station (a) and SRSI station (b).

Fig. 26 .
Fig. 26.Comparison of inversion results (red) and observations (blue) of the stacking signal receiver function at WKCM Station.

29.
Comparison of inversion results (red) and observations (blue) of the stacking signal receiver function at SRSI Station.

Fig. 31 .
Fig. 31.S-wave velocity profile at SRSI Station.Furthermore, to obtain depth information using Migration, the AK-135 velocity model is used to migrate receiver functions from the time domain to the depth domain.The migration results using the AK-135 velocity model indicate that the depth of the Mohorovicic Discontinuity layer beneath Station WKCM is approximately 32-37 km, and there is a negative polarity phase (Low Velocity Zone) at a depth of around 18-22 km (Figure 32(a)).Meanwhile, at Station SRSI, the Mohorovicic Discontinuity layer is located at a depth of about 27-32 km, with a negative polarity phase (Low Velocity Zone) at a depth of around 10-15 km (Figure 32(b)).

Fig. 32 .
Fig. 32.The results of the migration of receiver functions to depth at WKCM station (a) and SRSI station (b), which are located in the Metamorphic Zone.

Fig. 33 .
Fig. 33.The results of the Moho depth and Vp/Vs velocity using H-k stacking at WKCM station.

Fig. 34 .
Fig. 34.The results of the Moho depth and Vp/Vs velocity using H-k stacking at SRSI station.

4. 4
Ophiolites Zone The crustal structure in the Ophiolite Zone is represented by 2 three-component broadband seismic stations from BMKG.The earthquake event data used for analysis covers the period from March 2019 to March 2023.These stations are Stasiun Bubung, Banggai Regency, Central Celebes (LUWI), located at coordinates 1.04° S 122.77° E and an elevation of 6 meters, and Stasiun Katobengke, Bau-Bau City, Southeast Celebes (BBSI), located at coordinates 5.49° S 122.56° E and an elevation of 95 meters.(Figure 35(a)) displays 30 teleseismic receiver function events with radial component data from various back azimuths and the stacking results at station LUWI with a fit value ≥ 90%.(Figure 35(b)) shows 37 teleseismic receiver function events with radial component data from various back azimuths and the stacking results at station BBSI with a fit value ≥ 90%.The initial signals displayed at both stations exhibit the same pattern.

Fig. 35 .
Fig. 35.The results of the receiver function plotted based on back azimuth (left) and stacking results (right) at LUWI station (a) and BBSI station (b).

Fig. 36 .
Fig. 36.Comparison of inversion results (red) and observations (blue) of the stacking signal receiver function at LUWI Station.

Fig. 39 .
Fig. 39.Comparison of inversion results (red) and observations (blue) of the stacking signal receiver function at BBSI Station.

Fig. 41 .
Fig. 41.S-wave velocity profile at BBSI Station.To obtain depth using Migration, the AK-135 velocity model is employed to migrate the receiver function from the time domain to the depth domain.The migration results using the AK-135 velocity model indicate that the depth of the Mohorovičić Discontinuity layer beneath the LUWI Station is approximately 35-40 km, and there is a negative polarity phase (Low Velocity Zone) at a depth of around 15-30 km (Figure 42(a)).Meanwhile, at the BBSI Station, the Mohorovičić Discontinuity layer is situated at a depth of approximately 34-39 km, with a negative polarity phase (Low Velocity Zone) at a depth of around 18-24 km (Figure 42(b)).

Fig. 42 .
Fig. 42.The results of the migration of receiver functions to depth at LUWI station (a) and BBSI station (b), which are located in the Ophiolite Zone.

Fig. 43 .
Fig. 43.The results of the Moho depth and Vp/Vs velocity using H-k stacking at LUWI station.

Fig. 44 .
Fig. 44.The results of the Moho depth and Vp/Vs velocity using H-k stacking at BBSI station.

Table 1 .
Summary of crustal structure results from receiver function analysis at 8 BMKG seismic stations in Celebes.