Seismic evidence for a thermochemical mantle plume underplating the lithosphere of the Ontong Java Plateau

The Ontong Java Plateau in the western Pacific Ocean is the world’s largest oceanic plateau. It was formed 122 million years ago by a massive volcanic event that significantly affected Earth’s environment. The cause of the magmatic event remains controversial because the upper mantle structure beneath the plateau is poorly known. Here we use passive seismic data obtained through seafloor observations, alongside existing seismic data, to determine the three-dimensional radially anisotropic shear wave velocity to depths of up to 300 km. We find that the lithosphere–asthenosphere boundary is approximately 40 km deeper beneath the centre of the Ontong Java Plateau than beneath the surrounding seafloor. Based on our results and petrological and rheological constraints, we propose that the lithosphere–asthenosphere boundary has deepened as a result of underplating of dehydrated residual material beneath the pre-existing lithosphere during formation of the Ontong Java Plateau by a thermochemical mantle plume. The lithosphere beneath the Ontong Java Plateau is thickened by up to 40 km relative to the surrounding ocean floor which suggests it formed through the ascent of a mantle plume, according to analyses of seafloor seismic observations.

T he Ontong Java Plateau (OJP) in the western Pacific Ocean is the largest oceanic Large Igneous Province (LIP) on Earth 1 , with crust 30-40 km thick [2][3][4] . While the entire OJP is thought to have been emplaced at approximately 122 Ma by a single massive magmatic event, influencing the global climate, as well as causing oceanic anoxic events and the extinction of marine species 1 , the cause of the magmatism is still debated. The geochemistry of OJP basalts is extraordinarily homogeneous 5 , which differs from those of other LIPs and hotspots. The main OJP emplacement was submarine at depths greater than 800 m below sea level 6 , but partly subaerial 7 , and post-emplacement subsidence has been minimal compared to normal seafloor 8 .
Three major hypotheses have been proposed for the origin of the OJP to explain these characteristics: (1) the ascent of a hot mantle plume 1,9,10 with a potential temperature (T p ) higher than 1500°C; (2) passive mantle upwelling near a spreading ridge 11 ; and (3) a meteorite impact resulting in extensive melting 12 ; however, none of these hypotheses satisfy existing data or are conclusive. Identifying the origin of the OJP is not only essential to understand Cretaceous Earth geodynamics, but also to evaluate the magnitude and mechanism of a change in the climate that may be induced by a future massive magmatic event. The seismic structure beneath the OJP may be the key to resolving its origin. For example, a recent study 13 found strong radial anisotropy just beneath the Moho at the Hikurangi Plateau, which may have been a part of the single Ontong Java-Manihiki-Hikurangi Plateau 14 , with interpretations of its formation via a plume. Current knowledge on the structure beneath the OJP is sparse, particularly for the upper mantle, and seismic studies to date have yielded no direct evidence of a plume. Until now only two studies have focused on the upper mantle beneath the OJP, and the resulting models differ markedly. One study suggests the presence of a lowvelocity, rheologically strong mantle root extending as deep as 300 km beneath the entire OJP 15 . Shear wave splitting analysis around the OJP also supports this mantle root hypothesis 16 . More recently, a very high-velocity and compositionally different mantle has been proposed, extending to a depth of 100 km beneath the corner of the L-shaped OJP 17 , with no noticeable low-velocity anomalies. These opposing results may arise from the data used, which are solely from land stations well beyond the boundaries of the submarine OJP. To investigate the OJP's upper mantle more reliably, in situ seismic seafloor observations are required.
From late 2014 to early 2017, we conducted the first broadband seismic and electromagnetic investigation of the OJP using the 'OJP array' 18 on OJP and adjacent seafloor and neighbouring islands (Fig. 1). Here, we determine the three-dimensional radially anisotropic shear wave structure of the upper mantle in and around the OJP via surface wave tomography using seismograms recorded by the OJP array and data from existing land and seafloor stations ( Supplementary Fig. 1).

Results
Seismic structure. The shear velocity model obtained in this study shows that a high-velocity lid (4.45-4.6 km s −1 ) at a depth of~100 km in most regions overlies low-velocity regions at~150 km depth. This feature is seismic evidence of plate tectonics, i.e., the cold, rigid oceanic lithosphere overlies the hot, weak, viscous, and ductile asthenosphere. However, under the central OJP, the high shear velocity region (4.45-4.6 km s −1 , 2% high-velocity anomalies at a depth of 150 km) extends to 150 km, a remarkable feature of the upper mantle. Radial anisotropy is weak (ξ ≡ (β H / β V ) 2 ≈ 1.0) and strong (ξ > 1.1) in the lithosphere and asthenosphere, respectively (Figs. 2 and 3), which is similar to normal Pacific Ocean upper mantle 19,20 , but the depth of the lithosphere-asthenosphere boundary (LAB) is significantly greater beneath the OJP. The LAB depth, defined here as the depth of the negative peak in the vertical velocity gradient (see Methods, Fig. 3a), is~130 km, which is much deeper than the average LAB depth (80-90 km) beneath Pacific Ocean seafloor of a similar age in the study region as estimated from surface wave analysis 20 . The LAB depth beneath the OJP is deeper than the deepest LAB previously reported in the Pacific Ocean, i.e., 120 km in the northwest Pacific based on SS precursor analysis 21 .
The Nauru Basin, east of the OJP, has a LAB depth of 90 km, similar to the average depth of the LAB in the Pacific Ocean. North of the OJP is an east-west trending low-velocity region (4.35-4.45 km s −1 , 0.5-1% low-velocity anomalies, Fig. 3b) with depths of at least 250 km along the Caroline Islands. Radial anisotropy is strong (ξ > 1.1) at depths shallower than 150 km. The LAB depth is 70 km, which is consistent with the findings of a previous study that analysed SS precursors 22 . These features may be related to the Caroline plume, which is thought to be present below the easternmost Caroline Islands. The depth of the LAB gradually thins to the east, while thinning to the north and southeast of the OJP is abrupt (Fig. 3a, c, and e). As the crustal thickness of the central OJP is 30-40 km and that of the Nauru Basin is 10 km 23,24 , the lithospheric mantle is ≥20 km thicker beneath the OJP than beneath the Nauru Basin. However, owing to limited depth resolution and the definition of the LAB depth adopted in this study, the LAB depth beneath the OJP may be underestimated by 20 km or more (see Methods, Supplementary Note 1, and Supplementary Fig. 2). The observation of thick lithosphere is possible with multi-mode surface wave tomography, as the first higher-mode Rayleigh wave travelling through the OJP exhibits high phase velocities, which is sensitive to shear wave velocities at a depth of~100 km and insensitive to those at crustal depths (see Supplementary Note 2 and Supplementary  Figs. 3 and 4). Thick lithosphere beneath the OJP is qualitatively consistent with weak seismic attenuation in the OJP mantle 25,26 .

Discussion
Thermal versus compositional factors. What is the originthermal and/or compositional-of the thick lithosphere beneath the OJP? If thermal effects are the cause, the thick lithosphere, with high-velocity anomalies of~1-2%, should be accompanied by lower temperature anomalies of~100-200 K 20 compared to lithosphere beneath the Nauru Basin. One possibility is a difference in lithospheric ages between the OJP and Nauru Basin; according to a half-space cooling model 27 , the age of the OJP lithosphere should be~40 Ma older than that of the Nauru Basin (~130-160 Ma 23 ) to account for the temperature difference. However, interpreted seafloor ages underlying the southernmost OJP are~160 Ma 28 and of the OJP's uppermost crust~122 Ma 12 , which are unlikely to have caused such a large temperature difference. Other cooling models 29 provide even smaller temperature differences. Another possibility is the existence of small-scale downwelling flow 30 in the asthenosphere beneath the OJP. However, it appears unrealistic that a localised downwelling flow has moved with the OJP and maintained a cool lithosphere. Strong radial anisotropy beneath the OJP lithosphere at depths >150 km indicates the existence of horizontal flow rather than vertical flow. Therefore, we propose that the cause of the thick lithosphere is compositional, which agrees with studies of xenoliths sampled from the Malaita Island where OJP crust has been obducted 9,28,31-34 (Fig. 1). This xenolith population comprises diverse suites of peridotites and garnet pyroxenites that represent virtually the entire section beneath the southern margin of OJP lithosphere, extending to depths of ≤120 km. The upper part (<80 km) of OJP lithosphere is mainly composed of fertile lherzolites originating from pre-existing oceanic lithosphere with subordinate pyroxenitic dykes or sills, whereas the deeper part of the OJP lithosphere, at depths between 90 and 120 km, can be interpreted as melt residues of a mixed peridotite-pyroxenite source from a thermochemical mantle plume 9 . Based on petrological calculations, shear wave velocities of the peridotite beneath the OJP at depths of 60-120 km are estimated to range between 4.6 and 4.5 km s −1 , respectively 31 , which is consistent with our seismological results of this study ( Supplementary Fig. 5). The petrological model does not involve anomalous velocity structures, as proposed by previous studies, such as a broad lowvelocity anomaly (4.0-4.1 km s −1 ) at depths from 80 to 150 km 15 or very high anomalies (4.75 km s −1 ) extending to a depth of  17 . We consider the seismic structure described here to be more reliable than any structure previously published because our imaging is the first to be based on in situ data from the OJP. Therefore, we conclude that the seismically determined velocity model presented here is consistent with the petrologically determined velocity model, such that the thickened lithosphere was most likely caused by underplating of melt residues of a thermochemical mantle plume beneath pre-existing oceanic lithosphere at the time of OJP emplacement (Fig. 4).
The arrival of a hot and buoyant mantle plume comprising anhydrous peridotite at the base of the lithosphere should result in considerable subaerial volcanism 11 , which is inconsistent with the observation that the main OJP emplacement was submarine 6 . However, a thermochemical plume 35 , assuming the incorporation  of compositionally dense but fusible material such as eclogites or garnet pyroxenites, would decrease both the solidus temperature and net buoyancy flux of the plume, such that a large outpouring of basaltic magma could occur with minor topographic uplift of the OJP.
We speculate that the high potential mantle temperature (T p > 1500°C), derived from the geochemistry of OJP basalts 36,37 may be overestimated, as modelled 20-30% melting of homogeneous anhydrous peridotite is not a unique solution for the observed lava compositions. Evidence to date shows that OJP basalts are extraordinarily homogeneous 5 , depleted in incompatible trace elements including H 2 O and CO 2 5,35 , and slightly enriched in platinum-group elements (PGEs) 38 compared with other LIPs and hotspots thought to be generated by melting of heterogeneous plume sources. However, such observations in favour of the OJP magma source being homogeneous anhydrous peridotite remain controversial and inconclusive. Indeed, previous studies 11, 39 have pointed out that the high-Ni compositions of olivine from primitive OJP lavas are not compatible with the melting products of common peridotites. This led to the proposal of olivine-free lithology as a more appropriate source for the OJP lavas such as pure eclogite 11 or hybrid pyroxenite formed from reaction between eclogite-derived Si-rich melts and peridotite 39 .
We favour the latter because a series of hybrid pyroxenites is the main constituent of deeply derived xenoliths found in Malaita 9,32,33 . Moreover, it would be possible to reconcile the observed chemical and isotopic variations of OJP lavas if the hybridisation process during plume ascent created almost uniform pyroxenite as the main magma source 39 ; the melting of hybrid pyroxenite could enhance the release of compatible elements such as Ni and PGEs into the melt, resulting in elevated concentrations of Ni and PGEs in OJP lavas; and the olivine-free but pyroxene-rich residues left behind after the melting of hybrid pyroxenite could accommodate substantial amounts of incompatible trace elements that may account for incompatible element depletion in the primitive OJP magma.
Currently, geochemical modelling of trace element abundances and isotopic compositions of the OJP lavas is challenging because the compositional variability expected for recycled eclogites and resulting hybrid pyroxenites cannot be tightly constrained. This is partly because almost all OJP geochemistry is based on limited data obtained from the uppermost 3-4 km section of erupted lavas, presumably after effective mixing in large magma chambers in an already thickened crust. Nevertheless, a compositionally and isotopically enriched type of basalt was recognised as an overlying veneer of the main group of basalts at several localities, which seems consistent with the more fusible components in the heterogenous plume that was tapped during the waning phases of the OJP magmatism 36,40 . Future sampling of significantly deeper crustal sections 41 could clarify the extents of the heterogeneities in OJP lavas, which is required for viable trace element modelling and T p estimates for the thermochemical plume 35 .
Other hypotheses, such as the passive upwelling of dense material 11 and a meteorite impact 12 , which do not invoke deep mantle plume activity for the formation of the OJP, are consistent with the lack of convincing evidence of a post-OJP hotspot track (plume tail). The Louisville seamount trail is a potential plume trail 42 but remains controversial. However, the thickened lithosphere observed in this study is difficult to reconcile with these non-plume scenarios. In the passive upwelling hypothesis, basaltic magma would have been produced by almost complete melting of recycled crustal material that ascended via asthenospheric upwelling near a spreading ridge. This appears unlikely to have produced melt residues sufficient to thicken the OJP lithosphere as this essentially prevents the melting of anhydrous peridotite under high-pressure conditions. For impact-induced melting, basaltic magma would have been largely supplied first from the uppermost~150 km of pre-existing crust and mantle beneath the impact site where the most decompression and shock/release wave melting would have occurred 43 , although an extended secondary episode of additional melting would occur at greater depths 44 . However, if the melting region were sufficiently large to thicken the vast area of OJP lithosphere, Jurassic-aged oceanic crust and mantle identified in the Malaita xenolith population should have been disturbed 28 , which has not been observed. A lack of extraterrestrial signatures in coeval Pacific marine sediments would also suggest no impact event at the time of OJP formation 34 . We consider impact hypothesis to be unlikely, although absence of evidence does not equate to evidence of absence.
How can melt residues survive for over 120 Ma? During the formation of OJP crust, the volatile components of the plume, including H 2 O were incorporated into the basaltic melt, resulting in dry 45 and rheologically strong 46 residual plume material. The restite root hypothesis 47 , which assumes that a dry, hot, highly viscous restite root attached itself beneath the hotspot lithosphere, yielding a thick high-viscosity layer (lithosphere), was proposed to explain hotspot swell uplift. Recent receiver function analysis suggests that the LAB is as deep as 100 km beneath hotspot islands in the Pacific Ocean 48 , which is similar to the deep LAB beneath the OJP. Recent rheological and geophysical studies 49,50 indicate that the asthenosphere is wet and thus rheologically weak compared to the lithosphere, which suggests that the dry and rigid melt residues could have moved with Pacific lithosphere over long time intervals. A similar mechanism has been proposed to explain the extended life of continental cratons (>2 Gyr) 51 .
The area of >20 km thicker lithosphere beneath the OJP is 5 × 10 5 km 2 , so that the volume of the dry, dense residue that we have identified is~1-2 × 10 7 km 3 , which is significantly less than the mantle source volume estimated to form OJP crust (1-3 × 10 8 km 3 ) 35 . One possible explanation is that most residue moved laterally and did not underplate during emplacement and the last residue formed thick lithosphere during the waning phase of activity. An alternative explanation is that most dense residue has delaminated while the OJP has moved with Pacific lithosphere after emplacement. The latter explanation may be more plausible given that this might explain the less-than-predicted subsidence of OJP seafloor.

Methods
We employed a surface wave tomography method consisting of three independent stages 20,52,53 . First, we measured the path-specific multimode phase-velocity dispersion curves of Love and Rayleigh waves from the fundamental to the 4 th higher mode using a fully nonlinear waveform inversion method. The waveform-fitting process was fully automated, where the quality of the obtained dispersion curves was evaluated by a reliability parameter and determined by the degree of waveform fitting and relative power of each mode in the selected bandpass filtered time window. This empirical criterion is identical to that of our previous study 20 . We obtained, at most, 2000-4000 phase velocity dispersion curves for the fundamental Rayleigh waves, and 500-1500 for the higher modes in the 5-30 mHz frequency range. For the Love waves, the number of dispersion curves was approximately half the number obtained for the Rayleigh waves. Approximately 5-10% and 20-30% of the Love and Rayleigh wave dispersion curves, respectively, were obtained from broadband ocean bottom seismographs (BBOBSs). For the second stage, we inverted the path-specific multimode phase velocities for twodimensional phase-velocity mapping at each frequency, incorporating finite-frequency and ray-bending effects. Phase velocity maps were expanded in a set of B-spline functions with a rectangular grid, with the interval in each map varying from 2°to 6.5°d epending on the number of paths and the size of the Fresnel zone for each mode and frequency. Finally, we inverted the multimode phase velocity maps to generate a threedimensional radially anisotropic shear wave velocity model, employing iterative leastsquares inversion. The use of up to the 4 th higher modes of surface waves enables us to improve the depth resolution of the radial anisotropy ( Supplementary Fig. 6).
Following the protocol introduced in our previous study 20 , we used β V and β H as independent parameters, showing β iso , which is the Voigt-averaged isotropic shear wave velocity, where β iso 2 = (2/3)β V 2 + (1/3)β H 2 and ξ = (β H /β V ) 2 . In this inversion, two a priori parameters control the degree of perturbation and smoothness in the depth variations-the standard deviation, σ, and correlation length, L-which were respectively set tο 0.01 km s -1 and 3 km in the crust, 0.04 km s -1 and 15 km between the Moho and a depth of 400 km, and 0.03 km s -1 and 15 km between depths of 400 and 670 km. Here, σ decreases linearly from 0.03 km s -1 at 670 km to 0.02 km s -1 at 1,500 km, while L increases from 15 km at a depth of 670 km to 30 km at a depth of 1500 km. We confirmed L and σ through a series of synthetic tests, in which L and σ were varied to observe whether the retrieved models were consistent with various input models. Our one-dimensional initial velocity model was modified from PREM by smoothing the 220 km discontinuity with corrections for the crustal structure using the CRUST1.0 model and a recent result 4 from receiver function analysis beneath the OJP (Supplementary Fig. 7). We corrected for the anelastic effect using the attenuation model of PREM with a reference frequency of 1 Hz. The Moho depth of an initial model was fixed during the inversion. The mislocation of the Moho depth affected the inverted structure (Supplementary Note 3 and Supplementary Fig. 8); the obtained structure shallower than~50 km was less resolved because the phase velocities we used were longer than 30 mHz, which is less sensitive in shallow depths (See Supplementary  Fig. 4d). Checkerboard resolution tests suggested that the lateral resolution of the OJP is~4°for the isotropic structure and~7°for the anisotropic structure (Supplementary Fig. 9). We applied a spatial two-dimensional FFT low-pass Gaussian filter at 4°and 7°(grdfft in GMT 54 ). As the sensitivity of surface waves decrease with depth, we only show a three-dimensional radially anisotropic shear wave velocity structure for depths shallower than 300 km. We also performed a synthetic test to assess the influence of the subducting slab south of the OJP, which is indicated as a fast velocity anomaly. The recovered model suggested that the fast subducting slab anomalies have little effect on the mantle structure beneath the OJP ( Supplementary Fig. 10) A popular method for analysing the LAB is to use a reflected or converted seismic wave. In the northern part of the OJP, the depth of LAB is estimated to bẽ 70 km based on SS-precursor wave analysis 22 . However, SS-precursor waves do not sample the main part of the OJP due to the locations of the seismic stations and earthquakes. While S-wave receiver function analysis is also a common method for estimating the LAB depth 55 , a high level of horizontal component noise in the OJP array data prevents this type of analysis. Alternatively, we propose an approach that uses a negative gradient peak in the shear wave velocity as a proxy for the LAB depth. For this, we calculated the depth of the negative peak in the vertical gradient of the shear wave velocity obtained from the tomographic model to analyse the lateral variation of LAB depth. From a seismological perspective, the upper mantle structure can be described as a high-velocity lid (lithosphere) overlying the lowvelocity zone (asthenosphere), such that the negative peak in the vertical gradient can be used as a proxy for the LAB 20,53,56 . Thus, we defined the depth of the proxy as the depth of the LAB. To assess the depth resolution of the surface waves and the reliability of this proxy, we performed synthetic tests using a synthetic shear wave velocity model and by changing the thickness of a 3% high-velocity lid (Supplementary Note 1 and Supplementary Fig. 2). The recovered vertical profiles show a high-velocity lid at the top and a low-velocity body below. Although the sharp discontinuity was not recovered because surface waves are insensitive to this, the peak depth represents the depth of the discontinuity (i.e., LAB) for input LAB depths of 70-110 km within an accuracy of 15 km. The peak depths are greater than the input discontinuity depths for depths of 60-80 km and shallower for depths of 90-160 km, indicating that the deep LAB beneath the OJP is likely to be underestimated by 20 km or more whereas this is not the case for the surrounding regions (e.g., Nauru and Caroline).
Data. We collected seismic waveform data recorded by the OJP array for earthquakes with moment magnitudes (Mw) greater than 5.5 in the western Pacific Ocean to analyse the upper mantle shear wave structure beneath the OJP. As these data alone do not suffice for this purpose, we also collected seismic data recorded on land and on the seafloor in the western Pacific Ocean (Supplementary Fig. 1). These data are distributed by the Ocean Hemisphere network Project (OHP) and the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC), which were collected for earthquakes between 1990 and 2016. The study region in the western Pacific Ocean is similar to that of another recent study 17 in which only seismic data from land stations were used. Here, we used seafloor observation data recorded by BBOBSs as well as land data from a total of 161 stations, of which 115 are BBOBS stations (Supplementary Fig. 1). Vertical component seismograms recorded on the seafloor are contaminated by tilt noise caused by leaks of a high degree of noise in the horizontal components induced by seafloor currents 57 and compliance noise induced by the infragravity waves that are also recorded by seafloor pressure data 58 . We removed this noise from the BBOBS station data in the OJP array by applying a spectral transfer function between the vertical and horizontal components and between the vertical and pressure components, respectively.