{'fraction_non_alphanumeric': 0.032512461870396546, 'fraction_numerical': 0.02819730674800982, 'mean_word_length': 4.605504587155964, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Crustal Structure Study in Lampung Region Using Teleseismic Receiver Function Analysis

Lampung region is seismically and volcanic active because located in subduction zone of Indo-Australian and Eurasian plate. We applied receiver function and stacking H-k analysis to estimate the crustal structure in Lampung region. We used teleseismic earthquake data (epicenter distance 30°-90°) and M>6 recorded at 3 seismic broadband stations owned by Agency for Meteorology Climatology and Geophysics (BMKG). Those stations are PSLI (located on Sebesi Island approximately 20 km from Anak Krakatau) represented volcanic arc zone, KASI (located on Kota Agung, Lampung) represented Sumatran Fault Zone and KLI (located on Kotabumi, Lampung) represented back-arc basin. Crustal thickness estimated at PSLI station 32-36 km, KASI station 36-40 km, and KLI station 30-36 km. Furthermore, in 3 stations P wave velocity estimated 4.1-11 km/s, S wave velocity 2.2-6.2 km/s, while vp/vs value estimated 1.7-2.05. We estimated Anak Krakatau volcano’s magma chamber beneath PSLI station in depth 16-30 km, Great Sumatran Fault structure in depth about 8-14 km beneath KASI station, and thick sediment layer about 4 km near surface beneath KLI station. This study result is expected to explain more detail crustal of Lampung region and can be useful for developing of BMKG’s seismic monitoring systems and other geophysical fields in future.


Introduction
Lampung region located in the southern part of Sumatra Island has complex tectonic. This makes Lampung region seismically active that we can see on Lampung seismicity map at Figure 1. This condition is caused by Lampung located in subduction zone of Indo-Australian plate to Eurasian plate that move with velocity 7 cm/year [1]. Earthquake source in Lampung region divided in two, those come from subduction zone and from Great Sumatran Fault (GSF). There are 3 segments of Great Sumatran Fault crossed Lampung region, those are Sunda segment, Semangko segment, and Kumering segment [2]. This Kumering segment was recorded to caused earthquake in Liwa, West Lampung on 1994 with M6.6 and hipocenter depth 23 km. This earthquake caused 217 people died, more than 2,000 injured, more than 6,000 house, shops and buildings were permanently damaged [3]. Lampung region not only has an active seismically, but also has a very active volcanic condition. This condition is caused by the presence of Krakatau volcano which now being Anak Krakatau volcano in Sunda Strait. Anak Krakatau Volcano was recorded to have erupted in December 22 nd , 2018 and caused volcano tsunami with maximum run up 13.5 m and inundation of 330 m [6].
With the high level of seismicity and volcanic activity in Lampung region as well as the history of destructive earthquake and volcanic eruption have occurred, so it is very important to study about crustal structure in Lampung region. This study aim is to explain more detail crustal structure of Lampung region using receiver function analysis. Previous study using receiver function method such as [7][8][9][10][11]. We expect this study result can be useful for developing the capability of BMKG's seismic monitoring systems and research other geophysical fields in the future.

Method
In this study, we used time series earthquake seismogram data recorded at 3 broadband stations owned by Meteorology Climatology and Geophysics Agency (BMKG) which can be accessed in webdc BMKG earthquake catalog [4]. Those stations are PSLI represented volcanic arc zone (located on Sebesi Island, Sunda Strait approximately 20 km to Anak Krakatau volcano), KASI represented Great Sumatran Fault zone, and KLI stations represented back-arc basin. Earthquake data criteria used are teleseismic earthquake events with epicenter distance of 30°-90° from broadband stations and earthquake magnitude more than 6 (M>6) with high signal to noise ratio (SNR) and clear P arrival phase.

Receiver function method
Receiver function is a method to estimate earth structure using teleseismic earthquake recorded at 3 components broadband receiver station. This method introduced by Langston [12]. This method utilizes the phase conversion of seismic waves that occurs when the seismic waves pass through different impedances. By eliminating the influence of source effect and instrument response, the seismogram remains information about medium of seismic wave in its propagation from source to receiver. The diagram of receiver function can be seen on Figure 4.  [13] After data preparation, the processing start with removed instrument response to eliminate the influence of seismograph amplifier to get the original signal from the earthquake event. Then we rotated the seismogram to change seismogram orientation from ZNE coordinate system (vertical, north-south, and east-west) to ZRT coordinate system (vertical, radial, and transverse) to clarify receiver function signal with following equation: Where α is back-azimuth, Z is vertical component, T is transverse component, and R is radial component.  [14] is applied to produce receiver function with time window 5 second before and 25 second after P arrival. To produce receiver function signal with low frequency, we applied low pass gaussian filter with width of 1,5 which able to separate frequencies higher than 0.75 Hz. Furthermore, to modelling 1D P and S wave velocity, an inversion processing is applied to the stacked receiver function with the initial model using AK135 [15]. The process of removed instrument response, seismogram rotation, iterative time domain deconvolution, and inversion used program CPS330 [16].

Stacking H-k
To estimate crustal thickness and vp/vs value, we also used stacking H-k method in Matlab program. This method is developed by [17]. This method estimates crustal thickness and vp/vs value by calculating difference time arrival of receiver function phase Ps, PpPs, PpSS+PsPs to arrival of direct P phase. Stacking H-k method is matematically formulated as following equation.
Where is number of receiver function stacked, ∑ weighting factor with total 1, in this study we applied weight 0,7, 0,2, 0,1, ( ) is receiver function amplitude, and 1 , 2 , and 3 is the difference time arrival of Ps, PpPs, and PpSs phase against direct P phase.

PSLI Station
PSLI station also known as Rajabasa station is located in Sebesi Island, Sunda strait approximately 20 km from Anak Krakatau volcano. We got total 13 traces of receiver function fit>85% from earthquake event data recorded from 2020 to April 2021. Receiver function radial component of PSLI station is shown in the Figure below.  Figure 5, we observed direct P phase delayed and arrived in 1s. This delay time according to [18] this delay time is caused by thick sediment layer near surface and it is confirmed by [19] that identified 5 km sediment near surface of Anak Krakatau volcano. We also observed the strong negative phase consistent from all back-azimuth and according to [20] this negative phase indicated the low velocity zone associated with Anak Krakatau volcano's magma chamber. Then, PpPs observed in 10s.
To see the layering structure under the PSLI station, then we reprocessing the receiver function with high resolution (gaussian filter width 5) and migrated into the depth domain as we can see in Figure 6. The results obtained is fit with the model from [19] Figure 6. The first positive phase after direct P wave arrival was observed at depth of 5 km that showed discontinuity between the sedimentary layer and granitic basement. Then at depth of 10 km, the next positive phase is observed that showed a discontinuity between granitic basement layer and diorite/gabbro layer. Furthermore, at depth 29-40 km, a positive Ps phase was observed which indicated a discontinuity between the diorite/gabbro layer and Moho discontinuity in the upper mantle (shown by yellow dotted line). A strong negative phase is also observed which estimated to be a low velocity zone from Anak Krakatau's magma chamber at depth about 20 km (shown by black dotted line).

Figure 6. Receiver Function of PSLI Station Migration (left), and Schematic Illustration of Magma Plumbing System of Anak Krakatau Volcano [19]
The result of P and S wave velocity modelling from inversion processing are shown in the following Figure. We observed that the P and S wave velocity beneath PSLI is fluctuated. P wave velocity estimated in interval 4.1 km/s-11 km/s, while S wave velocity in interval 2.3 km/s-6.2 km/s. On the surface to a depth of 4 km, a low S wave velocity 2.3 km/s-3.2 km/s is observed and it is associated with a sediment layer. The P wave velocity is observed decrease at depth of 14-30 km from 11 km/s to 6,6 km/s, while S velocity also observed decrease at depth 16-30 km from 6.2 km/s to 3.6 km/s. This depth interval estimated as the low velocity of Anak Krakatau magma chamber. This result also fit with [21] using the seismic tomography method to identified the magma chamber at depth about 10-20 km, while [22] with local tomography identified magma chamber about down to 4 km. Furthermore, the most significant P and S wave velocity increation was observed at depth about 32-36 km/s which estimated to be Moho discontinuity.  26.59 km using stacking H-k method is different with crustal thickness identification by velocity modelling and receiver function migration. According to [23], the stacking H-k method application in anisotropic zone is less accurate, so crustal thickness estimation in PSLI is about 32-36 km by contrast increation of velocity modelling. This depth interval is not much different with several previous study such as, [24] identified Moho depth about 22 km beneath Anak Krakatau volcano using seismic method, and [25] using seismic reflection in Sunda strait identified Moho depth at 28 km.

KASI Station
KASI station also known as Kota Agung station is located in Kota Agung, Tanggamus district. We got total 44 trace receiver function fit>90% from earthquake event 2016-2018. Receiver function radial component of KASI station is show in the Figure 9.  Figure 9, we observed direct P phase at 0 s, followed by a small negative phase at 2s, Ps phase which is the conversion phase of P wave into the S wave at the Moho discontinuity is observed at 4s, and PpPs phase at 8s. Figure 10 shows P and S velocity model of KASI station. We observed that P and S wave velocity fluctuated, P wave velocity is in interval 4.1 km/s-8 km/s, while S wave velocity is in interval 2.2-4.4 km/s. P wave velocity is observed decrease from 6 km/s while S velocity decrease 3.2 km/s to 2.8 km/s at depth of 8-14 km. This depth interval estimated as the structure of Great Sumatran Fault. According to [26] who identified fracture zone of fault has relatively small density, and [27]

KLI Station
KLI station also known as Kotabumi station is located in Kotabumi, North Lampung. We got total 47 trace receiver function fit>90% from earthquake event 2016-2017. Receiver function radial component of KASI station is shown in the Figure 12. From Figure 12, we observed direct P phase delay and arrived in 1s. according to [18], this delay time is caused by thick sediment layer. Then Ps phase is observed at 3s, while PpPs in 5s.  Figure 12. Radial receiver function of KLI station Figure 13 shows P and S wave velocity model by inversion processing of KLI station. P and S wave velocity is observed fluctuated, P wave velocity is in interval about 4.6-9 km/s while S wave velocity is in interval about 2.6-5 km/s. A low S wave velocity observed from surface to a depth of 4 km associated with sediment layer. A significant P and S wave velocity increase was observed at depth about 30-36 km estimated as Moho discontinuity. Crustal thickness estimation in KLI station is about 28.14 km and vp/vs value 2.05 using stacking H-k method can be seen at Figure 14. This depth result by stacking is different with crustal thickness estimation by velocity modelling. According to [23], the stacking H-k method application in anisotropic zone is less accurate, so crustal thickness estimation in KLI station is about 30-36 km by contrast increation of velocity modelling. This depth interval is not much different with several previous study such as [28] and [29] in Southern Sumatran Basin zone using receiver function. The crustal thickness in this zone is about 27-33 km.  Figure 14. Stacking H-k of KLI station

Conclusion
We estimated crustal thickness, P and S velocity modelling, and vp/vs value in 3 seismic broadband station in Lampung. Crustal thickness estimated at PSLI station 32-36 km, at KASI station 36-40 km, and at KLI station is 30-36 km. Then P wave velocity in 3 seismic stations estimated in interval 4.1-11 km/s, S wave velocity estimated 2.2-6.2 km/s, and vp/vs value estimated 1.7-2.05. We observed high vp/vs value 2.05 and decrease of P and S wave velocity in depth 16-30 km beneath PSLI station associated with presence of Anak Krakatau volcano's magma chamber. Then we observed the decreasion of seismic velocity at KASI station in depth about 8-14 km associated with Great Sumatran Fault structure. Furthermore, in KLI station we also observed high vp/vs value 2.05 and low velocity of S wave near surface to depth 4 km association as sediment layer. We expect this study result can be useful for developing the capability of BMKG's seismic monitoring systems and research other geophysical fields in the future.