CASE STUDY ON OPTIMIZING DISC CUTTER SPACING BASED ON VERTICAL FORCE THRESHOLD IN COLORADO RED GRANITE CUTTING USING LS-DYNA SOFTWARE

3.1. Gravity data

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.

Figure 2. Research data (a) Topography (DEMNAS, 2023), (b) Complete Bouguer Anomaly (CBA), (c) Bouguer density estimation using Parasnis approach, where the gradient of the straight line represents the Bouguer density.

3.2. Upward continuation

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([TeX:] x_{0}$x_0$, [TeX:] y_{0}$y_0$, [TeX:] z_{0}$z_0$ = 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):

The source's gravity anomaly, [TeX:] \mathrm{\Delta}g_{0}$\Delta g_0$, at the mean sea level point, is given by the vertical component, [TeX:] dg_{z}$d{g}_z$.

At the upward continued surface, point A, the source's gravity anomaly, [TeX:] dg_{P}$d{g}_P$, equals

3.3. Edge detection technique

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:

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.

Figure 3. Illustration of the upward continuation technique (Kebede et al., 2020).

Where ∂G/∂x and ∂G/∂y are the first-order horizontal derivatives of the gravity field in the x and y directions.

An additional technique for enhancing source edges is presented by Oksum et al. (2021). This filter is known as fast sigmoid edge detection (FSED), which operates on gradient amplitudes. This filter can strengthen the edges of shallow and deep structures simultaneously. The mathematical equation is as follows:

Where:

3.4. Euler deconvolution depth solutions

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:

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.

4.1. Residual gravity anomaly

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

Figure 4. Separation of gravity anomalies a) regional anomaly data and b) residual gravity anomaly data.

Negative gravity anomalies shown in Figure 4b may indicate the extent of sedimentary layers from the surface to the sedimentary basin layer. Meanwhile, positive gravity anomalies, with a maximum value of 8.8 mGal, may indicate the presence of harder, thicker rocks surrounding the sedimentary rocks. This assumption will be discussed further in the inversion model of gravity data. In general, however, we can see the correlation between gravity anomalies and rock units representing the surface geology of the study area. A correlation of the gravity data with the geological map (see Figure 1b) shows that positive anomalies are generally associated with metamorphic, meta-sandstone, and meta-limestone units. Gravity anomalies ranging from -2.0 mGal to -0.2 mGal are generally associated with conglomerate rock units. Smaller gravity anomalies (i.e. less than -0.2 mGal to -7.7 mGal) are more prevalent in alluvial units. Hot spring manifestations are generally located in the high gravity anomaly zone, with gravity anomalies ranging from 1.8 to 8.8 mGal. Some coincide with the transition zone between positive and negative anomalies, such as hot springs M1 and M10.

4.2. Inversion of gravity data

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

Figure 5. Horizontal depiction of subsurface layers at various depths: a) 200 m, b) 500 m, c) 800 m, d) 1200 m; e) 1500 m, f) 2000 m, g) 2500 m, h) 3000 m, and i) 3500 m below sea level (bsl). Red indicates high density contrast and blue indicates low density contrast.

Figure 5 presents the horizontal slicing at depths of 200, 500, 800, 1200, 1500, 2000, 2500, 3000, and 3500 meters below sea level (bsl). The horizontal slicing reveals an intrusive rock pattern with a density contrast of 546 kg/m³. The peak of the intrusive rock (dark-red) occurs at a depth of roughly 1500 meters bsl and is consistent in deeper zones. The pattern of intrusive rock branches from west to east. In shallower zones (< 1200 m), density contrast varies from -20 to 366 kg/m³. These may demonstrate the diversity of rock types in the shallow zone.

The high-density contrast (red) in Figure 5 is associated with the high-gravity anomaly in Figure 4b, while the lower-density contrast (blue) is correlated with the lower-gravity anomaly. The lower-density mass anomaly in the shallow zone, extending to a depth of about 1,500 meters below sea level (see Figures 5a to 5e), contributes to the lower-gravity anomalies in the study area. Surface hot spring manifestations are generally located at the boundaries between areas of high and low-density contrasts. The presence of near-surface mass anomalies with large density contrast gaps around hot springs may indicate weak zones caused by geological structures, such as fractures and/or faults. The Boroboro Fault, Kaendi Fault, Landai Fault, and Lainea Fault, as shown in Figure 1, support the above interpretation. These faults are located around hot springs or the boundary zones between high- and lower-density contrasts.

The vertical slicing of the inversion model is presented in three variants based on direction: West-East, North-South, and oblique slices with certain azimuth angles. Figure 6 shows four vertical slices through the peaks of maximum anomalies, which actually coincide with the hot spring manifestations. The image in the top right of Figure 6 is included only to show where the four vertical slices are located. These slices are labelled (a), (b), (c), and (d). The first vertical slice runs West-East through latitude 9516500 UTM (see Figure 6a). Figure 6b illustrates the vertical slice's oblique direction, positioned at an azimuth angle of approximately[TeX:] \ 75^{o}${75}^o$. These two models clearly indicate the presence of the vertical intrusive rock model with a maximum density contrast of 546 kg/m³ in the West-East and Southwest-Northeast orientations. Furthermore, the layering system's strata, which are defined by density contrast differences ranging from -20 kg/m³ (dark-blue) to 546 kg/m³ (dark-red), are plainly evident. In addition, the layering model reveals the presence of basin-shaped and fault patterns in the environment surrounding the intrusive rock body. These basins are characterized by layers with varying density contrasts ranging from 237 kg/m³ (light blue) to 314 kg/m³ (light orange), where the layer with a density contrast of 263 kg/m³ (very light gray) complements the basin's central zone.

Figure 6. Vertical depiction of the subsurface layer in various directions: a) West-East direction at latitude 9516500 UTM, b) Southwest-Northeast direction with an azimuth of [TeX:] 75^{o}${75}^o$, c) North-South direction at longitude 455500 UTM, and d) North-South direction at longitude position 458700 UTM. Red indicates high density contrast, and blue indicates low density contrast. The position of the four vertical slices is shown in the uppermost right figure.

Figures 6c and 6d, respectively, represent the vertical slice model trending North-South at longitude positions 455500 UTM and 458700 UTM. Each of these slices passes through the top of an intrusive rock. The intrusive rock model, characterized by a density contrast value of 546 kg/m³ (a dark-red color), is modelled very clearly. The model also shows layering strata according to their density contrast and displays the presence of basin and fault patterns.

5.1. Geothermal system at LNVGF

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.

Figure 7. Geological model of the Lainea geothermal system: a) West-East direction at latitude 9516500; b) North-South direction at longitude 455500. The top image shows the positions of Lines (a) and (b), which represent the West-East and North-South positions of the geological model, respectively.

Figure 7b shows the geological model in the North-South direction, which cuts perpendicular to the West-East cross-sectional model (see Figure 7a) at the longitude position 455500 UTM. The shape of plutonic rocks in the North-South cross-section is much narrower than the shape in the West-East direction. This represents the orientation of the extension process, which tends to be west-east due to the Banda Rollback, which is located east of the Southeast Sulawesi arm. As previously explained, tectonic processes such as metamorphic exhumation, extension, crustal thinning, and magmatism centered in the central part of the study area cause the northern and southern basins to become thicker and deeper than the middle basin. The northern fault in Figure 7b is suspected to be the Kaendi Shear Fault, which is associated with the M1 and M2 hot spring manifestations. The southern fault is suspected to be the Boroboro Fault, which acts as a conduit for the M3 manifestation.

. The basement is the layer underlying the reservoir rock. Based on the geological model in Figure 7, the basement is the integration of rock layers with densities of 2948 ≤ ρ < 3045 kg/m³ and 3046 ≤ ρ ≤ 3122 kg/m³, which are metamorphic rock groups. Geological information (see Figure 1b) confirms the occurrence of the upper and bottom basement on the surface, which is represented by the same rock formation (i.e. metamorphic rock units). The identification of the basement is very important. It is a zone where plutonic intrusions occur, which are the heat sources of non-volcanic geothermal energy. Based on the model in Figure 7, the upper boundary of this basement is located at a depth of several meters asl to a depth of 1500 m bsl.

5.2. Subsurface fault lineament model

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.

Figure 8. Subsurface fault alignment model at LNVGF, Southeast Sulawesi, Indonesia: a) from the TAHG technique; b) from the FSED technique.

Figure 9. Euler deconvolution depth solution for the Lainea geothermal system, Southeast Sulawesi, Indonesia. Dots with various colors from dark blue to dark red represent variation depths of subsurface geological structure of the area.