XU Chong, XING Junhui, 2), *, GONG Wei, 2), ZHANG Hao, XU Haowei, and XU Xiaoyu

Density Structure of the Papua New Guinea-Solomon Arc Subduction System

XU Chong1), XING Junhui1), 2), *, GONG Wei1), 2), ZHANG Hao1), XU Haowei1), and XU Xiaoyu3)

1),,,266100,2),,266071,3),,250002,

The Papua New Guinea-Solomon (PN-SL) arc is one of the regions with active crustal motions and strong geological actions. Thus, its complex subduction system makes it an ideal laboratory for studying the initiation mechanism of plate subduction. However, the PN-SL subduction system has not yet been sufficiently studied, and its density structure has yet to be revealed. In this paper, we used the free-air gravity data, Parker-Oldenburg density surface inversion method, and the genetic algorithm density inversion method to obtain the density structure of an approximately 1000-km-long northwest-southeast line crossing the PN-SL subduction system under the constraints of the CRUST1.0 global crustal model, onshore seismic data, and the LLNL-G3Dv3 global P-wave velocity model. The density structure shows that density differences between the plates on the two sides of the trench could play a significant role in plate subduction.

Papua New Guinea-Solomon; plate subduction; gravity anomaly; density structure; genetic algorithm

1 Introduction

The PN-SL subduction system is located in the Southwest Pacific Ocean, in the convergence zone of the Pacific and the Indo-Australian Plates (Fig.1). It has been consi- dered one of the most active areas of crustal motion and geological action under the continuous effects of convergence, collision, subduction, and expansion among the Indo- Australian, the Pacific, and the Eurasian Plates since the Cenozoic (Cooper and Taylor, 1985; Holm and Richards, 2013; Yang, 2022). In the PN-SL subduction system, four trenches and three basins are distributed within a spatial area of less than 1000 km (Cooper and Taylor, 1985). Thus, the special geographical location and complex tectonic feature of the PN-SL subduction system make it a natural laboratory for revealing the initiation mechanism of plate subduction (Cooper and Taylor, 1985; Petterson, 1999; Chadwick, 2009).

Thus far, several geophysical experiments have been con- ducted in the PN-SL subduction system. Most of the seismic reflection surveys have focused on the structural features of the Solomon Basin, New Britain Island, the Bismarck Basin, the Woodlark Basin, Island Arc, and the On-tong Java Plateau (., Finlayson, 1972; Gerard and Debeglia, 1975; Inoue, 2008; Santos, 2015), while the structure of the subduction zones in the area have not yet been fully investigated. While seismic tomography studies have portrayed the large-scale plate subduction pat- terns in the study area, they have yet to explore the subduction mechanism and deep mantle motion (., Cordell and Henderson, 1968; Cooper and Taylor, 1985; Hall and Spakman, 2002). Taylor (1979) used magnetic anomaly data to characterize the tectonics and boundaries near the sea basin. Yang(2018) calculated the Bouguer gravity anomaly and the isostasy gravity anomaly in and around the New Britain Trench and Papua New Guinea, as well as analyzed the effect of plate subduction on the isostasy gravity anomaly in the region. There is less work on the density structure of the PN-SL subduction zone using gra- vity data compared with other geophysical methods. However, determining the density structure of the PN-SL subduction system by gravity data is significant, as it can help understand the role of plate density differences in the sub- duction activity of the PN-SL subduction system and promote the study of the initiation mechanism of plate subduction in the PN-SL subduction system.

This paper focuses on the density structure features of the PN-SL subduction system. Following the classical Parker- Oldenburg density surface inversion method, we construc- ted a genetic algorithm (GA) density structure inversion pro- cess under the constraint by combining the CRUST1.0 global crustal model, onshore seismic data, and the LLNL- G3Dv3 global P-wave velocity model. The Moho morpho- logy and density structure of an approximately 1000-km- long 2-D line crossing multiple basins and trenches in the PN-SL subduction system were obtained. The results indicate that the density differences between the plates on the two sides of the trench play a significant role in the plate subduction in this area.

Fig.1 (Inset) Location of the study area with large-scale faults. The yellow box indicates the area shown in the main figure. (Main figure) Layout of the gravity line overlain on a topographic map. The red solid line indicates the gravity line AC, which is divided into the AB and BC sections. The location of the onshore seismic station AU-MANU is indicated by the red-filled triangle. The white filled arrows indicate the plate movement direction. The legend shows different types of faults. The topography data are from SRTM15+ (Tozer et al., 2019). The figure was constructed using GMT 6 (Wessel et al., 2019).

2 Geological Setting

As shown in Fig.1, the PN-SL subduction system is located in the Southwest Pacific Ocean, specifically at the convergence boundary of the Indo-Australian and the Pacific Plates and the eastern end of the Neotethys tectonic domain, which is adjacent to the Caroline Plate to the north, the Eurasian tectonic domain to the west, the Australian Plate to the south, and the largest oceanic plateau, the Ontong Java Plateau, to the east (Wessel and Kroenke, 2000; Stotz., 2017). Since the late Cretaceous, the PN-SL subduction system has evolved through key tectonic periods of 50, 45, 25, 8, and 5 Myr (Schellart, 2006) into the cur- rent west-to-east distribution of the West Melanesian-North Solomon Trench, the Bismarck Basin, the New Britain-San Cristobal Trench, the Solomon Basin, the Trobriand Trench, the Woodlark Basin, and the Pocklington Trench within a spatial area of less than 1000 km.

At present, the Indo-Australian plate is moving at a speed of 90 – 110 mm yr?1(Wallace, 2004) (Fig.1). The convergence of high-speed multiple plates has caused multi-phase and complex magmatic activities in the PN-SL subduction system. Studies about the seismic spatial dis- tribution and tomographic imaging indicate that the Indo- Australian Plate and the Pacific Plate are subducted at high angles of over 70? and depths of over 500 km along the New Britain Trench and the West Melanesian-North Solomon Trench, respectively (Holm and Richards, 2013; Zhang, 2018).

Meanwhile, the Solomon Sea is located between the Indo- Australian and the Pacific Plates, and the presence of rich types of tectonic systems developed in this region makes it a focal area for deep density structure studies. As a post-arc spreading basin formed in the early Tertiary (Musgrave, 2013), the approximately 3 – 5-km-deep Solomon Basin is bounded by complex tectonic systems. On the north side, the Solomon Sea Plate is subducted northwest and northeast along the New Britain Trench, while on the south side, it is subducted beneath the Trobriand Trench. Its southeast boundary is separated from the Woodlark Rise by a transform fault. The New Britain Trench is around 6 – 8.5 km in depth and has a 70? change in trench strike at 153?E. The Trobriand Trench is around 5 – 5.3 km in depth and is around 300 km in length (Schellar, 1987). The Ontong Java Plateau, located on the north-eastern side of the Solomon Sea, was formed by two phases of the Labuan basaltic mag- matism of the mantle plume beneath the Pacific Plate at 122 and 90 Myr. The Ontong Java Plateau rises at about 2000 m above the seafloor and covers an area of 1.9 × 106km2, making this large igneous province the largest oceanic plateau in the world (Hanyu, 2017). To the north of the study area is the Bismarck Basin, which is a back-arc basin with Bismarck arc, a volcanic arc extending from New Britain to Schouten Island (Galewsky and Silver, 1997). The Manus Basin, located in the eastern part of the northern Bismarck Sea, was formed by asymmetric spreading 3.5 Myr ago. The West Melanesian Trench is a subduction zone of the Pacific Plate. The North Solomon Trench, which is also inactive, is connected to the West Melanesian Trench. Along with the subduction of the Ontong Java Plateau, the North Solomon Trench gradually stopped subducting in the Miocene, and in the late Miocene, the back-arc basin gra- dually started subducting along the San Cristobal Trench and the New Britain Trench.

3 Data and Methods

The data used in this paper include bathymetry, gravity, the CRUST1.0 global crustal model, the LLNL-G3Dv3 glo- bal P-wave velocity model, and the IRIS onshore seismic station data on the study area. In this paper, we determined the density structure of the PN-SL subduction system along the Line AC, which crosses the Pocklington Trench, Woodlark Basin, Trobriand Trench, Solomon Basin, New Britain Trench, Bismarck Basin, and West Melanesian Trench from southeast to northwest (Fig.1).

The bathymetry data were derived from the ETOPO1 (ETOPO1 Global Relief Model) of the National Geophy- sical Data Center and contained global absolute land elevation and seafloor topographic data with a data grid spacing of 1? × 1?.

The gravity data were derived from the WGM2012 mo-del (World Gravity Map Model) provided by the Interna-tional Gravimetric Bureau (BGI) with a data grid spacing of 2?× 2?. The topographic correction was conducted using the ETOPO1 model, and the Bouguer correction was calculated using the FA2BOUG code developed by Fullea. (2008). The average density of the Earth’s crust used in the Bourget correction was 2.67 g cm?3.

CRUST1.0 (Global Crustal Model) is a global 3D crustal model with a data grid spacing of 1? × 1? (Laske, 2013). Here, the grid points are defined at the center of the cell, and the data points are stored following the longitude priority criterion,., latitude is cycled first, followed by longitude. CRUST1.0 uses the bathymetry and topography data from the ETOPO1 model. It provides a unique eight- layer crustal profile containing water, ice, upper, middle, and lower sediments, and upper, middle, and lower crust. Moreover, the P-wave velocity, S-wave velocity, and density parameters are given explicitly for these eight layers and the mantle below the Moho.

Meanwhile, LLNL-G3Dv3 is a global P-wave velocity model for the crust and mantle with regional-scale details (Simmons, 2012). The model is parameterized using a spherical tessellation with node spacings of 1? and 2? in the upper and lower mantles, respectively. The model consists of a total of 57 layers from the surface to the core and 1.6 million nodes.

The number of onshore seismic stations with available data in the study area is limited. Fifty-two seismic events (Fig.2) were recorded at the AU-MANU station (Fig.1), a station on the Manus Island near the line AC, and were selected following the criteria of magnitude greater than 5.5, the inclusion of a complete set of three-component data, and the ability to effectively identify first arrivals. The Moho depth beneath station AU-MANU was obtained us- ing the receiver function method proposed by Zhu and Ka- namori (2000).

Considering the large extent and complex geological settings, we used the waveform data of a seismic station in the study area to calculate the Moho depth below it – a complement to the CRUST1.0 model – to further constrain the average depth of the Moho required for the inversion of the Moho morphology along the line AC. In addition, the WGM2012 Bouguer gravity anomaly data were used to invert the Moho morphology by the Parker-Oldenburg density interface inversion method. Due to the lack of publicly available seismic data in the study area, we determined the tectonic boundaries along the profile by combining the LLNL-G3Dv3 model with the hypocenter data and the calculated Moho depth and then used a constrained GA density inversion method to invert the density structures (Fig.3).

Fig.2 Selected 52 seismic event waveforms recorded at station AU-MANU and the stacked waveform (on the top).

Fig.3 Flowchart of the genetic algorithm density structure inversion with the density interface constraint.

The GA was used in this paper for the inversion of the density structure. Normally, for the forward problem, the complex geological bodies are dissected into several tiny rectangular prisms, and the sum of gravity anomalies of all rectangular prisms is then used as the gravity anomaly of the geological bodies. The GA requires at least one forward computation for each iteration, which places high de- mands on the efficiency of the forward computation (., Hinze, 2013). Thus, in the current paper, we introduced polygonal prisms in the calculation of the forward problem, which can accurately portray the shape of the geological anomalies and significantly reduce the time of each iteration.

4 Results

4.1 Moho Depth at Station AU-MANU

Fig.4 shows thestack result of the Moho depth and/ratio below the station AU-MANU. As can be seen, the best-fit Moho depth and/ratio are 19.5 km and 1.43, respectively. This best-fit Moho depth was then used to con- strain the mean Moho depth for the inversion of the Moho morphology along the line AC.

Fig.4 H-k stack result. Moho depth = 19.5 ± 0.0 km, Vp/Vs = 1.43 ± 0.01, Poisson’s ratio = 0.016 ± 0.013.

4.2 Moho Depth Along the Line AC

In the inversion of the Moho depth along the line AC, the density difference above and below the Moho was set to 0.6 g cm?3, with grid spacing set at 10 km. The mean Mo- ho depth for the AB part of the line was set to 16.9 km according to the CRUST1.0 model and the Moho depth at station AU-MANU (Fig.5b), while the mean Moho depth for the BC part was set at 18.3 km according to the CRU- ST1.0 model.

Along the line, the root mean square (RMS) misfit between the calculated and observed gravity anomaly is 8.61 mGal (Fig.5a). Fig.5c shows the comparison between the calculated Moho depth and the Moho depth from CRUST 1.0. In general, two Moho morphologies fit well. The Moho depth shallows from the Pocklington Trench to the West Melanesian Trench but vibrates significantly along the line. At the Pocklington Trench, the Trobriand Trench, the New Brain Trench, and the West Melanesian Trench, the Moho depth is relatively deeper (20 – 25 km) but gets shallower (12 – 15 km) at the Woodlark Basin, the Solomon Basin, and the Bismarck Basin. At 1460 km, where the line is closest to station AU-MANU, the Moho depth is 19.6 km, which is consistent with the Moho depth (19.5 km) below the sta- tion AU-MANU.

4.3 Density Structure Along the Line AC

Since 1900, there have been more than 1000 seismic events with magnitude ≥ 6.0 in the PN-SL region (USGS). The intense seismic activities make it possible to study the preliminary spatial and structural features of the PN-SL sub- duction system. Fig.6 shows the P-wave velocity structure derived from the LLNL-G3Dv3 model, with the distribution of the earthquake hypocenters along the line AC. The Moho depth along the line AC obtained by the Parker-Oldenburg density surface inversion method was used to constrain the Moho depth in Fig.6. The Indo-Australian Plate and the Pacific Plate are subducted along the New Britain Trench and the West Melanesian-North Solomon Trench at a high angle of more than 70?, respectively. The subduction depth of at New Brain Trench is more than 500 km (Cordell and Henderson, 1968).

Fig.7 shows the density inversion result along the line AC of the GA density inversion method with constraints. In particular, Fig.7a presents the calculated and observed gravity anomaly profile. The RMS misfit between the calculated and observed gravity anomalies along the line is 4.31 mGal, which is mainly due to the difficulty of the density blocks’ division considering the complexity of the tectonic structure within the PN-SL subduction system. Compared with the range of the observed anomalies (approximately 300 – 1550 mGal), the RMS misfit is much smaller, indicating the reliability of the results and the GA density inversion method.

Fig.7c shows that the crustal density varies from 2.60 to 3.09 g cm?3and increases near the subduction zone. There is a significant crustal density difference between the two sides of the Trobriand Trench and New Britain Trench sub- duction zone, while there is no such difference between the two sides of the Pocklington Trench and West Melanesian Trench subduction zone. Moreover, there is a density difference of 0.02 g cm?3at the top of the upper mantle between the two sides of the Trobriand Trench and New Britain Trench subduction zone. The density at the top of the upper mantle for the unsubducted plate is 3.26 – 3.34 g cm?3. For the subducted plate, its density gradually increases with the increase of the subduction depth. Furthermore, the subducted plate at the New Britain Trench has the lar- gest subduction depth, and its upper mantle density increases from 3.34 g cm?3at a shallow depth to 3.81 g cm?3at a deep depth.

5 Discussion

5.1 Density Structure Inversion for the PN-SL Subduction System

The PN-SL subduction system is an ideal natural laboratory for studies about plate subduction. However, previous studies mostly used natural seismic tomography to investigate the deep velocity structure of the region, resulting in the deep density structure of the region being relatively underinvestigated. In this paper, we determined the density structure along a 2-D line AC, which crosses the Pocklington Trench, Woodlark Basin, Trobriand Trench, So- lomon Basin, New Britain Trench, Bismarck Basin, and West Melanesian Trench from southeast to northwest. Con- sidering the complex tectonic setting and the lack of geophysical constraints, we established a density structure in- version method consisting of Moho depth inversion using Bouguer gravity anomalies and density structure inversion under the constraints of Moho morphology.

Fig.6 Tectonic boundary delineation along the line AC. Black solid dots indicate the locations of the earthquake hypocenters. The red-filled triangle indicates the location of station AU-MANU.

Fig.7 (a) Comparison between the observed and calculated gravity anomalies along the line AC. (b) Topography along the line AC. The red-filled triangle indicates the location of station AU-MANU. (c) Inversion result of the density structure along the line AC. Red lines indicate the Moho morphology. The numbers in the profile show the densities in g cm?3.

When performing the inversion of the Moho depth using the Parker-Oldenburg method, it is important to extract gra- vity anomalies related to Moho morphology. The multiscale wavelet analysis is a powerful tool for gravity anomaly ex- traction. However, due to the complex tectonic setting and widespread subduction zone, as well as the lack of a priori information in the study area, it is imprecise and unreliable to extract gravity anomalies only using wavelet multiscale analysis. In practice, we found that the original WGM2012 Bouguer gravity anomaly can be used for Moho depth inversion, which can provide detailed high-frequency information in the results and has better consistency with the Moho depth below the station AU-MANU compared with the CRUST1.0 model. Hence, it is feasible to obtain the Moho surface directly from Bouguer gravity anomaly inversion, especially in regions with complex tectonic settings where it is difficult to separate the gravity fields accurately.

In this paper, the GA was only applied to the inversion of the density structure along a 2-D line, and it showed good reliability and potential. Hence, constructing an inversion using GA that can be used for 3-D density structure inversion is highly recommended for future research. In addition, we also expect more seismic experiments to be conducted in the study area, the results of which can be used for the density structure inversion constraints.

5.2 Density Structural Features of the PN-SL Subduction System

Along the line AC, the vertical trend of the density of each subducted plate is similar (Fig.7); that is, as the plate subducts, the density of the plate gradually increases with the depth due to the effect of plate dehydration. For example, the density of the subducted Australian plate at the Pocklington Trench increases from 3.26 g cm?3at a shallow depth to 3.37 g cm?3at a deep depth, while the density of the subducted Solomon Basin at the New Britain Trench in- creases from 3.33 to 3.81 g cm?3. However, the upper mantle density differences between the two sides of the subduction zone at different trenches vary. At the West Melanesian Trench and the Pocklington Trench, where the initial subductions began in the late Cretaceous and in the Eocene, respectively, the corresponding upper mantle density differences between the plates on the two sides are both 0.01 g cm?3. This finding suggests that there is no gravity difference driving the plate subduction and that these two subduction zones tended to die out, which is in good agree- ment with the very low convergence rate and few seismic activities (., Weissel, 1982; Wallace, 2004). At the Trobriand Trench, where the initial subduction began in the late Miocene, the upper mantle density difference between the plates on the two sides is 0.03 g cm?3, which is relatively small and suggests weak subduction activities consistent with the lack of seismic activities. By contrast, at the New Britain Trench, the upper mantle density difference between the plates on the two sides reaches 0.07 g cm?3. The significant density difference between the two sides of the trench induces strong gravitational effects and contributes to the continuation of plate subduction, which is in good agreement with the strong and frequent seismic activities at the trench. Hence, we suggest that the density difference between the two sides of the trench plays an important role in plate subduction.

6 Conclusions

Using the free-air gravity data, Parker-Oldenburg density surface inversion method and GA density inversion method, we obtained the density structure of an approximately 1000- km-long northwest-southeast line crossing the PN-SL subduction system under the constraints of the CRUST1.0 glo- bal crustal model, onshore seismic data, and the LLNL- G3Dv3 global P-wave velocity model. Our major conclusions are as follows:

1) The density structure inversion method constructed on the basis of the GA showed good density structure inversion performance along a 2-D line crossing the PN-SL sub- duction system.

2) The density structure inversion results along the line AC reveal the vertical and horizontal density features of the Pocklington Trench, the Woodlark Basin, the Trobriand Trench, the Solomon Basin, the New Britain Trench, the Bismarck Basin, and the West Melanesian Trench, sugges- ting that the density difference between the two sides of the trench plays an important role in plate subduction.

Acknowledgements

This research was supported by the National Natural Sci- ence Foundation of China (Nos. 91858215, 42076224).

Chadwick, J., Perfit, M., McInnes, B., Kamenov, G., Plank, T., Jonasson, I.,., 2009. Arc lavas on both sides of a trench: Slab window effects at the Solomon Islands triple junction, SW Pacific., 279 (3-4): 293-302, DOI: 10.1016/j.epsl.2009.01.001.

Cooper, P. A., and Taylor, B., 1985. Polarity reversal in the Solomon Islands arc., 314 (6010): 428-430, DOI: 10.1038/ 314428a0.

Cordell, L., and Henderson, R. G., 1968. Iterative three-dimension- al solution of gravity anomaly data using a digital computer., 33 (4): 596-601, DOI: 10.1190/1.1439955.

Finlayson, D. M., Cull, J. P., Wiebenga, W. A., Furumoto, A. S., and Webb, J. P., 1972. New Britain – New Ireland crustal seismic refraction investigations 1967 and 1969., 29 (3): 245-253, DOI: 10.1111/j.1365-246X. 1972.tb06157.x.

Fullea, J., Fernandez, M., and Zeyen, H., 2008. FA2BOURG – A FORTRAN 90 code to compute Bouguer gravity anomalies from gridded free air anomalies: Application to the Atlantic- Mediterrannean transition zone., 34 (12): 1665-1681, DOI: 10.1016/j.cageo.2008.02.018.

Galewsky, J., and Silver, E. A., 1997. Tectonic controls on facies transitions in an oblique collision: The western Solomon Sea, Papua New Guinea., 109 (10): 1266-1278.

Gerard, A., and Debeglia, N., 1975. Automatic three-dimensional modeling for the interpretation of gravity or magnetic anomalies., 40 (6): 1014-1034, DOI: 10.1190/1.1440578.

Hall, R., and Spakman, W., 2002. Subducted slabs beneath the eastern Indonesia-Tonga region: Insights from tomography., 201 (2): 321-336, DOI: 10.1016/S0012-821X(02)00705-7.

Hanyu, T., Tejada, M. L. G., Shimizu, K., Ishizuka, O., Fujii, T., Kimura, J. I.,., 2017. Collision-induced post-plateau volcanism: Evidence from a seamount on Ontong Java Plateau., 294: 87-96, DOI: 10.1016/j.lithos.2017.09.029.

Hinze, W. J., Von Frese, R. R., and Saad, A. H., 2013.. Cambridge University Press. Cambridge, 525pp.

Holm, R. J., and Richards, S. W., 2013. A re-evaluation of arc- continent collision and along-arc variation in the Bismarck Sea region, Papua New Guinea., 60 (5): 605-619, DOI: 10.1080/08120099.2013.824505.

Honza, E., Davies, H. L., Keene, J. B., and Tiffin, D. L., 1987. Plate boundaries and evolution of the Solomon Sea region., 7 (3): 161-168, DOI: 10.1007/BF02238046.

Inoue, H., Coffin, M. F., Nakamura, Y., Mochizuki, K., and Kroenke, L. W., 2008. Intrabasement reflections of the Ontong Java Plateau: Implications for plateau construction: Ontong Java Pla- teau construction., 9 (4): Q04014, DOI: 10.1029/2007GC001780.

Laske, G., Masters, G., Ma, Z., and Pasyanos, M., 2013. Update on CRUST1.0 – A 1-degree global model of Earth’s crust.. Vienna, 2658.

Musgrave, R. J., 2013. Evidence for late Eocene emplacement of the Malaita Terrane, Solomon Islands: Implications for an even larger Ontong Java Nui oceanic plateau: Malaita Terrane and Ontong Java Nui., 118 (6): 2670-2686, DOI: 10.1002/jgrb.50153.

Petterson, M. G., Babbs, T., Neal, C. R., Mahoney, J. J., Saunders, A. D., Duncan, R. A.,., 1999. Geological-tectonic framework of Solomon Islands, SW Pacific: Crustal accretion and growth within an intra-oceanic setting., 301 (1-2): 35-60, DOI: 10.1016/S0040-1951(98)00214-5.

Santos, D. F., Silva, J. B. C., Martins, C. M., dos Santos, R. D. C. S., Ramos, L. C., and de Araújo, A. C. M., 2015. Efficient gravity inversion of discontinuous basement relief., 80 (4): G95-G106, DOI: 10.1190/geo2014-0513.1.

Schellart, W. P., Lister, G. S., and Toy, V. G., 2006. A late Cretaceous and Cenozoic reconstruction of the Southwest Pacific region: Tectonics controlled by subduction and slab rollback processes., 76 (3-4): 191-233, DOI: 10. 1016/j.earscirev.2006.01.002.

Simmons, N. A., Myers, S. C., Johannesson, G., and Matzel, E., 2012. LLNL-G3Dv3: Global P wave tomography model for improved regional and teleseismic travel time prediction., 117 (B10): B10302, DOI: 10.1029/2012JB009525.

Stotz, I. L., Iaffaldano, G., and Davies, D. R., 2017. Late Miocene Pacific Plate kinematic change explained with coupled global models of mantle and lithosphere dynamics., 44 (14): 7177-7186, DOI: 10.1002/2017GL07 3920.

Taylor, B., 1979. Bismarck Sea: Evolution of a back-arc basin., 7 (4): 171, DOI: 10.1130/0091-7613(1979)7<171:B SEOAB>2.0.CO;2.

Tozer, B., Sandwell, D. T., Smith, W. H. F., Olson, C., Beale, J. R., and Wessel, P., 2019. Global bathymetry and topography at 15 arc sec: SRTM15+., 6 (10): 1847- 1864, DOI: 10.1029/2019EA000658.

Wallace, L. M., Stevens, C., Silver, E., McCaffrey, R., Loratung, W., Hasiata, S.,., 2004. GPS and seismological constraints on active tectonics and arc-continent collision in Papua New Guinea: Implications for mechanics of microplate rotations in a plate boundary zone., 109 (B5): B05404, DOI: 10.1029/2003JB002481.

Weissel, J. K., Taylor, B., and Karner, G. D., 1982. The opening of the Woodlark Basin, subduction of the Woodlark spreading system, and the evolution of northern Melanesia since mid- Pliocene time., 87 (1-4): 253-277, DOI: 10.10 16/0040-1951(82)90229-3.

Wessel, P., and Kroenke, L. W., 2000. Ontong Java Plateau and late Neogene changes in Pacific Plate motion., 105 (B12): 28255-28277, DOI: 10.1029/2000JB900290.

Wessel, P., Luis, J., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.,., 2019. The generic mapping tools version 6., 20 (11): 5556-5564, DOI: 10.1029/2019GC008515.

Yang, G., Li, Y., Tong, L., Wang, Z., Si, G., Lindagato, P.,., 2022. Natural observations of subduction initiation: Implications for the geodynamic evolution of the paleo-Asian Ocean., 1 (1): 100009, DOI: 10.10 16/j.geogeo.2021.10.004.

Yang, G., Shen, C., Wang, J., Xuan, S., Wu, G., and Tan, H., 2018. Isostatic anomaly characteristics and tectonism of the New Bri- tain Trench and neighboring Papua New Guinea., 9 (5): 404-410, DOI: 10.1016/j.geog.2018.04. 006.

Zhang, Z., Li, S., Tian, J., Yang, P., Ding, D., Zang, Y.,., 2018. Formation mechanism of the moniliform seamounts outside the West Melanesian Trench., 53 (4): 1604-1610, DOI: 10.1002/gj.2979.

Zhu, L., and Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver functions., 105 (B2): 2969-2980, DOI: 10.1029/1999JB900322.

(May 11, 2022;

June 15, 2022;

June 21, 2022)

? Ocean University of China, Science Press and Springer-Verlag GmbH Germany 2023

. E-mail: junhuixing@ouc.edu.cn

(Edited by Chen Wenwen)