Wai Selabung geothermal reservoir analysis based on gravity method

Received: August 11, 2021 Accepted: October 16, 2021 Published: October 30, 2021 Research has been conducted using the gravity method in the Wai Selabung area, South Ogan Kemiring Ulu Regency, South Sumatra Province, correlated with geological data, magnetotellurics, and geochemical data. This research aims to get structural patterns, subsurface models and identify the heat source and reservoir areas of the Wai Selabung geothermal system. This study uses the gravity method to model the subsurface, which is correlated with magnetotelluric and geochemical data to identify reservoir prospect areas. The results obtained from this research include residual anomalies in the research area showing the presence of a northwest-southeast trending fault structure by the main fault structure of this area trending northwest-southeast and slightly southwest-northeast. Analysis of the Second Vertical Derivative value of zero indicates the boundaries of the geothermal reservoir in the middle of the research area. The results of the 3D inversion modeling of the research area show that low density (2 to 2.15 g/cm) indicates the location of the reservoir, medium-density values (2.2 to 2.4 g/cm) are tertiary sandstone sedimentary. The high-density distribution value (2.5 to 2.9 g/cm) indicates a potential heat source. And based on the analysis of the gravity method correlated with geological data, magnetotelluric, and geochemical data, the prospect area for the Wai Selabung geothermal reservoir, is around Teluk Agung, Perekan, and Talang Tebat.


INTRODUCTION
Indonesia has the potential for renewable energy, namely the most significant geothermal energy in the world. Indonesia's geothermal potential is 29,038 MW, and only 1,226 MW has been utilized or 4.2%. The National Energy Policy has targeted geothermal to support 5% of the national energy mix by 2025. However, geothermal has only contributed approximately 4% with slow development, so an investigation is needed to develop geothermal fields in Indonesia.
The Wai Selabung area is located on the Sumatran fault zone, where many hightemperature geothermal locations are found. Manifestations were found at the site in the form of hot springs on the banks of the Wai Selabung River, in a tributary of the Wai Selabung river, and alterations that appeared near Suban hot springs (banks of the Wai Selabung river) .
Several studies have been conducted in this area, whether related to geothermal activity or not, including Nasution et al. (1999), who made regional mapping using the Baturaja sheet system. Research related to geothermal activity has been carried out by the Center for Geological Resources from 2011 through an integrated survey of the Wai Selabung Geothermal Area . The last is the magneto-telluric (MT) Survey conducted in 2020 (Ilmi et al., 2020). The combined survey results show a potential geothermal resource of 68 MWe in this area, and the prospect area is in the southwest. Previous research revealed that the prospect area is only in the southwest, so further research is needed to identify the Wai Selabung geothermal area's structural pattern, heat source, and Wai Selabung geothermal reservoir (Hadi et al., 2011). Gravity surveys were conducted for 210 stations.
The gravity method is used to identify and describe subsurface geology based on variations in the earth's gravity field caused by differences in density between rocks. In the case of geothermal, the difference in rock density is a reference for gravity research. The heat source area and its accumulation below the earth's surface can cause density differences with the surrounding rock mass. Geological structures can form permeable zones that can escape to the surface so that manifestations appear on the surface. This study aims at a complete 3D inversion of the Bouguer anomaly, gradient analysis for structural analysis, and to determine the geochemical characteristics of the geothermal Wai Selabung area. Bouguer anomaly inversion modeling correlated with MT, geological, and geochemical data.

Geology
The structure of the study area is included in the volcanic environment of the magmatic pathway in southern Sumatra. The formation of the geothermal system in the Wai Selabung region is influenced by volcanic and tectonic activity in the direction of the Sumatran fault. In general, the rocks that make up the study area consist of volcanic rocks in the southwest and clastic rocks in the northeast of the Tertiary age (Hadi et al., 2011).
The geological structure in the Wai Selabung area is dominated by a northwestsoutheast direction cut by a fault in a southwest-northeast and a north-south. In general, the research area consists of volcanic rocks ( Figure 1). Based on the physical characteristics and age of the rocks, the rock units are obtained in order from old to young; Akar Jangkang Lava, Sandstone, Asadimana Lava, Pematang Gong Lava, Old Breccia, Ranau Pyroclastic Flow, Laai Lava, Bengkok Lava, Pandan Lava, Gedang Lava, Perean Lava, Tebat Gayat Lava, and Alluvium Deposits .

Surface Manifest
Based on the results of an integrated investigation in 2011, symptoms of geothermal appearances are characterized by the arrival of hot springs, which are divided into three groups of geothermal manifestations ( Figure  2), namely manifestations on the banks of the Wai Selabung River (2 groups), in Wai Selabung tributaries (2 groups). Alterations that occur on the banks of the Wai Selabung River appear near Suban hot springs (banks of the Wai Selabung river) . The physical characteristics of the hot springs are shown in Table 1. From the pH data and the temperature of the manifestations, it can be seen that the highest temperature of the hot springs is in Wai Selabung Hot Water I, which is 92.5℃, the lowest is in Wai Selabung III Hot Water, which is 40.2℃. Meanwhile, the highest pH is found in Wai Selabung II Hot Water, which is 9.43, and the lowest pH is in Selabung Damping Hot Water, 8.19. The degree of acidity (pH) is generally alkaline (8.19 -9.47). Strong electrolytes have strong electrical conductivity because they contain a more significant number of ions than weak electrolytes. In Table 1, it can be seen that the smaller the pH value will be, the greater the power of the electricity. Based on previous research, it was shown that the prospect area is only located in the southwest. Further research is needed to identify the Wai Selabung geothermal area's structural pattern, heat source, and reservoir Wai Selabung geothermal area using 3D inversion of the gravity anomaly. Furthermore, the gravity, geology, MT, and geochemical data compilation also supports identifying geothermal reservoir prospect areas.

METHODS
The gravity method was used to detect geothermal areas. The gravity method is a geophysical method based on the acceleration of gravity on the earth's surface. Changes influence variations in gravity values in the lateral density of the subsurface rock around the measurement point. This method identifies and describes subsurface geological structures based on variations in the earth's gravity field caused by differences in density between rocks.
In this study, a complete 3D inversion of the Bouguer anomaly was carried out, gradient analysis for structural analysis, and to determine the geochemical characteristics of the geothermal area of Wai Selabung. Bouguer anomaly inversion modeling correlated with MT, geological, and geochemical data.
The gravity data used in this study is secondary data from measurements carried out by the Center for Geological Resources in 2011 in the form of Complete Bouguer Anomaly. Gravity data was conducted griding data and contouring to produce the Bouguer Anomaly Map using Surfer 16 software. Spectrum analysis was carried out in several paths to determine the depth of regional and residual anomalies, so it is known that the window width is used for the anomaly Bouguer filter. Spectrum analysis was performed by Fourier transform (FFT) of trajectories representing the research area in a specific spatial domain ( ∆ ) (Blakely, 1996).
The separation of the Bouguer anomaly is done with a moving average filter. Moving average is done by averaging the anomaly value. This average is the regional anomaly, while the residual anomaly is obtained by subtracting the data from the gravity measurement with the regional anomaly.
, and N must be an odd number.
Second Vertical Derivative (SVD) analysis of the Bouguer anomaly was carried out to obtain fault structures, lithological boundaries, and sources of anomalies originating close to the surface. The zero contours flanked by high and low anomaly values indicate the presence of a fault structure in the area.
3D inversion modeling of residual Bouguer anomaly was carried out to obtain a subsurface density distribution model. 3D inversion modeling using Grablox 1.7 software.
Based on the analysis of gravity, MT, geochemistry, and geology of the research area, a compilation map is formed showing the reservoir prospect area of the study area. The research data processing process can be seen in the following flow chart.

RESULTS AND DISCUSSION Bouguer Anomaly
The Bouguer anomaly value of the Wai Selabung area is the result of data processing of the Geological Resources research group from field measurements with several long trajectories from northwest to southeast, namely from Sinamarga to Talang Najam and from southwest to the northeast, namely the Wai Selabung area which is then processed to obtain Bouguer anomaly map.
Based on the research of Wahyu et al. (2019), the estimation of rock density was carried out using qualitative analysis and quantitative analysis. Qualitative analysis was carried out by measuring the density of rock samples from the investigation site at the Geological Resources Center laboratory. In contrast, the quantitative analysis used was the Parasnis method. The laboratory analysis results show that the average value of the rock in the Wai Selabung geothermal research area is around 2.60 grams/cm 3 (Wahyu et al., 2019).
The Bouguer anomaly map and the distribution of gravity measurement points in Figure 4 show that the study area has varying anomaly values ranging from 4 to 46 mGal. The high Bouguer anomaly value ranges from 36 mGal to 46 mGal, which is in the northwest, north, and northeast of the study area. In comparison, the low Bouguer anomaly value ranges from 4 mGal to 22 mGal, which is in the middle of the study area flanked by moderate gravity values with contour values of 22 mGal to 36 mGal in the Sinamarga and Kotadalam areas.

Residual Anomaly
Residual anomaly is obtained from reducing the total anomaly value with the regional anomaly of the research area to bring up the residual anomaly value. The residual anomaly map shows a more complex pattern than the regional map and the Bouguer anomaly map. The residual anomaly map depicts an anomaly with a shorter wavelength that reflects the effect of a shallower object. The residual anomaly map of the study area has a residual anomaly value between -10 to 8 mGal. The residual anomaly map can be seen in Figure 5.
Low negative residual anomaly values range from -10 mGal to -4 mGal forming a closure in Kota Dalam, Teluk Agung, and Pematang Gong. The low anomaly is located in the center of the study area. Medium residual gravity values ranging from -3 mGal to 3 mGal dominate the investigation area in the north, south, east, and west directions.
Based on the research of Wahyu et al. (2019), the high anomaly was very clearly demarcating the low anomaly zone so that the high anomaly zone, which was near from the manifestation. This shows that deep structures probably cause this residual anomaly with high complexity (Wahyu et al., 2019).
High residual gravity values with values >4 mGal poles in the middle of the study area between Pematang Gong and Teluk Agung where geothermal manifestations are found and in the Pematang Petang area to the southeast of the study area.
The residual anomaly shows a contour lineation pattern in a northwest-southeast direction and shows positive and negative anomalous closures in the center of the study area. Where there are manifestations, indicates the presence of a northwestsoutheast trending fault structure by the main fault structure of this area trending northwest-southeast and slightly southwestnortheast.

Second Vertical Derivative (SVD)
The Elkins operator generates the Second Vertical Derivative (SVD) map from the anomaly filter process (Elkins, 1951). Black contour lines indicate anomaly contours with a value of zero. Zero contours flanked by high and low anomaly values indicate the existence of a fault structure in the area. From the SVD map, a fault structure line can be drawn against the anomaly contour with a value of zero. In the Second Vertical Derivative principle, the value of zero indicates the value of a boundary, so it can be assumed that the area has several faults.
Second Vertical Derivatives analysis was performed on residual anomalies. Figure 7 shows the SVD contour map of the study area. The filtered SVD contour map has anomaly values in the range of -8.5 to 5.5 mGal. The anomaly value is divided into three parts, namely low anomaly (-8.5 to -3.5 mGal), moderate gravity value (-2.5 to 2.5 mGal), and high gravity value (> 2.5 mGal).
The residual SVD contour map shows high SVD values in orange to red and low SVD values in purple to blue. In the center of the study area, the zero anomaly contours are identified as the boundaries of the geothermal reservoir. The moderate resistivity value pattern distribution with a value range of 45-380 Ωm is considered a reservoir because it is located in two zones with high and low contrast resistivity. A reservoir is a volume of a permeable rock layer that stores heat and a circulation fluid that extracts heat.
The distribution of high resistivity value (>380 Ωm) is thought to be a basement experiencing intrusion. It is suspected to be one of the heat sources given the high resistivity value compared to other zones.

3D Modeling of Gravity
The residual Bouguer anomaly 3D inversion modeling was carried out using Grablox software. This inversion modeling uses a box-shaped object model approach arranged in a grid with specific dimensions, as input is the residual Bouguer anomaly grid data. After running the data process, the output in the form of a density value is obtained, then displayed in a 3D model as shown in Figure 9. The inversion of the residual anomaly (seen from the southwest) produces a density distribution as shown in Figure 8 with a value range of 2 to 2.9 g/cm 3 with a maximum depth of 4000 meters. The distribution of low surface density values is in the north to a west area with a density value of 2 to 2.15 g/cm 3 . Low-density values indicate the location of the reservoir. This area is located in the geological setting of old Miocene breccia, which was overwritten by Pleistocene Tebat Gayat lava. The distribution value of medium density ranges from 2.2 to 2.4 g/cm 3 . The density distribution is in the middle to the ends of the study area; the north to the west has a depth of 225 to 4,000 meters. Medium-density values are probably tertiary-aged sandstone sediment deposits. The distribution value of high-density ranges from 2.5 to 2.9 g/cm 3 . High-density values are spread at a depth of 225 to 4,000 meters; high-density values are of interest in gravity surveys because they indicate potential heat sources and possible hydrothermal silicate deposits (Sarkowi & Wibowo, 2021). The density distribution model resulting from the 3D inversion of the residual Bouguer anomaly above is then sliced to obtain a 2D cross-section of the trajectory. Figure 10 shows the location of the reservoir in the study area based on the analysis of the density distribution of the study area. The A-A' trajectory is northwestsoutheast of the study area. The fault model is derived from a horizontal cross-section of order-1 (maximum or minimum) and order-2 (zero) gradient, while low density is assumed to be a reservoir in the area. The Pematanggong fault influences the fault in this trajectory with a southwest-northeast trending fault pattern and the Perean fault with a north-south trending fault pattern.

Figure 10. Density Distribution Cross-Sectional Model Trajectory A-A 'Section Results from the 3D Model
Density Inversion Modeling Figure 11 shows the location of the reservoir in the study area based on the analysis of the density distribution of the study area, which is correlated with the 1 MT trajectory. The B-B's trajectory is in the southwest-northeast direction of the study area. The fault model is derived from the cross-section of the first-order horizontal gradient (maximum or minimum) and the second-order gradient (zero). At the same time, the low density near the manifestation is assumed to be the reservoir in the area. The reservoir's existence is also supported by the MT model, which shows the existence of a reservoir in the area limited by faults. The fault in this path is supported by the emergence of hot springs. The Pematanggong fault influences the faults in this trajectory with a southwest-northeast trending fault pattern, the Wai Selabung fault and the Kotadalam fault, which control the manifestation of the Selabung Damping AP, and the Gistong fault with a southwestnortheast pattern which controls the occurrence of the Lubuk Suban AP manifestation.  Figure 12 shows the reservoir based on the density distribution analysis, which is correlated with the 2 MT model. The C-C' cross-section is southwest-northeast of the study area. The fault model is derived from the cross-section of the 1 st order horizontal gradient (maximum or minimum) and the 2nd order gradient (zero). At the same time, the low density near the manifestation is assumed to be the reservoir in the area. The reservoir's existence is also supported by the MT model, which shows the existence of a reservoir in the area limited by faults. The fault is one component of the geothermal system that functions as a fluid transport route from the reservoir. The number of fault zones will increase the amount of fluid that comes to the surface. This fault is reinforced by the presence of hot springs near track 3. The Akarjangkang fault influences the fault in this line with a north-south trending fault pattern that controls the emergence of the manifestations of AP Wai Selabung 1 and AP Wai Selabung 2.  Figure 13 shows the location of the reservoir in the study area based on the analysis of the density distribution of the study area. The D-D' trajectory is in the southwest-northeast direction of the study area. The fault model is derived from a horizontal cross-section of order-1 (maximum or minimum) and order-2 (zero) gradient, while low density is assumed to be a reservoir in the area.  Figure 14 shows the location of the reservoir in the study area based on the analysis of the density distribution of the study area. The E-E' trajectory is in the southwest-northeast direction of the study area. The fault model is derived from a horizontal cross-section of order-1 (maximum or minimum) and order-2 (zero) gradient, while low density is assumed to be a reservoir in the area.

Geochemical Analysis
Based on the results of an integrated investigation in 2011, the characteristics and types of hot water are based on the Cl-SO4-HCO3 triangle, which refers to Giggenbach (1988) with the type of chloride bicarbonate or chloride sulfate and type bicarbonate.
In determining the temperature of the geothermal reservoir in the research area, the Na-K-Mg triangle is used. Wai Selabung 1 hot spring and Wai Selabung 2 hot springs are in partial equilibrium and are thought to originate from the reservoir with a temperature of 160℃ to 180℃. The partial equilibrium zone indicates that the fluid reaction with the reservoir rock has reached partial equilibrium. The other hot springs samples are in the immature water zone, indicating that the hot water has been mixed with cold water on the surface in a high proportion.
To determine the origin of the recharge reservoir water using the stable isotopes Ddeuterium and 18 O. The plotting results of Wai Selabung 2 hot springs and Wai Selabung 3 hot springs are close to the meteoric water line, indicating that the hot springs are influenced by surface water. Meanwhile, the other hot springs, namely Wai Selabung 1, Lubuk Suban, Selabung Damping, Arumatai, and Kota Batu, show oxygen enrichment. So, that the position on the graph is to the right of the meteoric water line. An indication of the formation of Hot springs is related to the interaction between the hot fluid in the geothermal system and the rock that causes 18 O enrichment. This indication shows that Wai Selabung 1 Hot Water, Lubuk Suban Hot spring, Selabung Damping Hot spring, Arumatai Hot Spring, and Kota Batu Hot spring are more influenced by magmatic fluids.
In the 2011 integrated investigation, the bottom borehole temperature, surface temperature gradient, and surface heat flow were also carried out. The results obtained that the distribution of the anomaly zones of the three measurements is around the manifestation of the Wai Selabung hot spring to the Lubuk Suban hot spring.

Geothermal Reservoir Prospect Area Analysis
Previous research by the Center for Geological Resources in 2011 showed that the prospect area is only in the southwest . The analysis of gravity, MT, geochemistry, and geology of the study area is depicted in a compilation map showing the reservoir prospect area of the study area, as shown in Figure 15. The prospect area for the Wai Selabung geothermal reservoir is estimated to be located in old lava rock and sandstone sedimentary layers under the Ranau Pyroclastic layer, located around the Great Bay. This area is influenced by the Teluk Agung fault, which has a northwestsoutheast trend. The next reservoir prospect area is in the sandstone sediment layer under the Sapatuhu Pyroclastic layer and the Perean Lava rock around Perekan. The north-south Perean fault influences this area. Moreover, the last reservoir prospect area is in the Akarjangkang Lava Rock and the shale sandstone sediment layer around Talang Tebat. This area is influenced by the Pematanggong fault and the Gistong fault, which are trending southwest-northeast.

CONCLUSION AND SUGGESTION
The conclusions of this research can be drawn: the prospect area for the Wai Selabung geothermal reservoir based on the Bouguer anomaly inversion modeling, which is correlated with geological data, geochemical data, and magnetotelluric data, is estimated to be located around Teluk Agung, Perekan, and Talang Tebat. This reservoir prospect area can be exploited to meet energy needs in Indonesia.
It is necessary to conduct exploration drilling studies to support the results of subsurface research and to calculate the potential of geothermal energy in the research area so that a feasibility study can be carried out for the development of geothermal fields.

ACKNOWLEDGMENT
The author would like to thank all those who have helped the author in writing this journal.