Integration of Audio-Frequency Domain Magnetotelluric and Global Gravity Model Plus for Delineating Possible Karst Aquifer Zones of Southern Malang, Indonesia

2.1 Geological conditions

Southern Malang Regency is part of the Southern Mountain range of Java Island, Indonesia, which is dominated by karst rocks. Geologically, karst areas are dynamic and highly complex. Such zones are dominated by limestone, dolomite, and gypsum, which are relatively soluble in water and therefore can undergo weathering relatively easily [20]. This leads to the formation of geological features such as faults, sinkholes, cavities, and subsurface rivers. Eventually, these cavities and fractures serve as water storage media [21], forming subsurface flow systems. As a result, the surface of karst areas may appear dry, although they may store groundwater reserves at certain depths.

tu_pian_1.jpg

Figure 1. Geological map of the study area [22]

The research site is located in Sumberbening Village, Bantur District, Malang Regency, East Java Province, covering an area of approximately 2 km × 2.7 km. Based on the Turen geological map sheet (Figure 1), the study area belongs to the Wonosari Formation (Tmwl) and Nampol Formation (Tmn) [22]. The Wonosari Formation belongs to sedimentary rock unit composed of limestone, sandy marl, and claystone intercalations. While the Nampol Formation consists of tuffaceous or calcareous sandstone, black claystone, and sandy marl [23]. This combination of carbonate lithology and fine clastic material significantly influences the hydrogeological characteristics of the area, particularly by controlling porosity, permeability, groundwater storage mechanisms, and subsurface groundwater movement patterns within the karst aquifer system and valley deposits. Geomorphologically, the study area is dominated by hills and valleys, with relatively complex topography. Hilly zones generally develop on carbonate rocks that have undergone karstification, characterized by intensive limestone dissolution that produces fractures, cavities, and underground channels that serve as groundwater conduits. Meanwhile, the valley areas tend to be filled with deposits derived from weathering and the transport of sediment from higher areas.

2.2 Audio-Frequency Domain Magnetotelluric method

ADMT can provide high-resolution imaging of conductivity variations at different depths, enabling detailed analysis of complex subsurface structures [24]. ADMT can demonstrate strong performance for imaging electrical resistivity in near-surface aquifer systems by exploiting natural electromagnetic noise over a bandwidth of 1 Hz to 10 kHz. Because sand or soil saturated with groundwater can conduct electricity much better than dry sand or soil, ADMT surveys can produce continuous resistivity-depth profiles that correlate directly with water content.

Table 1. The parameters of Audio-Frequency Domain Magnetotelluric (ADMT)-3HT3 [25]

Data acquisition utilized the ADMT-3HT3 instrument, operating within a frequency range of 1-8000 Hz. The measurement duration was maintained at 420s per point to ensure stable signal-to-noise ratios. Before inversion modeling, rigorous data-quality control was executed to remove statistical outliers caused by cultural noise. The data processing was conducted using the instrument's proprietary software, which utilizes a fast 1D pseudo-inversion (Bostick transformation) to convert frequency data directly into apparent resistivity-depth profiles. Because this algorithm relies on direct depth-transformation rather than iterative 2D forward modeling, a traditional iterative root mean square (RMS) error value is not generated. However, data-quality control was rigorously maintained during acquisition by monitoring signal stability and removing isolated outlier measurements before the contouring process. Table 1 summarizes the acquisition and processing parameters for the ADMT method.

tu_pian_2.jpg

Figure 2. Audio-Frequency Domain Magnetotelluric (ADMT) measurement survey design

A total of 48 measurement points were distributed throughout the study area (Figure 2). The measurement points were designed to account for geomorphological variations, geological conditions, and field accessibility, enabling them to capture changes in subsurface characteristics across the different landform units in the study area. These randomly distributed points were projected onto three main 2D sections (Line L1, Line L2, and Line L3) using a nearest-neighbor interpolation method with a maximum projection distance of 150 meters. In addition to the geophysical measurements, field data from an existing local production well (90 meters deep) within the study area were incorporated to provide control and ground-truth validation for the survey design. Coordinate details for ADMT data points are provided in Table 2.

Table 2. Audio-Frequency Domain Magnetotelluric (ADMT) point coordinates

The 2D sections allow the identification of vertical and lateral changes in resistivity values; a good 2D model can be used to identify deep or shallow objects [26]. That model can also be used to interpret the presence of water-bearing layers, fracture zones, and boundaries between lithological units. The obtained ADMT data were then processed and modeled in 3D to provide a more comprehensive representation of aquifer distribution and groundwater flow patterns. The 3D model was generated by integrating all ADMT measurement points, enabling the spatial distribution of groundwater to be displayed within a broader subsurface volume. This approach provides advantages in understanding the relationships among geological conditions, subsurface structures, and aquifer distribution.

2.3 Global Gravity Model Plus

Global Gravity Model Plus (GGMplus) is a recent composite gravity data model with an ultra-high spatial resolution of 7.5 arc-seconds, or approximately 200 meters [27]. This sub-local spatial resolution can theoretically be developed through the integration of low-frequency satellite gravity data (such as GRACE and GOCE), the EGM2008 Earth gravity model, and modeling of the short-wavelength topographic gravity effects from Shuttle Radar Topography Mission (SRTM) data [27]. For this study, the data were acquired by downloading an anomaly-value grid through open access from the official Curtin University gravimetric data repository [28].

This composite method makes GGMplus data highly sensitive to the detection of regional anomalies [28]. In this integrated research framework, the GGMPlus satellite data functions strictly as a preliminary survey or reconnaissance tool. Due to its 200 m spatial resolution, it is utilized to delineate broad, macroscopic low-density zones that indicate regional intensive secondary porosity or local faults that act as main controls on underground river systems. The result of this regional gravity assessment is then used to guide the higher-resolution ADMT measurements for detailed, localized aquifer targeting. In addition, this satellite-based method does not require large-scale topographic land surveys in karst areas that are notoriously difficult to access [29]. The GGMPlus data used in this study consisted of 40 points spaced 200 meters apart, as shown in Figure 3.

tu_pian_3.jpg

Figure 3. Global Gravity Model Plus (GGMplus) survey design

In the data processing stage, the output obtained from the GGMPlus repository was still in the form of Free Air Anomaly (FAA) values [27]. To accurately interpret the karst aquifer system, these anomaly values were transformed into observed surface gravity values (Gobs) by extracting topographic elevation from a high-resolution 8 m Digital Elevation Model Nasional (DEMNAS). The theoretical gravity (Gtheo) was then applied based on latitude and station elevation correction [28]. The next step was to calculate the Complete Bouguer Anomaly (CBA) by applying Bouguer correction and terrain correction (using a terrain-correction radius of 5 km and a standard reduction density of 2.50 g/cm³, which is essential because karst morphology generally has rough and undulating topography [28].

The processed CBA data were then deconvolved with a Butterworth filter to focus on the target aquifer system and to remove basement-rock effects. A Butterworth pass-band filter (order 8, cutoff wavelength of 560 m) was selected because it gradually blocks signals outside the frequency range without causing a ringing effect in the visualization results [30]. The CBA values were then inverted into a map using Oasis Montaj software. Further processing focused on separating the CBA anomaly into a residual anomaly to isolate target anomalies indicating potential groundwater storage zones. To present the data more effectively, 2D cross-section modeling was then performed using ZondGM2D software. Through this complete sequence of topographic correction processing, lateral density-contrast variations in the subsurface can be accurately identified, allowing massive and solid limestone matrices to be distinguished from permeable cavity zones (filled with either water or air) above them as the main object of analysis [29].

4.1 Characterization of the karst conduit system through density contrast

Validation of the spatial dimensions and depth of dissolution zones was carried out using the 2-D inversion method on cross-section lines with ZondGM2D software. The cross-section model of the residual anomaly is generally shown in Figure 8. The model is presented as a 2-D vertical density-contrast section along four explicitly marked cross-section lines: A-A' and B-B' (South-North orientation), as well as C-C' and D-D' (West-East orientation), with a maximum observable depth of 80-85 meters below the ground surface (indicated by the depth scale in meters on the y-axis of each section). The inversion incorporated structural constraints based on smoothness criteria and depth weighting, achieving specific fitting errors ranging from 3.5% to 5% for the profiles. ZondGM2D software iteratively simulates the density distribution so that the surface observation data (residual anomaly) approaches the theoretically calculated anomaly curve. In the inversion process, the color range in the vertical section represents the density contrast of the rocks. The darker the blue color, the greater the negative density contrast, with values between -0.2 and -0.4 g/cm³.

tu_pian_8.jpg

Figure 8. Residual anomaly cross-section

It is important to acknowledge that its 200 m spatial resolution limits the ability to map small-scale individual karst conduits or localized dissolution voids directly. Therefore, this mass-deficiency zone is not interpreted as specific, isolated active conduits, but rather as macroscopic low-density karst features or bulk dissolution zones that function as a subsurface aquifer system [34]. The zones interpreted as potential macroscopic aquifers are marked by the red circles in Figure 8. In this study, these zones are defined as broader areas of low density, indicative of heavily dissolved limestone rather than distinct physical cavities. The pink-to-yellow colors indicate positive density anomalies, with values of 0.1 to 0.4 g/cm³, which are interpreted as massive, solid limestone matrices that have not undergone extensive dissolution. This integrative modeling result (Figure 8) shows that GGMPlus satellite data are sufficiently representative for mapping limestone dissolution zones and subsurface aquifer potential in a karst system [29].

4.2 Aquifer system heterogeneity

The most important finding was obtained by comparing the conceptual cavity model derived from GGMPlus data with the resistivity distribution from 3D ADMT modeling. The overlay of the 3D aquifer distribution model and Google Earth satellite imagery is shown in Figure 9. This overlay map clearly displays the spatial extent of the target aquifer system relative to the surface conditions in Sumberbening Village. This correlation provides a basic understanding of the hydrogeological conditions in Sumberbening Village. Groundwater conditions in this area cannot be assumed to be uniformly distributed, as in typical porous rock aquifers. Instead, karstification of the limestone formation over thousands of years has produced scattered and disconnected macroscopic dissolution pathways. Water entering the subsurface infiltrates and flows into these heterogeneous passages. This is highly consistent with a conduit-type aquifer system, which is characterized by very rapid, dynamic, and spatially nonuniform fluid flow [33]. Therefore, the overlay map in Figure 9 helps provide an understanding of subsurface connectivity associated with linear features and supports the determination of potential groundwater storage locations in this area.

tu_pian_9.jpg

Figure 9. Aquifer model and google earth overlay

The highly heterogeneous spatial macroscopic dissolution characteristics, as comprehensively visualized in the overlay map of the 3D model and satellite imagery (Figure 9), emphasize the advantage of the integrated approach proposed in this study. In hydrogeophysical studies in areas with karst morphology, Abir et al. [35] and Jacob et al. [36] concluded that conventional exploration using Electrical Resistivity Tomography (ERT) and land microgravity surveys is highly limited due to narrow area coverage, high cost, and access constraints in hilly and extreme terrain. In contrast, the use of GGMPlus data in this study required only a day for data acquisition and processing. Furthermore, the data were obtained at no cost, as GGMPlus provides open-access data. Although the accuracy of satellite data is not equivalent to direct field measurements, we use the ADMT method to validate the GGMPlus results.

As an efficient alternative, this study demonstrates that GGMPlus satellite data can be used for preliminary screening, and ADMT for field-scale follow-up. The very low residual anomaly values from GGMPlus (-0.2 to -0.4 g/cm³) show connectivity with the aquifer model generated from ADMT. The combination of these two methods can more efficiently visualize heterogeneous groundwater spatial domains.

4.3 Strategic drilling recommendations

Fulfilling the primary objective of this hydrogeophysical investigation, the integration of ADMT and GGMPlus data facilitates the delineation of strategic bore-well targets. While the ADMT pseudo-sections indicate broad potential aquifer zones at depths of 25–40 m for shallow aquifer and 70–100 m for recommended aquifer with details as in Figure 10, identifying optimal drilling points requires spatial convergence with the macroscopic dissolution zones derived from the residual gravity anomalies.

To explicitly demonstrate this spatial convergence, the 2D ADMT measurement lines were directly overlaid onto the GGMPlus residual anomaly map, as illustrated in Figure 10. This integrated visualization reveals that the very low residual anomaly values from GGMPlus (-0.2 to -0.4 g/cm³) spatially intersect with the deep low-resistivity zones identified by ADMT. This combination precisely pinpoints the optimal locations where massive secondary porosity coincides with high fluid saturation.

tu_pian_10.jpg

Figure 10. Audio-Frequency Domain Magnetotelluric (ADMT) and gravity Global Gravity Model Plus (GGMPlus) overlay

Table 3. Outlines recommended drilling locations

Table 3 outlines the recommended drilling locations aimed at mitigating the clean water shortage in Sumberbening Village. These strategic points were selected based on the strongest spatial overlap between low resistivity values (indicating fluid saturation) and significant mass deficiencies (indicating high secondary porosity). Due to the inherent non-uniqueness of geophysical data and the absence of prior local well-logging for calibration, an uncertainty assessment is provided for each point. These targets represent geophysical predictions that require careful validation during the initial drilling phase.