Spatial variations of shear-wave velocity anomaly derived from Love wave ambient noise seismic tomography along Lembang Fault (West Java, Indonesia)

1 Introduction

The Lembang Fault, located north of the city of Bandung in West Java, Indonesia, is an active fault that poses a significant earthquake hazard in the region [1,2,3,4,5,6,7,8]. The Lembang Fault extends for 29 km in an east-west direction, and the eastern part of the fault is estimated to have been formed about 200,000 years ago, while the western segment was formed about 27,000 years ago [4,7]. Based on surface geological observations and paleoseismological studies by Daryono et al. [7], the fault is capable of generating earthquakes with a magnitude of 6.5–7.0. The Lembang Fault is reported to be active at a rate of 6 mm per year based on geodetic studies, as shown by Meilano et al. [9]. The Lembang Fault zone has experienced several small earthquakes with a maximum magnitude of 3.4 [10,11,12,13]. This has caused moderate damage in neighborhoods as a result of insufficient structural strength of buildings.

Several studies have been conducted on the Lembang Fault, specifically regarding its activity. The Lembang Fault is a left-lateral strike-slip fault with a minor vertical component, as observed through surface geological analysis [1,7]. While there are differing opinions, some suggest it may be a normal fault [6,14,15]. The Lembang Fault is a left-lateral strike-slip fault with a displacement rate ranging from 3 to 14 mm/year, as observed through global positioning system (GPS) measurements from 2006 to 2009 [16]. This movement was later confirmed by subsequent observations from 2006 to 2011, which showed an identical displacement rate of 6 mm/year [9]. Seismic observations indicate that the earthquakes occurring along the Lembang Fault are caused by strike-slip mechanisms [10,11,12,13], the Lembang Fault seismicity, observed over 6 km deep, suggests that the fault may be thrusting [10]. Moreover, geophysical observations, such as gravity data [17,18] and Controlled Source Audio Magnetotelluric (CSAMT) measurements [19], strongly support the claim that the Lembang Fault is a normal fault.

Currently, there is no subsurface model supported by high-resolution data to accurately describe subsurface conditions around the fault. Therefore, to adequately describe the geometry of the Lembang Fault, a subsurface model with high-resolution data is necessary. This study presents a seismic velocity model that characterizes the distribution of low-shear wave velocity areas associated with fault planes in the Lembang Fault zone. The Love wave ambient seismic noise tomography method was used to develop the model.

Ambient seismic noise tomography is an effective method to visualize subsurface arrangements using seismic waves [20,21,22,23,24,25]. This method has been shown to have great potential for subsurface exploration and structural imaging. In ambient seismic noise tomography, seismometers record seismic waves from natural sources, such as ocean waves, wind, or human activity. This technique uses cross-correlations of ambient noise recordings from different stations to generate Green’s function, which quantifies the propagation characteristics of seismic waves across the established network. By analyzing the variations in the Green’s function over the region of interest, the subsurface structure can be inferred. This inexpensive and non-invasive method provides relatively high-resolution data, allowing continuous subsurface monitoring with minimal environmental impact, thus proving a valuable technique for fault zone characterization [26,27,28,29,30,31,32,33,34,35,36].

Love waves are a type of surface wave that propagates along the Earth’s surface and are particularly useful for imaging the uppermost layers of the Earth’s crust [37,38,39,40,41]. By analyzing the propagation of Love waves in a fault zone, we can gain insights into the characteristics of the fault, including its geometry, location, and orientation. This information is critical for understanding the behavior of the fault and for assessing the potential for seismic hazards in the region. In general, the shear wave (S-wave) velocity resulting from a love wave inversion tends to be lower in the fault zone compared to the surrounding rock due to the presence of fractures and faults that create pathways for fluids and weaken the rock [42,43]. The S-wave velocity can also vary within the fault zone itself, depending on the degree of fracturing and faulting. The presence of fluids in the fault zone can also affect the S-wave velocity [38,42,43,44,45,46,47,48].

The relatively dense network of data acquisition temporary stations allowed us to generate S-wave velocity models around the Lembang Fault zone with high resolution using the Love wave ambient seismic noise tomography method. The velocity contrast associated with the fault plane is shown in the model and is confirmed by the small earthquakes which occurred in this area.

2 Method

Seismic data were collected from 74 stations between November 2020 and May 2021, utilizing 15 portable seismometers. Three of the 15 portable seismometers were deployed as fixed stations, maintained at the same location throughout all measurement periods, while the other 12 were used as mobile stations, changing locations during each measurement period. Layout station locations are shown in Figure 1.

Figure 1

Data recording station layout. The blue triangle is the location of the mobile station and the red box is the location of the fixed station. The big red triangle is the crater of Mount Tangkuban Parahu. The continuous black lines are traces of faults (Lembang Fault) on the surface referred by Daryono et al. [7].

Data processing for this study consisted of several steps, such as preparation of individual data, Green’s function estimation using cross-correlation and stacking, dispersion curve extraction, tomographic inversion using the Love wave group velocity map, and one-dimensional (1D) depth domain S-wave velocity inversion, according to Bensen et al. [27]. The data from each station were segmented into daily intervals and trend corrected. Amplitude spectrum analysis was used to determine the period range of useful data. A consistent instrument response observed from a period range of 0.04–10 s amplitude spectrum has a value above 20 dB/Hz and is consistent at each location of the data, so we selected this as the filter operator for all data.

We use the cross-correlation of each pair of stations as Green’s function to obtain coherent noise, as described by Bensen et al. [27]. Green’s functions of surface waves, which estimate the Earth’s response between two stations using noise data from two stations, have been extracted using cross-correlation, a common interferometry technique used in ambient seismic noise tomography [49].

Love wave data processing is performed in this study by extracting the transverse components. This is accomplished by rotating the data collected from the horizontal components (X and Y) after stacking each component, as described by Igel et al. [50]. Then, the extraction of the dispersion curve is performed on the transverse components. A signal-to-noise ratio (SNR) cutoff of 10 is applied to the normalized amplitude of each cross-correlation result. We applied a velocity filter with a range of 1.5–3.2 km/s to the cross-correlation results after stacking as a waveguide to localize the seismic signal of interest [51,52,53,54,55]. Rotation of radial and transverse components is shown in Figure 2.

Figure 2

The results of the rotation of the radial and transverse components. The left side is the result of cross correlation of the radial component and the right side is the transverse component. The dotted red line in the left and right images is the result of cross correlation at an interstation distance of 5.9933 km and the blue line is dashed at an interstation distance of 33.1228 km, which are further shown by their respective dispersion curves in Figure 3. The relatively vertical dashed black lines on the left and right represent the velocity filter area with a value of 1.5–3.2 km/s.

Using the frequency–time analysis (FTAN) method [27,37], the group velocity dispersion curve of the Love wave was obtained for each period. FTAN is a multi-filter process based on the desired bandwidth and selected filter parameters. In this study, the frequency range was evaluated from 0.1 to 25 Hz with a frequency step of 0.1 Hz and a frequency width of 0.8 Hz. Figure 3 shows a sample dispersion curve extracted for distances between stations ranging from 5.9933 to 33.1228 km.

Figure 3

(a) and (d) The results of cross correlation of transverse component at an interstation distance of 5.9933 km and the blue line is dashed at an interstation distance of 33.1228 km. (b) The dispersion curve of the results of cross correlation trace a, and (e) the dispersion curve of the results of cross correlation trace d. (c) and (f) Love wave group velocities, respectively, of curves a and d. The color scale shows the energy value of the dispersion curve.

The Love wave group velocity tomography process uses only the dispersion curves of data pairs with wavelengths greater than one lambda. The Love wave group velocity model was inverted using fast-marching surface tomography (FMST) as described by Rawlinson and Sambridge in 2005. In order to determine the optimal parameterization of the model, we evaluated several parameters (Figure 4) and used the checkerboard (CKB) test to determine the grid parameters for the inversion process. The results of the CKB test were evaluated qualitatively. In tomography, the grid size varies due to the different number of rays at each period (Figure 5).

Figure 4

Variety of parameters were tested to determine the best model parameterization. (a) The effect of the damping factor on the root mean square error, (b) the effect of the subspace dimension on the root mean square error, and (c) the effect of the smoothing factor on the root mean square error.

Figure 5

The CKB test was used to test grid parameters for the inversion process, the size of the grid in the tomography process is varied due to the difference in ray number in each period. In the period from 0.04 to less than 1 s, the optimal grid size is 2.5 km. For a period of 1 to less than 5 s, a grid of 2 km is used, while for a period of 5–10 s a grid size of 4 km is used. The continuous black lines are traces of faults on the surface referred by Daryono et al. [7]. The thin dotted black line marks the boundary area’s which have fairly good of recovery CKB test.

To obtain 1D depth-domain S-wave velocity values from dispersion curves, we used the surf96 program, a component of the Computer Program in Seismology package, as described by Herrmann [56]. Seismic velocity and density tend to increase gradually as a result of overburden pressure. Therefore, we used an initial velocity model for the depth inversion process in this study that gradually increases as a function of depth. The surface S-wave velocity model value is 2.0 km/s, referring to the published S-wave velocity model around the Lembang Fault by Afnimar et al. [11]. At a depth of 10 km, the S-wave velocity value is 3.46 km/s, referring to the AK135 1D seismic velocity model. The Vp/Vs ratio employed is 1.73, and the density is obtained from compressional velocity conversion using the Gardner equation [57]. Figure 6 illustrates the 1D inversion process.

Figure 6

The illustration of the 1D inversion depth domain process. (a) The 1D S-wave velocity curve resulting from the inversion from the initial model up to the thousandth iteration. The initial model is indicated by a dotted black line. (b) Sensitivity analysis of the period 0.2, 1, 5, and 10 s. (c) The thousandth inversion dispersion curve, the black point is the data, and the continuous red curve is the modeling of the inversion result. (d) The effect of the number of iterations on the average residual velocity dispersion.

3 Results

The results of the cross-correlation of the signals, especially for the radial and transverse components, are quite satisfactory when the velocity fan filter is applied. In general, the arrival time of the Love waves shown in the cross-correlation results of the transverse component is faster than that of the radial component or the Rayleigh waves. As shown in Figure 3, the group velocity of Love waves has a higher value in longer periods than in shorter periods at both closer (5.9933 km) and farther (33.1228 km) interstation distances.

The optimal parameterization of the model was selected by trial and error evaluation of the root mean square error value. The damping factor had an optimal value of 1, while the subspace dimension and smoothing factor were also set to 1. The CKB test is employed to evaluate the inversion process’s grid parameters. The variation in grid size in tomography is influenced by the number of rays in each period. An optimal grid size of 2.5 km is selected for a range between 0.04 and less than 1 s. For periods of 1 to less than 5 s, a 2 km grid is utilized, and for periods of 5–10 s, a 4 km grid is used.

One-dimensional depth inversion was carried out for 1,000 iterations for each point, resulting in an average residual velocity of 0.01 km/s. The results of the period sensitivity analysis in the 1D depth inversion process show that the highest vertical resolution is obtained at the lowest period of 0.04 s around 200 m, and the maximum period of 10 s has a low normalized sensitivity value at depths greater than 8 km, so the depth value considered most optimal is only up to a depth of 8 km, as shown in Figure 6.

Figure 7 shows the Love group velocity map, and Figure 8 shows the ray paths at each period. It is important to note that the Love wave group velocity value for a period of 1 s decreases compared to the velocity map of less than 1 s and then increases again for higher periods. Figure 9 displays the quantities of rays and root mean square error, and the average Love wave group velocity in each period. The highest number of rays, 692, is observed in a period of 1 s with a lambda constraint of 1. The minimum root mean square error (3.37 s) is also observed in the 1-s period due to the large number of rays. A cross plot between the number of rays and the root mean square error value is presented in Figure 9c. In general, the number of rays has an impact on the root mean square error value. The root mean square error value decreases as the number of rays increases. Table 1 shows the values for each parameter.

Figure 7

Love wave group velocity maps produced by the FMST tomography process with the best model parameterization analysis results and grid size evaluation at periods of 0.04, 0.1, 0.5, 1, 5, and 10 s. The continuous black lines are traces of faults on the surface referred by Daryono et al. [7]. The thin dotted black line marks the boundary area’s which have fairly good of recovery CKB test.

Figure 8

Love wave group velocity rays were produced by the FMST tomography process at periods of 0.04, 0.1, 0.5, 1, 5, and 10 s.

Figure 9

(a) The number of rays in each period. (b) Root mean square error value in each period. (c) Cross plot number of ray and root mean square error value. (d) Love wave group velocity value as a function of period.

The 1D depth inversion process was carried out with a point spacing of 500 m over an area of 20 km × 40 km. The results of the inversion were used to generate a three-dimensional velocity model, which was then used to characterize the Lembang Fault zone. Figure 10 displays the horizontal S-wave velocity (VsH) map plotted at depths of 0–10 km. At depths shallower than 4 km, the Lembang Fault is characterized by low velocity values, especially in the central part of the Lembang Fault. Furthermore, in the southern part of the Lembang Fault, there are several areas with low velocity values, which are likely associated with the recent volcanic products of Mt. Tangkuban Parahu, which is located to the north of the Lembang Fault. Areas of low velocity also tend to be areas of relatively higher elevation than the surrounding areas. The lateral boundary plane of the Lembang Fault is challenging to identify in the resulting velocity model, possibly due to lateral limitations in the model. To identify the high velocity contrast and the low-velocity zone associated with the fault plane, a vertical cross-section of the velocity model was produced. Figure 11a shows the map of horizontal S-wave velocity at the surface and the positions of the cross sections at lines A–A′, B–B′, and C–C′ as illustrated in Figure 11b–d, respectively. The Lembang Fault plane is identified on a vertical cross-section by recognizing high velocity contrast and low-velocity zone, which is highlighted in black color on the cross section. These zones of high velocity contrast and low velocity are interpreted as fault planes. The three cross sections suggest the existence of three principal fault planes identified as F1, F2, and F3. The F1 faults are associated with fault traces that appear on the surface, several small earthquakes have occurred associated with this fault, while F2 indicates a velocity contrast that may suggest a fault located north of the surface fault trace. The next feature of note is F3, which indicates a fault south of the surface fault trace with a relatively low dip.

Figure 10

Shear wave velocity (VsH) map at depth 0, 1, 2, 3, 4, 5, 6, 7, 8, and 10 km. The blue triangle is the location of the mobile station and the red box is the location of the fixed station. The big red triangle is the crater of Mount Tangkuban Parahu. The continuous black lines are traces of Lembang Fault on the surface referred by Daryono et al. [7].

Figure 11

(a) The map of horizontal S-wave velocity at the surface and the positions of the cross sections at lines A–A′, B–B′, and C–C′ as illustrated in (b)–(d), (b) cross-section A–A′, (c) cross-section B–B′, and (d) cross-section C–C′. F1, F2, and F3 are fault interpreted based on high-velocity contrast and low-velocity anomaly, the blue triangle is the location of the mobile station and the red box is the location of the fixed station. The yellow star marks are the locations of the earthquakes that occurred, which may be associated with the Lembang Fault activity. The big red triangle is the crater of Mount Tangkuban Parahu. The continuous black lines are traces of Lembang Faults on the surface referred by Daryono et al. [7].

4 Discussion

To obtain accurate and reliable images of the subsurface, it is critical to consider survey design analysis and ray density in ambient seismic noise tomography. With denser ray density, tomographic inversion can capture more detailed velocity variations and accurately resolve small-scale structures [58,59,60,61,62,63,64]. To achieve the desired ray density for imaging, careful planning and design of the receiver layout are essential [65,66,67,68,69,70]. An acquisition survey was designed by simulating the placement of 12 mobile stations that switched locations every 4–6 weeks. The analysis was conducted by calculating the number of ray densities in each 1 km × 1 km cell. To achieve good resolution around the Lembang Fault, more than ten rays were needed to pass through each cell. The grid parameters in the inversion process were evaluated using the CKB test, and it was concluded that a grid of 2–4 km for varying periods of Love wave group velocity inversion produced the most optimal resolution outcomes.

The sensitivity and resolution of imaging results are significantly affected by the effects of wavelength and station spacing on ambient seismic noise tomography, when the inter-station spacing is smaller than the wavelength, the probability of noise between stations increases. This can lead to increased noise contamination, which can affect the reliability and accuracy of the tomography results [71,72,73,74,75]. In this study, we exclude data from the Love wave dispersion curve of pairs with an inter-station distance of less than one lambda as input for the Love wave group velocity tomography process. This is due to the shallow target and a data period of fewer than 10 s. Although a wavelength limitation can decrease the number of rays, it may also reduce spatial resolution for each period. Due to the variation in the number of rays in each period, we conducted a CKB test in each period to obtain the optimal grid.

Heterogeneity in S-wave velocity can result from the presence of different rock types, different sediment consolidation rates, fractures, faults, or other subsurface structures. In seismic studies of fault zones, the identification and characterization of high velocity contrast and low-velocity anomalies can provide valuable insight into the nature and behavior of faults. This has been demonstrated in several studies [42,43,44,45,46,47,48,76,77]. Low-velocity anomalies associated with fault planes can be caused by several factors, such as the presence of fluid-filled fractures, velocity decrease due to increased porosity and fluid saturation, weakened mechanical properties, which in this case can help us identify fault planes in order to mitigate seismic hazards.

The S-wave is a type of body wave, and it is possible to estimate the velocity of S-waves from surface wave data because surface waves are the result of the superposition of body waves, both compressional waves and S-waves propagate along the surface [78,79,80]. The S-wave velocity obtained from surface wave especially Love wave analysis in ambient seismic noise tomography primarily reflects the horizontal shear wave velocity. Love waves have horizontal particle motion, providing details on S-wave velocity in the horizontal direction. However, they are less sensitive to the vertical S-wave velocity component, commonly known as the vertical S-wave velocity. Combining Love wave ambient seismic noise tomography with Rayleigh wave ambient seismic noise tomography can be beneficial for fault zone characterization. Love and Rayleigh waves have different characteristics and sensitivities to subsurface structures. Combining these waves provides a comprehensive understanding of subsurface properties, including faults. Integrating data from different wave types provides accurate and reliable estimates of subsurface velocity, seismic anisotropy, and other relevant parameters. This can help to quantify fault zone characteristics and understand how they affect seismic hazard assessment.

In Figure 11b–d, three north-south cross sections display contrasting high-velocity and low-velocity zones that correspond to fault planes, with fault F1, which shows a relatively vertical plane, being the main fault and associated with the surface fault trace. Although not significant in terms of velocity contrast, the presence of the F2 fault is confirmed by earthquakes in its vicinity. The F3 fault dips fairly low and indicates that it extends from west to east toward the center of the Lembang Fault. However, it is not continuous to the east of the Lembang Fault.

Quaternary volcanic products generally dominate the surface around the Lembang Fault Zone and enclose the underlying Tertiary sedimentary rocks [7,15,81,82,83,84]. The Lembang Fault zone lies at the geological boundary between the Bandung and Bogor zones. The Bogor zone represents the central depression zone of West Java, while the Bandung zone is relatively elevated. When there is a discrepancy in the altitudes of the two physiographic zones, the Lembang Fault is triggered [14]. A fault slope near the Lembang Fault (F2) is thought to have developed due to landslides triggered by volcanic eruptions from the Sunda region around the youthful Tangkuban Parahu volcano [14,83]. Some experts assert that the Lembang Fault is a fault with a left-lateral strike-slip [1] and a minor vertical component (thrust fault) [7], geological observations made with a geomorphological approach show a shift in river patterns indicating an active shear fault in the area [7], while others argue that the fault is a normal fault [6,14,15].

Seismic observations suggest that earthquakes on the Lembang Fault are caused by a strike-slip faulting mechanism [11,12,13]; however, the 6–7 km deep seismicity is associated with thrust faulting, while the 3–4 km shallow seismicity is associated with strike-slip faulting [10]. Due to the small magnitude of the earthquakes that occur and the uncertain focal mechanism, it is difficult to determine the mechanism of faulting in this region from seismic activity data. The claim that the Lembang Fault is a normal fault is supported by gravity measurements [17,18] and CSAMT measurements [19]. In this region, GPS geodetic observations provide compelling evidence for tectonic movement due to left-lateral fault strike-slip [16,9]. Four earthquakes can support the interpretation of our velocity model out of several earthquakes that can be associated with the Lembang Fault zone, as shown in Figure 11. Earthquake 1 [13] may be associated with the F2 fault, while earthquakes 2–4 may be associated with the F1 fault [11,13,85]. Earthquake 1 is on the F2 fault at 2 km depth, it has NE orientation and is close to the Burangrang Hill.

To gain a better understanding of tectonic activity in this region, it is necessary to review research studies conducted around the Lembang Fault. Some studies, based on surface geological observations [86] and audio-magnetotelluric measurements [87], have indicated the presence of a thrust system on the Cimandiri Fault, which is located to the southwest of the Lembang Fault (Figure 1) and is believed to be a continuation of it. Our interpretation of the model shows similarities to the model of the Cimandiri Fault, and this finding suggests the presence of a thrust fault (F3) that is located to the south of the surface trace of the Lembang Fault. In addition, a regional tectonic study showed the existence of an active back-arc thrust in northern west Java [88], the study was conducted on the Baribis Fault, north of the Lembang Fault. Our model shows that the Lembang Fault is part of the Java Island fault system and is still influenced by a larger thrust fault system.

Fault characterization is essential for earthquake hazard mitigation. The analysis is based on previous geological, geophysical, and geodetic studies. The Lembang Fault (F1) is the main fault acting in this area and is confirmed by the presence of several small earthquakes associated with this fault, in addition to geodetic observations also confirm that the Lembang Fault is a left-lateral slip fault. F3 is an interesting feature as it indicates a thrust fault located south of the fault trace. This is shown by the presence of a large velocity contrast and a low-velocity anomaly. The fault has the potential to produce earthquakes as it is located in a seismogenic layer, which typically ranges from 3 to 20 km in depth [89]. Thrust faults are caused by high stress levels. In the case of thrust faults, the rock mass moves upward in the direction of the least principal stress, resulting in the lifting of the rock mass. Thus, the overburden has the lowest principal tension. Regardless of the dimensions of the fault, variations in the stress regime in different fault environments can contribute to differences in the stress drop during faulting, which in turn affects the radiant energy. In comparison, earthquakes on thrust faults tend to produce stronger shaking that affects a larger geographic area when compared to earthquakes on strike-slip faults that are of comparable magnitude [76,77,90]. Thrust faults can cause more severe shaking and damage to structures due to the vertical movement of the ground. Such faults can also trigger large landslides by rupturing rock and soil layers through vertical movement. These landslides can pose additional risks to communities, infrastructure, and transportation networks. They can result in increased destruction and potential loss of life. To gain a more complete understanding, the combination of body wave tomography (Vp and Vs), and Love and Rayleigh waves tomography, as well as seismic anisotropy analyses (azimuthal and radial anisotropy), is recommended for the identification of areas of high strain. Alternatively, data from magnetotelluric and gravity methods, which have strong depth penetration capabilities, could also improve model accuracy.

5 Conclusion

Using the Love wave ambient seismic noise tomography method, we generated S-wave velocity models (Vs horizontal or VsH) with high vertical and lateral resolution around the fault zone due to the relatively dense network of data acquisition stations. To better interpret the subsurface conditions of the Lembang Fault zone, we constructed a 3D S-wave velocity model by compiling multiple 1D S-wave velocity depth inversion results.

The analysis showed significant velocity contrasts and low-velocity anomalies, suggesting that the vertical planes may correspond to fault planes. The three cross sections are interpreted as representing three possibilities major fault planes, labeled F1, F2, and F3. F2 is new blind active fault found at 2 km depth with NE orientation associated with hypocenter of an earthquake event.

Lembang Fault has two fault planes: One is a vertical fault and the other is a south dipping fault of thrust fault. The interpretation that the Lembang Fault is a left-lateral slip fault is supported by the F1 fault, which is consistent with the location of the Lembang Fault trace on the surface. There is an indication of the presence of a thrust fault (F3) in the southern part of the Lembang fault. The existence of a thrust fault (F3) needs to be a special concern because the potential earthquake generated by this thrust fault can be a potential disaster threat with a considerable impact.