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

In recent years


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 (V S ) 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 2 [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.

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].

International Journal of Geophysics
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].

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 (g ϕ ) was calculated using the 1967 international gravity formula and was adjusted with free air (g FA ), the Bouguer (g BC ), and terrain (g TC ) corrections.The Bouguer anomaly (g BA ) was calculated by using the following well-known equation [36]: where h is the elevation of the gravity station and ρ is the average rock density of 2.67 g/cm 3 used to calculate the Bouguer correction (g BC ).The terrain correction (g TC ) 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].

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 (E X and E Y ) and magnetic (H X , H Y , and H Z ) 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 Z relating the horizontal components of the electric field E to magnetic field H 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 Z contains information about the earth's resistivity at a given location as function of depth.As a complex number, each component of Z = ½Z i,j consists of real and imaginary parts and thus can be expressed in terms of apparent resistivity ρ a related to the amplitude jZj and phase φ as follows [38]: 3 International Journal of Geophysics where ω is the angular frequency, μ 0 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 4 Hz down to 10 -1 Hz at all stations, with only 11 stations among them that were complemented with lower-frequency (10 2 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  [13]).The surface geology was from [19,20], while seismicity was compiled from [1,34].International Journal of Geophysics 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 −3 °< β < 3 °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 −3 °< β < 3 °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 com-plex 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.
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 Figure 3: The Bouguer anomaly map of southern Garut and its surroundings from [35] covering a larger area beyond the interest area (dashed blue outline) to estimate the regional anomaly component.Active faults (full red line) and possible active faults (dashed red line) are from [3].Black and red triangles represent dormant and active volcanoes, respectively.International Journal of Geophysics 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 (Z XX and Z YY ) and maximize the off-diagonal components (Z XY and Z YX ) in the twodimensional (2D) modeling.Rotation to the geoelectric strike results in an apparent resistivity component obtained with the electric field parallel to the geoelectric strike (Z XY ) as the transverse electric (TE) mode while the other one obtained with the magnetic field parallel to the geoelectric strike (Z YX ) as the transverse magnetic (TM) mode [38].

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).

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    International Journal of Geophysics 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 cal-culated 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.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.

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.
C' C   13 International Journal of Geophysics 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 highresistivity 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   14 International Journal of Geophysics 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 Eq. 2017    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 lowresistivity 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 highlow 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.

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.

Figure 1 :
Figure 1: Location of the study area in southern Garut, West Java, Indonesia, indicated by the blue box, along with the relocated shallow earthquakes from 2009 to 2015 [1], destructive earthquakes represented by focal mechanisms [2], and active fault traces [3].

Figure 4 :Figure 5 :
Figure4: Dimensionality analysis based on skew angle (β) from the phase tensor analysis[39] for MT data in line C-C′, D-D′, and E-E′ (see Figure2for MT profile lines and sounding sites).

Figure 6 :
Figure 6: The value of regularization parameter (tau) was obtained from L-curve analysis for MT data line C-C′ (left), D-D′ (middle), and E-E′ (right).

Figure 8 :
Figure 8: The results of 3D gravity inversion modeling, measured data (a), calculated data or inverse model response (b), misfit (c), and density block model sliced at A-A′ and B-B′ cross-sections (d).

9
International Journal of Geophysics 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 threedimensional 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.

Figure 9 :
Figure 9: Apparent resistivity and phase pseudosection of MT model for line C-C′ (northernmost line); observed (obs) and calculated (cal).

Figure 10 :
Figure 10: Apparent resistivity and phase pseudosection of MT model for line D-D′ (central line); observed (obs) and calculated (cal).

Figure 11 :
Figure 11: Apparent resistivity and phase pseudosection of MT model for line E-E′ (southernmost line); observed (obs) and calculated (cal).

Figure 12 :
Figure 12: Resistivity model from 2D MT inversion of line C-C′ (a), D-D′ (b), and E-E′ (c).R1, R2, C1, and C2 are resistive and conductive bodies tested in the sensitivity test (see text for details).

Figure 14 :
Figure 14: Interpretation of the density (top) and resistivity (bottom) models of line C-C′.Interesting feature is observed at Mt. Malang, which is characterized by high-density and low-resistivity values.

Figure 16 :
Figure 16: Interpretation of the density (top) and resistivity (bottom) models of line E-E′.The resistivity model shows a basement discontinuity in the western part of the cross-section and a conductive body beneath Mt.Kencana.MEQ (grey shaded area) cluster observed 4 km south Mt.Kencana [45].

Figure 17 :
Figure 17: Basement depth map based on the 3D density model using a threshold of 2.525 g/cm 3 as the minimum basement value.