This study uses elevation data from the Indonesia Digital Elevation Model (DEMNAS), as shown in Figure 2a, as well as complete Bouguer gravity anomaly data (see Figure 2b). The gravity data we utilize is a combination of field data we measured and data from the Indonesian Geological Resources Centre (PSDG). The total number of data points is 344. A total of 216 data points have regular grids with intervals of 250 m spanning from west to east, while the remaining 128 data points have intervals ranging from 250 m to 550 m. Hot springs, represented by red squares, are located on the hillside.

We calculate the complete Bouguer anomaly (CBA) data by applying several corrections to the observed gravity data. This correction begins with tidal correction, instrument height correction, latitude correction, free-air correction, Bouguer correction, and terrain correction. The last two corrections are very dependent on the Bouguer density, namely the density of the rock layer between the mean sea level (MSL) and the measurement points. The Bouguer density used was 2.6281 g/cm³. The estimated Bouguer density was calculated using the Parasnis method, which is represented by the gradient of the elevation scale (0.04192 h) versus the free-air gravity anomaly (see Figure 2c). This estimated density is comparable to those presented by the PSDG (2010), which used an average density of 2.65 g/cm³ for field rock samples. Figure 2b shows the CBA contour, which varies in the range of 39 mGal to 55 mGal. The distribution of CBA forms a trend in the northwest-southeast and northeast-southwest directions. Geothermal manifestation points are located around high gravity anomalies with a value range of 50 mGal to 55 mGal. The presence of high-density rocks or the position of the basement near the surface may be associated with high gravity anomalies.

This step involves the process of separating residual and regional anomalies. For this purpose, we use the upward continuation technique on CBA data. This technique applies a mathematical filter that amplifies data at low frequencies associated with regional anomalies. In contrast, this upward continuation technique weakens high-frequency data associated with residual anomalies (Chisenga et al., 2019; Kebede et al., 2020). Therefore, this technique emphasizes regional scale anomalies.

Figure 3 shows the upward continuation technique. The gravitational attraction per unit mass, dm, of an anomalous body at mean sea level, B(x0, y0, z0 = 0), and a distance R from the source point, B(x,y, z), may be determined using the following equation (Jacobsen, 1987; Kebede et al., 2020):

d g = Gdm R 2 = Gdm ( x − x 0 ) 2 + ( y − y 0 ) 2 + ( z − z 0 ) 2 (1)

The source's gravity anomaly, Δg0, at the mean sea level point, is given by the vertical component, dgz.

d g z = Δ g 0 = Gdm R 2 z R = G d m z R 3 = G d m z [ ( x − x 0 ) 2 + ( y − y 0 ) 2 + ( z ) 2 ] 3 / 2 (2)

At the upward continued surface, point A, the source's gravity anomaly, dgP, equals

d g p = Gdm S 2 z + h S = G d m z + h S 3 = G d m z + h [ ( x − x 0 ) 2 + ( y − y 0 ) 2 + ( z + h ) 2 ] 3 / 2 (3)

Several recent edge detection techniques are often used to analyze the Bouguer gravity map and refine the edges/boundaries of the subsurface causative bodies. Numerous mathematical strategies are available for enhancing the accuracy of gravity gradient data in edge detection. As a result, distinct lineaments appear at slightly different spatial coordinates. This study used two edge detection approaches, which will be discussed below.

According to Ferreira et al. (2013), the TAHG (tilt angle of the gradient amplitude) is the arctangent of the gradient amplitude's derivative ratio. The TAHG filter highlights both shallow and deep geological contacts that cause gravity anomalies at the surface (Pham et al., 2023; Safani et al., 2023). It can be expressed as follows:

TAHG = atan ( ∂HG ∂z ( ∂HG ∂x ) 2 + ( ∂HG ∂y ) 2 ) (4)

Where ∂HG/∂x and ∂HG/∂y are the first-order horizontal derivatives of the horizontal gradient magnitude (HG), and ∂HG/∂z the first-order vertical derivative.

Recently, the Euler deconvolution (ED) methodology has emerged as a highly effective method for semi-quantitative analysis of potential field data (Hinze et al., 2013; Chen et al., 2022; Safani et al., 2023). The ED approach has been applied to gravity data to enhance the characterization of geological structures and ascertain their depth. Marson and Klingele (1993) devised the Euler deconvolution for vertical gravity data to enable this method to function with grid gravity data sets, as illustrated in the equation below:

( x - x 0 ) ∂ 2 G ∂z∂x + ( y - y 0 ) ∂ 2 G ∂y ( z - z 0 ) ∂ 2 G ∂ z 2 = N ∂ G ∂ z (8)

Where (x₀, y₀, z₀) denotes the coordinates of the anomaly's source in the (x, y, z) axis. G represents the gravitational force at the coordinates (x, y, z). N represents the degree of homogeneity as a structural indicator (SI). The index structure delineates the geological conditions that produce the anomaly (Marson and Klingele, 1993; Pham et al., 2023). ∂G/∂x and ∂G/∂y represent horizontal derivatives in the x and y planes, respectively, whereas ∂G/∂z denotes a vertical derivative.

The complete Bouguer anomaly (CBA) combines both local and regional anomalies. We use the upward continuation technique to separate the two anomalies. The continuation process involves several elevation levels, starting from an elevation of 100 m to 700 m above mean sea level (MSL) to obtain representative regional anomalies. We chose a continuation result, 500 m above MSL, to represent the regional gravity anomaly (see Figure 4a). This is because the contour patterns in these results tend to be stable, and the influence of local anomalies tends to diminish. Subtracting the regional gravity anomaly from the complete Bouguer anomaly produces a local Bouguer anomaly with varying values ranging from -7.7 mGal to 8.8 mGal (see Figure 4b).

The study aims to estimate the three-dimensional structure of the subsurface layers of the Lainea geothermal region, utilizing gravity data and the GROWTH 3.0 tool developed by Camacho et al. (2021). This package offers both default values and user-specified options. This coding package utilizes all default values, except for a few specific ones. Through trial and error, we selected a balance factor of 800 rather than the package's default of 80, as an increased balance factor results in a greater percentage of cell filling in the subsurface layer. We used a flattening coefficient of 10 and a lateral smoothing value of 0.3 (selected after trying various options) to enhance cell filling stability.

The inversion results of the residual gravity data describe the subsurface layer in two ways: horizontally (see Figure 5) and vertically (see Figure 6).

In this section, we discuss the geothermal system at the LNVGF, which includes reservoir basins, basements, heat sources, and subsurface faults that control geothermal activity. For this purpose, we created two geological models (namely, Figures 7a and 7b), utilizing the gravity inversion model from Figures 6a and 6c. The choice of these two inversion models is considered appropriate because they represent two different cutting directions and cross the area where hot water manifestations occur. This model improves the interpretation models of PSDG (2010), Sugianto (2011), and Zarkasyi et al. (2014), especially in describing the patterns and geometry of basins, basements, and heat sources.

The two geological models show the underlying layout in the research area, which determines the mechanisms of the geothermal system. We interpreted six rock units in the LNVGF based on rock density values. The first rock units have a density range of 2608 ≤ ρ < 2762 kg/m³. These density values are the sum of density contrasts (-20 ≤ ∆ρ < 134 kg/m³) and Bouguer density (2628 kg/m³). The second rock units have a density range of 2762 ≤ ρ < 2865 kg/m³, after summation of the Bouguer density and the density contrast 134 ≤ ∆ρ < 237 kg/m³. The density range of the third rock units is 2865 ≤ ρ < 2942 kg/m³ (with the range of density contrast of 237–314 kg/m³). The fourth rock units have a density range of 2948 ≤ ρ < 3045 kg/m³ (the density contrast in this stratum ranges from 315 to 417 kg/m³). The fifth rocks have a density of 3046 ≤ ρ ≤ 3122 kg/m³ and a density contrast 418 ≤ ∆ρ < 494 kg/m³. The sixth rock unit, which is the bottom layer, has a density of ρ = 3174 kg/m³ or a density contrast 546 kg/m³, indicating the presence of an intrusive rock shape.

Figure 7a shows subsurface basin formed at depths ranging from 500 m to about 1200 m represented by the rock units with the density range of 2865 ≤ ρ < 2942 kg/m³. This basin appears deeper on the West and East sides of the cross-section than the one in the middle. The maximum thickness of this layer is about 300 m. The formation of the basin and its flanking faults is associated with metamorphic exhumation that began in the Middle Miocene and continued into the Pliocene (Mawaleda et al., 2018; Muzani et al., 2023). Furthermore, the development of the basin and its associated normal and strike-slip faults was predominantly driven by persistent extension caused by Banda rollback (Nugraha and Hall, 2022). This extension caused crustal thinning and magmatic activity, creating an environment conducive to pull-apart basin formation. These basins created optimal conditions for geothermal activity and hot spring formation. Based on these characteristics, the central basin zone likely acts as a reservoir basin in the LNVGF. The central basin zone has depths ranging from approximately 600 to 900 meters. Using magnetotelluric and time domain electromagnetic methods, Zarkasyi et al. (2014) estimated the existence of a reservoir at a depth of 750-1000 m, which is comparable to the reservoir estimate from this study. This reservoir layer is likely composed of metamorphic rocks (PSDG, 2010; Zarkasyi et al., 2014) with densities ranging from 2865 kg/m³ to 2942 kg/m³.

Above the reservoir basin there are several rock blocks with various densities that act as cap rock. The cap rock layer is located at a depth of several tens of meters above sea level (asl) to a depth of about 1000 m bsl (see Figures 7a and 7b). PSDG (2010) and Zarkasyi et al. (2014) suspect that the cap rock comes from altered metamorphic/sedimentary/limestone rocks. Our study also suspects the potential of the altered metamorphic, meta-sandstone and meta-limestone rocks to function as cap rock in the LNVGF. Rocks with density ranges of 2608 ≤ ρ < 2762 kg/m³ (interpreted as meta-limestone) and 2762 ≤ ρ < 2865 kg/m³ (interpreted as meta-sandstone), as well as altered metamorphic rocks (with densities equivalent to meta-sandstone), located in the western to southern zones of the Boroboro Fault, act as caprock. In the northern to eastern zones, a combination of altered metamorphic rocks and meta-limestone acts as caprock.

Under the reservoir layer, there are two high-density rock units with the density range of 2948 ≤ ρ < 3045 kg/m³ and 3046 ≤ ρ ≤ 3122 kg/m³, respectively. These rock outcrops appear on the surface and represent metamorphic rocks (see Figure 1b). The difference in density between these two metamorphic rocks may be due to the different pressures and mineral deposits they experienced during their alteration. The higher density of the metamorphic rocks at depths of 2-4 km is considered to result from a combination of intrusive bodies located directly below it and mineral deposition from thermal water in that depth zone. This condition resembles the geothermal field depicted in the Imperial Valley, California (Gupta and Roy, 2010). These two metamorphic layers located on the left and right sides of the cross-section can be thousands of meters thick. However, the exhumation and extension during Miocene-Pliocene unconformity (Nugraha et al., 2022) cause the metamorphic rocks thinning in the center of the research area.

The interpretation of the faults in Figure 7 is based on the geometric pattern of the resulting density model (see Figures 6a and 6c), as well as the geological information in the area (see Figure 1b). The interpretation of the Boroboro Fault as a normal fault in Figures 7a and 7b aligns with the geological data (see Figure 1b), which indicates that the northwest-southeast trending Boroboro Fault is a normal fault. Figure 7a shows the location of the Boroboro fault on the west side of the M3 hot spring manifestation. However, this visual cannot explain the western Boroboro fault as a conduit for the M3 hot spring. The Boroboro Fault also lies beneath the M7 and M8 manifestations, and it is suspected that this fault acts as a conduit for both. The faults beneath the M5 and M6 manifestations may be part of the Boroboro Fault system due to their similar fault patterns and proximity to the hot spring manifestations. The Boroboro Fault and the Landai Fault are suspected to act as conduits for the M5, M6, M7, and M8 hot springs due to the confluence of the two faults around the four manifestations. The fault beneath the M3 hot spring in Figure 7a is suspected to be associated with the Boroboro Fault but requires confirmation with additional evidence. Meanwhile, the suspected fault located east of the M7 and M8 hot springs is based solely on gravity model interpretation; there is no geological evidence to confirm it. These estimated faults are located at depths ranging from a few hundred meters to about 2000 m bsl. The Boroboro Fault extends into the upper basement at a depth of approximately 2000 m bsl.

Figure 7a also demonstrates the existence of intrusive/plutonic rocks at depths ranging from 1500 m to more than 3000 m. The form is quite vast, widening from west to east, with the peak located in the center of the track. This intrusive rock contains branches that produce two dome summits. This rock is the result of magmatic activity after metamorphic exhumation and crustal thinning resulting from extension driven by the Banda rollback in the eastern part of the SE arm (Nugraha and Hall, 2022). The intrusive plutonic rock serves as a heat source in the Lainea geothermal system. In line with this study, Zarkasyi and Widodo (2014) demonstrated the presence of the plutonic rock model at the depth of 1500 m using a magnetotelluric method. The gravity-based inversion model in this study has the advantage of presenting a model of underlying layered strata based on rock density contrast, which makes interpretation easier.

Fault systems play a crucial role in geothermal systems, serving as permeable conduits for the ascent of hot fluids or steam from the reservoir to the surface, as well as facilitating heat transfer from thermal sources into shallower strata. Advanced gravity data processing approaches, including tilt angle horizontal gradient (TAHG), fast sigmoid edge detection (FSED), and Euler deconvolution, demonstrate dependability in estimating subsurface faults (Safani et al., 2023).

Figure 8 illustrates the application of TAHG and FSED techniques to residual gravity anomaly data. The TAHG method shows a very clear contrast between the strong (red) and weak (blue) amplitudes of the gravity anomaly (see Figure 8a). Meanwhile, the FSED method emphasizes positive amplitudes by degrading negative amplitudes (see Figure 8b). Both edge detection methods exhibit analogous positive amplitude patterns. The highest amplitude reveals the edges of anomalous sources that appear to depict subsurface fault characteristics (thick black lines) inside the research area. This is because the alignment of the identified subsurface faults exhibits two directional trends, northwest-southeast and southwest-northeast, which typically correlate with the pattern of surface faults observed using the geological approach (yellow). The orientation of the Boroboro Fault, identified as a normal fault, indicates that its northwestern and southeastern ends were once aligned. Additionally, the Windo and Kaendi faults, which are strike-slip faults, consistently bend in a southwest-northeast orientation. This type of fault aligns with the geological information presented in Figure 1b. The core of these evolving faults aggregates all geothermal manifestations. The activity of these faults and the fractures they generate significantly influence the migration of hot fluids to the surface. Further discussion of the identified subsurface faults using other geophysical approaches is needed. There is very minimal evidence of earthquake activity around the faults during the period 1972-2024. Only two earthquakes (yellow asterisks in Figure 9) are recorded throughout that period.

Euler deconvolution (available in Oasis Montaj software) was applied to the gravity anomaly data in order to estimate the depth of the anomaly source. The solution was found through trial and error while considering the deconvolution window size. After trying various options, the optimal deconvolution window size for estimating the depth of anomalous sources was found to be 10 × 10. This study used a structural index (N) of 0 to represent the geological contact response formed in the study area. This index is suitable for estimating geological structures, such as faults (Pham et al., 2023; Safani et al., 2023). Figure 9 depicts the solution of Euler Deconvolution represented by dots with colors varying from dark blue to dark red. This figure illustrates depth variations of subsurface geological structures ranging from 168.6 to 1998.4 m. These depth variations align with the fault depth interpretation shown in Figure 7. The hot springs in the central part of the study area are surrounded by shallow geological structures (blue dots), while deeper structures (light pink to dark red dots) surround the hot springs at more distant positions. This pattern of structure depths also illustrates the uplift pattern that occurred in the area. Figure 9 shows the locations of geological structures associated with strike-slip fault lineaments in the study area. These include the Windo Fault, Kaendi Fault, Lainea Fault, and Damba Fault. These strike-slip faults have shallow depths, as indicated by the blue dots.

The solution of the Euler deconvolution reveals subsurface anomaly direction trends that are analogous to those estimated with TAHG and FSED. This good agreement may be due to selecting an appropriate index structure (N = 0) to describe the geological contact in the form of a fault. Figure 8 and Figure 9 illustrate the northwest-southeast (NW-SE) oriented Boroboro Normal Fault, which developed during the extensional process, followed by crustal thinning. Strike-slip faults trending NE-SW cut the Boroboro Fault relatively shallowly after crustal thinning. These strike-slip faults formed as secondary faults after the main fault developed during the Miocene-Pliocene deformation phase. The rose diagram of the resultant data constructed from the FSED’s predicted faults indicates that the predominant orientation of the geological structures is southwest-northeast, succeeded by northwest- southeast. The presence of two earthquake occurrences (yellow asterisks) in the center of Figure 9 indicates the dynamic nature of the geological features in the studied area.