Imaging crustal features and Moho depths through enhancements and inversion of gravity data from the Philippine island arc system

The Philippine archipelago is a complex island arc system, where many regions still lack geopotential field studies. The high-resolution isostatic anomaly and free-air anomaly digital grids from the World Gravity Map (WGM) were processed and analyzed to present a general discussion of the Philippines’ gravity signatures and contribute to understanding its regional geology and tectonics. The isostatic anomaly map was continued upward to investigate the high-density ophiolitic basement rocks and low-gravity sedimentary basins at depth. The first vertical derivative (1VD) filter was applied to the free-air anomaly grid map to locate regional structures represented by density contrast boundaries. The depth to the top of the Moho and basement rock over the Sulu Sea was computed using the two-dimensional (2-D) radially averaged power spectrum analysis. Three-dimensional (3-D) gravity inversion was applied to some major sedimentary basins in the Philippines to present 3-D subsurface density contrast models. The interpreted gravity maps highlighted prominent geologic features (e.g., trench manifestation, ophiolite distribution, basin thickness). The negative isostatic anomalies (< 0 mGal) represent the thick sedimentary basins, while the moderate signatures (0 to 80 mGal) correspond to the metamorphic belts. The distinct very high-gravity anomalies (> 80 mGal) typify the ophiolitic basement rocks. The gravity data’s upward continuation revealed contrasting deep gravity signatures; the central Philippines with continental affinity (with 20–35 mGal) was distinguished from the regions with oceanic affinity (with 45–200 mGal). The 1VD map over the Sulu Sea showed anomalies associated with shallow features dominantly related to the Cagayan Ridge. The 2-D radially averaged power spectrum analysis exposed gravity anomalies with tectonic significance (e.g., basement characterization, Moho depth estimation). The estimated average Moho depth in the Sulu Sea is from 12 to 22 km, while the average basement depth is within the range of 5 to 11 km. Lastly, the 3-D subsurface density contrast models characterized the very low-density zones representing the deep (> 7 km) sedimentary basins in the northern Cagayan Valley and southern Central Luzon basins. Furthermore, thin (~ 3.5 km) sedimentary formations are inferred for the low-density areas in northern Agusan-Davao and eastern Cotabato basins.


Introduction
Gravity data are fundamental in understanding and modeling Earth's interior (e.g., subsurface, crust), especially in studying its relationship to geology and structures. With the advancement of technology, high-resolution satellite gravity data are being utilized for geologic exploration and tectonic studies. Satellite gravity data were processed and interpreted for bathymetry prediction (Majumdar and Bhattacharyya 2005), lineament investigation (Braitenberg et al. 2011), crust-mantle boundary study (Steffen et al. 2011), sediment basin survey (Vaish and Pal 2015), and geologic mapping (Pal et al. 2016). This area of research was made possible by acquiring a more precise Earth gravitational model. The Earth Gravitational Model 2008 (EGM2008) is a recent Earth's geopotential field model. This model integrates satellite gravimetry, satellite altimetry, and surface gravity measurements (Pavlis et al. 2008). Several studies already assessed and validated the accuracy of EGM2008 (Arabelos and Tscherning 2010; Pavlis et al. 2012). The gravity field data used in generating the high-resolution World Gravity Map 2012 (WGM2012) are derived from the EGM2008.
In the Philippines, regional gravity exploration began in the twentieth century when Teodoro (1970) compiled Luzon Island gravity surveys. A simple Bouguer anomaly map could not have been generated in those years due to the lack of detailed topographic maps (Teodoro 1970). The Philippines' earliest regional gravity anomaly maps (i.e., free-air, bouger, isostatic) were presented, and different gravity anomalies were discussed relative to various geologic factors by Sonido (1981). Gravity surveys have undergone continuous development during the past twenty years. Ground and marine gravity surveys were employed by several studies focusing on specific regions-e.g., the crustal structure and tectonic evolution along Manila Trench (Hayes and Lewis 1984), the emplacement of Bohol ophiolite (Barretto et al. 2000); the regional tectonics of northern Luzon (Milsom et al. 2006), the arc-continent collision in the central Philippines (Dimalanta et al. 2009), the crustal thickness of Central Philippines , the upper crustal structure beneath Zambales Ophiolite Complex (Salapare et al. 2015), and the terrane boundary in northwest Panay (Gabo et al. 2015).
The Philippines' historical overview of gravity surveys presents a wide range of gravity survey scales and applicability. Earlier studies generated and presented gravity maps based on limited point data from local to regional surveys (e.g., ground, marine). With the advent of satellite-derived gravity data and global gravity dataset, geologic studies' scope is no longer limited to the previously available point data. This study used the isostatic and free-air anomaly grids from WGM to investigate the gravity anomalies in the Philippine island arc system and to identify regional geologic and tectonic features (e.g., structures, sedimentary basins, and basement rocks) including those which may have not been previously delineated on gravity anomaly maps. Two-dimensional (2-D) radially averaged power spectrum analysis was employed to calculate new depth estimates to the top of the Moho and basement rock over the Sulu Sea. Lastly, three-dimensional (3-D) gravity inversion was applied to characterize the density contrast of subsurface geology of major sedimentary basins in the Philippines. This work offers an efficient and effective means of understanding the Philippines' regional geology and tectonics through gravity signatures.

Tectonic and geologic setting
The Philippine island arc system is a complex and tectonically active region. It is characterized by ophiolite accretion, arc magmatism, ocean basin closure, and other tectonic processes (Mitchell et al. 1986;Rangin 1991;Yumul et al. 2008a;Aurelio et al. 2013). The Philippine archipelago consists of two general terranes: the Palawan-Mindoro Microcontinental Block and the Philippine Mobile Belt (PMB) (Mitchell et al. 1986;Rangin 1991;Yumul et al. 2008a) (Fig. 1). The Palawan-Mindoro Microcontinental Block was once part of mainland Asia, while the PMB originated from the sub-equatorial regions (Taylor and Hayes 1980;Holloway 1981;Rangin et al. 1990;MGB 2010). The PMB is an actively deforming zone between two oppositely dipping subduction systems (Mitchell et al. 1986;Rangin 1991;Yumul et al. 2008a;Aurelio et al. 2013) (Fig. 1). The eastern side of the PMB is bounded by the west-dipping East Luzon Trough and the Philippine Trench (e.g., Hamburger et al. 1983;Yumul et al. 2008a). The archipelago's western side is marked by southeast-dipping subduction zones: Manila Trench, Negros Trench, Sulu Trench, and Cotabato Trench (e.g., Mitchell et al. 1986;Yumul et al. 2008a). The bathymetric map of the Philippines shows very deep features (~ 8 km to 9 km) corresponding to the mentioned subduction zones (Fig. 2). Generally, the eastern side of the Philippines has a deeper bathymetry than the western flank. The left-lateral strike-slip Philippine Fault, which traverses the entire island arc system, accommodates the oblique convergence between the Philippine Sea Plate and Eurasian Plate (Barrier et al. 1991;Aurelio 2000). The amalgamation of different terranes paved the way to form tectonic collages with diverse lithologic characteristics categorized into ophiolitic rocks, metamorphic rocks, magmatic arcs, and sedimentary basins (MGB 2010). Ophiolitic and metamorphic basement rocks overprinted by relatively younger volcanic series and thick sedimentary basins define the Philippines' present geology.

Gravity signatures over the Philippine island arc system
The isostatic gravity anomaly map extracted from the WGM reveals the Philippines' regional geologic and tectonic features (Fig. 3). The PMB (red to pink; ~ 10 to 100 mGal) is generally surrounded by very low negative anomalies (blue; ~ 0 to − 70 mGal) attributed to the deep trenches and troughs that bound the archipelago. Areas underlain by denser materials reflect more positive anomalies, while lower density zones generate more negative signatures (e.g., Lowrie and Fichtner 2019).
The broad north-trending negative anomaly characterizes Luzon Island's eastern offshore area; it represents the East Luzon Trough (ELT) forearc basin. The very low negative anomaly zones (< − 70 mGal) correspond to the thick accumulation of sediments, as confirmed by previous seismic and bathymetric surveys (Fig. 3b). Lewis and Hayes Page 5 of 20 Casulla et al. Progress in Earth and Planetary Science (2022) 9:16 (1983) defined the plate boundary along eastern Luzon as a young active zone that decreases its activity toward the north. They also noted that gravity signatures do not follow the ELT's trend, except the low anomalies south of 17° latitude. This observation is consistent with the trace of the ELT gravity anomaly (which propagates to the northeast) that was delineated in the isostatic anomaly map (Fig. 3a).
The discrepancy between the active tectonic zone and the ELT suggested that the ELT exemplifies a portion of past subduction episodes (Lewis and Hayes 1983). The northern, very low-gravity zone (< − 70 mGal) was identified as the Sierra Madre Basin (SMB), with a maximum sediment thickness of 4.5 km (Lewis and Hayes 1983). The remnant of the Oligocene subduction zone, the Isabella Ridge, was also delineated in the isostatic gravity anomaly map (Lewis and Hayes 1983). Section A-A' on the map presents the pattern of the negative Sierra Madre Basin (forearc basin), positive Isabella Ridge (subduction complex), and moderate gravity anomaly East Luzon Trough (trench) (Fig. 3b). The Benham Rise, which is interpreted as a large igneous province (Barretto et al. 2020), is exemplified by the circular high-gravity region on the gravity anomaly map (BR in Fig. 3a). To the south of the ELT, very low negative gravity anomaly zones may indicate very thick sediment accumulations, which defines the active tectonism along eastern Luzon (e.g., Lewis and Hayes 1983). The gravity anomalies along the ELT are in agreement with the known geology and tectonics delineated by previous studies (e.g., Lewis and Hayes 1983;Barretto et al. 2020). The Philippine Trench is described as a young subduction system with an accretionary prism that disappears toward the Mindanao area (Cardwell et al. 1980;Karig et al. 1986). The isostatic anomaly map showed varying gravity signatures along the Philippine Trench; the northern part (P1) has a higher-gravity anomaly than the southern portion (P2). The very low negative gravity zone (< − 70 mGal) along the southern part of the Philippine Trench system may correspond to very thick sediment accumulation along the forearc basin. The positive low (green) anomaly zones (0 to 20 mGal) (P1 in Fig. 3a) represent the elevated higher-density mantle rocks (seaward) and thinning of sedimentary deposits (landward) (e.g., Lewis and Hayes 1983;Lowrie and Fichtner 2019). The bathymetric anomalies defining the Philippine Trench do not coincide well with the gravity low anomalies indicating a heterogenous subduction zone morphology. Furthermore, a discontinuous feature (L1) was also delineated on the gravity anomaly map, which may represent a structural boundary between the P1 and P2 (Fig. 3a). This gap can also be recognized as linear bathymetric low in Fig. 2.
The east-dipping Manila Trench shows non-uniform negative gravity anomalies that generally correspond to the variable thicknesses of sedimentary deposits overlying basement rocks. Hayes and Lewis (1984) reported that the Manila trench's forearc basins have a maximum sediment thickness of 4.5 km. They also suggested that the thickness variation in the forearc basin is due to sediment accumulation and the accretionary prism's local uplift rate. The distinct negative gravity values (< − 40 mGal) on the northern (M1) and southern (M3) portions of the Manila Trench represent a balance between the local accumulation of sediments and the uplift rate of accretionary prisms (Fig. 3a). In contrast, the absence of very low negative gravity anomaly values in the central part (M2) corresponds to the lower rate of local sediment accumulation relative to the rate of accretionary prism uplift (complex forearc) (Hayes and Lewis 1984). The very low negative gravity anomalies (< − 70 mGal) at the northern and southern portions of the Manila Trench correspond to the very thick sediment deposits; high sediment supply comes from the collision zones of Taiwan-Eurasia (north) and the Palawan-Mindoro Microcontinental Block (south) (Hayes and Lewis 1984). Recent seismic reflection and bathymetry data confirm the deformational patterns in the forearc region of the Manila Trench (Armada et al. 2020). The gravity anomalies along the Manila Trench (i.e., M1, M2, M3) ( Fig. 3a) correlate well with the bathymetric and seismic data from earlier studies (e.g., Hayes and Lewis 1984;Armada et al. 2020). However, additional regional geologic and tectonic features can be recognized in the gravity anomaly map. The very high and contiguous gravity anomaly along the offshore western Luzon Island (Fig. 3a) was interpreted as the extension of the Zambales Ophiolite Complex. Furthermore, the gravity anomaly map shows a gap (L2) that divides the M3 (Fig. 3a). This west-northwest linear feature can also be recognized in the bathymetric map ( Fig. 2) and may correspond to a regional tectonic feature related to the Philippine Fault Zone. The Negros, Sulu, and Cotabato Trenches have prominent gravity lows generally associated with thick lowdensity sediments. Based on the generally expected gravity anomaly signatures of trenches and the previously defined correlation between the processed isostatic gravity anomaly map and detailed ground surveys, the complex forearc basin system of these three trenches (i.e., Negros, Sulu, Cotabato) can be understood. The peculiar very low-gravity zones were noted at the intersection of Negros and Sulu Trenches (NS) and the southern part of the Cotabato Trench (C) (Fig. 2a). Since there are no detailed studies about these three trenches, we can deduce that the very low-gravity zones may suggest a very thick accumulation of sediments; these may indicate active local tectonics along the negative zones.
The isostatic gravity anomalies correlate with the subsurface geology, sedimentary basins, and basement rocks of the Philippines. The map reflects the variations of gravity fields caused by density differences of materials in the upper crust. Based on the lithologic classification from published works (e.g., Aurelio and Peña 2002;MGB 2010), regional groupings were overlaid on the gravity anomaly map (Fig. 4). The summary of the regional geologic groupings concerning the gravity anomaly map is presented in Table 1. Generally, negative gravity signatures represent the sedimentary basins (< 0 mGal), moderate gravity anomalies correspond to the metamorphic rocks (0 to 80 mGal), and very high-gravity anomalies typify ophiolitic basement rocks (> 80 mGal).
The gravity anomaly map also correlates well with the three major basins of the Philippines, namely, Ilocos-Central Luzon Basin (ICL), Cagayan Valley Basin (CV), and Agusan-Davao Basin (AD) (Fig. 4a). These sedimentary basins have distinct and defined north-trending negative anomalies (< − 20 mGal). The isostatic gravity anomaly map only shows negative gravity anomalies on Circular gravity lows were also delineated across the Bohol Sea (BS), signifying a very thick sediment accumulation. This feature was previously interpreted as the proto-Southeast Bohol Trench that bound the Western Visayan Block (Yumul et al. 2008b). The distribution of metamorphic rocks generally coincides with moderate gravity anomaly values (0 to 80 mGal) (Fig. 4a). MGB (2010) classified metamorphic rocks into Pre-Cretaceous (continental) and Cretaceous (island arc) metamorphic zones. Pre-Cretaceous metamorphic zones in the east-central Philippines (i.e., Northern Palawan-Mindoro, Antique Range) are represented by lower gravity anomalies (0 to 30 mGal). The Cretaceous metamorphic rocks, which are sparsely distributed in eastern Luzon, Visayas, and Mindanao islands, have higher-gravity signatures (30 to 60 mGal). The Cretaceous zones are characterized by mafic-to-ultramafic rocks (Aurelio and Peña 2002;MGB 2010). The exception to the positive correlation between the moderate gravity signatures and metamorphic rocks are those areas that are dominantly underlain by ophiolitic rocks. The very high-gravity anomaly signature of ophiolitic rocks masks the gravity lows that represent the metamorphic regions. The documented metamorphosed ophiolitic rocks along the eastern Luzon (Geary et al. 1988;Billedo 1994) and eastern Mindanao (Pubellier et al. 1991;Quebral 1994) are consistent with these observations.
The regional groupings of ophiolitic rocks, delineated by previous works (e.g., Aurelio and Peña 2002;MGB 2010), coincide well with areas having very highgravity anomalies (> 70mGal) (Fig. 4b). The occurrence of ophiolitic rocks, which serve as basement rocks of most islands, is extensive within the Philippines. Lowergravity anomalies are due to metamorphism in some ophiolitic zones (e.g., southeastern Luzon). Among the identified ophiolitic regions, the gravity anomaly map presents cluster of very high-gravity zones. These clustered regions have distinguishable massive outcrops of ultramafic rocks: northern Luzon (Ilocos Norte Ophiolite), western Luzon (Zambales Ophiolite Complex), eastern Luzon (Isabela-Aurora Ophiolite), southern Palawan (Palawan Ophiolite), Samar-eastern Mindanao (NE Leyte, Samar, SW Leyte, Dinagat, Surigao, Pujada ophiolites), Central Mindanao (Central Mindanao ophiolites), and western Mindanao (Zamboanga Ophiolite). These regions were described in McCabe et al. (1982), Schweller et al. (1984), Rangin et al. (1985), Mitchell et al. (1986) and MGB (2010). This gravity information is essential in constraining geologic boundaries in areas difficult to map due to inaccessibility. The massive occurrence of ultramafic rocks resulted in highly positive anomalies in these regions. This has previously been reported for some Philippine ophiolites such as the Southeast Bohol Ophiolite (Barretto et al. 2000) and the Zambales Ophiolite Complex (Salapare et al. 2015). Gravity-high anomalies related to ophiolites were also observed in some studies in the USA (Godfrey et al. 1997) and Oman (Manghnani and Coleman 1981). Complete ophiolite suites were also reported in some areas (i.e., Zambales, Isabela, southern Palawan, Pujada). These wide-ranging and regional gravity signatures provide a better picture of the complex Philippine island arc system in correlation with available ground data.

Data used
The Philippines' isostatic anomaly digital grid was acquired from the World Gravity Map (WGM) of the Bureau Gravimetrique International (BGI) (2012). The BGI produced global gravity anomaly maps and digital grids considering an Earth model that accounts for the influence of most surface masses (e.g., atmosphere, land, oceans, lakes) . Different corrections were applied to the gravity data to remove the non-geologic effects; three WGM anomaly maps were produced (i.e., surface free air, Bouguer, isostatic) by BGI taking into account the elevation data from ETOPO1 Global relief . The gravity anomalies were computed based on the spherical geometry of the isostatic equilibrium (Airy-Heiskanen model). In this computation, the effects of deep isostatic roots and anti-roots were removed ). Thus, the isostatic anomaly map shows the gravity anomalies that correspond to the geologic features in the upper crust (Simpson et al. 1985;Lowrie and Fichtner 2019). The isostatic anomaly grid has a gravity dataset with a 1-min by 1-min (1′ by 1′) spatial resolution . The high-resolution isostatic anomaly digital grid of WGM was processed to reveal the Philippines' geologic structures and features from surface to upper crustal depths. Nonlinear color zoning was used to generate the Philippines' isostatic anomaly map to represent the wide range of grid values (-280 to 200 mGal) efficiently. The elevation data were downloaded from the 30-m Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM). The bathymetric map of the Philippines was modified from the General Bathymetric Chart of the Oceans (GEBCO) (2020). The distribution of the general geologic groupings was outlined from the 'Geology of the Philippines' (e.g., Aurelio and Peña 2002;MGB 2010). ArcGIS software was used to register and overlay secondary data and visualize the features related to gravity anomaly. Geosoft Oasis Montaj software was utilized to process, filter, analyze, and generate gravity anomaly maps.

Grid enhancement processes
Upward continuation: In understanding deeper largescale crustal features, gravity anomalies due to smaller local small structures are less important than the regional anomalies. The regional gravity signals can be enhanced to focus on deeper features (Lowrie and Fichtner 2019). Since deep and large bodies produced long-wavelength and broad anomalies, upward continuation was applied to smooth out near-surface effects (e.g., Nabighian et al. 2005). The isostatic gravity anomaly was continued upward to investigate the Philippines' density distribution according to depths. The digital grid was processed by applying an upward continuation filter at 5 and 20 km distances. The upward continuation estimates and emphasizes the gravity anomaly at a minimum depth of half of the input filter (e.g., 5 km filter = 2.5 km minimum depth) (Jacobsen 1987). The upward continuation was implemented to investigate the high-density ophiolitic basement rocks and low-gravity sedimentary basins at depth (Figs. 5 and 6). First vertical derivative: The first vertical derivative (1VD) filter was applied to the FAA grid map to highlight the edges of gravity anomalies in the Sulu Sea (Fig. 7b). This grid enhancement process helped locate regional structures represented by density contrast boundaries. The first vertical derivative presents the rate of change of the gravity field in a vertical direction. The resolution of the short-wavelength anomalies is significantly enhanced, and the regional (long-wavelength) gravity field signal is attenuated when the first vertical derivative is applied (Nabighian et al. 2005). The shallow near-vertical contacts of the subsurface bodies in the Sulu Sea are represented by the zones that correspond to the zero value. The resulting map was correlated with the Sulu Sea upper crustal structures. The 1VD was only performed over the Sulu Sea to present an example of a supplementary way of deducing large-scale geologic structures and features, especially in remote areas with complex tectonics and scarce ground data.

Two-dimensional (2-D) radially averaged power spectrum analysis
The depth to the top of the interfaces with varying regional density contrast (e.g., Moho, basement rock) in the Sulu Sea was estimated using the 2-D radially averaged power spectrum technique. The power spectrum radial value of potential field data can be obtained using a Fourier transform (Spector and Grant 1970;Mishra and Pedersen 1982;Tselentis et al. 1988 Kunnummal and Anand 2019). In this study, the radial power spectra were prepared using the forward fast Fourier transform (FFT) function of the Geosoft software. The natural logarithms of the spectral values were plotted against the wavenumbers. The slope of each segment gives a proportional average depth of the density contrast interface (e.g., Tselentis et al. 1988). The WGM2012 FAA grid in the Sulu Sea region was used in this process to delineate the estimated depths of the abrupt density discontinuities at deeper levels. The Sulu Sea was chosen as the representative area to apply the 2-D radially averaged power spectrum analysis for characterizing the crustal features and estimating the Moho depth because of its complex tectonics (Murauchi et al. 1973;Kudrass et al. 1986;Rangin 1989;Schlüter et al. 1996;Liu et al. 2014). Firstly, the optimum sub-window size was identified by analyzing power spectra of different grid cell sizes of 0.5° by 0.5° (~ 55 km by 55 km), 1° by 1° (~ 110 km by 110 km), and 2° by 2° (~ 220 km by 220 km). Since the geophysical data available in the Sulu Sea are scarce, an initial estimate of the depth to Moho was obtained using the Crust 1.0 model (Laske et al. 2013) and some limited seismic reflection and refraction data (Murauchi et al. 1973). The Crust 1.0 model presents a range of 18 to 23 km Moho depth over the Sulu Sea, while the seismic refraction measurements estimate a crustal thickness of around 11 km along the NW Sulu Basin. By comparing the calculated depths of the first segment (at different sub-window sizes) using the 2-D spectral analysis and the initial estimate of Moho depths, a good correlation was obtained for a sub-window dimension of 1° by 1°. A total of 60 sub-windows having 50% overlap in all directions were generated. The spectral graph can be generally subdivided into three different slopes that correspond to the average depth of the interfaces. The first segment (lowest wavenumber) represents the deeper depth. The second segment corresponds to the shallower depth; the slope with the highest wavenumber is usually due to the shallow noise.

Three-dimensional (3-D) gravity inversion
Several studies about the geology of the sedimentary basins of the Philippines have been published focusing on drilling data and seismic profiles (e.g., Caagusan 1978;Tamesis 1981;Saldivar-Sali et al. 1986). These previous studies presented results and geologic cross sections with higher resolution (e.g., basin thickness). However, since not all of the portions of the basins were assessed using seismic and drilling surveys, the wide-ranging coverage of the WGM2012 dataset can fill the gaps or areas without available ground data. The 3-D gravity inversion program was originally developed by the co-author, Dr. Hideki Mizunaga. The program is based on the theory of inverting surface gravity data to retrieve a 3-D distribution of density contrast (Li and Oldenburg 1998). The  et al. Progress in Earth and Planetary Science (2022) 9:16 coarseness of the 3-D inversion result presents a limitation for local and detailed interpretation. Yet, these 3-D models provide supplementary information for a largescale density estimate of the upper crustal zones. The 3-D density contrast models of four large representative basins were estimated by applying a linear solution approach. The solution was determined using a linear relationship between the shape/dimension variation and density contrast. A large rectangular superblock was discretized into small rectangular blocks to enable the numerical solution of the inverse problem of the gravity data. A total of 638 gravity data points from the isostatic anomaly map of the WGM 2012 were used to generate a synthetic model of 8,000 blocks. The subsurface was divided into 50 by 40 by 4 in X, Y, and Y directions. The program employed iterations using a least-square solution to minimize the objective function (Haáz 1953). An L2-norm regularization was used with varying smoothness factors to minimize the non-unique solution. Several inversion processes were performed using different smoothness factors. The available geologic cross sections were used to constrain the 3-D models. The smoothness factor that generated the best match with the geologic cross sections (Tamesis 1981) based on drill hole and seismic reflection data was used for the 3-D gravity inversion. A parameter that minimizes the sum of the differences between the observed and calculated gravity anomalies (residual function) should be obtained to solve the inverse problem. The best values of the parameters were calculated by minimizing the sum of squares of the residuals between the observed and calculated anomalies.

Basement rocks and basins
The Philippines' upward continuation maps show that the very high-gravity anomalies (> 75 mGal), associated with the dense features, are distributed in Luzon, southern Visayas islands, Mindanao, and southern Palawan (Fig. 5). The 5 km upward continuation delineates areas underlain by very dense ophiolite rocks or may indicate the occurrence of massive magmatic arcs (e.g., Negros, Daguma Range). Very high (> 90 mGal) gravity anomaly signatures coincide with the well-known massive ophiolitic outcrops (e.g., Tamayo et al. 2004;Yumul 2007). The 5 km upward continuation of gravity anomaly can be clustered into four regions: western Luzon, eastern Visayas-Mindanao, western Mindanao, and southern Palawan (Fig. 5a). Very high-gravity anomalies were recognized in western Luzon-representing the Zambales Ophiolite Complex (e.g., Abrajano and Pasteris 1989;Yumul and Dimalanta 1997), and offshore extension of Ilocos Norte Ophiolite northeast of Luzon (e.g., Arai et al. 1997;Pasco et al. 2019). These high-gravity anomalies characterize the dense ultramafic rocks separated by the thick Ilocos-Central Luzon basin. In southern Palawan, the very high-gravity anomaly corresponds to the Palawan Ophiolite (e.g., Rammlmair et al. 1987;Aurelio et al. 2014) perceivable at the eastern offshore of central Palawan. High-gravity signatures of Zamboanga Ophiolite (i.e., Polanco, Titay) ) are apparent in western Mindanao. Finally, the continuous very high-gravity anomalies along the Leyte and Samar islands due to the Tacloban and Samar ophiolites (e.g., Balmater et al. 2015;Guotana et al. 2017) are very prominent on the 5 km upward continuation map (Fig. 5a). The same anomalies are also remarkable along the easternmost Mindanao due to the Dinagat and Surigao ophiolites (Yumul 2007;MGB 2010). Similar to the case in northern Luzon, the signatures of the very high anomaly zones in western Mindanao and northcentral Mindanao (Central Mindanao Ophiolite) are separated by the negative anomaly signature of the ~ 4.5 km thick Agusan Davao Basin (Ranneft et al. 1960). High-gravity signatures of the massive Negros and Daguma magmatic arcs that persist at deeper levels may indicate dense ophiolitic basement rocks. Limited regional studies of southern Mindanao mentioned the occurrence of serpentinized peridotite as part of the Basement Complex of western Mindanao (e.g., Ranneft et al. 1960). After applying the 20 km upward continuation (Fig. 5b), the exceptionally high anomalies (> 90 mGal) are only recognizable in western Luzon and southwest Mindanao. These anomalies indicate that the sources of the signal are located at a deeper level. The persistence of the very high-gravity anomaly values in southern Mindanao may suggest a very massive and dense ophiolitic basement complex. It is also interesting to note that the central Philippines (CP) has generally lower-gravity signatures (20-35 mGal) compared to the distinct very high-gravity values  in Luzon and Mindanao (Fig. 5b). This is a significant indication of dissimilar major basement rocks (i.e., continental and oceanic origins) of the Philippine archipelago, revealed by their characteristic gravity signatures (e.g., Manalo et al. 2015).
In contrast to the high anomaly zones of dense and massive basement rocks, the sedimentary basins manifest a strong negative anomaly due to the mass deficiency of the underlying thick sedimentary rocks and Quaternary alluvium. The gravity anomaly data from Luzon Island were separately presented to understand the range of gravity anomaly values corresponding to the sedimentary basin. Figure 6a shows Luzon Island's gravity signatures after the 5 km upward continuation filtering; negative gravity anomalies characterize the Philippines' two major basins, i.e., Ilocos-Central Luzon Basin (ICL), Cagayan Page 13 of 20 Casulla et al. Progress in Earth and Planetary Science (2022) 9:16 Valley Basin (CV). The Oligo-Miocene magmatic belts generally separate the two basins along Central Cordillera (CC) (MGB 2010) (Fig. 6b). The main north-trending negative anomalies (− 15 to − 37 mGal) are still present until the 20 km upward continued depth (Fig. 6c). The Ilocos-Central Luzon Basin (west) exemplified a larger negative anomaly zone than the Cagayan Valley Basin (east). The maximum thickness of the Oligocene to Pleistocene sedimentary deposits underlying the Ilocos-Central Luzon Basin (14 km) is thicker than the sedimentary sequence in Cagayan Valley Basin (10 km) (Tamesis 1976;Bachman and Lewis 1983). The portions dominated by very low anomalies (< 5 mGal) gave a general indication of where sediment accumulation is thickest within the basins. At a 10 km upward continuation, the lowest gravity anomalies were delineated in the central portion of the Cagayan Valley Basin and the southern part of the Ilocos-Central Luzon Basin. The 20 km upward continuation map shows the negative gravity anomaly diminished in the northern part of the Ilocos-Central Luzon Basin; it implies that the dense basement rock is shallower in the northern Ilocos Region than in south-central Luzon. These new regional processed data provided supplemental knowledge in understanding the Philippines' basins and basements in relation to their gravity signatures.

Crustal features and Moho depths in the Sulu Sea
Large-scale features in the Sulu Sea were characterized using the FAA grid data of the WGM 2012. The FAA gravity map of the Sulu Sea shows values that range from − 86 to 85 mGal (Fig. 7a). This map reveals regional geologic features (with varying density) representing shallow and deep sources. The Sulu Sea has three distinct gravity features, which will be the focus of this section. The very high-gravity anomaly (> 60 mGal) of Cagayan Ridge (CR) subdivides the basin into two: NW Basin (NWB) and SE Basin (SEB) (Fig. 7a). This eastnortheast-trending feature represents a submerged ridge of an extinct volcanic arc, which is generally composed of volcanic rocks (e.g., basalt, andesite) overlain by thick carbonate sediments (Kudrass et al. 1986;Rangin 1989). Gravity low (< − 24 mGal) represents the NWB, while moderate gravity anomalies (20 to 40 mGal) characterize the SEB. The NWB represents a suture zone that marks the collision between the Palawan Continental Block and Cagayan Ridge (Rangin and Silver 1991). There are two prominent hypotheses for the origin of the SEB. The first hypothesis implies that the subduction of the Proto South China Sea (PSCS) triggered the Cagayan arc volcanism and 'SE Basin' back-arc formation. The counter hypothesis suggests that the SEB is a marginal basin opened before the Cagayan arc volcanism (Rangin and Silver 1991).
The first vertical derivative (1VD) map was prepared to delineate upper crustal geologic features in the Sulu Sea (Fig. 7b). The 1VD map generally presents how much the gravity potential changes in the vertical direction (e.g., Nabighian et al. 2005). Steep and semi-vertical features, where potential does not change, are represented by zero or near-zero values. Prominent linear features and some large-scale structures (e.g., vertical fault, fold axis) identified by the previous studies (e.g., Rangin 1989;Schlüter et al. 1996) were overlaid on the map. The 1VD map shows that anomalies are associated with shallow features that are dominantly related to the Cagayan Ridge.
Gravity anomaly features at the central part have a general NE-SW orientation. These structures are related to normal faults as delineated from the seismic and Seabeam data (Rangin 1989). The southeast extension of the CR fracture zones was also traced to the map's southwesternmost portion. These features imply the extension of CR to the southwest as represented by the similar rock types identified in the area (Rangin 1989). The NWB and SEB are relatively quiet in terms of 1VD features. However, some large-scale features bound the northern and southern NWB. The ENE-trending features in the northern part of the NWB are associated with the set of folds (Rangin 1989). In contrast, NE and NNW-trending features are attributed to the normal faults along Balabac Basin (BB) and Banggi Ridge (BR) (Rangin 1989). Two NW-trending features were identified in the northern part SEB. The previous studies did not recognize these features; these areas can be sites of interest for future detailed studies. Generally, the good correlations present us with a supplementary way of deducing large-scale geologic structures and features, especially in remote areas without available ground data.
The result of the 2-D spectral analysis presented an estimated depth of Moho and basement rocks along the three prominent gravity high and low zones. Figure 8 shows representative examples of 2-D radially averaged power spectral curves in the Sulu Sea. The results obtained using the 2-D radially averaged power spectrum analysis revealed that the average depth to the top of the deeper segment ranges from 12 to 23 km (Fig. 9a). On the other hand, the average depth to the top of the shallow segment is within the range of 5 to 11 km (Fig. 9b). The deeper segment is associated with the Moho depth based on the 11 to 23 km depth range obtained from the Crust 1.0 model (Laske et al. 2013) and some available seismic reflection and refraction data (Murauchi et al. 1973). The gravity zone associated with the CR can be attributed to a relatively shallow Moho (~ 14 to 16 km). The NWB is generally represented by deep Moho (~ 19 to 22 km), while the SEB is dominantly characterized by very shallow Moho (~ 12 to 14 km) (Fig. 9a). The shallower segment is interpreted as the depth to the basement rock (5 to 11 km). This estimate was constrained by previous works (Murauchi et al. 1973;Rangin 1989). Based on the seismic refraction and reflection study of Murauchi et al. (1973), the layer at a depth of 6 km (SEB) to 10 km (NWB) may represent a volcanic basement or metasediments. The gravity zone associated with the NEB has a deep basement rock (north: 9 to 11 km, south: 7 to 9 km)   . 9b). This occurrence can be attributed to the thickening of the sedimentary deposits (1 to 7.5 km) toward Palawan (Murauchi et al. 1973). The SEB is dominantly characterized by basement rock with shallow depth to the north (5 to 7 km) and relatively deep (7 to 9 km) to the south. The shallow basement rock is due to a basement rock overlain by thin sedimentary deposits (Murauchi et al. 1973;Rangin 1989). Generally, the southeastern side of the CR has ~ 7.5 km depth to basement rock. However, varying depths to the basement were noted on the northwestern side. Some portions at the north have very shallow basement rock (5 to 7 km), while some of the southern parts have deeper basement rocks (8 to 9 km). That may be attributed to the difference in sedimentary thickness along these sides (Murauchi et al. 1973;Rangin 1989). The three prominent FAA gravity anomalies discussed in this section are attributed to the depths of basement rock and Moho. However, based on the pattern and orientation of the depth maps, the FAA gravity anomaly associated with CR and the two sub-basins: NWB and SEB, are predominantly controlled by Moho undulation. The new estimate of Moho depth and pattern, based on the 2-D spectral analysis along the Sulu Sea, supports the back-arc basin evolution hypothesis of the SE Basin (Rangin and Silver 1991). The deep Moho and basement rock (thicker sedimentary deposits) in NWB can be associated with the thickened crust due to the subduction/ collision of the Proto South China Sea (PSCS) (Rangin and Silver 1991;Liu et al. 2014). This event triggered the volcanism of the Cagayan arc and the opening of the 'SE Basin' oceanic back-arc basin (shallow Moho and basement rock). Although the 2-D spectral analysis scope is limited within the Sulu Sea, the results still present the usefulness and applicability of the WGM 2012 dataset for large-scale mapping and reconnaissance surveys.

Three-dimensional (3-D) gravity inversion in sedimentary basins
This section presents a general discussion of the results of 3-D gravity inversion in representative areas of major sedimentary basins in the Philippines. The result of the 3-D gravity inversion provides additional knowledge by presenting new 3-D density contrast models of the subsurface with extensive coverage (e.g., estimated thickness, morphology). The 3-D density contrast model of the Cagayan Valley Basin (CV) shows a low-density anomaly belt, which has a general trend NNW direction. This low-density zone corresponds to the low-gravity anomaly of sedimentary basin features in CV (Fig. 10a). Figure 10b shows a low-density zone characterized by deep extension at the northern part of the map (7 km depth). These features were not presented in the southern part of the map. The thick low-density zones to the north were correlated to the upper Neogene to lower Neogene sedimentary deposits underlain by older and denser Cretaceous-Paleogene sedimentary rocks (Tamesis 1981) (Fig. 10e). On the contrary, the 3-D density of the Central Luzon Basin (CL) presents deeper and more extensive low-density zones, especially to the southern part of the map (Fig. 10c). Relatively thin low-density zones, representing thin sedimentary deposits, are recognized in the central part of the map (Fig. 10d). The southern part reveals the very thick Quaternary and Neogene deposits (> 7 km) (Fig. 10f ). The very dense zones for both the CV and CL 3-D models represent the basement rocks of the sedimentary basins (e.g., ultramafic complexes, ophiolitic rocks) (e.g., Tamesis 1981). The basin morphologies in the Agusan-Davao Basin (AD) and Cotabato Basin (C) were also mapped using the 3-D gravity inversion results. Figure 11a and b shows relatively thin and narrow low-density zones representing the Pliocene to Miocene sedimentary deposits. The AD is generally oriented NNW and thickens from south to north. The sedimentary basin thickness can be estimated at 3 to 5 km. Figure 11c and d shows maps and sections of Cotabato Basin that are dominantly influenced by positive density contrast. Higher-density zones correspond to the very shallow basement rocks in the area. The very thin (2 to 3 km) low-density zones are interpreted as Neogene sedimentary deposits distributed in the central and eastern parts of the maps. The very low negative density contrast at the eastern part is attributed to the Quaternary pyroclastics in the area. The model's resolution is coarse due to the coarseness of the dataset and the computational capacity of the program. Nonetheless, the large-scale sedimentary thickness distribution and basin morphology can be mapped using the 3-D gravity inversion results. The use of the WGM2012 dataset for basin characterization and 3-D visualization presents a fast and Page 18 of 20 Casulla et al. Progress in Earth and Planetary Science (2022) 9:16 efficient way of understanding the regional geology of vast areas of interest.

Conclusions
Many parts of the Philippine island arc system still lack geopotential field studies; earlier gravity anomaly maps are only based on limited data points. The wide-ranging and complete gravity datasets from the WGM 2012 are very useful in understanding the regional geologic and tectonic features of the Philippine archipelago by providing significant constraints in evaluating the structures from upper to lower crustal depths (e.g., basin, basement, Moho). The investigation of gravity anomalies in the Philippine island arc system presents negative gravity anomaly zones that correspond to the surrounding trenches that bound the PMB and thick sedimentary basins. Areas with moderate gravity anomalies are associated with metamorphic belts and the very high-gravity anomaly regions define the ophiolitic basement rocks. The upward continuation of the gravity data reveals a relatively low gravity (20-35 mGal) in continental central Philippines compared to the gravity highs (45-200 mGal) of the island arc PMB. This study also presents a new estimate of the average Moho (12 to 22 km) and basement (5 to 11 km) depths along the Sulu Sea using the 2-D radially averaged power spectrum analysis. The application of the 3-D gravity inversion in some major sedimentary basins in the Philippines provides an efficient and effective way of characterizing and visualizing the regional subsurface density contrasts. With the availability and proven efficiency of the WGM 2012 data, these techniques can be applied for future structural and geologic explorations in the Philippines and areas with similar geologic and tectonic settings.