1 Introduction

The New Guinea-Solomon arc system (PN-SL) exists in the convergent zone between the Pacific Plate and Indo-Australian Plate, surrounded by the Caroline Plate, Australian Plate, and Ontong Java Plateau, which is the largest oceanic plateau on Earth. This complex tectonic and dynamic background causes the PN-SL to exhibit a special structural system characterized by 'four trenches and three basins', where the West Melanesian-North Solomon Trench, Bismarck Sea, New Britain-San Cristobal Trench, Solomon Sea, Trobriand Trench, Woodlark Basin, and Pocklington Trough are distributed in a narrow space within a width of less than 1000 km. The PN-SL is a region with one of the most active crustal movements and most intense geological processes worldwide (Gong et al., 2019a). Its unique geographic location and geological structure make it a natural laboratory for subduction dynamic studies. Since the 1970s, scientists have paid considerable attention to this area to reveal the deep structure and tectonic dynamics of the PN-SL (Taylor, 1979; Crook and Taylor, 1994; Taylor et al., 1995; Mann and Taira, 2004; Wallace et al., 2004, 2014).

Since 1900, more than 1000 earthquakes with magnitudes greater than 6.0 have occurred in the PN-SL (https://earthquake.usgs.gov/earthquakes/search/), and this intense seismic activity enabled the determination of the spatial and structural characteristics of the subduction zone. The work of Taylor (1979) provides an overview of the structural characteristics of the Bismarck Sea, including the sea-floor topography, focal mechanisms, and seafloor magnetic anomaly data, demarcating the boundary between its southern and northern blocks. The pattern of the Wadati-Benioff zone gradually changes northward from the double subduction on the southern side of the Solomon Islands arc to the northward subduction of the Solomon Sea Plate along the New Britain Trench at a high angle (Cooper and Taylor, 1985). The natural epicenter and focal mechanism data demonstrate that the subduction angle of the Solomon Basin along the New Britain Trench reaches 70°, whereas the dip angle of the slab is less than 40° at the Trobriand Trench (Abers and Roecker, 1991; Pegler et al., 1995). A slab tear exists in the subducted Solomon Sea slab between the New Britain and Solomon Islands, as indicated by the epicenters (Holm and Richards, 2013). Moreover, seismic reflection data and tomography have been used to explore the structure of the crust and mantle in the PN-SL. The crustal velocities calculated by the deep seismic reflection data present that the Moho surface gradually deepens from the Bismarck Sea toward the New Britain Island (Finlayson et al., 1972). According to Miura et al. (2002) and Mann and Taira (2004), the seismic reflection data collected in the North Solomon Trench showed that the crustal thickness of the Ontong Java Plateau can reach 33 km. Richardson et al. (2000) presented the Rayleigh-wave tomography of the Ontong Java Plateau using the waveforms of more than 140 earthquakes and observed that the P-wave velocity within the crust-mantle beneath the Ontong Java Plateau is slightly lower than that in the surrounding Pacific area. Zelt et al. (2001) and Abers et al. (2002) described the crustal and deep mantle structure of the western side of the Woodlark Basin using the active and passive source seismic data observed by fixed stations and ocean bottom seismometer array, confirming the westward expansion of the Woodlark Basin in terms of the deep dynamic characteristics of the crust and mantle. The large-scale seismic tomography study in the Indonesia-Tonga area in the study of Hall and Spakman(2002, 2003) demonstrated the rough and deep geometry of the northward subduction of the Australian Plate along the Pocklington Trough and the southward subduction of the Pacific Plate along the West Melanesian Trench.

From the above information, although previous studies have explored the subduction geometry of the PN-SL using seismic information, the detailed structure and movement of the slab in the deep mantle still have not been obtained. In addition, the specific relationship between the subducted Solomon Sea Plate and the early subducted Pacific and Australian Plates is still unclear. To counter these problems, we used the recorded P-wave seismic data to analyze the velocity structure of the PN-SL system using the teleseismic tomography method. We used the abundantly available epicenters and P-wave traveltime data to explore the deep subduction geometry of the PN-SL and provide key insights into the subduction dynamics of the PN-SL.

2 Geological Setting

The PN-SL is located at the intersection of the Indo-Australian Plate and Pacific Plate. In the Cenozoic period, the Pacific Plate subducted beneath the Australian Plate along the West Melanesian and North Solomon Trenches (Schellart and Spakman, 2015). Induced by the subduction of the Pacific Plate, the Australian Plate is characterized by a strong tension deformation, forming a series of back-arc basins, such as Pocklington and Coral Sea Basins (Crawford et al., 2003). In the middle - late Eocene (45 Myr), as a new back-arc spreading basin, the Pocklington Basin began to subduct northward under the Pacific Plate. In the late Eocene (40 - 35 Myr), the Solomon Sea Basin expanded (Schellart et al., 2006). During the Oligocene, the Pock-lington Basin completely subducted and disappeared, and the New Guinea Island collided with the passive continental margin on the northern edge of the Australian Plate. Then, the expansion of the Solomon Sea Basin stopped. From the late Oligocene to the early Miocene, the Ontong Java Plateau and Solomon Islands collided and converged, and the subduction movement of the Pacific Plate along the West Melanesian and North Solomon Trenches gradually decreased (Petterson et al., 1997). At 8 Myr, the subduction activity migrated to the back-arc region, and the subduction of the Solomon Sea Basin along the New Britain Trench initiated (Holm et al., 2016). At 5 Myr, the Woodlark and Bismark Basins expanded, influenced by the subduction of the Solomon Sea Basin (e.g., Schellart et al., 2006; Holm et al., 2016). Affected by the Cenozoic convergence between these two plates, throughout the critical tectonic periods (50, 45, 25, 8, and 5 Myr), the PN-SL gradually formed a complicated subduction system containing a series of subduction zones (Fig.1). The PN-SL includes not only the mature subduction system (over 40 Myr) in the West Melanesian and North Solomon Trenches (Schellart et al., 2006) but also the subduction system (< 8 Myr) of the young Solomon Sea Plate (< 35 Myr) along the New Britain-San Cristobal and Trobriand Trenches (Holm et al., 2016).

The Solomon Sea is a diamond-shaped basin located at the center of the PN-SL and an Eocene-to-early Oligocene back-arc basin (Joshima et al., 1987; Hall and Spakman, 2003). The average water depth of the Solomon Sea is 4000 - 5000 m. It is one of many microplates in the convergent boundary system of the Indo-Australian and Pacific Plates (Li et al., 2018). At 25 Myr ago, the Ontong Java Plateau collided with the Solomon Islands; this event terminated the southwest-dipping subduction of the Pacific Plate. In addition, the Solomon Sea became bounded by the New Britain Trench with a depth of 6000 - 8500 m at its northern margin and the Trobriand Trench with a depth of 5000 m at its southern margin in the late Miocene (Holm et al., 2016). Compared with the high convergence rate of 9 - 11 cm yr^-1 at the New Britain Trench, the oblique subduction convergence rate of the Trobriand Trench was slower (2 cm yr^-1) (Honza et al., 1989; Tregoning et al., 1998; Wallace et al., 2004; DeMets, 2010). With the complete subduction of the Solomon Sea Plate in the western part of the Solomon Sea, the Trobriand and New Britain Trenches gradually merged, resulting in the convergence of New Guinea Island, Finistere Massif, and New Britain Island (Cooper and Taylor, 1987; Pegler et al., 1995; Whitmore et al., 1997; Holm and Richards, 2013). Moreover, the convergence of the New Britain and Trobriand Trenches presents an evident trend of eastward migration (Crook, 1989; Silver et al., 1991). In the eastern part of the Solomon Sea, the New Britain Trench parallels the distribution direction of the Solomon Islands in the southeast direction, whereas the Trobriand Trench merges with the Nubra fault, which is characterized by a right-lateral strike slip (Cameron, 2014). The Woodlark Basin on the southern side of the Nubara fault, which is characterized by an average water depth of approximately 3000 m, is a young sea basin, and it began to spread 6 Myr ago (Little et al., 2007). According to seismological studies (Yoneshima et al., 2005), the Woodlark Basin subducted beneath the Solomon Islands along the San Cristobal Trench, and the spreading center of the basin extended into the region at the northeast side of the Papuan Peninsula. The Woodlark Basin is a special area where the young sea basin exhibits a strong subduction activity with rapid expansion (Weissel et al., 1982). The Pocklington Trough is on the southern side of the Woodlark Basin and represents a relict trench that accommodated the north-dipping subduction of the Australian Plate beneath New Guinea (Holm et al., 2016). During the Oligocene, with the complete subduction of the Pocklington Basin (Schellart et al., 2006), the collision between New Guinea and the passive continental margin of the Australian Plate was initiated. The subduction along the Pocklington Trough gradually became inactive during the mid-late Miocene (12 Myr) (Holm et al., 2015).

The Ontong Java Plateau is located on the eastern side of the PN-SL, and the plateau summit is approximately 2000 - 2500 m above the ocean floor. This plateau was created by two stages of magma activity from 122 Myr to 90 Myr ago and is the world's largest oceanic plateau, with an area of 1.9×10⁶ km² (Mahoney et al., 1993; Taylor, 2006; Hanyu et al., 2017; Zhang et al., 2020). According to the seismic reflection data collected from the area of the Ontong Java Plateau and the area southwest across the North Solomon Trench, the average crustal thickness of the plateau can reach 30 km (a maximum of 38 km), which is 3 - 4 times thicker than the normal oceanic crust (Richardson et al., 2000; Miura et al., 2004; Inoue et al., 2008). The thick oceanic plateau is difficult to submerge in the trench, and the strong buoyancy induced by the thick crust-mantle of the plateau has a significant retarding effect on the plate subduction process, which can promote the collision between the plateau and the PN-SL.

The Bismarck Sea, with an average water depth of nearly 2000 m, is the back-arc basin corresponding to the New Britain Trench. The record of submarine magnetic anomalies in the Bismarck Sea indicates a seafloor spreading that has been occurring since the Pliocene (Holm et al., 2016). The Bismarck Sea back-arc basin is composed of two microplates (North Bismarck and South Bismarck). The active center of the basin is now concentrated in the Manus Basin, which has expanded asymmetrically since 3.5 Myr ago (Taylor, 1979; Both et al., 1986; Binns and Scott, 1993; Zhao et al., 2014; Ma et al., 2017). The West Melanesian Trench is on the north side of the Bismarck Sea and connected to the North Solomon Trench, which is the subduction area of the Pacific Plate during the early Cenozoic. The West Melanesian-North Solomon Trench is currently inactive or slightly active and closely related to the subduction of the Ontong Java Plateau (Schellart et al., 2006; Gong et al., 2019a).

3 Data and Methods

The P-wave traveltime data used in this study were obtained from the seismic catalog of the International Seismological Centre (ISC) from January 1960 to January 2017 and corrected using the Engdahl-van der Hilst-Buland method (Engdahl et al., 1998). The range of the inversion area is 138° - 162°E and 0° - 12°S (black rectangle in Fig.2). As published on the ISC website, this region has 82 seismic stations, and their distribution is shown in Fig.2. Given the absence of stations in the ocean, the spatial distribution of the seismic stations is limited by the distribution of islands in the study area.

The teleseismic P-wave travel times used in this study were those from the seismic traveltime catalog published on the ISC website. The criteria for earthquake selection were as follows: 1) the epicentral distance from the center of the study area is between 30° and 90° to avoid influences on the travel time of the seismic rays in the core-mantle boundary and the upper mantle outside the study region (Zhao et al., 1994, 2013); 2) the magnitude of the earthquake is greater than 5.0; 3) the number of seismic stations in the study region recording each seismic event is greater than or equal to five. In accordance with the above criteria, 2011 teleseismic events and 15009 P-wave arrival times were selected. Fig.3 shows the distribution of these seismic events, which were distributed around the study region, which has a good back azimuth coverage.

The one-dimensional (1-D) IASP91 velocity model was used to calculate the theoretical travel time of the teleseismic P-waves (Fig.4) and obtain the absolute traveltime residuals of each ray. In teleseismic tomography, relative traveltime residuals are generally used to avoid the influences of hypocenter mislocation, origin times, and complex structures outside the study region (Zhao et al., 1994), that is, the absolute traveltime residuals with the average traveltime residuals of all stations for each event in the area are removed. The traveltime residuals were corrected for crustal heterogeneity and topography. The crust model CRUST 1.0 model (http://igppweb.ucsd.edu/~gabi/crust1.html) was used in this process.

The teleseismic tomography used in this study applied the shooting ray tracing method (Julian and Gubbins, 1977; Rawlinson et al., 2006; Zhang et al., 2014, 2018) to calculate the travel time and ray paths; this process is a traditional ray tracing method and is commonly used in seismic tomography. The 3-D ray tracing was conducted following Snell's law, and by constantly adjusting and modifying the initial values of the ray parameters, relatively accurate ray paths and travel times can ultimately be obtained. The inversion problem of teleseismic tomography was solved by iteratively applying the damped least-squares algorithm (LSQR; Paige and Saunders, 1982), which can solve huge and sparse matrix equations with relatively less computation.

4 Reliability Analysis of Inversion Results

The 1-D IASP91 (Kennett and Engdahl, 1991) velocity model was used as the initial velocity model (Fig.4). For the initial grid division, considering that the coverage density of the seismic stations in the area is relatively low, after numerous tests, we decided that the optimum grid interval in the longitudinal and latitudinal directions was 1° for our dataset. In the vertical direction, the grid interval was 100 km, and the maximum depth was 900 km.

Given the multi-solution problem of geophysical inversion, the reliability of inversion needs to be verified and evaluated. To test the resolution of the observation system composed of the stations and seismic events used in this study, we conducted a checkerboard resolution test (CRT). In the CRT, the seismic stations, teleseismic events, and ray path distribution were completely consistent with the actual inversion. A 3-D perturbed velocity model was generated by adding the positive and negative velocity perturbations to an initial velocity model, and it used the stations and event distribution described above to calculate the theoretical relative travel time residuals of each ray for the inversion. By comparing the perturbed model and its inversion results, the scale of the smallest velocity anomaly that can be distinguished in the inversion results is the resolution of the model space.

The CRT velocity model used in this study was based on the initial IASP91 model with added ±3% perturbations and 1% random noise. Via repeated experiments on the scale of the velocity anomaly, when the scale of the abnormal block was 2°×2°×180 km, the results of the CRT became ideal. Fig.5 shows the horizontal slices of the CRT results at different depths from 100 km to 700 km. Compared with the output and input models, the resolution in the shallow (< 100 km) study region was relatively low, and only in the middle part of the study region, near the New Britain Island and New Guinea Peninsula, was the velocity perturbation recovered due to the seismic station distribution, which severely limited the shallow ray path. With the increase in depth, most of the study region exhibited a better resolution, especially at depths greater than 300 km.

The seismic stations are primarily distributed on the islands in the middle western part of the study region (141°- 154°E), whereas the eastern part of the study region is mostly sea and therefore lacks seismic stations, which resulted in an improved ray coverage in the middle and western parts of the study region. Therefore, the resolution of the middle and western parts (> 154°E) of the study region was better than that of the eastern part at all depths.

The inversion results of the LSQR algorithm were affected by the values of the damping parameter, which controls the convergence speed of the solution and smoothness of the model. By comparing the inversion results of different damping factor values, we obtained the relationship between the root mean square of the traveltime residuals and the smoothness of the model (Fig.6). When the damping factor was 20, the model smoothness and root mean square traveltime residuals were well balanced. Therefore, 20 was selected as the optimal value of the damping factor. Fig.7 shows the statistical results of the relative traveltime residuals before and after the inversion with a damping factor value of 20. Most of the residuals prior to the inversion were between -1.5 and 1.5 s. After the inversion, the residual distribution shrunk to the middle of the distribution, with the main body concentrated between -1 and 1 s. The distribution pattern conformed to the characteristics of a normal distribution, and the root mean square of the traveltime residual decreased by approximately 20%, indicating that the velocity model after the inversion can match the observed relative traveltime residuals.

5 Results

With the 1-D IASP91 velocity model as the initial reference model, the velocity structure beneath the PN-SL was obtained using teleseismic traveltime tomography. Fig.8 shows the imaging results at various depth levels.

In the depth range of 100 - 700 km, the northwestern and northeastern sides of the Solomon Sea were characterized by high-velocity anomaly zones. The spatial distribution of the high-velocity anomaly zones presented a clear correspondence with the New Britain Trench. With the increase in depth, the high-velocity anomaly zone at the New Britain Trench gradually migrated northward. The high-velocity anomaly zone along the New Britain Trench became wider at a depth of 400 km, and the southern boundary of the anomaly zone gradually expanded toward the Solomon Sea. At a depth of 700 km, the linear-shaped highvelocity anomaly zone at the New Britain Trench was no longer complete, but the PN-SL presented a mixture of highand low-velocity anomalies dominated by the highvelocity anomalies.

The Bismarck Sea is an evident low-velocity anomaly area, and this low-velocity anomaly can be traced to a depth of 600 km. At the West Melanesian Trench north of the Bismarck Sea, remarkable high-velocity anomalies were distributed in the depth range of 100 - 200 km; these were gradually dominated by low-velocity anomalies as the depth increased. At the Trobriand Trench, a high-velocity anomaly zone can be observed in the velocity slices at 100 - 200 km. Although the high-velocity anomaly was still observed near the Trobriand Trench as the depth increased, its spatial distribution pattern was not limited to the region along the Trobriand Trench. Evident low-velocity anomalies were recorded within the area of the Woodlark Basin. In particular, the western part of the basin showed significant characteristics of low-velocity anomalies in the range of 400 - 600 km. In the area of New Guinea Island, a remarkable high-velocity anomaly can be observed at a depth of 100 km and can be traced down to a depth of 700 km. Although this high-velocity anomaly within the depth range of 400 - 500 km is divided by a low-velocity anomaly, the New Guinea Island still exhibits velocity characteristics dominated by high-velocity anomalies.

6 Discussion

Based on the above analysis of seismic inversion results, we extracted four vertical sections along the red lines in Fig.1 to determine the deep shape of the subducted slab beneath the New Britain Trench and the eastern part of the New Guinea Island using the abundant distribution of the epicenters and good resolution (Fig.5). On the basis of the analysis of the deep-velocity anomalies beneath the PN-SL, we discussed the subduction pattern of the Solomon Sea Plate and its spatial geometric relationship with the subducted Pacific and Australian Plates.

As shown in Fig.9, as the most active subduction system of the PN-SL, the velocity structure of the subducted Solomon Sea Plate presented significant spatial differences. At the New Britain Trench, the subduction depth of the Solomon Sea Plate increased gradually from west to east. The subduction depth was 300 km on the P1 profile, which is similar to the findings of Pegler and Woodhouse (1995), and it increased significantly to 600 km on the P2 and P3 profiles, which is also supported by the spatial distribution of the epicenters (Holm and Richards, 2013). The detached slab also extended to a depth of 800 km on the P4 profile. However, the shape of the detached slab was unclear under the constraints of the number of seismic events. This finding indicated that the width of the Solomon Sea Basin prior to subduction could have been at least 1000 km. If the Solomon Sea had started being consumed by the New Britain Trench 8 Myr ago (Holm et al., 2016), the subduction rate of the slab would have been 10 cm yr^-1. As indicated by the high-velocity anomalies, the subduction angle of the Solomon Sea Plate shared a similar spatial variation trend with the subduction depth of the slab and gradually increased from west to east. The subduction angle of the slab increased from 45° on the P1 profile to 75° on the P2 and P3 profiles. On the P4 profile, the subducted slab was nearly vertical, which is similar to the subduction pattern of the old Pacific Plate at the Mariana Trench (e.g., Jaxybulatov et al., 2013; Barklage et al., 2015; Wu et al., 2016; Gong et al., 2019b). Confusingly, the subduction pattern at the New Britain Trench was completely different from that at the Peru-Chile Trench (e.g., Anderson et al., 2007; Hu et al., 2016), which is characterized by a similar seafloor age and subduction rate as the New Britain Trench (Müller et al., 2008; Schellart et al., 2011). In previous studies, the geometry of the subducting slab was interpreted to be controlled by its own gravitational pull (Pegler and Woodhouse, 1995), diachronous subduction initiation, and the eastern Solomon Sea remaining partly coupled to the Pacific (Hall and Spakman, 2002). For the New Britain Trench, with its young subducted plate and rapid subduction rate, we proposed that the high angle of the slab may be attributed to the migration of the trench and mantle convection (e.g., Faccenna et al., 2007; Boutelier and Cruden, 2008).

In contrast to the subduction at the New Britain Trench, the subduction depth of the Solomon Sea Plate at the Trobriand Trench changed from 200 km along the P1 profile to 100 km along the P4 profile. The eastward weakening of the subduction activity was supported by spatial changes in the seismic intensity, with a decrease in the seismic events from the P1 to P4 profiles. In addition, the subduction angle of the Solomon Sea Plate at the Trobriand Trench was 45°. Overall, the subduction activity of the Solomon Sea Plate at the Trobriand Trench was weaker than that at the New Britain Trench. According to Lock et al. (1987), the multichannel seismic reflection data showed a gentle onlap of the uppermost trench sediment onto the frontal anticline at the toe of the accretionary wedge. Moreover, the small earthquakes, shallow subduction depth, and lack of arc magmatism suggest that the subduction zone at the Trobriand Trench is presently inactive and does not accommodate a major component of the regional convergence (Davies et al., 1984; Wallace et al., 2004; Holm et al., 2016). Along with the collision between the Ontong Java Plateau, the subduction of the Solomon Sea Basin along the New Britain and Trobriand Trenches is the main adjustment style to accommodating the convergence of the plates. Different from the active and strong subduction along the New Britain Trench, the subduction activity at the Trobriand Trench weakened at the early Pliocene and may be stagnant at present. The different times of subduction activity led to the heterogeneous slab patterns in the two trenches. Further, as shown on the P1 profile, the Solomon Sea Plate formed an 'inverted U-shaped' double-dipping structural pattern at the Trobriand Trench and New Britain Trench (Pegler et al., 1995). As indicated by the focal mechanism and stress field analysis (Woodhead et al., 2010; Holm and Richards, 2013), the double subduction of the Solomon Sea Plate may be completely decoupled from the overlying Australian Plate. Thus, the Solomon Sea Plate at this location has been completely consumed beneath the convergence zone.

As shown in Fig.9, the deep structure of the early subduction system and its geometric relationship with the present active subducted Solomon Sea Plate is slightly more complicated. The Pacific Plate subducted strongly in the Cretaceous period along the West Melanesian and North Solomon Trenches, and the subduction of the Pacific Plate became extremely weak following the collision between the Ontong Java Plateau and Solomon Islands during the Oligocene period (Petterson et al., 1997). According to the inverted velocity model, the early subduction structure of the Pacific Plate showed spatial differences under the influence of later subduction activities. On the P1 profile, the Pacific slab extended to a depth of 200 km, and no slab was present in the deep mantle. According to the P2, P3, and P4 profiles, the late subduction of the Solomon Sea Plate had been intercalated within the early subduction of the Pacific Plate. Given the influence of preexisting faults within the slab and the late subduction activity, the Pacific slab in the deep mantle is discontinuous. However, the Pacific slab subducted downward at a large angle of > 60° and is stagnant in the mantle transition zone. The P2 profile showed that the detached Pacific Plate extended to a depth of 800 km. This finding is consistent with that of the inverted velocity model from the work of Hall and Spakman (2002), and our inversion results showed a considerably clearer shape of the subducted slab. As indicated by the P2, P3, and P4 profiles, at least 1000 km of the Pacific Plate has been consumed in the deep mantle at the West Melanesian Trench.

The Pocklington Trough in the southern part of the PNSL was active during the mid-late Eocene (45 Myr ago) and was inactive from 12 Myr ago (Schellart et al., 2006; Holm et al., 2015). On the P1 profile, we observed that the Australian Plate has a similar subduction scale to the eastern Pacific Plate. The Australian Plate extended to a depth of 600 km at an average dip of 60°, and 800 km of the slab has been consumed in the mantle. Although no similar interpenetrating relationship exists between the Pacific Plate and Solomon Sea Plate, the subducted Australian Plate collided with the subducted Solomon Sea Plate at the Trobriand Trench at a depth of 200 km. Further, the Australian slab is characterized by a break-off within the depth range of 300 - 400 km, which is evidenced by outcropped igneous rocks (Holm et al., 2015). Based on the tomographic inversion result and petrological records, we proposed that the Australian slab break-off may be correlated to a subduction polarity reversal caused by the southward subduction of the Solomon Sea Plate at the Trobriand Trench. This reversal would have removed the mantle wedge beneath the arc and culminated in ophiolite obduction in the Papuan Peninsula (Fitz and Mann, 2013).

In addition, the velocity profiles showed evident lowvelocity anomalies. According to our comprehensive analysis, three primary low-velocity anomalies, namely L1, L2, and L3, were observed (Fig.9). The L1 low-velocity anomaly was located between the subducted Solomon Sea Plate and the Pacific Plate and primarily distributed within the mantle above the Solomon Sea Plate. This anomaly was likely caused by the dehydration of the subducted plate, which led to the melting of the mantle, and should represent the deep dynamic factor of the Pliocene expansion of the Bismarck Sea back-arc basin with respect to the New Britain Trench (Baldwin et al., 2012; Holm et al., 2016).

The L2 low-velocity anomaly was located beneath the Solomon Sea and above the subducted Pacific Plate and primarily distributed at depths of 200 - 400 km. In terms of the current spatial distribution, the anomaly may be primarily caused by the melting of mantle induced by the dehydration of the subducted Pacific Plate, which in turn led to the expansion of the Solomon Sea. Thus, the expansion of the Solomon Sea is a 'back-arc expansion basin of a single subduction system' induced by the Pacific Plate subduction (Honza et al., 1987; Hall and Sparkman, 2003) rather than a back-arc basin of a double subduction system (Schellart et al., 2006). In addition, the expansion of the Solomon Sea may not be closely related to the subduction activity of the Australian Plate along the Pocklington Trough. Notably, a key issue was neglected in the above analysis: the early subducted Australian Plate might have been modified by later tectonic activities. In particular, the spread of the Woodlark Basin since the Pliocene induced the migration of the western Pocklington Trough and the deep subduction of the Australian Plate, which can be observed in the velocity structure shown on the P1 profile. As indicated by the P1 profile, which is on the western edge of the Solomon Sea and does not cross the Solomon Sea, the subducted Australian Plate lies beneath the double-subducted Solomon Sea Plate. In addition, a significant low-velocity anomaly body exists between the subducted Australian Plate and the double-subducted Solomon Sea Plate. From the tectonic location, this structure may be the same as that in the L2 low-velocity anomaly in the east, which dominated the back-arc expansion of the Solomon Sea. The subduction of the Australian Plate also played a role in the expansion of the Solomon Sea Basin. As for the coupling mode of the Solomon Sea spreading with the subduction of the Pacific and Australian Plate, further research is needed, in particular on the subduction dynamic simulations in the PN-SL.

The L3 low-velocity anomaly was located in the deep zone south of the Trobriand Trench and shown on the P4 profile. From the P4 profile, the L3 low-velocity anomaly extended from at least 100 km down to a depth of 600 km, with a scale larger than that of the L2 low-velocity anomaly. Using the seismic data recorded by the fixed stations, Abers et al. (2002) observed a significant low-velocity anomaly within the depth range of 0 - 240 km south of the Trobriand Trench; its position is consistent with the L3 low-velocity anomaly inverted in this paper. To a certain extent, the study is a supplement to the work of Abers et al. (2002), revealing the deep distribution of the L3 low-velocity anomaly. Combined with the study of normal focal mechanisms, young magnetic anomalies, and rapid exhumation (> 1 cm yr^-1) of Pliocene eclogite at the D'Entrecasteaux Islands (Weissel et al., 1982; Crook and Taylor, 1994; Taylor et al., 1999; Baldwin et al., 2004; Little et al., 2007), the northward extension of the Woodlark Basin is closely related to the generation of the largescale L3 low-velocity anomaly. The generation of the L3 low-velocity anomaly provided a dynamic factor for the expansion of the Woodlark Basin. As to whether the asthenospheric mantle corresponding to the L3 low-velocity anomaly body is an active or passive upwelling, or whether the Woodlark Basin was caused by back-arc expansion (Fitz and Mann, 2013), the activities of a tectonic block beneath the Goodenough Bay (Wallace et al., 2014), or the subduction and rotation of the Solomon Sea Plate (Benes et al., 1994; Baldwin et al., 2012; Webb et al., 2014), this study alone may be insufficient to provide a definite answer. In addition, the L3 low-velocity anomaly was found on the southern side of the subducted Pacific Plate with respect to the P4 profile. Thus, the L3 lowvelocity anomaly may not be related to the subduction of the Pacific Plate. Combining the horizontal slices of the tomography results and the P4 profile (Figs.8 and 9), the formation of the L3 low-velocity anomaly may be dominated by the subduction of the Australian Plate at the Pocklington Trough.

7 Conclusions

Through the teleseismic P-wave tomography of the PNSL, multiple and diverse types of subduction systems induced the complex deep structure beneath the PN-SL. As the most active subduction system of the PN-SL, the young Solomon Sea Plate is characterized by typical and distinct spatial differences in its subduction patterns, which cannot be revealed by previous studies. As the subduction angle of the Solomon Sea Plate increased from west to east along the New Britain Trench, the subduction depth also deepened, extending to a depth of 800 km. At the Trobriand Trench, the subduction depth of the Solomon Sea Plate is less than 200 km, and the subduction scale is significantly smaller than that of the New Britain Trench on the northern side of the Solomon Sea. Under the background of strong convergence of the Australian and Pacific Plates, the subduction rate of the Solomon Sea slab is close to or more than 10 cm yr^-1 since the late Miocene. During its large-scale subduction, the Solomon Sea Plate strongly changed the subduction structures of the early subducted Pacific and Australian Plates, and the subduction plates intersected each other. In addition, the subduction of the Solomon Sea Plate induced mantle melting and led to the formation of low-velocity anomaly bodies, which directly became a deep dynamic source for the expansion of the back-arc sea basin. As for the genetic mechanism of the expansion of the Solomon Sea back-arc basin, the inverted velocity structure indicates that the early subducted Pacific Plate played a key role in its expansion. However, constrained by the influence of the late tectonic transformation and general resolution of the velocity structure within the deep mantle, an important role of the Australian Plate in the expansion of the Solomon Sea cannot be completely excluded. Thus, the double subduction systems possibly dominated the formation of the Solomon Sea back-arc basin. Further, the subduction of the Australian Plate may have induced the formation of the lowvelocity anomaly below the Woodlark Basin. As for whether the expansion of the Woodlark Basin was caused by the upwelling of mantle melt or the movement of the microplate, further studies in combination with regional structural numerical simulations and other aspects are necessary.

Acknowledgements

We thank Wessel & Smith (https://www.soest.hawaii.edu/gmt/) for the free use of GMT software, which was used to produce most of the figures used in the study. We are grateful for the algorithmic support of the traveltime tomography study from Dr. Fengxue Zhang at the Institute of Geophysics, China Earthquake Administration. This work was supported by the National Natural Science Foundation of China (Nos. 91858215 and 41906048), the Fundamental Research Funds for the Central Universities (No. 201964015), and Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology (No. MMRZZ201801).