Application Power Spectrum Analysis of Gravity Data to Investigate Subsurface Structures of the North Banggai-Sula Microcontinent in Maluku Sea Area

Abstract


Introduction
The measured gravity acceleration on the surface is composed of the predicted value of the Earth's gravity fields and the contribution stemming from anomalous sources below the surface.The direction of gravity field varies vertically, depending on the position of the anomaly source.The gravitational field in theory is a field that results from non-geological causes and the pricing is determined using theoretically defined formulas.Mark Latitude influences this terrain.height and the surrounding topographic mass of the point.In comparison to the Earth's gravity field, the magnitude of the gravity anomaly is very small (Grant, 1966).
The source of the anomaly in the northern Banggai-Sula microcontinent area in the Maluku Sea Regional (Fig. 1) is dominated by horizontal fault structures associated with the separation of several pieces of the Australian continent (Hall, 2002(Hall, , 1987;;Socquet et al., 2006).The tectonic framework of this regional is characterized by the existence of a very thick collision zone (about 14 km), and its width in the narrowest part reaches 150 km (McCaffrey et al., 1980;Socquet et al., 2006).This area has a large negative free-air anomaly value of -230 mGal (McCaffrey et al., 1980).The movement of the Banggai-Sula microcontinent westward was restrained by a collision with the East Sulawesi Arm, believed to have occurred during the Middle Miocene to the Pliocene (Satyana, 2011;Simanjuntak, 1986).Limited studies have been conducted in the northern area of the Banggai-Sula Islands, specifically the northern part of Taliabu Island (Watkinson et al., 2011).Several fault zones, such as the Sula (Hamilton, 1973;Pigram and Supandjono, 1985;Silver and Smith, 1983) and Balantak Faults, have been identified in the northern Taliabu.Most of Sulawesi's current seismicity was associated with the subduction of the Celebes and Mollucas Seas under the North Arm and in the north of the Banggai-Sula microcontinent, respectively (Cardwell et al., 1980).The free-air anomaly data may be acquired via Topex (https://topex.ucsd.edu/cgi-bin/get_data.cgi)using satellite altimetry measurement.The retrieved information was then processed to derive terrain corrections and Complete Bougeuer Anomaly (CBA) data (Maulana and Prasetyo, 2019).CBA data encapsulates total gravity anomaly information that represents the complete subsurface geological structure, for both residual and regional.Furthermore, the anomaly measured below the surface in gravity exploration activities is a superposition of various sources of detected anomalies.
The majority of studies conducted in the Banggai-Sula microcontinent region mostly focus on terrestrial areas, with only a limited number of studies giving consideration to the northern region of Banggai-Sula that shares a border with the Maluku Sea.A multitude of seismic observations were recently conducted in the vicinity using a multibeam system.Nevertheless, there has been limited scientific investigation of the gravity method in this particular region.This article provides a structural analysis of interpretations that heavily rely on newly acquired data from the northern area of the Banggai-Sula microcontinent.The southern Maluku Sea serves as the foundation enhanced comprehension of the significance of structure within the regional framework.
Satellite serves as an excellent alternative to data acquisition directly at the ocean locations.To estimate the rock structure and sediment thickness, this study employed 3D gravity modeling and spectrum analysis which was recognized for its ability to show significant depth contrast within the depth of the crust (Saibi et al., 2006).

Regional Stratigraphy
The northern Banggai-Sula area has a relatively simple stratigraphy (Garrard et al., 1988;Pigram and Supandjono, 1985).The lower layers are characterized by Paleozoic or older metamorphic rocks, intruded by Permo-Triassic granites associated with acid volcanic rocks.Towards the island's western part, limestone formations span from Eocene to Miocene ages, extending to the more recent Neogene period (Ferdian et al., 2011;Socquet et al., 2006;Watkinson et al., 2011).Some locations have Triassic limestone and dolomite, while terrestrial conglomerates and sandstones, usually rich in quartz, could originate from the Bobong Formation (Late Jurassic) exposed at the Taliabu, Peleng, and Banggai fault boundaries.Shale coal originates from the undated Kabauw Formation, believed to be equivalent to the Bobong Formation in Mangole and Sulabesi (Kusnama, 2008).
The Taliabu Liassic marine strata, characterized by mollusk and echinoid fragments, transition into the Bobong Formation which is covered by black shales from the Buya Formation (Middle Jurassic) to the Late Cretaceous (Ferdian et al., 2011).The Tanamu Formation was deposited at a depth of more than 200 m and is dominated by carbonate mudstone with planktonic foraminifera.From the results of a seismic survey by Watkinson et al. (2011), it was discovered that these mudstones were aligned in the Buya Formation, as shown in Fig. 2. Additionally, Neozoic carbonates were almost absent in the Sula archipelago (Taliabu, Mangole, and Sulabesi) (Ferdian et al., 2011;Watkinson et al., 2011).The seismic package at the lowest depth exhibits a predominantly featureless composition, suggesting its probable identification as the crystalline basement.The cultivation of this crop primarily occurs in the central region and the southern periphery of Taliabu and Banggai.The geological composition of the area consists of various folded metasediments, marbles, schists, gneisses, and amphibolites, which are believed to have formed during the Permo-Carboniferous period (Watkinson et al., 2011).To the northeast of the Sula Archipelago is the Maluku Sea collision complex, arising from the collision between the Halmahera and Sangihe arcs (Hall, 1987).Despite being predominantly submerged, studies of the small islands within the Central Moluccas Sea, North Sulawesi, and Halmahera suggest the likelihood of ophiolitic basements and volcanic arcs overlain by deformed volcaniclastic sediments from the Miocene to the present (Hall, 2002;Silver and Smith, 1983;Watkinson et al., 2011).

Materials and Methods
This study was conducted in the northern Banggai-Sula Islands, situated within the Maluku Sea at coordinates 0° 11-1° 32 S and 123° 52 -125° 43 E, with an area of 30,689 km 2 as shown in Fig. 1.The Bouguer anomaly encompasses both local and regional variations in gravity, which can be further classified as shallow and deep gravity anomalies (Novianto et al., 2022).Residual anomalies are caused by their shallow counterpart, with high frequencies and short wavelengths, while regional anomalies arise from deep sources with low frequencies and long wavelengths (Telford et al., 1990).Additionally, noise predominantly lies shallower than residual anomalies.To interpret subsurface structures effectively, the separation of regional and residual anomalies was essential and can be achieved through various methods.In geological studies, gravity explorations are commonly applied through various techniques to separate regional and residual anomalies (Karunianto et al., 2017;Purnomo et al., 2016;Shandini and Tadjou, 2012;Supriyadi et al., 2017).
This study used CBA data to estimate the boundaries of the shallow and deep discontinuities.Power spectrum analysis, a technique commonly used to estimate the subsurface depth limits of the anomaly source, was also adopted (Kumar et al., 2017).Spectrum analysis uses the principle of the Fourier transforms by changing data from space to frequency domain.
The gradient of spectrum graph was proportional to the depth of the anomaly source field.Areas of regional (deep) and residual (shallow) anomaly, were represented by higher and lower gradients, respectively (Apriani et al., 2018).The study relied on Topex gravity data derived from altimetry satellites.The primary information obtained from these satellites pertains to sea level topography, representing deviations in sea level from the geoid surface (Abidin, 2001).Using geoid undulation data, a free-air gravity anomaly can be generated (Balmino et al., 1987).The equation used in its calculation was based on the negative radial derivative of the disturbing potential.
With ∆  is the free-air anomaly value (mGal),  0 is the standard gravity value (mGal), r is the radius of the earth (m), N is the geoid undulation (m), and ∅ is a disturbing potential (mGal.m).The free air anomaly map is shown in Fig. 3.The free-air anomaly, in actuality, is insufficient for describing the underlying geological state.Therefore, it necessitates many standardized adjustments in order to provide gravitational data that is connected with the rock composition beneath the surface (Yanis et al., 2023).

Complete Bouguer Anomaly (CBA)
The Bouguer correction can be applied when both the density and depth of the ocean are known.To facilitate the correction on land, the concept of topographical mass was used (illustrated as filling the ocean with concrete).In an oceanic context, the Bouguer correction in the ocean can be imagined as an infinite slab equivalent to the water depth and the density difference between water and rock. = (0,0419)( − ℎ  ) = 0,0419(  −   )ℎ  (2) with  as the Bouguer correction at sea (mGal) and ℎ  as water depth from the observation point (m).Meanwhile,   =1,03 g/cm 3 and   = 2,67g/cm 3 , are seawater and seafloor densities, respectively.The thickness of the slab was equal to the depth of sea (Lillie, 1999).By using equation (2) above, the Bouguer correction at sea was subtracted from the free air anomaly to obtain the following Bouguer gravity anomaly at sea (∆  ): ∆  = ∆  +  (3) With ∆  is the anomaly-free air value.

Spectrum Analysis
Spectrum analysis was used to examine the phenomenon of harmonic oscillators in nature.Furthermore, it aimed to obtain spectrum distribution of the oscillator phenomenon and present its statistical characteristics.In one-dimensional spectrum analysis, gravity field anomaly data was examined within a single section.
The distribution of gravity acceleration anomaly was represented in the form of a gravity anomaly map termed CBA.The value of this map was a combination of regional and residual anomalies, which can be separated in the frequency domain.This process is mathematically feasible because subsurface objects generate gravity responses with both high and low frequencies.Since gravity response occurs in the distance domain, conversion techniques are employed to shift data to the frequency domain.Gravity field anomaly is a superposition of regional and local (residual) gravity field anomalies.The power spectrum approach involves applying a high pass filter on gravity data to obtain the residual, then using a low pass filter to obtain the regional field (Al-Khalidi et al., 2023).The following equation presents the frequency domain relationship between the gravity field anomaly and the density distribution: ∆() = 2∆() − (4) Where ∆() is the frequency response of gravity anomaly field,  is the frequency response of the density contrast, d is the depth of field boundary of the spheroid reference field, and G is the universal gravity constant.When the density distribution is random and has no relationship with each gravity value, the response frequency, ∆() becomes 1.  =  −2|| (5) Where  is a constant,  = 2πk is the corner frequency, k is the wave number (cycles/meter), and d is the depth of field below the reference spheroid.The following represents equation ( 6): =   − 2|| (6) The two logarithmic values are the difference between the two power spectra in equation ( 7) below, where E1 and E2 are power spectrum, k1 and k2 is the number of waves, and ∅ is the angle of inclination of power spectrum curve.Equation ( 7) showed that the mean depth of the discontinuity field was proportional to the slope gradient of power spectrum curve.

3D Inversion Modeling
Inversion modeling is a technique through which model parameters are directly obtained from observational data (Zarkasyi and Suhanto, 2013).According to Grandis (2009), often referred to it as data fitting because the process involves the pursuit of model parameters to produce responses that match the observed data.The objective was to attain a high degree of concordance between the model response and the observational data, leading to an optimal model representation (Supriyanto, 2007).
Inversion modeling was the focus of most or almost all fields of geophysics because there is a need to estimate models or model parameters based on observations or field data measurements.An example of its application in this study involves the estimation of the subsurface structure model, specifically in the context of density value distributions derived from gravity method measurements.
3-dimensional (3D) inversion technique was used for subsurface structure modeling.To achieve this, the CBA gravity anomaly data were inverted using Grablox 1.6 software by Markku Pirttijärvi, enabling the production of a 3D cross-sectional density model.Furthermore, the outputted model can be visualized as a 2D or 3D section.The Grablox 1.6 software combined both Singular Value Decomposition (SVD) and Occam inversions (Hjelt, 1992), which were processed sequentially.

Results and Discussion
Terrain correction was applied based on the bathymetric map, to correct the effect of irregular mass distribution around the measurement point.Subsequently, the Bouguer anomaly was obtained at each measurement point, while its data was presented on a map using the Geosoft Oasis Montaj software, as shown in Fig. 4.
Fig. 4 shows the CBA map with the distribution ranging from -59.8 to 340 mGal.These anomalies can be grouped into four categories, namely negative, low-positive, medium, and high.
The negative anomaly range encompasses values ranging from -59.8 to -3.0 mGal, extending primarily to the western regional of the study area, as indicated by dark blue coloration.From this distribution, the anomaly showed a relatively wide area.A low positive anomaly is characterized by values ranging from 7.8 to 45.7 mGal.This anomaly is predominant, covering almost all of the total study area in the middle, with a light blue to green regional.Moderate positive anomaly, with values between 72.5 to 249.6 mGal, clusters around the north and south of the study area.
Lastly, areas with high positive anomaly values, ranging from 240 to 340 mGal are indicated with red and white to the north and south of the study area.Low anomaly indicated a fault collapse zone caused by oceanic and continental crust subduction.The moderate anomaly was attributed to thick sedimentary rocks covering the center of the study location.Additionally, the high anomaly, in the southern regional, is attributed to the prevalence of highdensity igneous rocks within a sedimentary rock environment characterized by lower densities.The gravity anomaly measured in the field was a combination of regional and residual anomalies.

Power Spectrum Analysis
Power spectrum analysis distinguishes the distribution of CBA anomaly values in regional and residual areas.Eleven cross-sections are analyzed according to the direction of the cross-sections in the 3D modeling.This section is a sampling of data from the CBA map (Fig. 4) to be transformed by the Fourier.
Power spectrum analysis serves to differentiate the distribution of CBA values within both regional and residual zones.Eleven cross-sections undergo examination based on their orientation within the three-dimensional model.This segment involves extracting data samples from the CBA map (refer to Fig. 5) that are subsequently subjected to Fourier transformation.Spectrum analysis result was presented as a scatter diagram of the amplitude and wavenumber, which were approximated by a linear function, as shown in Fig. 6.The intersection of the linear function of regional and residual areas, represented by the value of k was a reference in separating the anomalous regional.

Sub-surface Modeling
The results of the gravity modeling were presented in 3D, employing data extracted from the CBA.The inversion process is performed by optimizing three crucial variables: base, density, and block geometry.This optimization is based on the SVD and Occam methods.An optimization process was conducted to maximize the coefficients of the second-order polynomial that are linked to regional anomalies.The implementation of this improvement leads to an RMS inaccuracy of 23%.The density and block shape have been optimized to achieve an even distribution of contrast density, thickness, and depth of subsurface layers in the region.This optimization results in a 12% reduction in the root mean square (RMS) error compared to the baseline optimization.The SVD inversion approach is used to minimize the discrepancy between observation data and computation response by optimizing the contrast density and block design.This optimization is also performed using the Occam approach to generate a more refined and precise model by using all the restricted inputs in the reconstruction of the initial model.The root mean square (RMS) error resulting from Occam inversion for both the data and model of density contrast and block geometry is 4% and 3% respectively.
The outcome demonstrates a strong correlation between the observed data (Fig. 7a) and the computed data (Fig. 7b), with a minor difference (Fig. 7c).The density model resulting from the inversion process had intersections in two directions.As a result, the 2D cross-section model can only be generated vertically (intersecting the x and y axis) and horizontally (intersecting the z-axis).The direction of the section to be analyzed is not arbitrary but should follow the cross-section that the software had determined.In making the initial model, the x-axis was divided into twenty blocks, such that twenty 2D cross-sections are produced perpendicular to the yaxis in the final model.The y-axis (divided by twenty blocks) also generated twenty 2D cross-sections.However, not all sections are interpreted, as only a few were considered more representative of the subsurface conditions.
Several considerations were made before determining the orientation of the 2D sectional crosssection, as presented in Fig. 9.Some include the direction and distribution of faults and the morphology of the study location.In light of these factors, eleven vertical paths were chosen in the direction of the x-axis, which can represent the situation under the surface Table 1 showed the analysis results of the cross section in 3D modeling and power spectrum conducted on the CBA map.From Table 1, the deep and shallow discontinuities obtained from the 3D modeling method and power spectra have similar relative values.The average for the shallow and deep areas ranged between 1.23 km and 1.39 km, as well as 14.52 km and 15.14 km, respectively.The sediment thickness value becomes thicker towards the center of the research area.
The models are analyzed by considering geological data to determine the position of the study target.15.14The model picture depicts blue-colored rocks with density values ranging from 2 g/cm 3 to 2.7 g/cm 3 scattered at a depth of 2 km.The density of these rocks dominates the study area.Meanwhile, those with green, yellow, orange, and red distributions are scattered and clustered in several places starting at a depth of 10-30 km.High-density rock with a value of 3-4 g/cm 3 in the east and west of the study area (Fig. 10), originated from below to a depth of 20 km.They are interpreted as intrusive igneous rocks.In contrast, those with a low-density range of 2.0-2.3 g/cm 3 were discovered at a depth of 2 km to 25 km.Fig. 11 indicated that the shallow discontinuity area (< 2 km) was dominated by sediment rocks with a 2.6 g/cm 3 to 2.8 g/cm 3 density, the Buya Formation superimposed by the Tanamu Formation.The deep discontinuity area with a depth of up to 16 km, indicates the presence of a fault.In more detail, areas of high density (density value is >3 g/cm 3 ) were identified by orange, and red adjacent to the dark blue areas, which presented low rock density (density value is 2 g/cm 3 to 2.3 g/cm 3 ).From the description of the geological map, the area was a cross fault of the Sula and Balantak faults as presented in Fig. 12.According to the publication by (Watkinson et al., 2011), the existence of the Sula fault was estimated to be at a depth of >20 km.However, in 2000, the Peleng tsunami was generated at a focal depth of 18.6 km (Baeda and Husain, 2012).In this study, the rock structure around the fault area began to appear at a depth of 16 km and was more clearly visible at a depth of >20 km The observed above the basement are interpreted as fault-bounded formations consisting of terrestrial conglomerates and sandstones from the Bobong Formation.These formations are located within a half-graben and are found to rest unconformably on the basement towards the southern onshore region.
According to the results of the 3D modeling analysis, the shallow discontinuity area was interpreted to be composed of sediment rocks with a density of 2 g/cm 3 to 2.75 g/cm 3 .Meanwhile, the deep discontinuity field was composed of rock with a density of 2 g/cm 3 to 3.47 g/cm 3 .A rock density of more than 3 g/m 3 was suspected to be a basement (Ferdian et al., 2011).

Conclusions
In conclusion, the analysis results based on power spectrum method and 3D modeling of the CBA showed that the average thickness of the shallow discontinuity area ranges between 1.23 km and 1.39 km.Meanwhile, the average thickness of the deep discontinuity area ranges between 14.52 km and 15.14 km.The shallow discontinuity area (< 2 km) was dominated by sediment rocks with a 2.6 g/cm 3 to 2.8 g/cm 3 density.The area is interpreted as the Buya Formation superimposed by the Tanamu Formation.The deep discontinuity area with a depth of up to 16 km.We found areas of high density (density value is >3 g/cm 3 ) adjacent to low rock density (density value is 2 g/cm 3 to 2.3 g/cm 3 ).The observed are interpreted as fault-bounded formations of terrestrial conglomerates and sandstones from the Bobong Formation.These formations are located within a half-graben and are found to rest unconformably on the basement towards the southern onshore region

Fig. 1 .
Fig. 1.The northern Banggai-Sula Archipelago is restrained to the west by the eastern arm of Sulawesi, the study location is shown in the box, (the source of the faults line refers to the Indonesian earthquake source and hazard map published by the National Earthquake Center (Pusgen) in 2017).

Fig. 2 .
Fig. 2. Brief stratigraphy of a stratigraphy interpretation from the seismic data (red line on the map insert) and study area (black box) in Banggai-Sula north of the Maluku Sea regional (modified from Watkinson et al. (2011))

Fig. 3 .
Fig. 3. Free air anomaly map at the study area.The red area has a high value for Free Air anomaly, whereas the blue area has a low value for Free Air anomaly

Fig. 4 .
Fig. 4. The CBA map of the study area with fault lines.The red area has a high Bouguer anomaly value, whereas the blue area has a low Bouguer anomaly value

Fig. 7 .
Fig.7.The research area's 3D inversion modeling findings were obtained using the Singular Value Decomposition (SVD) and Occam methods.The final results displayed are from the Occam h procedure.(a) observation, (b) computational, (c) difference response

Fig. 8 .Fig. 9 .
Fig. 8.The 3D model of the density distribution with the x-axis direction of the cross-section in 3D modeling

Fig. 12 .
Fig. 12. Interpretation of rocks structure around the fault area in the deep layers by 3D modeling

Table 1 .
Estimation of deep and shallow depths for each line from power spectrum and 3D modeling