Low-velocity anomaly in the Coral Sea associated with subducting slabs and the Woodlark rift

Introduction

Plate tectonics provides insights into Earth’s surface processes, including the formation, migration, and subduction of tectonic plates, as well as the concentration of seismic activity and volcanism predominantly at plate boundaries. Complementing this framework, mantle plumes, originating from the core-mantle boundary, contribute to the formation of intraplate volcanism with seamount chains. These two processes are intertwined, generating distinctive tectonic features. For instance, interactions between plumes and ridges can result in the creation of large oceanic plateaus¹, asymmetric seafloor spreading ², or ridge jumps³. Similarly, when upwelling plumes intersect with subducting slabs, the slabs may stagnate, altering the direction of plume flow^4,5. Hence, elucidating enigmatic tectonic features may benefit from examining the interplay between mantle plumes and plate tectonics, rather than solely relying on either plate tectonics or plume models. However, classical thermal plume models often depict cylindrical plumes with plume heads measuring approximately 1000 km in diameter, followed by plume tails spanning about 100 km⁶. Such simplified models sometimes hinder explanations for puzzling surface tectonic features.

The southwestern Pacific (Fig. 1) presents a rare opportunity to capture the interaction of the three tectonic entities, as ridges and subduction zones are in close proximity owing to the complex tectonics involving microplates between the Australian and Pacific plates. Moreover, there are indications of potential plume upwelling⁷ in the Coral Sea, evidenced by positive heat flow anomalies⁸, positive dynamic topography⁹ and tomographic results^7,10,11,12. This area is characterized by a multitude of distinctive tectonic features. For instance, the Woodlark basin exhibits a seafloor spreading rate of 20–40 mm/yr, extending westward to the Woodlark continental rifting, one of the youngest (<8.3 Ma) and most rapidly extending continental rifts in the world¹³, which experiences 10–15 mm/yr of extension¹⁴. It hosts the youngest ultra-high pressure metamorphic rocks in the D’Entrecasteaux Islands (DEI), where domal uplifts are formed at a rate of approximately 20 mm/yr¹⁵. Adjacent to the DEI, a geochemical signature of enriched mantle like ocean island basalt (OIB) is observed along the Woodlark seafloor spreading centers¹⁶. Additionally, Pliocene-Quaternary intraplate volcanism is evident in Queensland, northeastern Australia⁷. However, the mechanisms driving these geological events remain unclear.

**Fig. 1: Tectonic setting of the study region with residual topography⁹.**

Here, we conduct multimode waveform tomography utilizing dense seismic networks in the southwestern Pacific, which incorporate ocean-bottom seismometers (Supplementary Fig. 1 and the Methods section). We present S-wave velocity and radial anisotropy models, reaching down to the bottom of the mantle transition zone (MTZ). Our analysis reveals a prominent vertical low-velocity anomaly beneath the Coral Sea, extending laterally at depths of 270–410 km toward the Woodlark rift and northeastern Australia. The widespread propagation of the low-velocity anomaly, dividing into several directions, shows complex behaviors of a thermochemical plume, which has been elusive in seismic tomography studies. The association of the Coral Sea plume with subducting slabs and rifts offers insights into the enigmatic geological features and the behaviors of thermochemical plumes.

Results and discussion

Low-velocity anomaly beneath the Coral Sea attributed to a thermochemical plume

In comparison to other back-arc basins in the western Pacific, the Coral Sea stands out with its unique combination of positive heat flow anomalies, measuring 23 ± 18 mWm^-2 relative to the predicted heat flow derived from the heat flow-age relationship⁸, and notably high dynamic topography, reaching up to ~1.7 km⁹ (Fig. 1). These phenomena are attributed to the presence of a mantle plume⁷. While some global tomographic models^10,11,12 have identified a cylindrical low-velocity anomaly in the lower mantle beneath the Coral Sea, rooted at the core-mantle boundary (Supplementary Fig. 2), its extension into the upper mantle remains elusive. Therefore, it is suggested that subducting slab remnants in the MTZ¹⁷ may impede the upward movement of the mantle plume, with a portion of it navigating through these obstacles⁷.

Current global^10,11,12 or local tomography models^7,18,19 offer too coarse resolution or have limited spatial and depth extents, making them inadequate to fully address these issues. Therefore, there is a pressing need for a new high-resolution upper-mantle velocity model, focusing on this area with denser data coverage. We conducted multimode waveform tomography²⁰ (Methods) by fitting S waves and multimode surface waves to obtain detailed upper mantle S-wave velocity and radial anisotropy models down to a depth of 660 km for the southwestern Pacific. We gathered three-component waveform data from the Incorporated Research Institutions for Seismology (IRIS) and the Ocean Hemisphere network Project (OHP) Data Managing Centers spanning the years 2004 to 2021. Our dataset comprised 171 land stations and 239 ocean-bottom seismometers. Further details on data acquisition, methodology, resolution test results, radial anisotropy results, and comparisons with global and regional tomographic models are provided in the Methods section and Supplementary Material.

At the bottom of the MTZ (Fig. 2f), high-velocity anomalies are widespread throughout the region, interpreted as slab remnants accumulated since the early Cenozoic¹⁷. However, a notable low-velocity anomaly is observed just south of the Woodlark basin, shifting southeastward with decreasing depth (Figs. 2d–f and 3f). We attribute the origin of the low-velocity anomaly to the Coral Sea plume arising from the lower mantle, as indicated by global tomographic results^10,11,12. At 650 km depth, the low-velocity anomaly overlaps with the location of a mantle plume in the lower mantle reported by French and Romanowicz¹⁰ and the region with thin MTZ reported by Yu et al.²¹. (Fig. 1) that indicates hot mantle upwelling from the lower mantle. Additionally, its position at 400 km depth well aligns with the peak of the positive residual topography (Fig. 1). Alternative explanations, such as mantle convection or upwelling associated with slab downwelling can be excluded, as these typically occurs due to dehydration on the front side of subducting slabs^22,23. In contrast, the low-velocity anomaly is observed behind the Solomon Sea slab.

**Fig. 2: Depth slices of Vs perturbations at various depths.**

**Fig. 3: Vertical cross-sections of Vs and ξs perturbations across the transects in Fig. 2a.**

From depths shallower than 400 km depth, the low-velocity anomaly extends radially, with northward and westward extensions reaching the Woodlark basin and the east coast of northeastern Australia, respectively. While the connection to the Woodlark basin is apparent (Fig. 3f), the linkage to the low-velocity anomaly beneath northeastern Australia is less pronounced (Fig. 3b). Compared with recent global^10,11,12 and regional^24,25 models (Supplementary Figs. 3–7), our model shows a clearer detection of the upwelling mantle plume in the upper mantle that aligns more closely with the peak of the positive residual topography. This improvement may be attributed to the use of dense networks, including ocean-bottom seismometers, and body and multimode surface waves with periods as short as 10 s.

Yu et al.²¹. found a thinner average MTZ thickness of 242 ± 7 km (minimum 229 km) to the east of 152°E (orange shaded region in Fig. 1) using receiver function analysis. Conversely, a thicker MTZ averaging 289 ± 13 km is observed to the west of 152°E beneath the Woodlark rift, indicating the presence of slab remnants within the MTZ. With Clapeyron slopes of +2.9–4.0 MPa K^-1 and −2.5 – −1.3 MPa K^-1 for the 410 and 660 km discontinuities^26,27,28,29, respectively, the average 8 km and the maximum 21 km thinner thickness than the normal average of 250 km thickness from the AK135 model³⁰ found east of 152°E corresponds to a mean temperature anomaly of 44–68 K with a maximum anomaly of 115–178 K. The center of the Coral Sea plume is located off the easternmost part of the MTZ thickness study area²¹ (Fig. 1), suggesting it may be hotter. Converting the lowest seismic velocity perturbation in section B-b (Fig. 3b) to temperature, following the approach of Bao et al.³¹, yields a maximum excess temperature of approximately 220 K (Supplementary Fig. 8a).

The origin of the Coral Sea plume appears to be thermochemical rather than purely thermal, with a more basaltic composition above 660 km depth, based on the following evidence. Segregation of harzburgite and basalt may occur at the bottom of the MTZ due to density contrasts, resulting in a higher proportion of harzburgite being positioned below 660 km and basalt above 660 km³². Thermochemical plumes with a dense basaltic composition in the upper mantle tend to pond at depths of 270–410 km³³(Supplementary Fig. 9) due to the density contrasts between basalt (denser) and pyrolite (lighter), bounded by the coesite-stishovite and olivine-wadsleyite phase transitions, as predicted in geodynamical modeling³⁴.

This phenomenon is observed in our S-velocity model, as the Coral Sea plume spreads laterally in this depth range (Fig. 3b, f). The influence of the Moresby-Aure slab on the lateral extent of the Coral Sea plume is limited, because the slab has already detached from the surface, sinking vertically without rollback, resulting in minimal surrounding mantle flow. Furthermore, the lateral extent of the plume is situated well away from the Solomon Sea slab (Fig. 3f), ruling out any substantial influence of the slab on the plume’s migration. The lack of a surface expression for the Coral Sea plume and its lateral flow at depths of 80–160 km just beneath the lithosphere, observed from lateral low-velocity anomalies (Fig. 3b, f) and positive radial anisotropy (Fig. 3d, h, Supplementary Fig. 10), also support its lower buoyancy.

Wamba et al.^35,36. identified horizontal ponding zones of mantle plumes at depths of 100–250 km and 660–1000 km beneath the Indian Ocean, attributed to low viscosity and the negative Clapeyron slope of the 660-km discontinuity. However, cross-sections in Fig. 1 of Wamba et al.³⁶. shows broad horizontal low-velocity anomalies extending from beneath the lithosphere to 410 km depth, which are better resolved into two distinct horizontal ponding zones at 80–160 km and 270–410 km depths in our model. This finer resolution to identify separate horizontal ponding zones caused by thermochemical characteristics of the mantle plume adds a new perspective to previous studies. The radius of ~200 km and the excessive temperature of 220 K of the Coral Sea plume, along with its nearly neutral buoyancy, may correspond to a basaltic composition of 10–15% (Supplementary Fig. 8b). Here, we consider composition, temperature, and plume size only for buoyancy estimation following the approach of Bao et al.³¹, while excluding other factors such as water content. This analysis further supports the interpretation of the Coral Sea plume as thermochemical in origin.

The lateral extent of the Coral Sea plume exhibits a notable shift in radial anisotropy strength with depth: strong positive radial anisotropy is evident at 80–160 km depth (Fig. 3d, h), while substantially weaker anisotropy appears within the lateral extent of the plume at 270–410 km depth. This variation may be attributed to changes in olivine types as the plume ascends³⁷. As the water-rich plume rises, dehydration likely occurs at >100 km depth due to high temperatures³⁸, prompting a transition from C-type olivine, which dominates within the water-rich plume at deeper depths, to A-type olivine in drier conditions at shallower depths. Because C-type olivine produces radial anisotropy oppositely to A-type olivine³⁷, C-type olivine could explain the strong positive radial anisotropy associated with vertical low-velocity anomalies indicating upwelling, and weaker anisotropy aligned with the lateral extent of the low-velocity anomaly at 270–410 km depth.

Association of the Coral Sea plume with slab remnants

Currently, active subduction of the Solomon Sea plate is observed at the New Britain and San Cristobal Trenches, with seismicity along the subducting slab. However, there is no definitive evidence for the presence of subducting slabs at the inactive troughs bounding the Woodlark basin, including the Trobriand Trough to the north and the Moresby-Aure-Pocklington Trough to the south (Fig. 1), leading to ongoing debate^14,39.

In the cross-sections C-c and D-d (Fig. 3e, f), we identified a steepened high-velocity anomaly from the New Britain Trench to the MTZ, accompanied by seismicity, which is interpreted as a subducting Solomon Sea slab undergoing rollback. Additionally, another steepened high-velocity anomaly is observed from the Moresby-Aure Trough south of the Woodlark rift (Fig. 3e), appearing to be detached from the surface. In contrast, no subducting slab is observed from the Pocklington Trough in cross-section D-d (Fig. 3f), which is connected to the Moresby-Aure Trough at the surface (Fig. 1). The high-velocity anomalies, interpreted as subducting slabs, are associated with negative radial anisotropy in our anisotropic model. This is consistent with, but more clearly observed than previous research¹⁹, suggesting vertical sinking of the slabs (Fig. 3g, h).

The slab morphology along the Moresby-Aure-Pocklington Troughs provides key insights into the unique characteristics of the Coral Sea plume. As the Coral Sea plume ascends in the upper mantle, it exhibits a southeastward shift toward the Pocklington Trough where no subducting slab is present. However, its southwestward movement is hindered by the presence of the Moresby-Aure slab, acting as a barrier. The steepening of the Moresby-Aure slab offers a compelling explanation for the orogenesis and the migration of arc volcanism in Papua New Guinea. Approximately 12 million years ago, the initial uplift of the New Guinea Orogen occurred due to continental collision with the Australian plate, coinciding with the southward migration of arc volcanism³⁹. This phenomenon may have been caused by the steepening and rollback of the Moresby-Aure slab (Fig. 3e). Additionally, the break-off of the Moresby-Aure slab around 7 million years ago likely contributed to a second phase of orogenesis, potentially explaining Pliocene-Quaternary volcanism in Papua New Guinea⁷. Following the slab break-off, the Moresby-Aure slab may have shifted ~400 km southward due to the northward Australian plate motion at a speed of ~60 mm/yr^13,40. In contrast, our model does not image the subducting slab at the Trobriand Trench, consistent with previous studies indicating no evidence of active subduction in that region^41,42.

The Coral Sea plume, transitioning southeastward to the Coral Sea with decreasing depth, exhibit a westward extension toward the east coast of Australia at ~400 km depth, attributed to the weak buoyancy of the thermochemical plume (Fig. 3b). This extended plume material may have generated an additional low-velocity anomaly, likely upwelling toward the surface and contributing to intraplate volcanism in northeastern Australia⁷. High attenuation is observed beneath the intraplate volcanism at depths shallower than 120 km⁴³ (Fig. 1), extending into the Coral Sea at greater depths, consistent with the observed plume morphology (Fig. 3b).

Influence of the Coral Sea plume on the Woodlark rift/basin

The Woodlark rift/basin has been extensively studied as a model of the transition from seafloor spreading to continental rifting, believed to have been predominantly shaped by slab-pull forces originating from the New Britain and San Cristobal Trenches at the northern boundary of the Solomon Sea plate^14,44. Additionally, it is proposed that dehydration from slab remnants subducted at the Moresby-Aure Trough weakened the upper mantle, facilitating localized extension^19,21, which is supported by the presence of a subduction-enriched source found in volcanoes along seafloor spreading centers⁴⁵. The Woodlark rift hosts the exhumation of youngest ultra-high pressure rocks in the DEI, with uplifts of 1500–2500 m above sea level¹⁵, despite anticipated subsidence associated with observed crustal thinning of 10–15 km⁴¹. Abers et al. ⁴². demonstrated that Bouguer gravity anomalies are indicative of isostatic compensation, and estimated temperatures from velocity inversion are higher than predicted from pure-shear models of lithospheric thinning, implying additional buoyancy sources supporting the topography. Therefore, it is suggested that the high topography of the DEI may be compensated by a hot, low-density, low-velocity upper mantle^42,45. However, the role of a hot mantle plume in the rifting has been discounted due to several observations: azimuthal anisotropy parallel to the extensional direction as observed in typical mid-ocean ridges⁴⁶, positive radial anisotropy indicative of passive rifting¹⁹, the absence of a deep-seated low-velocity anomaly below 250 km depth¹⁹, and MTZ thickening found west of 152°E below the rift suggesting the presence of slab remnants rather than mantle plumes that would cause MTZ thinning ²¹.

In our model, a strong low-velocity anomaly is observed at shallow depths along the Woodlark rift/basin, reaching a maximum of −8% (Fig. 3a). The anomaly extends down to ~200 km depth without deep-seated low-velocity anomalies, a characteristic typically associated with mid-ocean ridges. However, our model clearly indicates an upwelling low-velocity anomaly from the base of the MTZ, diverging into the Woodlark basin and the Coral Sea after spreading widely at depths between 270 and 410 km (Fig. 2b–d and Fig. 3f). At depths shallower than ~60 km depth, the density contrast between basalt and pyrolite is reversed (lighter basalt, denser pyrolite) due to de-eclogitization, potentially increasing the upwelling speed of the Coral Sea plume beneath the Woodlark basin, thereby inducing positive buoyancy. In our model, negative radial anisotropy is observed in the eastern part of the Woodlark basin near the surface, where the upwelling Coral Sea plume contacts the preexisting spreading centers (a dotted red circle in Fig. 3c), indicating vertical mantle flow. In contrast, positive radial anisotropy predominates in other parts of the Woodlark rift/basin, consistent with previous studies¹⁸. The plume material, fed by the vertical flow in the eastern part of the basin, could have migrated westward along the pre-existing Woodlark seafloor spreading centers, contributing to the OIB signature observed near the DEI¹⁶. This influx of heat could concentrate crustal thinning and buoyancy, partially contributing to the uplift of the DEI and potentially facilitating continental rifting. Hence, the formation of the Woodlark rift/basin likely involves a complex interplay, including slab-pull force from the New Britain Trench, dehydration from slab relics^19,21, and plume activity.

The southwestern Pacific exhibits a plethora of unique geological and tectonic characteristics, including positive heat flow anomalies⁸, positive dynamic topography ⁹, most rapidly extending continental rifting^13,14, exhumation of the youngest ultra-high pressure rocks¹⁵, OIB signatures¹⁶, and intraplate volcanism⁷ in northeastern Australia. These features present challenges to classical plate tectonics and thermal mantle plume frameworks. Our study demonstrates that the enigmatic nature of these features finds resolution in the association of a concealed thermochemical mantle plume, devoid of surface expression, with subducting slabs and ridges (Fig. 4).

**Fig. 4: 3-D isosurfaces from our S-velocity model.**

Methods

Data

In this study, we used three-component seismograms collected from the IRIS and the OHP Data Managing Centers during years 2004–2021. The IRIS dataset comprises 171 temporary and permanent land stations, along with 190 ocean-bottom seismometers (OBSs) in the western Pacific. Additionally, 49 OBS data were obtained from the OHP Data Managing Center. To ensure robust waveform data quality, we selected earthquakes with magnitudes ranging from 5.5 to 7.5 (Supplementary Fig. 11) and epicentral distances of (4^{circ} -65^{circ}). This yielded a total of 410 stations and 3373 events, resulting in a total of 53,222 ray paths (Supplementary Fig. 1). The hit-count map of ray paths used in this study demonstrates extensive ray coverage, particularly in the Solomon Sea and the eastern Coral Sea regions.

Waveform tomography

Witek et al.²⁰ proposed an efficient and automated algorithm based on the partitioned waveform inversion method developed by Nolet⁴⁷ to analyze large datasets and calculate radial anisotropy. The initial step in waveform inversion involves nonlinear waveform fitting to establish linear constraints. This fitting process yields perturbations in the path-averaged S wave velocity and radial anisotropy for each ray path.

In the first stage of waveform fitting, fundamental-mode waveforms are extracted using a phase-matched filter based on the 3-D reference model consisting of CRUST1.0⁴⁸ and AK135³⁰, along with source information from a published GCMT solution⁴⁹. After extracting the fundamental mode, we determine a minimum fitting frequency, from which we apply a series of overlapping Gaussian bandpass filters, (Gleft(fright)=exp left[-alpha {left{left(f-{f}_{0}right)/{f}_{0}right}}^{2}right]), where ({f}_{0}) is the center frequency of the filter, and α (set to 12.5) is the filter parameter. To identify the minimum frequency, we require at least three wavelengths between the source and receiver to satisfy the far-field approximation. We then evaluate the radiation pattern and discard frequencies where the amplitude at the station azimuth drops below 50% of the maximum. A signal-to-noise ratio of at least 10 is required, defined as the ratio of the maximum amplitude of the filtered signal to the root-mean-square of the pre-P-wave noise. If these criteria are not satisfied, the minimum frequency is slightly increased, and the process to find a suitable minimum frequency is repeated.

Once the minimum frequency is determined, we calculate center frequencies such that the Gaussian filters overlap at 90 % of their amplitudes. The maximum fundamental-mode fitting frequency is set to 62.5 mHz. The observed waveform is divided into several time-frequency windows, and the misfit ({F}_{i}left({{{boldsymbol{gamma }}}}right)) for i-th window is defined as

$${F}_{i}left({{{boldsymbol{gamma }}}}right)=frac{{int }_{{t}_{i}}{w}_{i}left(tright){left[dleft(tright)-uleft({{{boldsymbol{gamma }}}},tright)right]}^{2}{dt}}{{int }_{{t}_{i}}{w}_{i}left(tright)d{left(tright)}^{2}{dt}}$$

(1)

where (dleft(tright)) and (uleft({{{boldsymbol{gamma }}}},tright)) represent the observed and synthetic waveforms, respectively, ({w}_{i}left(tright)) is a weighting function that downweights large-amplitude wave groups within a window by using the inverse of the envelope, and γ is a model parameter vector. To account for amplitude differences from uncertainties in attenuation structure or the source mechanism, both the observed and synthetic waveforms are normalized by the energy of the signal within each time-frequency window,

$$dleft(tright)leftarrow frac{d(t)}{sqrt{{int }_{{t}_{i}}{d}^{2}left(tright){dt}}},,uleft({{{boldsymbol{gamma }}}},tright)leftarrow frac{u({{{boldsymbol{gamma }}}},t)}{sqrt{{int }_{{t}_{i}}{u}^{2}left(tright){dt}}}.$$

(2)

Nonlinear waveform fitting proceeds from the lowest to the highest frequency window. During the process, the 3-D reference model consisting of CRUST1.0⁴⁸ and AK135³⁰ was used to account for crustal effects. Optimization of the results from the nonlinear waveform fitting was performed to minimize misfits with the observed waveforms, employing the BFGS algorithm (after Broyden-Fletcher-Goldfarb-Shanno)⁵⁰, as shown in Supplementary Fig. 12. Synthetics were calculated using normal mode codes based on MINEOS⁵¹.

If the fundamental-mode waveform fit was successful, the resulting model perturbations serve as the initial model for the second stage of full waveform fitting. Here, the dominant misfit arises from overtones, as fundamental mode was already fitted. The fitting frequency of overtones extends to 100 mHz with a series of Gaussian bandpass filters. The mode sum includes modes n = 0, 1, 2, …, 20 to ensure accurate reconstruction of the S waves and multiple S waves. Fitting procedure in the second stage is the same as for the fundamental-mode waveform fitting, except that we fit the full waveform directly rather than an extracted fundamental-mode waveform using a phase-matched filter. An example of the full waveform fitting is shown in Supplementary Fig. 13. Histograms of misfits between observed and synthetic waveforms are presented for both fundamental-mode and full waveform fitting in Supplementary Fig. 14, respectively. Significant misfit reductions are achieved with the final models.

After acquiring the path-averaged linear constraints from the waveform fitting, we conducted an iterative inversion to derive 3-D perturbations of Vs and radial anisotropy relative to the 3-D reference model using the LSQR algorithm as detailed in Paige and Saunders⁵². Optimal damping and flattening parameters were determined through a trial and error process, balancing variance reduction and model roughness. To mitigate the influence of dominant ray paths, we applied inverse weighting based on the number of similar ray paths, defined as those where the events and stations are within 2% of the path length. Further details can be found in Witek et al.²⁰.

Model parameterization

We parameterize our model using isotropic S wave velocity and radial anisotropy as follows:

$${V}_{S}^{2}=frac{1}{2}({V}_{{SH}}^{2}+{V}_{{SV}}^{2})$$

(3)

$${zeta }_{S}=frac{{V}_{{SH}}^{2}-{V}_{{SV}}^{2}}{2{V}_{S}^{2}}$$

(4)

We link variations in density and isotropic P wave velocity to S wave velocity according to Anderson et al.⁵³ and Robertson & Woodhouse⁵⁴:

$$frac{{{{rm{dln}}}}rho }{{{{rm{dln}}}}{V}_{s}}=0.4$$

(5)

$$frac{{{{rm{dln}}}}{V}_{P}}{{{{rm{dln}}}}{V}_{s}}=0.5$$

(6)

For anisotropy, we convert radial anisotropy ({zeta }_{S}) to the conventional radial anisotropy parameter ({xi }_{S}={({V}_{{SH}}/{V}_{{SV}})}^{2}) for the comparison purpose. We calculate the perturbations of the S wave velocity and radial anisotropy and converted to percentages as follows:

$${{{rm{dln}}}}{V}_{S}=frac{{V}_{S}-{V}_{S}^{0}}{{V}_{S}^{0}}times 100( % )$$

(7)

$${{{rm{dln}}}}{xi }_{S}=frac{{xi }_{S}-{xi }_{S}^{0}}{{xi }_{S}^{0}}times 100( % )$$

(8)

where ({V}_{S}^{0}) and ({xi }_{S}^{0}) represent the reference S wave velocity and radial anisotropy calculated via the 3-D reference model consisting of CRUST1.0⁴⁸ and AK135³⁰, respectively.

We use the spherical tessellation method of Wang & Dahlen⁵⁵ to parameterize model perturbations on nested shell grids. The grid is centered at 0°N, 168°E and extends 55° in all directions. Each depth layer has 42,556 points, with an inter-grid spacing of ~50 km at the surface. The 28 depth layers are located at 0, 5, 10, 20, 35, 55, 75, 95, 120, 145, 170, 200, 230, 260, 290, 320, 350, 380, 410, 450, 490, 530, 570, 610, 660, 710, 780 and 900 km depths.

Resolution tests

To validate our tomographic model, we performed several checkerboard resolution tests based on the methodology of Rawlinson & Spakman⁵⁶, incorporating sparse 5% sinusoidal variations to avoid misleading results. We generated synthetic data vectors by multiplying the sparse checkerboard model with the sensitivity kernel matrix, using the same regularization parameters as for the inversion of observed data. Gaussian noise was added to the synthetic data at a level proportional to the ratio between the observed data norm and the estimated noise norm. We set input models with either isotropic or anisotropic structures of sizes of 2°, 3°, and 4° and with amplitudes of ±5%.

Supplementary Figs. 15 and 16 show the results of the sparse checkerboard resolution tests for isotropic ({V}_{S}) and (zeta), respectively. Both tests show good recoveries down to 600 km depth, except for some smearing observed in the western Coral Sea. Anomalies with sizes of 2° are well resolved down to 600 km depth, despite some smearing effects beneath the western Coral Sea. Anomalies of 3° and 4° sizes are also well resolved down to 600 km depth, with improved constraint on the locations and amplitudes of peak anomalies compared to the 2° sized anomalies. Therefore, our model demonstrates a fair resolution down to the MTZ beneath the study region.

We also conducted vertical checkerboard tests for both isotropic ({V}_{S}) and (zeta) checkers with sizes of 2° and amplitudes of ±5% along latitudes of 6°S, 10°S, 14°S, and 18°S, and along longitudes of 150°E, 155°E, and 160°E (Supplementary Figs. 17 and 18). For each transect, checkers were positioned at depths of 100 km, 300 km, and 500 km, with a lateral spacing of 200 km and a thickness of 50 km. The results showed that both isotropic and anisotropic checkers are well resolved, except for some smearing beneath the western Coral Sea. Leakage between isotropy and anisotropy is less than 1% for both horizontal and vertical checkerboard tests, indicating that errors from leakage are well constrained in our model.

We then performed structural resolution tests to validate the anomalies identified in our model. Supplementary Figs. 19 and 20 present the results of structural resolution tests for isotropic ({V}_{S}) and (zeta) anomalies, respectively. The isotropic ({V}_{S}) input models were designed to resemble realistic structures, featuring low-velocity anomalies extending to 200 km depth beneath the Woodlark rift/basin, low-velocity anomalies reaching the MTZ beneath the Coral Sea, and high-velocity anomalies down to the MTZ beneath the New Britain Trench and the Moresby-Aure Trough. For radial anisotropy, input models included negative radial anisotropy along the subducting slabs and beneath the eastern Woodlark basin, and positive radial anisotropy beneath the remaining Woodlark rift/basin and the Coral Sea. Some horizontal smearing effects were observed in Sections B-b, C-c, and D-d that cross the western Coral Sea. However, the results confirm that main anomalies in our model are well recovered, with only weak smearing effects. This validates the robustness of our model in capturing the primary features of the mantle structure beneath the study region.