Magnetotellurics (MT) and Gravity Study of a Possible Active Fault in Southern Garut, West Java, Indonesia

Arisbaya, Ilham; Wijanarko, Edy; Warsa, undefined; Sumintadireja, Prihadi; Sudrajat, Yayat; Handayani, Lina; Mukti, M. Ma'ruf; Grandis, Hendra

doi:https://doi.org/10.1155/2023/4482074

International Journal of Geophysics

On this page

Abstract Introduction Results Discussion Conclusion Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2023 | Article ID 4482074 | https://doi.org/10.1155/2023/4482074

Magnetotellurics (MT) and Gravity Study of a Possible Active Fault in Southern Garut, West Java, Indonesia

Ilham Arisbaya,^1,2Edy Wijanarko,^1,3 Warsa,¹Prihadi Sumintadireja,⁴Yayat Sudrajat,²Lina Handayani,²M. Ma'ruf Mukti,⁵and Hendra Grandis¹

Academic Editor: Salvatore Gambino

Received21 Jan 2023

Revised08 Mar 2023

Accepted11 Mar 2023

Published31 Mar 2023

Abstract

In recent years, dozen low-intensity earthquakes occurred in southern Garut, West Java Indonesia; two of them were reported destructive. However, those shallow earthquake clusters are hardly associated with well-known active faults in the area. Hence, we conducted 3D gravity combined with 2D magnetotellurics (MT) inversions to study the subsurface. Gravity and MT modeling confirm a basin with around 5 km depth consisting of two subbasins separated by a NE-SW trending local-high ridge. The local high coincides with the magmatic intrusion in geothermal fields and aligns with a series of volcanic bodies’ lineament observed on the surface. We interpret this structural high as a preexisting fault that serves as a magma pathway in the tectonomagmatic interaction. Shallow low-magnitude seismicity in the southern Garut area tends to occur in the resistive bodies. We interpret that heat from the cooling magmatic intrusion may decrease the effective fault-normal stress of the rocks, leading to a decrease in fault failure resistance and may initiate rupture. The resistivity structure around the initial rupture may affect whether or not the nucleation will end up as a large-magnitude earthquake. Furthermore, the unconsolidated young volcanic cover in this area could amplify ground shaking when earthquake occurs that might lead to more extensive damage.

1. Introduction

In recent years, the Indonesian Meteorology, Climatology, and Geophysical Agency (BMKG) recorded the occurrence of low-magnitude earthquakes in southern Garut, West Java, Indonesia (Figure 1). At least two, one in 2016 and the other in 2017, of these low-magnitude earthquakes were reported destructive [1, 2]. These earthquakes appear relatively local yet destructive due to their shallow hypocenter depth. These two events were associated with an active fault, identified only recently as the Garsela Fault [2, 3]. The lineament of recent seismicity was used in defining the Garsela Fault which consists of two segments, the Kencana and Rakutak segments [3]. The latter is also known as the Kendang Fault in geothermal exploration studies [4] but was not considered an active fault.

Southern Garut area was characterized as a bowl-like feature of low shear-wave velocity () anomaly that was interpreted as a sedimentary basin with a depth of up to 5 km [5]. However, this information is resulted from a regional-scale study that must be further constrained. Our limited knowledge of the subsurface condition in this area raises further concern since unconsolidated basin can increase ground motion when earthquakes occur [6] and could cause much more damage. Furthermore, southern Garut hosts three geothermal fields, i.e., Darajat, Kamojang, and Wayang-Windu, producing more than 700 MW of electricity. Together with Patuha and Karaha-Talaga Bodas in the west and east, they form a cluster of large geothermal fields in a relatively small area of 100 by 50 km² [7]. Hence, in addition to being an area with disaster potential with possible extensive damage, it is also interesting to study further its remarkable subsurface condition.

Identifying and mapping the geometry of an active fault are crucial and considered as one of the basic knowledge needed in assessing seismic hazard analysis [3]. However, the surface trace of an active fault is often masked and covered by various factors, from thick sediment or soil to young volcanic products in an active volcanic area. The condition becomes more challenging in the low slip-rate regime, where fault lines may only be recognized when the next major event occurs [8]. Since a fault deforms and cuts brittle rocks in response to stress, rocks’ physical parameters will change and differ from the surrounding. In this study, the magnetotelluric (MT) method is used because of its sensitivity to conductive (low resistivity) bodies that may indicate the water-filled weak zone or fault [9]. Furthermore, previous MT studies show the possible association of earthquake’s magnitude with subsurface resistivity distribution. Low-magnitude earthquakes tend to occur within conductive bodies, while higher-magnitude earthquakes are commonly found near the high-low resistivity boundary on the resistive side [10, 11].

Fracturing will also decrease the bulk density of rock mass. However, gravity method is often used to delineate geological structures bounding the faulted blocks with different densities, rather than the fault directly. The gravity method is widely used, especially when traces are not observed on the surface [12]. In this study, we complement the MT analysis with the gravity data from our previous study in southern Garut, West Java, Indonesia [13]. Although the combination of MT and gravity methods has been well established in depicting the subsurface structures in many places, these methods cannot directly distinguish between active and inactive faults. However, MT and gravity methods may assist seismological analysis in explaining seismic phenomena in an area, e.g., Aizawa et al. [14]. We report our attempt to image the subsurface geometry of the southern Garut zone using two-dimensional (2D) MT modeling combined with three-dimensional (3D) gravity modeling to analyze the possible relationship with shallow earthquake events in the southern Garut.

2. Geological Setting

The geological configuration of Java Island is often assumed as a product of orthogonal subduction in the south of the island with a relatively simple structure. The physiographic zone extends along the island axis (W-E), parallel to the subduction line in the south of Java [15]. This W-E trending trench was thought to be a new subduction path that has been active since the Middle Eocene [16]. The former subduction path was interpreted based on Cretaceous ophiolite found in Ciletuh, Karangsambung, and Bayat [17] and then turned NE-SW towards the Meratus Complex in Borneo. This Cretaceous subduction was supposed to have ceased when the fragment of Gondwana microcontinent collided in the Late Cretaceous [16].

The subduction in southern Java was then restarted in a new W-E trench direction, with the main product of Middle Eocene to Miocene volcanic rocks known as old andesite formations (OAF) [16]. OAF, also known as the southern mountain zone, is dominant in the southern part of West Java and East Java [18]. Volcanism shifted northward in the Plio-Pleistocene, and the modern active volcanoes take the same place as the former volcanism period. Figure 2 shows that Quaternary volcanic products dominate the surface geology in the southern Garut area, covering the OAF [19, 20]. There are several prominent volcanoes in this area. Some volcanoes are dormant (e.g., Mt. Cikuray, Mt. Malabar, Mt. Kencana, and Mt. Kendang), while the other two (Mt. Papandayan and Mt. Guntur) are active [21, 22].

Current subduction in the south of Java Island is almost trench-perpendicular with a convergence rate from 58 to 65 mm/year from west to east [23]. The neotectonic process generates earthquakes in the subduction megathrust and onshore faults, and most of them rupture in the relatively W-E direction, such as the the Cimandiri [8, 24], Lembang [25, 26], and Baribis Faults [27, 28] in the western part of Java island. Several other seismic activity clusters have been reported with a relative NE-SW direction in the southern part of Java Island in recent years, one of which is the Garsela Fault [29]. Little is known about the Garsela Fault, and the interpretation was mainly based on the focal mechanism analysis of the 2016 and 2017 earthquakes [2]. The 2016 earthquake showed a strike-slip mechanism, while the 2017 earthquake showed a normal fault mechanism. The Garsela or Kendang Fault is interpreted as a preexisting fault that acts as a path for magma rise for Mt. Kendang that is active in the Plio-Pleistocene [4].

NE-SW is reported as dominant effective fracture direction in the reservoir of Wayang-Windu, Darajat, and Kamojang geothermal fields [4, 30, 31]. The NE trending fault was considered a regional shear fault that had been reactivated as a normal fault [30]. Resistivity model based on MT data in the southern Garut shows a subsurface horsts and graben structure, indicating an extensional tectonic regime [32]. Furthermore, the Java Island deformation model based on regional GPS observation data also shows that the southern region of Garut is currently in an extensional tectonic regime [33].

3. Methodology

3.1. Gravity

The gravity method measures the earth’s gravitational field variations due to the density inhomogeneity of rocks that compose the earth’s crust. In practice, the subsurface density variation may be interpreted as basement geometry of a basin, lateral lithological changes, and geological structures. We utilized the results of our previous study in southern Garut using the gravity method [13] to complement the MT analysis. Here, we highlight some key points of the study as an overview.

We digitized the Java Bouguer anomaly map produced from the collaboration between the Geological Survey of Indonesia and the Geological Survey of Japan [35]. The original map was compiled based on terrestrial measurements with station interval ranging from 2 to 5 km. The resulting Bouguer anomaly map is shown in Figure 3. Raw field measurement data were processed using standard gravity data workflow including corrections for tidal and instrumental drift. The theoretical gravity () was calculated using the 1967 international gravity formula and was adjusted with free air (), the Bouguer (), and terrain () corrections. The Bouguer anomaly () was calculated by using the following well-known equation [36]: where is the elevation of the gravity station and is the average rock density of 2.67 g/cm³ used to calculate the Bouguer correction (). The terrain correction () was calculated up to 10 km radius from each gravity station or zone K according to Hammer’s chart scheme. These corrections aim at removing the effect of the nongeological component observed in the gravity measurement data.

We performed 3D gravity modeling to obtain subsurface density distribution for a limited local subsurface area, relative to overall gravity data coverage and to a maximum depth of 10 km. In this context, gravity modeling requires assuming the shape of the base anomaly, i.e., the regional field [37]. We used the result from the 7500-meter upward continuation transformation for the regional-residual decomposition. Detailed procedures and parameters for the 3D gravity inversion can be found in Arisbaya et al. [13].

3.2. Magnetotellurics (MT)

3.2.1. Data Acquisition and Analysis

The magnetotelluric (MT) method employs the natural electromagnetic field to study the conductivity properties of the rocks. The temporal variations of the electric ( and ) and magnetic (, , and ) fields are measured using two pairs of electrodes and induction coils, respectively. Time series data from the field measurement are transformed into frequency domain data. Then, the MT data are commonly expressed as 2 by 2 impedance tensor relating the horizontal components of the electric field to magnetic field estimated from their respective spectra following the standard method [38]. We used the MT method to depict the subsurface electrical resistivity structure down to 10 km depth.

The complex impedance tensor contains information about the earth’s resistivity at a given location as function of depth. As a complex number, each component of consists of real and imaginary parts and thus can be expressed in terms of apparent resistivity related to the amplitude and phase as follows [38]: where is the angular frequency, is the magnetic permeability of free space (4π 10^-7 H/m). The apparent resistivity can be considered as the volume (bulk resistivity) at a certain frequency depicting the Earth’s subsurface structure at a certain depth [38].

In this study, we used two MT datasets acquired in 2011 and 2018 (see Figure 2). The 2011 dataset [32] consists of 41 MT sounding sites, forming two lines in a nearly W-E trending direction (lines C-C and D-D). The 2018 survey consists of 32 sounding stations forming a WNW-ESE trending line along 22 km (line E-E). The interval between stations was around 750 m. We have audio-frequency MT (AMT) data covering the frequency range from 10⁴ Hz down to 10^-1 Hz at all stations, with only 11 stations among them that were complemented with lower-frequency (10² Hz down to 10^-3 Hz) MT data. The interval between stations with MT data was approximately 2 km. The frequency range reflects the depth of investigation, with the lower frequency (or longer period) associated with a deeper depth. Both surveys (2011 and 2018) used two sets of Phoenix MTU-5A instruments. We have a total of 73 sounding data in three lines, but 4 stations’ data were excluded in further analysis due to poor data quality. We processed raw field MT data, including instrument calibration, transformation from time series data to power spectra, and MT transfer function estimation. Transient EM (TEM) data were not available for static shift correction for MT data. However, we assumed that the static shift is very limited since it appears only in a few stations and can be accounted for in the modeling.

We used the phase tensor analysis [39] for data dimensionality and regional geoelectric strike estimations. A small skew angle value is necessary (but not sufficient) for quasi-2D, and the value of is recommended as the threshold for the conductivity to be considered 2D [40]. Figure 4 shows the results of the phase tensor analysis of the three MT lines, and the value dominates almost all measurement stations up to a frequency range of 10^-1 Hz. Large skew angle values were observed at lower frequencies at all observation stations, implying a more complex structure at depth. Large skew angle values at high frequencies were observed locally at some stations, indicating shallow inhomogeneity.

By assuming the same regional geological conditions for the entire investigation area, the geoelectric strike in this area was estimated based on the principal axis direction (-) of the phase tensor ellipse of all stations from the three lines. Figure 5 shows that the geoelectric strike has almost N-S trend (N0°E to N15°E) for the whole frequencies, with relatively more well-defined direction for higher frequency range (>10 Hz). The geoelectric strike is more varying in the medium (10-0.1 Hz) and low (<0.1 Hz) frequency ranges.

(a)

(b)

(c)

(d)

The direction of the N15°E geoelectric strike correlates well with the trend of geological structures in the Wayang-Windu and Darajat [4, 30], basement structure patterns in West Java [41], and the dominant pattern found in the residual anomaly of the gravity data. Therefore, we deduced that the geoelectric strike in the study area is NNE-SSW and choose the N15°E direction. The geoelectric strike value was used to rotate the impedance tensor to minimize the diagonal components ( and ) and maximize the off-diagonal components ( and ) in the two-dimensional (2D) modeling. Rotation to the geoelectric strike results in an apparent resistivity component obtained with the electric field parallel to the geoelectric strike () as the transverse electric (TE) mode while the other one obtained with the magnetic field parallel to the geoelectric strike () as the transverse magnetic (TM) mode [38].

3.2.2. 2D MT Inversion

The data components used in the inversion were TE and TM modes in the frequency range of 0.1-10,000 Hz (five decades). 2D MT data inversion was done using the nonlinear conjugate gradients (NLCG) algorithm [42]. The initial model used was a 100 Ohm.m homogeneous half-space with a horizontal mesh size of approximately 500 m wide while the vertical mesh size increasing with a factor of 1.3. The inversion was carried out using an error floor value of 15% and 10% for the apparent resistivity of TE and TM modes, respectively, and 10% for phases of both TE and TM modes. We incorporated the static shift as model parameters in the inversion, as it is allowed by the software. For the weighting function and the smoothing parameter, we used the default values of the software, 1 for alpha and 1.5 for the beta [42]. The regularization parameter (tau), which is used for trade-off between the misfit and the model smoothness, was determined by performing several inversions using different tau values. The RMS misfit versus roughness plot produces an L-curve [43] to determine the optimum value of tau (Figure 6).

4. Results

4.1. Density Model

The gravity data shows that the southern Garut zone is located in the low Bouguer anomaly zone with a geometry of approximately 25 km by 35 km elongated in the W-E axis (see Figure 3). The low gravity anomaly in southern Garut forms two subareas separated by a NE-SW trending local-high ridge, whereas Darajat and Kamojang geothermal fields are located in this local-high ridge. The eastern gravity anomaly appears to have shifted towards NE, or conversely, the western low area anomaly is shifted towards SW, supporting the regional strike-slip fault hypothesis proposed by Bogie et al. [30]. Figure 7 shows the residual anomaly map in the study area. It is dominated by segmented high anomalies in the N-S direction, which match the surficial regional geological structure information [41]. The Bouguer anomaly map in the interest area shows that an anomalous pattern matches the NE-SW structure, while several other structural patterns (N-S, W-E, and NW-SE) correspond to residual anomalies.

Figure 8 are maps and model associated with 3D gravity inversion, i.e., the Bouguer anomaly as observed data, calculated model response, difference between observed and calculated data or misfit, and the perspective view of the 3D final density model. It shows that the gravity data have been modeled relatively well by the 3D inversion. A prominent misfit is observed in the NE area. However, the density model in the high misfit area can be ignored and should be avoided in interpretation since the study focuses on the low anomaly in the center of the area. The A-A section shows a basin consisting of two subbasins, bounded by a local high approximately 10 km wide. The B-B section shows a low-density model to a depth of around 6 km below the active volcano Mt. Papandayan. However, with spatial data resolution ranging 2-5 km between measuring points, attempts to describe the detailed features of the active volcano were not feasible.

(a)

(b)

(c)

(d)

4.2. Resistivity Model

Figures 9–11 present the pseudosections of apparent resistivity and phase in TE and TM modes for the observed data and calculated responses from 2D MT model at C-C, D-D, and E-E lines, respectively. The pseudosection can help visually identify subsurface resistivity distribution before modeling. All three pseudosections show the dominance of low to medium apparent resistivity near the surface, with thickening observed in the eastern part of the line C-C. Significant low apparent resistivity thickening was observed in the western part of line D-D, whereas in line E-E, the low apparent resistivity was distributed more evenly with slight thickening in the middle part of the line. However, such a qualitative interpretation is only a simplistic description, especially in a complex geological environment.

Pseudosections can be helpful to roughly compare the measured data with the resulting model response. In general, pseudosections of the three MT lines show a good agreement between the model response and the measured data. However, slight differences are observed at low frequencies. Since MT data at longer periods typically have lower quality compared to the short periods, indicated by larger error bars, modeling the low-frequency data becomes more challenging. Furthermore, the low-frequency data indicates three-dimensional conductivity. It is possible that the 2D inversion failed to model the data properly.

The resistivity model from 2D MT inversion of line C-C, D-D, and E-E is shown in Figure 12. The global misfits between observed and calculated data (RMS error) for MT models at line C-C, D-D, and E-E are 2.823, 2.450, and 2.956, respectively. All three cross-sections show superficial low (less than 15 Ohm.m) to moderate (15-15 Ohm.m) resistivity. The thinnest low-resistivity layer is observed in the C-C line, while the thickest is spotted in the D-D line. High resistivity (more than 500 Ohm.m) at depth characterizes the basement overlain by a thin transition layer of 150 to 500 Ohm.m.

(a)

(b)

(c)

Basement disturbance in the form of narrow resistivity gaps of about 30 Ohm.m is observed in the western part of the C-C (R1) and in the central part of D-D (R2). Since the MT method is sensitive to conductive objects, a sensitivity test was needed to test the robustness of the resistive bodies (R1 and R2). The high-resistivity basement appears to be cut by a medium resistivity (15-150 Ohm.m) in a subvertical geometry (C2) at the western part of E-E cross-section. A small, localized conductive body is also observed in the middle part of the same section to a depth level of about 3 km (C1).

A model sensitivity test was conducted to test whether C1, C2, R1, and R2 are part of a robust model or modeling artifacts. The sensitivity test was carried out by changing the resistivity value of the block being tested to its background resistivity. The block is then locked, while the remaining blocks are left free and subject to change when the inversion was rerun. The fitting curve of apparent resistivity and phase at the closest station to the tested block is compared, before and after changing the resistivity value.

The sensitivity test result shows that the station G36 and G37 misfits on the line C-C tend to remain the same after the basement discontinuity is removed (Figure 13). The sensitivity test also shows that the global misfit on the line C-C tends to slightly decrease, from 2.823 to 2.798. The same thing was observed at stations G08 and G09 on the line D-D, i.e., the global misfit only changed slightly from 2.450 to 2.449. The sensitivity test results show that the basement discontinuity feature in the lines C-C and D-D is inconclusive. Such phenomena may be generated by relatively small resistivity contrasts or modeling artifacts generated in the inversion. In contrast to the results of the C-C and D-D lines, the sensitivity test of line E-E showed an increase in global misfit from 2.956 to 3.254. Misfit at stations ST01 and ST19 also increased (especially in the apparent resistivity curve). Thus, the conductive bodies C1 and C2 are indeed needed in modeling and are robust features.

5. Discussion

Figures 14–16 show the collocated resistivity and density models of line C-C, D-D, and E-E, respectively. The near-surface part of all three cross-sections is dominated by low to medium resistivity and density. Surface geology shows that the low to medium resistivity and density values are associated with young volcanic products (Pleistocene-Holocene) which may not have been well compacted. On the other hand, the underlying basement was interpreted based on high-resistivity and high-density values.

An interesting near-surface feature is observed in C-C at a distance of 35 to 40 km from the west end of the line (Figure 14). In this position, the gravity anomaly increases while the resistivity model shows a low-resistivity body. The surface geological map at the location of the conductive body shows a hill (Mt. Malang) which is also the intersection of the NE-SW fault with the WNW-ESE fault. There are two possible interpretations for this phenomenon. The first possibility is its parallel location with Darajat and Kamojang geothermal fields, which raises a suspicion that the conductive body may be a hydrothermal area. However, to date, no surface hydrothermal manifestation has been reported in Mt. Malang. Another possibility is that the intersection of two faults causes an intensive weak zone and when filled with fluid may decrease rock resistivity values. A similar result was conveyed earlier in the Wayang-Windu area based on an analysis of SAR satellite imagery [44]. The intersection area of two fault directions shows a low rock strength index value, which is interpreted as an intensive weak zone.

A conductive rock body was observed in the center of the cross-section at a depth of about 4 to 5 km that coincides with the location of Mt. Kencana that was active in the Pleistocene (Figure 16). The 2016 earthquake hypocenter is located under the conductive body at a depth of about 17 km. The microearthquake (MEQ) cluster occurred at about 4 km south of Mt. Kencana and down to a depth of about 10 km and is thought to be caused by tectonic activity [45]. The gravity anomaly decreases in this location and can be interpreted as pyroclastic products of Mt. Kencana or a fault that failed to be imaged in the resistivity model at greater depths. The failure of the resistivity model to image the weak zone may be due to its too small geometry and lack of fluid content. A microgravity study reported that the NW-SE trending faults in the Wayang-Windu area have a sealing property [46] which inhibits fluid circulation. Another possibility is that the weak zone has been overprinted by the magmatic intrusion, as shown at Darajat [4].

The shallowest basement, interpreted based on high-resistivity and high-density values, is observed in the C-C section, and the deepest is spotted in the D-D section at about 5 km deep. A disturbed basement is observed in the western part of the C-C cross-section, which coincides with Mt. Malabar’s opening to NW. Thus, the disturbed basement may be related to the graben in the magma pocket of old Mt. Malabar. The basement discontinuity in the western part of E-E can be interpreted as a fault, but the limited MT data distribution makes it difficult to determine the lateral continuity.

Basement on the D-D cross-section is observed at a depth of 4 to 6 km, shallowing to a depth of 2 to 3 km in the middle of the cross-section forming a local high of about 10 km wide. Based on the resistivity and density model, as well as borehole data in Darajat geothermal field, the Garsela/Kendang Fault is interpreted to be overprinted by igneous rock intrusion at a depth of about 1500 m [4]. Meanwhile, a high-resistivity anomaly in Kamojang is confirmed as an igneous rock intrusion at a depth of about 2000 m from the surface [31]. Another basement undulation was observed in the western part of the D-D cross-section, both in the density and the resistivity model. Based on the location, this undulation may represent microdiorite intrusion, which acts as the heat source of Wayang-Windu geothermal field [47].

The position of the basement undulations coincides with the magmatic intrusion in geothermal fields. It indicates that the magmatic intrusion may be part of the NE-SW trending local-high ridge which is observed as a young volcanic body lineament on the surface [29]. The lineament of intrusive bodies can be interpreted as a series of dyke intruded through the weak zone of fault plane. A previous study reported MEQ clusters around Darajat and Wayang-Windu related to geothermal activity [45]. However, the earthquake relocation indicates that the depth of the 2017 earthquake is difficult to be interpreted as related to geothermal activity.

Gravity modeling in the Darajat field indicates that the upper part of the microdiorite intrusion was identified at about 1500 m deep and is associated with a density contour of 2.5-2.53 g/cm³ [48]. Here, we use the value of 2.525 g/cm³ as the basement-sediment threshold to generate a basement depth map presented in Figure 17. The deepest part of the basement is at a depth of 5-6 km, supporting the result from ambient noise tomography [5]. The basement depth map is dominated by local deep with NE-SW direction flanking the local high where Darajat, Kamojang, and Mt. Guntur are situated. In addition, an NW-SE trending local deep extended from Mt. Papandayan towards Wayang-Windu matches the low-velocity anomaly reported previously [45].

The epicenter of the 2017 earthquake, located between two active volcanoes, Mt. Papandayan and Mt. Cikuray, raises the question whether the earthquake has a relationship with the volcanic activity. However, Supendi et al. [1] stated that the 2017 earthquake waveforms could be detected clearly by the BMKG seismic stations located far from the area, while volcanic earthquakes are usually very local and are only detected by seismometers around the volcano. Thus, the authors argue that the 2017 earthquake is most likely a tectonic earthquake. However, volcano unrest may reactivate preexisting fault nearby as reported in Mt. Agung [49] and Mt. Etna [50].

The cases mentioned above highlight the interaction between tectonic and magmatic activity, which in the southern Garut area the magmatic activity occurs in the form of a modern active volcano and/or paleomagmatic activity that currently hosts geothermal fields. Recent studies have shown that earthquake nucleation is not only influenced by rheology and stress accumulation but also factors such as temperature, fluid pore pressure, and permeability that also play an important role in tectonomagmatic interactions [11, 51]. It is thought that the Garsela/Kendang Fault formerly acted as a magma pathway for Mt. Kendang and currently appears to overprint the fault in the depths [48]. Rigid solidified cooling intrusive magmatic rock was proposed to be able to weld two adjacent sliding rock blocks of a fault, increasing the bonding between the two blocks and possibly decreasing the chance of an earthquake [51]. On the other hand, the same authors also argue that the upward force and heat carried by the magma body may decrease the effective fault-normal stress and lead to a decrease in fault failure resistance. A similar interpretation was also proposed to explain the relationship between seismicity at the center of the Mosha Fault and the Damavand active volcano [52]. We interpret that the shallow low-magnitude earthquakes in southern Garut also resulted from tectonomagmatic interactions.

Low-magnitude earthquakes often occur within conductive bodies, whereas higher-magnitude earthquakes are frequently located close to the high-low resistivity boundary on the resistive side. Stress can accumulate significantly on the rheologically rigid resistive block and produce large earthquakes while on the conductive side stress is released in the form of microearthquakes due to the decreasing effective strength caused by fluid-rich content in the rock mass [10]. A recent MT study highlighted that significant earthquakes in the 2016 Kumamoto earthquake sequence are more likely to start near the high-low resistivity boundary than those initiated far or within the fluid-rich low-resistivity zone [11]. They propose that fluid pore pressure plays a significant role in the evolution of crustal earthquake rupture. If the rupture nucleates near the high-low resistivity anomaly boundary, the high prefailure pressure/temperature (PT) gradient in the pore fluids may promote propagation and increase the size of the rupture, resulting in macroscopic rupture. In contrast, if rupture nucleates far from the high-low resistivity boundary, the low PT gradient in the pore fluid may be less likely to promote rupture growth, resulting in a low-magnitude earthquake.

Two possible weak zones in this study are observed in the eastern part of line C-C and the central part of line E-E. Interestingly, low-magnitude earthquakes in this study tend to occur in resistive bodies rather than the conductive part. Since they are located relatively distal to the high-low resistivity boundary, it is understandable that the earthquakes in southern Garut tend to show low magnitudes. However, large-magnitude earthquakes could still occur in southern Garut, and further discussion remains open. Since sediments in the shallow basin may amplify ground motion during earthquakes, the interpreted basin in southern Garut should raise a concern to mitigate the potential damage that may occur in the future.

Several MT surveys have imaged low-resistivity zone in the midcrust of several island arc settings, often spatially related to earthquakes. The resistivity model presented in this study tends to show high-resistivity homogeneity from 5 km depth down. Based on slight differences at long periods that indicate 3D conductivity, we suspect that this is due to the incompatibility of the 2D inversion we used here. Thus, we plan to implement 3D inversion as well as acquire more MT data in our near future projects.

6. Conclusion

Active fault earthquakes usually occur within a few tens of kilometers of hypocenter depth in the upper crust. Geophysical methods that can describe the subsurface up to that depth order are important for studying the possible condition or structures related to earthquakes. However, each geophysical method has its advantages and disadvantages and is complementary to one another. The gravity method has been widely used in fault research, with relatively simple data acquisition allowing for rapid coverage of vast regions with a specific spatial resolution. However, the drawback of the gravity method is its ambiguity especially in the vertical direction, such that gravity modeling requires extra information as constraint. The MT method, with better vertical resolution, can provide complementary information. Nevertheless, since MT data collecting time is rather long, research with MT is usually limited for a particular selected area defined by results from gravity survey with broad regional or subregional coverage.

We have utilized 2D MT modeling combined with 3D gravity modeling to image the subsurface geometry of the southern Garut region. Our gravity and MT modeling confirmed the existence of a basin in the study area that consists of two subbasins separated by a NE-SW trending local-high ridge. The local high coincides with magmatic intrusion in the geothermal fields and a series of volcanic bodies’ lineament observed on the surface. We interpret this feature as a preexisting fault that serves as magma pathways in the tectonomagmatic interaction. It implies that a significant interplay of tectonic and magmatic activity has shaped the geology of southern Garut, including the recently observed shallow seismicity. We interpret that heat from the cooling magmatic intrusion decreases the effective fault-normal stress of the rocks, which in turn leading to earthquakes of low magnitude as observed recently.

This study supports the interpretation suggesting that the pore fluid’s PT gradient contributes to earthquake initiation and propagation. Thus, the subsurface resistivity structure may provide hints in investigating the nucleation of crustal earthquakes. Furthermore, our study demonstrates an opportunity to explain the schematic mechanism of earthquakes within tectonomagmatic framework widely observed in various places such as Indonesia. Similar studies can be conducted in other earthquake-prone areas with thick sediments and young volcanic products covering the surface, as is common in many areas, particularly on the heavily populated island of Java.

Data Availability

The dataset analyzed during the current study is available from the corresponding author on reasonable request.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Authors’ Contributions

IA, HG, W, and PS conceived the initial ideas of the study. IA and YS participated in the field survey. IA, EW, HG, and YS processed and analyzed the data. IA, HG, W, PS, LH, and MMM contributed to the data interpretation. IA, EW, HG, LH, and MMM contributed to the preparation of the manuscript. All authors read and approved the final manuscript.

Acknowledgments

H. Grandis, W. Warsa and P. Sumintadireja are grateful to Institute for Research and Community Services (LPPM) ITB for continuous support. National Research and Innovation Agency of Indonesia (BRIN) provides the Saintek scholarship awarded to I. Arisbaya.

References

P. Supendi, A. D. Nugraha, N. T. Puspito, S. Widiyantoro, and D. Daryono, “Identification of active faults in West Java, Indonesia, based on earthquake hypocenter determination, relocation, and focal mechanism analysis,” Geoscience Letters, vol. 5, no. 1, pp. 1–10, 2018.
View at: Publisher Site | Google Scholar
P. Supendi, A. D. Nugraha, and S. Widiyantoro, “Recent destructive earthquakes around Garut area, West Java, Indonesia: an unidentified fault?” AIP Conference Proceedings, vol. 1987, article 020077, no. 1, 2018.
View at: Google Scholar
M. Irsyam, S. Widiyantoro, D. H. Natawidjaja et al., Map of Indonesia’s Earthquake Sources and Hazards 2017, Center for Housing and Settlements Research and Development, Ministry of Public Works and Public Housing, Bandung, 2017.
R. G. Intani, G. U. Golla, Y. Syaffitri et al., “Improving the conceptual understanding of the Darajat geothermal field,” Geothermics, vol. 83, article 101716, 2020.
View at: Publisher Site | Google Scholar
S. Rosalia, S. Widiyantoro, P. R. Cummins et al., “Upper crustal shear-wave velocity structure beneath Western Java, Indonesia from seismic ambient noise tomography,” Geoscience Letters, vol. 9, no. 1, 2022.
View at: Publisher Site | Google Scholar
A. Cipta, P. Cummins, M. Irsyam, and S. Hidayati, “Basin resonance and seismic hazard in Jakarta, Indonesia,” Geosciences, vol. 8, no. 4, p. 128, 2018.
View at: Publisher Site | Google Scholar
I. B. Raharjo, R. G. Allis, and D. S. Chapman, “Volcano-hosted vapor-dominated geothermal systems in permeability space,” Geothermics, vol. 62, pp. 22–32, 2016.
View at: Publisher Site | Google Scholar
G. I. Marliyani, J. R. Arrowsmith, and K. X. Whipple, “Characterization of slow slip rate faults in humid areas: Cimandiri fault zone, Indonesia,” Journal of Geophysical Research - Earth Surface, vol. 121, no. 12, pp. 2287–2308, 2016.
View at: Publisher Site | Google Scholar
M. Becken, O. Ritter, P. A. Bedrosian, and U. Weckmann, “Correlation between deep fluids, tremor and creep along the central San Andreas Fault,” Nature, vol. 480, no. 7375, pp. 87–90, 2011.
View at: Publisher Site | Google Scholar
A. Gürer and M. Bayrak, “Relation between electrical resistivity and earthquake generation in the crust of West Anatolia, Turkey,” Tectonophysics, vol. 445, no. 1-2, pp. 49–65, 2007.
View at: Publisher Site | Google Scholar
K. Aizawa, S. Takakura, H. Asaue et al., “Electrical conductive fluid-rich zones and their influence on the earthquake initiation, growth, and arrest processes: observations from the 2016 Kumamoto earthquake sequence, Kyushu Island, Japan,” Earth, Planets and Space, vol. 73, no. 1, p. 12, 2021.
View at: Publisher Site | Google Scholar
D. Ombati, J. Githiri, M. K’orowe, and E. Nyakundi, “Delineation of subsurface structures using gravity data of the shallow offshore, Lamu Basin, Kenya,” International Journal of Geophysics, vol. 2022, Article ID 3024977, 12 pages, 2022.
View at: Publisher Site | Google Scholar
I. Arisbaya, E. Wijanarko, H. Grandis, W. Warsa, and P. Sumintadireja, “An attempt to depict deformation under young volcanic covered area using 3D gravity data inversion in southern Garut,” in 1st ITB Graduate School Conference, pp. 1–15, Bandung, 2021.
View at: Google Scholar
K. Aizawa, H. Asaue, K. Koike et al., “Seismicity controlled by resistivity structure: the 2016 Kumamoto earthquakes, Kyushu Island, Japan,” Earth, Planets and Space, vol. 69, no. 1, p. 4, 2017.
View at: Publisher Site | Google Scholar
R. W. van Bemmelen, The Geology of Indonesia, Govt. Printing Office, The Hague, 1949.
H. R. Smyth, R. Hall, and G. J. Nichols, “Cenozoic volcanic arc history of East Java, Indonesia: the stratigraphic record of eruptions on an active continental margin,” in Formation and Applications of the Sedimentary Record in Arc Collision Zones, A. E. Draut, P. D. Clift, and D. W. Scholl, Eds., vol. 436, no. 10, pp. 199–222, Geological Society of America, 2008.
View at: Google Scholar
K. Wakita and W. Bambang, “Cretaceous radiolarians from the Luk-Ulo Melange Complex in the Karangsambung area, Central Java, Indonesia,” Journal of Asian Earth Sciences, vol. 9, no. 1-2, pp. 29–43, 1994.
View at: Publisher Site | Google Scholar
B. Clements, R. Hall, H. R. Smyth, and M. A. Cottam, “Thrusting of a volcanic arc: a new structural model for Java,” Petroleum Geoscience, vol. 15, no. 2, pp. 159–174, 2009.
View at: Publisher Site | Google Scholar
M. Alzwar, S. Bachri, and N. Akbar, Geological Map of the Garut and Pameungpeuk Quadrangle, Java, Scale 1:100.000, Geological Survey of Indonesia, Bandung, 1992.
M. Koesmono, Kusnama, and N. Suwarna, Geological Map of the Sindangbarang and Bandarwaru Quadrangle, Java, Geological Survey of Indonesia, Bandung, 2nd edition, 1996.
S. Byrdina, H. Grandis, P. Sumintadireja et al., “Structure of the acid hydrothermal system of Papandayan volcano, Indonesia, investigated by geophysical methods,” Journal of Volcanology and Geothermal Research, vol. 358, pp. 77–86, 2018.
View at: Publisher Site | Google Scholar
D. K. Syahbana, C. Caudron, P. Jousset et al., “Fluid dynamics inside a "wet" volcano inferred from the complex frequencies of long-period (LP) events: an example from Papandayan volcano, West Java, Indonesia, during the 2011 seismic unrest,” Journal of Volcanology and Geothermal Research, vol. 280, pp. 76–89, 2014.
View at: Publisher Site | Google Scholar
A. Koulali, S. McClusky, S. Susilo et al., “The kinematics of crustal deformation in Java from GPS observations: implications for fault slip partitioning,” Earth and Planetary Science Letters, vol. 458, pp. 69–79, 2017.
View at: Publisher Site | Google Scholar
I. Arisbaya, L. Handayani, M. M. Mukti, Y. Sudrajat, H. Grandis, and P. Sumintadireja, “Imaging the geometry of Cimandiri fault zone based on 2D audio-magnetotelluric (AMT) model in Nyalindung, Sukabumi-Indonesia,” Pure and Applied Geophysics, vol. 176, no. 11, pp. 4833–4845, 2019.
View at: Publisher Site | Google Scholar
A. Afnimar, E. Yulianto, and R. Rasmid, “Geological and tectonic implications obtained from first seismic activity investigation around Lembang Fault,” Geoscience Letters, vol. 2, no. 1, p. 4, 2015.
View at: Publisher Site | Google Scholar
M. R. Daryono, D. H. Natawidjaja, B. Sapiie, and P. Cummins, “Earthquake geology of the Lembang Fault, West Java, Indonesia,” Tectonophysics, vol. 751, pp. 180–191, 2019.
View at: Publisher Site | Google Scholar
S. Aribowo, L. Husson, D. H. Natawidjaja et al., “Active back-arc thrust in North West Java, Indonesia,” Tectonics, vol. 41, no. 7, 2022.
View at: Publisher Site | Google Scholar
D. S. Widarto, F. Febriani, B. Nurdiyanto et al., “Magnetotelluric resistivity imaging of the Baribis fault zone’s Majalengka segment in West Java, Indonesia,” Acta Geodaetica et Geophysica, vol. 57, no. 1, pp. 177–194, 2022.
View at: Publisher Site | Google Scholar
I. Arisbaya, H. Lestiana, M. M. Mukti, L. Handayani, H. Grandis, and P. Sumintadireja, “Garsela Fault and other NE-SW active faults along the southern part of Java Island,” IOP Conference Series: Earth and Environmental Science, vol. 789, article 012065, 2021.
View at: Google Scholar
I. Bogie, Y. I. Kusumah, and M. C. Wisnandary, “Overview of the Wayang Windu geothermal field, West Java, Indonesia,” Geothermics, vol. 37, no. 3, pp. 347–365, 2008.
View at: Publisher Site | Google Scholar
M. Y. Kamah, A. Armando, D. L. Rahmani, and S. Paramitha, “Enhancement of subsurface geologic structure model based on gravity, magnetotelluric, and well log data in Kamojang geothermal field,” in 6th ITB International Geothermal Workshop, vol. 103, p. 012013, 2017.
View at: Publisher Site | Google Scholar
L. Handayani, K. Kamtono, and D. D. Wardhana, “Extensional tectonic regime of Garut Basin based on magnetotelluric analysis,” Indonesian Journal on Geoscience, vol. 8, no. 3, pp. 127–133, 2014.
View at: Publisher Site | Google Scholar
E. Gunawan and S. Widiyantoro, “Active tectonic deformation in Java, Indonesia inferred from a GPS-derived strain rate,” Journal of Geodynamics, vol. 123, pp. 49–54, 2019.
View at: Publisher Site | Google Scholar
A. D. Nugraha, H. A. Shiddiqi, S. Widiyantoro et al., “Hypocenter relocation along the Sunda Arc in Indonesia, using a 3D seismic-velocity model,” Seismological Research Letters, vol. 89, no. 2A, pp. 603–612, 2018.
View at: Publisher Site | Google Scholar
M. Untung and Y. Sato, Gravity and Geological Studies in Java, Indonesia, vol. 6, Geological Survey of Indonesia and Geological Survey of Japan, 1978.
R. J. Blakely, Potential Theory in Gravity and Magnetic, Cambridge University Press, Cambridge, 2009.
M. Pirttijarvi, “GRABLOX: gravity interpretation and modelling software based on a 3-D block model,” User’s guide - Version 1.6b, 2008.
View at: Google Scholar
F. Simpson and K. Bahr, Practical Magnetotellurics, Cambridge University Press, Cambridge, 2009.
T. G. Caldwell, H. M. Bibby, and C. Brown, “The magnetotelluric phase tensor,” Geophysical Journal International, vol. 158, no. 2, pp. 457–469, 2004.
View at: Publisher Site | Google Scholar
J. R. Booker, “The magnetotelluric phase tensor: a critical review,” Surveys in Geophysics, vol. 35, no. 1, pp. 7–40, 2014.
View at: Publisher Site | Google Scholar
Sribudiyani, N. Muchsin, R. Ryacudu et al., “The collision of East Java microplate and its implication for hydrocarbon occurences in the East Java Basin,” in 29th Annual Convention Proceedings of Indonesian Petroleum Association, 2003.
View at: Google Scholar
W. L. Rodi and R. L. Mackie, “Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion,” Geophysics, vol. 66, no. 1, pp. 174–187, 2001.
View at: Publisher Site | Google Scholar
P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Review, vol. 34, no. 4, pp. 561–580, 1992.
View at: Publisher Site | Google Scholar
A. Saepuloh, R. H. Saputro, M. N. Heriawan, and D. Malik, “Integration of thermal infrared and synthetic aperture radar images to identify geothermal steam spots under thick vegetation cover,” Natural Resources Research, vol. 30, no. 1, pp. 245–258, 2021.
View at: Publisher Site | Google Scholar
P. Jousset, R. Sule, W. Diningrat et al., “Seismic tomography and monitoring of magmatic geothermal and natural hydrothermal systems in the south of Bandung, Indonesia,” World Geothermal Congress 2015, 2015.
View at: Google Scholar
Setianingsih, E. J. Wahyudi, W. G. Kadir, S. Alawiyah, R. A. Wicaksono, and Y. I. Kusumah, “Case study of time-lapse microgravity data on the Wayang Windu geothermal field, Indonesia: quality control of data acquisition and data application for reservoir analysis,” in 10th SEGJ International Symposium, pp. 465–470, Society of Exploration Geophysicists of Japan, 2011.
View at: Google Scholar
Y. I. Kusumah and H. H. Suryantini, “Horizontal derivative from gravity data as a tool for drilling target guide in Wayang-Windu geothermal field, Indonesia,” in Proceedings World Geothermal Congress, Bandung, 2010.
View at: Google Scholar
R. G. Intani, S. Wicaksono, and G. Golla, “Updated geologic structures and stratigraphy of the Darajat geothermal field in Indonesia,” IOP Conference Series: Earth and Environmental Science, vol. 732, article 012012, 2021.
View at: Google Scholar
F. Albino, J. Biggs, and D. K. Syahbana, “Dyke intrusion between neighbouring arc volcanoes responsible for 2017 pre- eruptive seismic swarm at Agung,” Nature Communications, vol. 10, no. 1, p. 748, 2019.
View at: Publisher Site | Google Scholar
S. Gambino and G. Distefano, “Intrusive seismic swarms as possible precursors of destructive earthquakes on Mt. Etna’s eastern flank,” International Journal of Geophysics, vol. 2022, Article ID 8565536, 10 pages, 2022.
View at: Publisher Site | Google Scholar
W. Huang, W. Gao, and G. Ding, “Neogene volcanism and Holocene earthquakes in the Tanlu fault zone, eastern China,” Tectonics, vol. 260, no. 4, pp. 259–270, 1996.
View at: Publisher Site | Google Scholar
S. Momeni and R. Madariaga, “Long-term triggered seismicity on the Mosha Fault by Damavand volcano, Iran: implications on the seismic hazard of Tehran metropolis,” Frontiers in Earth Science, vol. 10, article 945297, 2022.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2023 Ilham Arisbaya et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1031

Downloads

825

Citations