The igneous rock intrusion beneath Ambon and Seram islands, eastern Indonesia, based on the integration of gravity and magnetic inversion: its implications for geothermal energy resources

Richard LEWERISSA, Sismanto SISMANTO*, Ari SETIAWAN, Subagyo PRAMUMIJOYO Department of Physics, Faculty of Mathematics and Natural Sciences, Papua University, Manokwari, Indonesia Department of Physics, Faculty of Mathematics and Natural Sciences, Universitas Gadjah Mada, Yogyakarta, Indonesia Department of Geological Engineering, Faculty of Engineering , Universitas Gadjah Mada, Yogyakarta, Indonesia


Introduction
Eastern Indonesia is in the convergence zone of the three major world plates, including Eurasia, Pacific, and Australia. This region is a complex island arc of plate boundaries, collision zones, and active subduction associated with the Banda arc (Hamilton, 1979;Audley-Charles et al., 1988;Honthaas et al., 1999;Hinschberger et al., 2005). In general, geothermal systems in Indonesia are characterized by the manifestation of fluid at the surface with boiling temperatures, which occur in the area of Quaternary and active volcanoes along the volcanic arcs. All Quaternary volcanic regions are related to cooling magma or igneous rock intrusion, which in turn is a source of heat from the geothermal system in the active arc (Hochstein and Sudarman, 2008). Ambon is in the southern part of Seram Island, which belongs to the outer Banda arc (Figure 1). The Banda arc is mostly composed of mélange rocks as a result of Tertiary subduction (Marini and Susangkyono, 1999). From this process, Ambon and the surrounding islands have many geothermal manifestations that appear on the surface (Poorter et al., 1989).
Several studies on the Banda arc in eastern Indonesia have been carried out previously. Milsom et al. (2001) qualitatively reviewed the short-wavelength high gravity anomalies around the Banda Sea caused by Ophiolite rock dominance. Widiwijayanti et al. (2004) conducted studies on the evolution of geodynamics in the Maluku sea based on three-dimensional (3D) inversion of the Earth's gravitational field. Fichtner et al. (2010) studied lithosphere layer subduction in the Seram Sea using a combination of full-waveform tomography and isotope ratios. Widiyantoro et al. (2011) determined a complex structural model of the lithosphere slab under the Banda arc using a seismic tomographic model. Pownall et al. (2017) researched the evolution of the tectonometamorphics of the islands of Ambon and Seram, eastern Indonesia, based on the geochronology of 40Ar /39Ar. From this information, it is known that there has not been a specific study of subsurface structures that can explain the physical phenomena of the appearance of geothermal manifestations in Ambon and its surroundings based on the integration of gravity and magnetic methods.
The aims of our research are to determine the intrusion model of igneous rocks beneath Ambon and Seram by utilizing gravity data from the World Gravity Map (WGM) in 2012 and magnetic data based on the Earth Magnetic Anomaly Grid (EMAG2-V3) in 2017. The WGM is a world gravity anomaly map derived from the global earth gravity model or EGM 2008 and the terrain correction using the ETOPO1 topographic model Bonvalot et al., 2012). EMAG2-V3 is the latest version of the global earth magnetic anomaly model developed by Meyer et al. (2017), which is an update of the EMAG2-V2, with the addition of the global magnetic field anomaly on geoid and upward continuation at an altitude of 4 km from the geoid. Maus et al. (2009) first estimated the coefficients of a spherical harmonic expansion of the magnetic potential to degree 150 from the merged grid at 4 km above the geoid.
This study was conducted through the integration of 3D inversion modeling of gravity and magnetic anomalies by adding a constraint model from the CRUST 1.0 global earth crust model of Ambon and Seram (https://igppweb. ucsd.edu/~gabi/crust1.html#matlab). CRUST 1.0 is a global earth crust model and a substitute for CRUST 2.0 and CRUST 5.1 with a spatial resolution of 1° × 1°. This model provides the latest information in the form of global surface topography data, seabed bathymetry, refraction seismic, and also thickness data of ice sheets, sediments, and worldwide crust (Laske et al., 2013). We also determine the Curie point depth (CPD) of magnetic data from the EMAG2-V3 to estimate the bottom source of magnetic rocks where the ferromagnetic mineral changes to paramagnetic properties due to high temperature influence. It is expected that the results of gravity and magnetic inversion will be able to  Poorter et al., 1989;Widodo et al., 2007;Ryan et al., 2009;Villeneuve et al., 2010). explain and determine the geodynamic mechanisms under Ambon and Seram, eastern Indonesia.

Geology and tectonic setting
In the southern part of Seram Island, which belongs to the outer Banda arc, Ambon and the surrounding islands are situated in the plate subduction zone with an extreme curvature of up to 180° (Pownall et al., 2016). This arc is a nonvolcanic arc that mostly consists of mélange during tertiary subduction, including old continental crust slices (Marini and Susangkyono, 1999). Overall, there are similarities in the geology of Ambon, Haruku, Saparua, Nusalaut, and western Seram (Figure 2). This region is dominated by felsic to basaltic volcanic rocks with glassy characteristics (Poorter et al., 1989). Ambon Island consists of two main geological sections, namely Hitu in the north and Leitimor in the south. Hitu is dominated by volcanic rocks composed of cordierite and dacite pads (ambonites), locally covered by a coral reef and quarter alluvium deposits (Honthaas et al., 1999).
Leitimor consists of upper Triassic rocks and metamorphic limestone overlaid by Ophiolitic rocks (Van Bemmelen, 1949;Honthaas et al., 1999). Haruku and Saparua islands are also generally composed of felsic volcanic rocks, while southwestern Seram is dominated by metamorphic rocks of Paleozoic time, which are tectonically overlain by ophiolitic rocks (Monnier et al., 2003).

Materials and methods
The study area covers Ambon and western Seram in Maluku province, eastern Indonesia, with an area of 2° × 2°, located at coordinates 2°S to 4°S and 127°E to 129°E. The study area showed subduction activity in the southern part of Seram, which is considered to be a trigger for potential geothermal resources in or around Ambon (Figure 3). Modeling of gravity inversion was conducted using Bouguer anomaly data of Ambon and Seram from WGM 2012 (Bonvalot et al., 2012).
WGM 2012 is a high-resolution gravity anomaly map based on global-scale digital computing of the reference Figure 2. Regional tectonic and geological settings of Ambon and Seram islands, eastern Indonesia (Tjokrosapoetro et al., 1993). earth gravity and elevation models. WGM 2012 consists of three types of anomaly maps: free air, Bouguer, and isostatic anomalies derived from the EGM geopotential system and global relief model ETOPO1 with a grid resolution of ±3.7 km Bonvalot et al., 2011). We used the physical parameter data of the CRUST 1.0 global earth crust developed by Laske et al. (2013) as the constraint in gravity modeling. CRUST 1.0 is a replacement for global crust models with a spatial resolution of 1° × 1°, including CRUST 2.0 and CRUST 5.1. CRUST 1.0 provides the latest data covering global surface topography, seabed bathymetry, seismic refraction, and ice sheet, sediment, and earth crust thickness data (Laske et al., 2013).  (Table).
Our work also used EMAG2-V3 magnetic anomaly data from Ambon and Seram for subsurface inversion modeling (Meyer et al., 2017). The magnetic anomaly data have the same spatial grid as the WGM 2012 gravity data, which is 3.7 km. It is focused on more than 50 years of analysis of magnetic field anomaly measurements from satellites, ships, and air measurements. The EMAG2-V3 model consists of geoid or average sea level magnetic anomalies and global magnetic anomalies starting upward at an elevation of 4 km from the geoid. Maus et al. (2009) developed the Earth Magnetic Anomaly Grid (EMAG2) with 2 arc min resolution, a new grid with more data, improved spatial resolution, and 4 km reduced altitude, while EMAG2-V3 combines marine and airborne observations of track lines, satellite data, and magnetic observatory data to map the location, intensity, and extent of lithospheric magnetic anomalies (Meyer et al., 2017).

Horizontal gradient and analytic signal
The research uses a combination of horizontal gradient evaluation and analytical signals from gravity and magnetic data in the Ambon and Seram islands, eastern Indonesia. Horizontal gradients are used to locate the boundaries of density or susceptibility contrast of the potential field data. However, this approach is powerful to delineate either shallow or deep sources relative to the vertical gradient, which is only effective for the shallower structures (Aderbi et al., 2017). Analytical signals of gravity and magnetic  data were obtained to delineate the edge of the subsurface in the study area, such as faults and intrusive body contact. The horizontal gradient magnitude (HGM) of gravity data in the x and y directions can be given as follows: (1) where g is the gravity anomaly.
It is also possible to give the HGM of magnetic data in the x and y directions by: where T is the magnetic anomaly.
The maximum amplitude of the analytic signal is a combination of horizontal and vertical gradients (Hsu et al., 1996). The basic concept of an analytic signal has been described extensively using 2D magnetic data (Nabighian, 1972). Analytical signals from gravity data have also been developed by Hansen et al. (1987) using Poisson's relation between magnetic and gravitational fields. In this section, the depths and boundaries of magnetic data were located (Nabighian, 1972(Nabighian, , 1974. The analytic signal amplitude is defined as the square root of the sum of the squares of the vertical and horizontal derivatives (Roest et al., 1992). Marson and Klingele (1993) in Saibi et al. (2006) define the analytical signals from observed gravitational field data generated by 3D sources as follows:  (Blakely, 1996). Analytical signals from magnetic fields at a location (x,y) can be derived from three orthogonal total magnetic field gradients using the following equation (Roest et al., 1992;MacLeod et al., 1993;Yadav et al., 2018): where A(x,y) is an absolute value from the analytic signal at position (x,y) and T is an observed magnetic field. 3.2. Regional and residual anomaly separation Regional and residual anomaly separations of gravity and magnetic data have been determined to obtain the main input data for gravity and magnetic inversion. Gravity and magnetic anomalies are a combination of regional influences related to deep and wide structures, as well as local effects by shallow structures near the surface (Nishijima and Naritomi, 2017). We used the upward continuation method for anomaly separation, which is an operation to shift data at a certain height above the Earth's surface that is used to estimate regional-scale trends from potential field data . Upward continuation is the opposite of downward continuation to eliminate high frequency using low-pass filters. The upward and downward continuations of a potential field can be written as (Blakely, 1996;: is a radial wave number.

Depth of magnetic sources
CPD analysis of the upward continuation of magnetic data at an altitude of 4 km was performed to predict the depth of magnetic sources in Ambon and Seram. CPD is a lower magnetic source theoretical surface with a temperature of 580 °C (Hsieh et al., 2014). Techniques for determining a magnetic source's depth rate can be divided into two categories: determining the isolated anomaly shape and defining the magnetic anomaly statistical characteristic. Both of these methods relate the magnetic anomaly spectrum to the magnetic source depth by transforming spatial data into the frequency domain (Tanaka et al., 1999;Dolmaz et al., 2005). Depending on the power spectrum of magnetic anomalies, the top bound and the centroid of a magnetic source Z t and Z 0 , separately, are used to determine the basal depth of a magnetic source Z b , which is known as the CPD. The CPD can be calculated by the following equation (Bektas et al., 2007): The top bound and centroid of a magnetic surface composed of a horizontal (equivalent) sheet are determined from the slope of the power spectrum. The benefit of a magnetic anomaly spectrum analysis is that estimates of the top bound and the centroid of a magnetic source can be obtained by simply assuming the magnetic source dimension and the magnetization vector (Tanaka et al., 1999). We use Fourpot software  to determine the spectrum of magnetic anomalies in the frequency domain and then estimate the depth of the upper bound and centroid of magnetic anomalies.

Gravity and magnetic inversion modeling
Gravity inversion modeling was done to obtain an Ambon and Seram subsurface model that could clarify the presence of igneous rock intrusion as a source of geothermal energy in the study area. Gravity forward and inversion processes were previously conducted using GRABLOX 2.1 ) and BLOXER software for editing and visualizing the block model (Pirttijarvi, 2012). In particular, GRABLOX is based on a 3D block model, where large rectangular blocks represent the volume below the gravity observation regions, subsequently divided into a small minor block. In addition, the inversion modeling objective is to minimize the difference (e) between measured gravity data (d i ) and estimated gravity model responses (y i ). The use of GRABLOX was based on two methods: singular value decomposition (SVD), with adaptive damping, used within unconstrained inversions, and Occam, adopted with two important parameters in the process of reducing the roughness of the resulting model. These include the Lagrange multiplier adopted to declare discontinuities, which is an important control to minimize data errors and roughness . In comparison to SVD, Occam inversion appears to produce more detailed models, and the RMS error between response and observations is written as follows : where Δd = maximum value of the gravity anomaly (d max )minimum value of gravity anomaly (d min ) used for scaling data. The inversion process was done repeatedly until a certain number of iterations to get a small difference. The main problems of the inversion are its nonunique nature and instability (Widiwijayanti et al., 2004), so we reconstructed the initial model ( Figure 5) in the form of 3D minor blocks using CRUST 1.0 global earth crust information (Table).
The initial model was made in five layers without seawater by adding a margin of 1 km 2 from the Bouguer anomaly data area to reduce the edge effect during the inversion process. The inversion on GRABLOX was done by optimization of the base anomaly, density, and block height.
The modeling of magnetic inversion on Ambon and Seram islands was performed using Oasis Montaj from Geosoft Inc. technology. Magnetic anomalies on the surface of the Earth that correspond to the susceptibility properties of the subsurface can be represented as follows in a linear equation (Li and Oldenburg, 1996):  that has the element of g ij , which is quantitatively the contribution of the susceptibility unit to the jth cell to the ith (datum) point of measurement. The g ij function is a projection in a particular direction from the magnetic field produced by a square cell, so Eq. (8) is validly used to calculate various magnetic anomalies. Magnetic inversion modeling is defined by the misfit between measurement and computational data calculated using the following equation: (9) where W d is a diagonal matrix for the standard deviation of noise associated with the whole dataset and d represents predictive data from the magnetic susceptibility model. The inversion problem can be solved by finding the model parameter m that minimizes the misfit of data with a specified amount. The susceptibility model can be recovered through an optimization problem solved by: φ = φ d + μφ m (10) where φ is the objective function and μ is the trade-off parameter to balance the low data misfit that is important to the model's objective function.
Inversion modeling began with the reconstruction of a square block model in the form of a mesh grid based on inversion equations developed by Li and Oldenburg (1996) and Macleod and Ellis (2013).
The main input parameters of the model consist of a geometry model, inclination = 1.19°, declination = -22.21°. The total magnetic intensity values should be about 41430 nT. The angles of inclination and declination used for pole reduction (RTP) are 90° and 0°. Spatial discretization involved 46 W-E directed grid blocks, 46 N-S directed grid blocks, and 8 grid blocks in the z direction ( Figure  6). Magnetic anomaly data are in the form of residual RTP anomaly as the main input as well as regional RTP anomaly and topographic data of the study area.

Bouguer anomaly and analytic signal
Based on the WGM 2012 gravity model, a Bouguer anomaly in Ambon and Seram has a positive value ranging from 73.25 mGal to 493.57 mGal divided into three parts, specifically with a high anomaly in the northwest and extending from east to west in the southern part of the study area (Figure 7a).
Most of this anomaly is in the thin-thick oceanic crust and is thought to be caused by high-density rocks lifted by bedrock or mantle. CRUST 1.0 model information ( Figure  4) supports the thin thickness of the crust, which usually has a thin ocean in the south. Medium anomalies are dominant in Ambon and Seram, whereas low anomalies are found in the north, usually associated with the Seram trench (Figure 7a). In the Seram trench, which is an area of destruction or transformation due to plate subduction, lowgravity anomalies are thought to be related to low-density deep-sea sedimentary rocks. Apart from sedimentary rocks, a strong negative anomaly also causes the water mass deficit in the trench (Lowrie, 2007). Marotta et al. (2006) stated that high positive density anomalies are correlated with interrelated negative temperature anomalies, whereas negative density anomalies may be caused by the negative topography of the surface associated with trenches or crust material trapped in subduction zones. Milsom et al. (2001) showed that the islands of the Banda Sea are generally made up of igneous ophiolite rocks composed of deep-sea sediments, basalt, gabbro, and peridotite. Ophiolite rocks have a high density of rock associated with high variations of gravity that are assumed to be a layer of oceanic crust driven by continental crust obduction. Tjoktrosapoetro et al. (1993) showed the volcanic rocks dominated by andesite, dacite, breccia, and tuff formed at the end of the Pliocene to earlier Pleistocene in the regional geological map of Ambon and its surrounding islands. Granite was also found in some locations in the western part of Ambon as intrusive rocks. Subsequently, horizontal gradient and analytical signal analysis of the Bouguer anomaly was performed to determine the edge of the subsurface structure. On the Ambon and Seram islands, the horizontal gradients of gravity anomaly data are positive, ranging from 0.001 mGal/m to 0.025 mGal/m (Figure 7b), while the amplitudes of the analytical signal of gravity anomaly varied from 5.37 mGal/m to 5.41 mGal/m (Figure 7c). The horizontal gradient map of the study area shows the highest value, reflecting the lithological edge or boundaries. In regard to that, the maximum amplitude of horizontal gradient provides a clear boundary of the "graben" structures that separate Ambon in the southwest to the northeast. The maximum amplitude of the analytical signal of gravity anomaly shows the geological contact from the subsurface structure. The lineament pattern of geological structures found in Ambon and Seram is predominantly influenced by Banda tectonic activity in the eastern part of Indonesia. Audley-Charles et al. (1979) suggested for Seram a tectonostratigraphic scheme involving allochthonous units of Asian similarity thrusting on Australian para-autochthonous sequences. Ultramafic rocks over western and central Seram and Ambon were originally ascribed to an allochthonous thrust sheet and were thus interpreted as comprising part of an extensive ophiolite thrust from the Banda Sea.

Regional and residual anomaly of gravity data
As the main input data in inversion processing, our study uses residual gravity anomaly, so it needs to be separated from the regional anomaly. The separation is implemented using the upward continuation method at an altitude of 4 km from the topographic surface with the continuation height being adjusted to the information of EMAG2-V3 regional magnetic data of the study area. Regional anomalies in Ambon and Seram are positive with a range of 130.36-401.35 mGal (Figure 8a), while residual anomalies are negative to positive, ranging from -48.64 mGal to 41.81 mGal (Figure 8b). Regional anomalies are generally influenced by wide and deep basement rock structures, while residual anomalies describe a local anomaly associated with shallow structures (Khamies and El-Tarras, 2010). Regional anomalies show high values in the northwest and in the southern part of the study area from east to west. The middle value is mostly at Ambon and Seram, while low values were found in the southern part correlated with the Seram trench.
Regional anomaly patterns have similarities with the Bouguer anomalies so that in the oceanic crust, which has a high rock density, the contribution of high values is also present. Anomalies of medium to low gravity tend to be associated with low density of rock, found in the Seram trench. The residual anomaly in the study area is more complex than regional anomalies, with high values in the northwest and south of the study area from the southwest to northeast. In the mountain area, high-gravity anomalies also dominate, whereas low anomalies occur in the basin. All of the variations are caused by changes in the density of rocks in the subsurface associated with lithological or structural interactions.

EMAG2-V3 magnetic anomaly and Curie point depth (CPD)
The latest version of the EMAG2-V3 magnetic anomaly data has been used in our research. In several areas of the world, many geodynamic, tectonic, and geological studies have used data from EMAG2 (Maus et al., 2009;Meyer et al. 2017). The magnetic anomalies consist of total magnetic intensity (TMI) at sea level ( Figure 9a) and an upward continuation of 4 km altitude (Figure 9b). The sea-level TMI extends from -90.14 nT to 127.99 nT, while the upward continuation ranges from -81.21 nT to 98.86 nT. The residual anomaly is created from the difference between sea-level and upward continuation anomalies ranging from -19.61 nT to 32.69 nT (Figure 9c).
There are similarities between sea-level and upward continuation anomalies with the highest circular pattern at Ambon, while a low anomaly is seen in the northwest and south of the study area. High magnetic anomalies could describe the existence of intrusion and outcrops of ultramafic igneous rocks during the lifting process (placement tectonics) under Ambon and Seram, while the low anomaly is connected to sediment rocks in the Seram trough and the northwest. The location of the study area in the southern part of the hemisphere is at low to medium Figure 8. Gravity anomalies after upward continuation at an altitude of 4 km from the topography: (a) regional anomaly, (b) residual anomaly. latitude, so the TMI data need to be reduced to pole (RTP). For TMI data at sea level (Figure 10a), upward continuation (Figure 10b), and residual anomaly (Figure 10c), the RTP process was carried out using Magick V 3.25 software developed by Tchernychev (www.geometrics.com).
Overall, the magnetic anomalies of RTP in the southern part of Ambon are negative to positive with the highest value. At sea level, the RTP magnetic anomalies range from -147.71 nT to 294.79 nT, the upward continuation RTP ranges from -108.47 nT to 198.44 nT, and the residual anomaly values range between -35.50 nT and 74.20 nT. We also evaluated the horizontal gradient and analytical signal of the RTP residual anomaly to determine the geological edge or boundaries in Ambon and Seram. The horizontal gradient of the RTP residual map ranges from 0.001 nT/m to 0.005 nT/m (Figure 11a), while signal amplitude values range from 0.007 nT/m to 0.020 nT/m (Figure 11b).
Horizontal gradient and analytic signal maps have similar patterns that illustrate the contact limits of structures in the study area. Both of these maps show full magnetic interactions over the magnetic sources, indicating strong magnetic susceptibility in the south of Ambon with a circular pattern, which may be correlated with the intrusion of igneous rocks. Magnetic anomalies are also associated with short and long wavelengths due to the deep and shallow magnetic source (Thébault et al., 2010;Saada, 2016). We use the CPD method of RTP regional anomaly to estimate the depth of the magnetic source. The CPD is widely used to determine the thickness of shallow structures, complex basements, and geological subsurface structures (Araffa et al., 2018). The CPD is calculated through the 2D radial power spectrum of the magnetic anomaly using Fourpot software .
The range consists of two components based on the average value of the slope: z 0 shows the centroid depth of the magnetic source, which reaches 17.30 km ( Figure  12a), while z t shows the upper bound of the magnetic source, reaching a depth of 8.15 km (Figure 12b). For the depth of the bottom of the source (z b ), estimated at 27 km, if subtracted by the height of the upward continuation to 4 km, the depth of the Curie point reaches 23 km. These results are consistent with CPD analysis in East and Southeast Asia by Tanaka et al. (1999), in which the Indonesian archipelago's CPD reaches 25 km.

3D gravity inversion
The aim of the inversion of gravity data is to determine the density distribution and shape of the intrusion of igneous rocks under Ambon and Seram, East Indonesia. The inversion process is conducted by optimizing three variables including base anomaly, density, and height of the minor block. The optimization process is carried out on residual and regional gravity data to generate the density distribution influenced by shallow and deep structures. 3D density inversion results with the CRUST 1.0 constraint are indicated by the comparison between observations and calculated gravity data (Figure 13a and 13b), as well as the difference value indicating the error of the inversion (Figure 13c). Figures 13a and 13b show the results of the inversion, giving the similarity of shapes and values between observations and computational data. The root mean square (RMS) error generated from the data is 0.02, while that of the model is 0.003. Figure 13c shows deviation from the inversion process of only about 0.3% (distribution of light blue and yellow). These results indicate that there is a high match between the observations and calculated data. The distribution of rock density below Ambon and Seram is shown in Figure 14.
The average density of the whole model is 2.58 g cm -

3
, with a minimum average density of 1.86 g cm -3 at the surface and a maximum average density of 3.17 gr cm -3 in the basement (Figure 14a). Figure 14b shows the igneous intrusion model under the islands of Ambon and Seram, which is dominated by high-density rocks (yellow to orange), which are thought to be heat-source rocks in the study area. From rock and mineral density tables (Doo et al., 2009;Hunt et al., 2013), bedrock dominating the study area is estimated to be composed of basalt, gabbro, and peridotite to pyroxenite, which are fragments of oceanic crust and upper mantle (Milsom et al., 2001). To reinforce the igneous intrusion model as a source of geothermal energy in the study area, we have created four crosssection profiles each of two west-east trending sections (A-A' and B-B') and two south-north-trending sections (C-C' and D-D'). Cross-section A-A' intersects Buru, Ambon, Haruku, Saparua, and Nusalaut islands ( Figure  15a), B-B' cuts Seram island (Figure 15b), C-C '( Figure  15c) intersects Seram and Ambon islands, and D-D' cuts across Seram and Nusalaut islands (Figure 15d). Each cross-section profile of the inversion shows a fit between the observed data and computation data. The A-A' profile with a west-east direction shows an ultramafic igneous rock intrusion under the Buru, Ambon, Haruku, Saparua, and Nusalaut islands with the highest intrusion just below Ambon at a depth of less than 10 km. This can explain and support the hypothesis about the existence of geothermal potential in the study area. Profile B-B' shows the existence of ultramafic igneous intrusion patterns under the island of Seram at a depth of less than 10 km, as well as the C-C' profile, which cuts across the islands of Ambon and Seram from south to north, allegedly occurring during the lifting process. Profile D-D' crossing the islands of Saparua and Seram from south to north also shows a pattern of ultramafic igneous rock intrusion. The C-C' and D-D' cross-sections also show subduction in the northern part of the associated study area in the Seram trench.
Based on the inversion results, the average depth of bedrock is estimated to range between 20 km and 25 km. This is supported by the results of the CPD analysis of EMAG2-V3 magnetic data reaching 22 km. The whole cross-section clearly shows the intrusion of igneous rocks on Ambon and the surrounding islands. In addition, this result is also supported by the cross-section profile of topography-bathymetry data from GEBCO 2019.

3D magnetic inversion
Magnetic data inversion modeling has also been carried out to support the results of subsurface models from gravity data inversion in Ambon and Seram, eastern Indonesia. The inversion process produces contrast values of rock susceptibility ranging from -0.0034 SI to 0.0058 SI. The subsurface model was constructed in four horizontal layers at depths of 5 km, 10 km, 15 km, and 20 km (Figure 16). The first layer at a depth of 5 km is dominated by rocks with negative to moderate positive contrast, having values ranging from -0.00252 SI to 0.00298 SI. This layer looks more complex with high contrast in the north and low contrast in the south; on the island of Ambon and in its surroundings, there is a low negative susceptibility contrast (yellow circle), which is thought to be the result of intrusion of igneous rocks with high susceptibility contrast, which is rich in magnetic minerals (Figure 16a).
The second layer at a depth of 10 km shows a more homogeneous model with contrast susceptibility ranging  from -0.0023 SI to 0.00347 SI (Figure 16b). High susceptibility contrast dominates in Ambon, northern Seram, and western Buru, which also has many geothermal manifestations on the surface. Figures 16c and 16d show the contrast distribution of the susceptibility of the third and fourth layers at depths of 15 km and 20 km, which is dominated by rock layers with high negative to high positive contrasts, having values ranging from -0.00317 to 0.00483 SI and -0.00339 to 0.00575 SI.
The distribution of susceptibility contrasts is similar to the second layer, where the contrast is high in Ambon and Seram. The third layer and the place also show the existence of a massive igneous intrusion under the islands of Ambon and Seram.
To confirm the igneous rock intrusion under Ambon and the surrounding islands, we have made a cross-section from west to east cutting across the islands of Ambon, Haruku, Saparua, and Nusalaut with the coordinates of the crossing point in the Tulehu geothermal field of 422461.97 E, 9601865.7 S (Figure 17). It can be seen that the prediction of massive igneous rock intrusions occurred just below the island of Ambon and its surroundings, specifically in the regions of Suli and Tulehu.
When associated with the inversion results from gravity data, the results obtained through magnetic inversion modeling are mutually supporting and compatible, in which the gravity results show intrusion patterns of igneous rocks under Ambon and Seram at depths reaching 10 km.

Conclusion
High gravity and magnetic anomalies are thought to have contributed to the intrusion of igneous rock with high density and rich in magnetic minerals on the Ambon and Seram islands, eastern Indonesia. Low gravity and magnetic anomalies are usually associated with low rock densities that have changed magnetic properties of susceptibility and are present in the Seram trench in the destroyed zone. The edges or boundaries of the geological structure in the study area are clearly visible and correlate with the maximum value of the horizontal gradient and analytical signal from gravity and magnetic data. CPD analysis of the RTP regional magnetic anomaly shows the magnetic source depth reaching to 25 km, with the depth of the upper bound of the magnetic source reaching 6 km, while the centroid depth reaches 17 km. The 3D inversion modeling of gravity and magnetic anomalies constrained by the CRUST 1.0 model clearly shows that beneath the Ambon and Seram islands there is a massive igneous rock intrusion reaching depths of less than 10 km. This pattern of intrusion is clearly seen through the cross-section model of the directions of W-E and N-S that cuts the islands of Ambon and Seram perpendicularly. This has implications for the discovery of geothermal energy in the form of hot springs on the surface, where it is assumed that the igneous rock is part of thin oceanic crustal fragments and upper mantle.