Moho Geometry of the Okinawa Trough Based on Gravity Inversion and Its Implications on the Crustal Nature and Tectonic Evolution

The Okinawa Trough (OT) is an incipient back-arc basin, but its crustal nature is still controversial. Gravity inversion along with sediment and lithospheric mantle density modeling are used to map the regional Moho depth and crustal thickness variations of the OT and its adjacent areas. The gravity inversion result shows that the crustal thicknesses are 17–22 km at the northern OT, 11–19 km at the central OT, and 7–19 km at the southern OT. Because of the crust with a thickness larger than 17 km, the slow southward arc movement, and scarce contemporaneous volcanisms, the northern OT should be in the stage of early back-arc extension. All of the moderate crustal thickness, high heat flow, and intense volcanism at the central OT indicate that this region is probably in the transitional stage from the back-arc rifting to the oceanic spreading. A crust that is only 7 km thick, lithosphere strength as low as the mid-ocean ridge, and MORB-similar basalts at the southern OT demonstrate that the southern OT is at the early stage of seafloor spreading.


INTRODUCTION
Crust is the outermost solid shell of the Earth and the crust types vary significantly among different tectonic units, such as continental margins, oceanic basins, and island arcs. The variation of crustal thickness is a critical factor for understanding the processes of continental rifting and breakup, and determining the crust nature for the incipient back-arc basin (Sutra and Manatschal, 2012).
The Okinawa Trough (OT) is a back-arc basin developed under the area of East China Sea. Although a number of studies have been performed to examine its crustal structure (Iwasaki et al., 1990;Nakamura et al., 2003;Gungor et al., 2012;Klingelhoefer et al., 2012;Shang et al., 2017;Qi et al., 2020), it is still controversial if the nature of its crust is continental, transitional, or oceanic   (Figure 1). Iwasaki et al. (1990) and Nakahigashi et al. (2004) suggested a thinned continental crust at a rifting stage for northern OT based on OBS data. Sibuet et al. (1998) had the same suggestion for southern OT. Han et al. (2007), however, suggested that the central and southern OT are at a transitional stage based on the extremely high heat flow and intense volcanic activities. Arai et al. (2017) and some earlier researchers (Lee et al., 1980;Sibuet et al., 1987) suggested that oceanic spreading have already occurred in OT. Liu et al. (2016) and other earlier researchers (Sibuet et al., 1987;Sibuet et al., 1998) reported the founding of the linear magnetic anomalies in central OT. Also, regarding the oceanic spreading, Kimura (1985) gave an average half spreading rate of 2 cm/year to the southern OT since early Pleistocene. The time of initial rifting of the northern and central OT is considered to be at the Middle Miocene (Gungor et al., 2012;Xu et al., 2014), while that of the southern OT is believed to be just in the Quaternary (Wu et al., 2007;Shang et al., 2017). Until now, the references have no answer to why there is so long a delay time between northern and southern OT (Shang et al., 2017).
Many researchers have mapped the Moho of the OT by gravity methods, but their results vary considerably from each other (Hao et al., 2006;Ding et al., 2017;Xuan et al., 2020). In addition, the crustal thickness predicted by gravity method is quite different from the reflection and wide-angle seismic methods in southern OT (Klingelhoefer et al., 2012;Liu et al., 2016). The reason for this difference might be that the low-density anomaly in the mantle due to magma upwelling was not considered in the previous gravity study (Zhou et al., 2001).
Recently, a new alternative gravity inversion method was provided by Bai et al. (2019b), which incorporates sediment and lithospheric mantle density corrections to map the regional Moho topography. Considering the low coverage of deep seismic refraction and broadband seismogram data (Aitken, 2010), here we adopt Bai's method to map new regional Moho topography of OT by combining it with heat flow, OBSs, and lithospheric strength data, to further examine and discuss the crustal nature of the OT and provide new insight into the ongoing processes of the back-arc basin.

GEOLOGIC SETTING
The Philippine Sea Plate is characterized by three seafloor highs, known as the Amami Plateau, the Datio Ridge, and the Oki-Daito Ridge to the north (Nishizawa et al., 2014) and relatively flat topography to the south (Taylor and Andrew, 2004). The Philippine Sea Plate is subducting beneath Ryukyu Arc along Ryukyu Trench (Sibuet et al., 1998) with a current subduction rate from ∼8 to ∼13 cm/year from north to south progressively (Argus et al., 2013). The earthquake data indicate that the dip of the Wadat-Benioff zone of the Philippine subduction plate is 25°-27°at the north and 55°-75°at its south.
OT extends ∼1,200 km in the NE-SW direction along and at the back of Ryukyu Arc. It can be divided into northern, central, and southern OT by the Tokara Fault in the north and the Kerama Fault in the south (Fabbri et al., 2004;Gungor et al., 2012) ( Figure 1). The depth of seafloor of the southern OT is deeper than that of the central and northern OT, suggesting more rapid subsidence in the south. Seismic data suggest that the southern OT is characterized by well-developed symmetric deep faults, while more diffuse rifting occurred in the north (Gungor et al., 2012). It indicates that the forces along the OT are uneven (Doo et al., 2018). The OT ends at the collision zone between the Luzon Arc and the Taiwan Island since 3.5-4.0 Ma, and this collision probably induced the opening of southern OT (Letouzey and Kimura, 1986;Rateb et al., 2017).
OT has an abnormally high heat flow and intense magmatic activity (Zhang et al., 2019). The average heat flow of the OT based on the 348 available measurements from the global heat flow database of the International Heat Flow Commission is 458 mW/m 2 , which is much higher than the global average value of 86 mW/m 2 (Davies, 2013). The extremely high heat flows (>1,000 mW/m 2 ) in the OT are almost distributed in the central axis of OT, where active volcanoes and hydrothermal vents are developed considerably (Figure 1) . The geochemistry studies indicate that magmatic activity is affected obviously by fluids derived from the subducting slab dehydration (Guo et al., 2017). Active magmatism at the back-arc and nonspreading central and northern OT may be induced by subducting of the Datio Ridge and the Amami Plateau (Sibuet et al., 1998;Rateb et al., 2017).
The evolution of the OT can be divided into two (Letouzey and Kimura, 1986;Sibuet et al., 1987;Shang et al., 2017) or three stages (Kimura, 1985;Gungor et al., 2012;Liu et al., 2016). The first stage is the initial rifting, triggered by the back-arc normal faulting. The second stage is the passive extension, triggered by the strike-slip pull-apart process or transtensional NNE-trending faults. The third stage is the initial divergence indicated by the newborn oceanic crustal spreading .

METHODOLOGY AND DATA
The Moho defined as the boundary between crust and mantle is one of the largest density boundaries with the lithosphere (Lai et al., 2016), so gravity method can be used to image the Moho geometry.
The free air gravity anomaly (Δg f ) is given by: where g m is the mantle residual gravity anomaly that reflects Moho undulations; g b and g s are the gravity anomaly induced by seawater and sediment, respectively, when taking continental crust density (2.7 g/cm 3 ) as the background; g t is the anomaly that originates from mantle density perturbations because of the thermal expansion when taking the normal mantle density (3.3 g/ cm 3 ) as the background; g o is the gravity anomaly caused by the other sources, and this anomaly is relatively small and can be ignored (Bai et al., 2014;Kusznir et al., 2018). After the isolation of the mantle residual gravity anomaly, the Moho burial depth can be mapped based on this mantle residual gravity anomaly (g m ) in the frequency domain (Oldenburg, 1974). Seawater density is taken as a constant of 1.03 g/cm 3 , but the density variations of sediment and lithospheric mantle should be modeled in detail. The advantage of this method is that it can remove the gravity effect of density variations of sediments due to compaction and those of lithospheric mantle due to thermal expansion from the observed free air gravity anomaly by density modeling, which is essential for the gravity inversion in the backarc basin with hot mantle upwelling (Bai et al., 2019b).

Estimating Gravity Effect of the Sediment Layer
The sediment gravity effect is unavoidable for the Moho inversion when the sediment layer is thick. The key issues for estimating gravity effect of the sediment layer are the sedimentary thickness and the density variation. The data and method for mapping sediment thickness variations will be explained in Methodology and Data. The relationship between sediment burial depth, porosity, and density can be used to calculate the sediment density (Sawyer, 1985;López-Coto et al., 2013;Bai et al., 2019a). The porosity variation is a function of buried depth (z) as where Φ 0 is the initial sediment porosity and c is an empirically determined constant with a unit of 1/depth. The parameters Φ 0 and c vary with lithology. Based on the drilling data in the Xihu Sag ( Figure 1), Φ 0 0.55 and c 0.45 × 10 −4 m −1 (Zhang et al., 2009). When sedimentary pore is filled by seawater, the sediment density (ρ z ) varying with depth z can be modeled via: where ρ w is the seawater density with a value of 1.03 × 10 3 kg/m 3 and ρ g is grain density with a value of 2.65 × 10 3 kg/m 3 (Sawyer, 1985;López-Coto et al., 2013).

Estimating Gravity Effect of the Lithospheric Mantle
The density perturbations of the lithospheric mantle caused by thermal expansion can generate a large gravity anomaly at young oceanic basins and rifted continental margins (Chappell and Kusznir, 2008). There are different lithospheric mantle temperature modeling methods (McKenzie, 1978;Stein and Stein, 1992;Afonso et al., 2008). The gravity effect of density perturbations due to thermal expansion can be modeled based on the thermal expansion coefficient and the temperature structure (McKenzie et al., 2005). The pure shear model by McKenzie (1978) is adopted here for modeling lithospheric temperature field as the work by Chappell and Kusznir (2008). Crustal age is an important parameter for modeling the temperature structure of lithospheric mantle (Cowie and Kusznir, 2012). The crust of the West Philippine Sea Basin is oceanic, and its age has been interpreted from marine magnetic lineation (Müller et al., 2019). The continental crust is usually much older than that of the oceanic crust, so we assigned the continental crust to have a thermal age of 300 Ma (Currie and Hyndman, 2006). Our tests show that the lithospheric mantle FIGURE 2 | The thickness map of the sedimentary layer in our study area. The grid is created via Kriging interpolation method based on the global marine sediment thickness grid (Straume et al., 2019) and recent interpretation results of seismic reflection data (Fang et al., 2020). The black dashed lines represent the boundaries of the Okinawa Trough.
FIGURE 3 | The free-air gravity anomalies (Sandwell et al., 2014) and the heat flow stations from the database of the International Heat Flow Commission.
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 752488 temperature almost does not change when the crustal age varies larger than 300 Ma. Since the crustal nature is unclear at the central and southern OT, a series of thermal ages will be assigned to examine which can yield the Moho inversion result that can best fit with the seismic interpretation.

Mapping Sediment Thickness Variations
The sediment thicknesses data used in this study are partly from the global sediment thickness model for oceans and marginal basins (Straume et al., 2019) and from recent seismic survey by CGS and reported by Fang et al. (2020). Kriging interpolation method is used to merge these two grids. Figure 2 shows the sediment thickness data created and used in this study. The depocenters distribute in NE trending, such as the Xihu and Jilong sags. The maximum sediment thicknesses of the two sags exceed 11 km. In the OT, the maximum sediment thickness decreased from 8 km at the central OT to only 4-5 km at the northern and southern OT.

Input Data
The free-air gravity anomalies ( Figure 3) are from the 1-minresolution global free-air gravity model based on the altimetry data by the satellite named Geosat, ERS-1, Envisat, GryoSat-2, and Jason-1 (Sandwell et al., 2014). The bathymetric data used in gravimetric correction are from the ETOPO1, which is also a 1 arc-minute global relief model (Amante and Eakins, 2009). The Moho interpretation results of 11 geophysical (mainly OBS) profiles at the OT and the Ryukyu Arc ( Figure 1; Table 1) are collected in order to estimate the gravity inversion uncertainty. The heat flow (Figure 3) can be used to recover the thermal structure of the Earth and the heat flow data are extracted from the global heat flow database of the International Heat Flow Commission.

GRAVITY INVERSION RESULT
According to our tests, when the ages of the southern and central OT are 1 and 10 Ma, respectively, the Moho inversion result can fit with the seismic interpretation best and the RMS between the gravity and seismic results is 2.28 km. The gravity effect of the sediment layer and lithospheric mantle are highly correlated with the sediment thickness and crustal age variations, respectively ( Figure 4). Figure 5 shows our final Moho inversion result and the crustal thickness map, and the later one is based on the Moho inversion, bathymetry, and sediment thickness data. The Moho and crustal thickness are characterized by NE trending lineation, consistent with characteristic of the Ryukyu trench-arc-basin system. The Moho depths are 25-30 km in the East China Sea Shelf Basin and 15-25 km in the OT, respectively ( Figure 5A). The southern OT holds shallowest Moho and thinnest crust in the whole OT. The crustal thickness of the northern OT is greater than that of the central OT ( Figure 5B). One belt with highly thinned crust locates along the Xihu and Jilong sags, extending in the NE direction in the East China Sea Shelf Basin even though there is no Moho shallowing in this region.

UNCERTAINTY ANALYSIS
We evaluate the reliability of our Moho depth estimates derived from inversion of gravity data against previous seismic imaging studies, and discuss the effect of lithospheric mantle temperature on Moho inversion.

Comparison With Seismic Profiles
The Moho depths interpreted from seismic studies are considered to have a better accuracy than gravity inversion. So, seismic interpretation results are always taken as the reference for estimating gravity inversion uncertainty (Chappell and Kusznir, 2008;Bai et al., 2019b). However, note that the seismic interpretation itself also contains uncertainty. For example, at the intersection point between line 7 and line 8 (Figure 1), the Moho depth from line 7 is 19.0 km (Wu et al., 2020), but that from line 8 is 17.9 km (Gao et al., 2006); thus, the difference between them reaches 1.1 km. Figure 6 shows the comparison between the Moho depths from the collected seismic interpretation and our gravity inversion. The root mean square (RMS) between the Moho depths from interpreted seismic profiles and those from our final gravity inversion at the same sampling stations is 2.28 km. However, the RMS is 4.12 km when ignoring gravity effect induced by both sediment and mantle, and that is 3.16 km when ignoring only the mantle gravity effect.

The Effect of Lithospheric Mantle Temperature on Moho Inversion
We set a 0-300 Ma OT crustal age span for testing. At first, each possible age pair for central and southern OT combination will be applied to temperature modeling and gravity inversion by setting age varying step as 10 Ma. When the crustal age of the central and southern OT is 10 and 0 Ma, the gravity inversion can fit with the seismic interpretation best. Then, the age step is reduced to 1 Ma; age varying range for the central OT is narrowed to 20-0 Ma, and that for the southern OT is narrowed to 10-0 Ma. Finally, when the crustal age for the central OT is 10 Ma and that for the southern OT is 1 Ma, the smallest RMS of 2.28 km is derived. Some of the results are listed in Table 2. Please note that, when the lithospheric mantle density perturbations due to temperature variations are ignored for gravity inversion, the RMS will be increased to 2.67 km. Therefore, the lithospheric mantle density modeling based on the age setting for the central and southern OT can improve the Moho inversion accuracy.

TECTONIC IMPLICATIONS
The large-scale Moho and crustal thickness variation trend is an important implication for the tectonic characteristics. There is no corresponding mantle upwelling beneath the Xihu and Jilong sags judging from the Moho depth ( Figure 5A and Figure 7), but the Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 752488 crust here has been thinned considerably ( Figure 5B and Figure 7). It indicates that the attenuation of the Xihu and the Jilong sag is mainly in the upper crust. However, the mantle has upwelled under the OT and there is also obviously high heat flow at the OT. It means that the crustal thinning at the OT is not contributed by the upper crustal necking predominantly. Therefore, the crustal thinning mechanisms are different between the OT and the Xihu-Jilong sags. The crustal thicknesses are 17-22 km at the northern OT, 11-19 km at the central OT, and 7-19 km at the southern OT ( Figure 5B). The crustal thickness variations along the OT is quite similar to the Lau-Havre-Taupo back-arc basins (Kimura, 1985;Sibuet et al., 1987;Yan and Shi, 2014), so probably multi-evolution stages from rifting to spreading have been developed along the OT.

The Northern Okinawa Trough
The Moho depth variations at the northern OT agree with the seismic interpretation results of the OBS line 1, line 2, line 3, and line 4 (references are list in Table 1). Even though the velocity structures obtained from the OBSs data also indicate that the northern OT is a thinned continental crust in an arc rifting domain (Arai et al., 2017), the average crustal thickness here is larger than those of the central and southern OT. We suggest that the subduction of the seafloor highs, such as the Kyushu-Palau Ridge and the Amami Plateau, beneath the northern Ryukyu Arc had hampered the back-arc extension at the northern OT. The conversion from subduction to collision causes fore-arc rotation and also results in plate boundary curvature. In addition, GPS measurements show that the southward movement of the northeast Ryukyu Arc is much slower than that in the southwest (Nakamura et al., 2003). Contemporaneous volcanisms are concentrated on the northern Ryukyu Arc and scarcely occur in the northern OT ( Figure 1). Therefore, the northern OT should be in the stage of early back-arc extension.

The Central Okinawa Trough
Compared with the northern OT, the crust of the central OT has been highly thinned, with the present thickness of 11-19 km and the Moho depths of 16-23 km. The seismic interpretation along FIGURE 6 | Moho depth from gravity inversion vs. that from seismic interpretation. The root mean square (RMS) between them is 2.28 km.  profile location is shown in Figure 5A. The bathymetry is from ETOPO1, the sediment thickness is extracted from the data shown in Figure 2, and the Moho geometry is from our gravity inversion. The heat flow data are extracted from the global heat flow database of the International Heat Flow Commission. The dashed part of the heat flow curve represents that the heat flow here is much higher than 350 mW/m 2 .
FIGURE 8 | The crustal structure and the heat flow along the profile BB′; profile location is shown in Figure 5A. The data sources are the same as those in Figure 7.
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 752488 line 7 (Wu et al., 2020) and our gravity inversion result (Figure 8) show an obvious mantle upwelling beneath the axis of the central OT. In addition, the crust here has many similar characters with the oceanic crust, such as extremely high heat flow (Figure 8), hydrothermal fields, and intense volcanism (Figure 1) (Hao et al., 2004). Unlike only arc volcanism at the northern part of the Ryukyu trench-arc-basin system, the volcanism migrated from the arc to the back-arc region at the central OT. This is called the volcanic arc-rift migration phenomenon (VAMP) (Sibuet et al., 1987). A series of NE-trending faults occurred in the region with VAMP (Gungor et al., 2012). It indicates that the magma upwelling is correlated to the extensional faults. Active rifting structures have been observed from seismic reflection data in the Iheya Graben and the adjacent area at the central OT (Ikegami et al., 2015). However, the crust here is still not thin enough for crust break and further seafloor spreading. In addition, the volcanism here is still dispersed. Similar to the northern OT, the subduction of the Datio Ridge and the Oki-Daito Ridge beneath the central Ryukyu Arc also had hampered the spreading of the central OT. Therefore, the central OT is probably in the transitional stage from back-arc rifting to oceanic spreading, characterized by moderate crustal thinning, high heat flow, and intense magmatic activity. Figure 7B shows that the southern OT has the thinnest average crustal thickness among the three OT sections. However, solely based on the crust thickness, it is difficult to determine whether it is highly thinned continental crust or oceanic crust. We will further discuss this issue from the following three aspects.

The Southern Okinawa Trough
1) From the thermal state. In general, the age of the oceanic crust is obviously younger than that of the continental crust. Our lithospheric mantle density modeling result indicates that the lithosphere thermal state of the southern OT is quite similar to the thermal state of the mid-ocean ridge, and so the crust here is apt to be oceanic. 2) From the plate strength. The effective elastic thickness (Te) of the lithosphere at the southern OT estimated from topography and gravity data is 5-7 km. The minimum effective elastic thickness occurred around the Yaeyama Graben with a thickness of 2.5-3.0 km (Fu et al., 2002), which is similar to the Te value of the mid-ocean ridge (Cochran, 1979). 3) From the rock type. Fresh basalts were collected by TV grab on the western end of the Yaeyama Graben, and this graben holds the thinnest crust of the whole OT (Lai et al., 2016). Furthermore, MORB-similar basalts at the southern OT is considered to be an important evidence for seafloor spreading here (Zong et al., 2016).
Therefore, we suggest that the southern OT is at the early stage of seafloor spreading, especially at the Yaeyama Graben.

Comparisons to Other Similar Tectonic Regions
The progressive variations of deformation style along the rift axis as occurred in the OT can be found in other basins, such as the East Gakkel Ridge-Laptev Sea Margin area in the Arctic Ocean (Franke et al., 2001) and the Woodlark Basin off Papua New Guinea in the western Pacific (Benes et al., 1994). The rifting to drifting transitions has generally been suggested as the result of differentiated lithospheric strength and the transfer or shear zones are always the boundaries between the rifting and drifting regions (Dunbar and Sawyer, 1996;Van Wijk and Blackman, 2005). Since the NW-SE-trending Kerama and Tokara faults, which separate the northern-central-southern OT, are thought as the result of the subduction of the high and buoyant topography in the Philippines Sea Plate (Sibuet et al., 1998;Gungor et al., 2012) and no evidence can demonstrate that pre-rifting OT has varying lithospheric strength, we suggest that the diffuse rifting in the OT is also due to the topographic high subduction.

CONCLUSION
The acoustic sedimentary basement geometry was mapped based on the published sediment thickness grid and recent seismic interpretations. The density variations in the sediment layer and the lithospheric mantle have been modeled for gravity inversion. The crustal thickness of the Ryukyu trench-arc-basin systems was estimated from the high-resolution free-air gravity anomaly.
The variations of the crustal thickness, along with the heat flow, the fault pattern, and the petrology data, indicate that the three sections of the OT are at different back-arc extension stages. The back-arc extension of the northern and central OT had been hampered by the subduction of bathymetry highs in the West Philippine Basin, but the southern section had not.
1) The northern OT holds the thickest crust among the three sections, the slow southward arc movement, and scarce contemporaneous volcanisms. These three facts indicate that the northern OT is in the stage of the early back-arc extension.
2) The central OT holds moderately thinned crustal thickness, extremely high heat flow, and intense volcanism. Therefore, the central OT is probably in the transitional stage from backarc rifting to oceanic spreading. 3) The thinnest crust at the southern OT is only with 7 km thick. In addition, both the lithosphere thermal state and the lithospheric strength of the southern OT are quite similar to those of the mid-ocean ridge; MORB-similar basalts have been found in the southern OT. Therefore, we suggest that the southern OT is at the early stage of seafloor spreading.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.