A Mid‐Lithospheric Discontinuity Detected Beneath 155 Ma Western Pacific Seafloor Using Sp Receiver Functions

This study probes the lithosphere‐asthenosphere system beneath 155 Ma Pacific seafloor using teleseismic S‐to‐p receiver functions at the Pacific Lithosphere Anisotropy and Thickness Experiment project ocean‐bottom‐seismometers. Within the lithosphere, a significant velocity decrease at 33–50 km depth is observed. This mid‐lithospheric discontinuity is consistent with the velocity contrast between the background mantle and thin, trapped layers of crystallized partial melt, in the form of either dolomite or garnet granulite. These melts possibly originated from deeper asthenospheric melting beneath the flanks of spreading centers, and were transported within the cooling lithosphere. A positive velocity increase of 3%–6% is observed at 130–155 km depth and is consistent with the base of a layer with partial melt in the asthenosphere. A shear velocity decrease associated with the lithosphere‐asthenosphere boundary at 95–115 km depth is permitted by the data, but is not required.


Introduction
The seismic structure of the upper mantle beneath oceanic seafloor has the potential to constrain models of the thermal evolution of the oceanic lithosphere, melt emplacement within the lithosphere, and the presence of partial melt in the asthenosphere.We examine these processes beneath the PLATE (Pacific Lithosphere Anisotropy and Thickness Experiment) array of ocean-bottom seismometers (OBS).These stations lie in an area southwest of the Shatsky Rise that has some of the oldest identifiable magnetic anomalies (150-160 Ma old) and is largely devoid of large seamounts (Figure 1).
The cooling and thickening of the lithosphere as it moves away from spreading centers has been described by two conductive models: the half-space cooling (HSC) model in which lithosphere thickness, seafloor depth, and inverse heat flow increase in proportion to the square root of seafloor age; and the cooling plate model (CPM), in which these quantities asymptotically approach fixed values at infinite time due to convective supply of heat to the base of the lithosphere (Parsons & Sclater, 1977;Stein & Stein, 1992).Although some seismic studies support the HSC model, others conclude that the CPM or a more complex thermal evolution is needed at older seafloor ages (e.g., Beghein et al., 2019;Fischer et al., 2020;Kawakatsu et al., 2009;Ma & Dalton, 2019;Maggi et al., 2006;Ritzwoller et al., 2004).
A significant number of receiver function (RF) and passive and active source reflected wave studies have detected velocity gradients beneath the oceans that have been interpreted as representing the lithosphere-asthenosphere boundary (LAB) (e.g., Kawakastu et al., 2009;Mehouachi & Singh, 2018;Rychert & Harmon, 2018;Rychert & Shearer, 2011;Rychert et al., 2018;Schmerr, 2012;Stern et al., 2015;Tharimena et al., 2017).However, conductive thermal models produce a LAB thermal gradient that is too gradual to produce seismic phases that are comparable to the data (Fischer et al., 2010(Fischer et al., , 2020;;Rychert & Harmon, 2018), suggesting an additional mechanism at the LAB to sharpen the velocity gradient, such as partial melt from hydrated and/or carbonated asthenosphere (Kawakastu et al., 2009;Takeo et al., 2013).This view is supported by RF studies that detect a positive velocity gradient at ∼150 km depth beneath some oceanic regions, which can be explained as the base of a partially molten asthenospheric layer (e.g., Hua et al., 2023;Rychert et al., 2013Rychert et al., , 2018)).
Localized velocity gradients within the lithosphere, often called mid-lithospheric discontinuities (MLDs), have also been observed in converted/reflected wave studies.For example, an active-source experiment conducted beneath 128-148 Ma northwestern Pacific seafloor, east of the Shatsky rise, imaged wide-angle reflectors in the depth range of 37-59 km, and interpreted them as frozen-in melt pockets embedded within the lithosphere (Ohira et al., 2017).RF studies have also detected phases within a similar depth range (Hannemann et al., 2017;Kawakatsu et al., 2009).Across the Pacific seafloor, SS precursors have found structures at depths of 40-70 km depth (Rychert & Shearer, 2011;Schmerr, 2012;Tharimena et al., 2017) that are consistent with widely distributed MLDs.Numerical models show that volatile-rich partial melt can be emplaced in the mid-lithosphere at the flanks of mid-ocean ridges, and transported away from the spreading center along with the cooling plate (Keller et al., 2017).This study provides a unique opportunity to probe these processes in an old end-member of relatively undisturbed oceanic lithosphere and asthenosphere.

Data and Method
S-to-p (Sp) converted waves, which are generated at localized velocity gradients beneath the station, arrive earlier than the direct S arrival, avoiding the overprinting of crust and mantle phases by direct phase reverberations in water, sediment, and crustal layers which is typical with P-to-s (Ps) converted waves.In the analyzed PLATE array, the western and eastern sub-arrays each consisted of eight stations, but five OBS were lost or provided no data, and only seven returned all the three components of motion required for RF studies (Figure 1).Teleseismic events were chosen with magnitudes ≥5.7 and epicentral distances between 55°and 80°in the USGS National Earthquake Information Center (NEIC) global catalog.
In the data processing, we first removed the instrument response from the waveforms and integrated to displacement.We reoriented the horizontal components based on Rayleigh wave polarization (Stachnik et al., 2012), and then removed tilt noise and compliance noise from the vertical component (Bell et al., 2015) with transfer functions calculated from 5-hr noise windows prior to the P arrival (Rychert et al., 2018).The waveforms were bandpass filtered (7-25 s period) to avoid the long periods dominated by ocean noise while retaining most of the energy of the direct S arrival and converted phases.The choice of this filter, similar to Rychert et al. (2013), also circumvents sediment reverberations which are predominantly at shorter periods (Text S1 in Supporting  (Nakanishi et al., 1992).Information S1). S arrivals with less than a 10 s difference relative to the S arrival prediction from the TauP toolkit (Crotwell et al., 1999;Kennett et al., 1995;Rychert et al., 2018) on the radial component were picked.
We calculated Sp RFs by deconvolving the radial from the vertical component.Using the vertical components to represent the relatively small amplitude converted phases, rather than reconstructing the converted phases from radial and vertical components, avoids contamination by the higher noise levels on the OBS horizontals and the stronger P-wave water column reverberations on the radial component (Text S1 in Supporting Information S1).Sp RFs were obtained using iterative time-domain deconvolution (e.g., Ligorría & Ammon, 1999) that deconvolved a ∼13 s radial window from a 50 s vertical window using up to 100 iterations (Figure S3 in Supporting Information S1, upper panel), but stopping when the improvement in fit between iterations was <1% (Chen et al., 2020).RFs that passed quality control checks (Text S2 in Supporting Information S1) were bootstrapped and stacked (Abt et al., 2010), weighting individual RFs by waveform-specific signal-to-noise ratios.

The PLATE Sp RFs
The PLATE stations yielded 16 Sp RFs (Figure 2a) that passed quality control from five OBSs, showing broadly consistent phases between 0 and 20 s.The stack of the RFs (Figure 2b), which we interpret as representing the PLATE region, contains two positive phases at ∼1.3 and ∼18 s and one negative phase at ∼6.5 s.A synthetic test of the effects of background noise (Text S3 in Supporting Information S1) demonstrates that the observed phases that exceed two standard deviation uncertainties represent real structure rather than random noise.We also verified that simple azimuthal anisotropy has a minimal impact on the observed RF phases (Text S4 in Supporting Information S1), and thus we focused on vertical isotropic velocity gradients when modeling the data.

A Simple Structure With Sharp Boundaries
Based on the stack of multiple stations (Figure 2b), we considered a range of possible models for the velocity gradients beneath the PLATE array.To test these models, we compared their predicted synthetic RF stacks (using ray parameters corresponding to the data) to the observed RF stack.The synthetic seismograms were generated using Telewavesim (Audet, 2016) which is capable of handling reverberations from thin, low-velocity sedimentary layers and the water column, including those generated by the converted phases that appear in the Sp RFs.To match the frequency content of the OBS data, the phase peaks from the synthetic seismograms were identified and convolved with a Gaussian function (half-width of 1 s).We then calculated synthetic RF stacks using processing that was identical to that applied to the real data.
We first tested simple models that include a ∼6 km water column, a 50-m-thick sedimentary layer with Vs of 0.1 km/s and Vp of 1.8 km/s (Kim et al., 2022(Kim et al., , 2023)), a crust and uppermost mantle model from ambient noise cross-correlation (Takeo et al., 2014) and two step-function velocity boundaries in the mantle.We determined the depths and velocity contrasts of the two mantle boundaries that minimize misfit with the observed stack through forward modeling (Figure 3a), finding values of 44 km ( 3%) and 146 km (3%).The synthetic RF stacks overpredict the amplitude of the shallowest positive phase, which represents the combined effects of the water layer, sediments and Moho, suggesting that velocities in the lower crust of our model are too low.However, in the individual RFs with high signal-to-noise ratios, this phase amplitude is higher, suggesting its amplitude in the stack may also be artificially low due to noise.In any case, at the relatively long periods represented by the RFs in this study, adjusting shallow layer parameters to better match the amplitude of the shallowest phase does not significantly impact the interpretation of the deeper mantle phases, and we chose to retain the shallow layer parameters that are constrained by independent studies.Although this simple model matches the observed RF stack, a single velocity decrease of 3% at ∼44 km is difficult to reconcile with surface wave constraints on absolute shear-wave velocities in old Pacific mantle (e.g., Kawano Between 130 and 155 km depth, a 6.2% gradual velocity increase provides a better fit to the deeper positive phase.(c) Similar to panel (b), except using a dolomite layer from crystallized carbonated silicate melts.A 4.5% gradual velocity decrease between 95 and 115 km depth is possible at the lithosphere-asthenosphere boundary.At greater depth, a ∼5% gradual velocity contrast is needed to match the deepest observed phase.et al., 2023;Nishimura & Forsyth, 1989;Takeo et al., 2018).We therefore turn to models where the MLD is created by layers of limited thickness.

Crystallized Partial Melts as the Origin of the Oceanic MLD
On the eastern flank of the Shatsky rise, multiple MLDs were inferred to be frozen partial melt pockets that had originated at the flank of a mid-ocean ridge as basalt (MORB) and horizontally aligned in response to shear deformation (Ohira et al., 2017).Trapped basaltic melts, which would be crystallized in old oceanic lithosphere given its pressure-temperature conditions (Hirschmann, 2010), are the first source of low velocity material that we considered.Based on the phase diagram for MORB (Hacker et al., 2003), these crystallized melts would be garnet granulite at ∼50 km depth, which would create a velocity reduction of ∼9.6% relative to the background oceanic lithosphere (Abers & Hacker, 2016).
To evaluate how garnet granulite would affect the RFs, we conducted waveform modeling for a velocity model in which the lithosphere contains three 1-km-thick layers of garnet granulite at 50-55 km depth.RF synthetics for this model agree well with the observations (Figure 3b).The predicted negative amplitude MLD phase appears at a slightly shallower depth due to the cumulative interfering effects of the upper and lower boundaries of the garnet-granulite layers.At 100 km depth, a sharp boundary with a 3% velocity decrease was included to compensate for a small positive phase created by water reverberations (Figure S6 in Supporting Information S1).The data do not require this velocity decrease (taking into account two standard deviation uncertainties), but they do permit it.Between 130 and 155 km depth, a 6.2% gradual velocity increase provides a better fit to the deeper positive phase in the RF stack than the step-function increase in the first model (Figure 3a).
Alternative models involving volatile-rich melts in the oceanic lithosphere-asthenosphere system have been predicted by experimental data (Dasgupta, 2018;Hirschmann, 2010) and numerical simulations (Keller et al., 2017).Low-degree partial melt containing CO 2 and H 2 O can ascend into and be trapped within young oceanic mantle lithosphere at depths of 20-60 km, and then be transported within the cooling mantle lithosphere as it ages (Keller et al., 2017).In these models, the uppermost mantle lithosphere is largely free of trapped melts as they migrate upward along the sloping LAB to form the crust.CO 2 -rich melts would crystallize to a form of dolomite (Dasgupta, 2018;Hirschmann, 2010).
We tested models that contain one 1-km-thick layer of dolomite embedded within the oceanic mantle lithosphere.The elastic properties of dolomite at the pressure-temperature conditions of the mid-lithosphere beneath the old oceanic seafloor have not been directly measured.Instead, dolomite properties at low temperatures and pressures (Bakri & Zaoui, 2011;Humbert & Plicque, 1972;Mavko et al., 2009;Peng & Mookherjee, 2020) were extrapolated to mid-lithospheric conditions using the depth-dependent variations of two similar species, calcite and magnesite (Abers & Hacker, 2016), producing a shear velocity drop of 3.1%-4.2%relative to shallow, cold dolomite properties, and a 12.7% velocity reduction relative to the surrounding mantle at 50 km depth.Synthetic RFs (Figure 3c) show that one 1-km-thick dolomite layer matches the data well.In this case, a 4.5% velocity decrease between 95 and 115 km depth that agrees with surface wave constraints in old oceanic mantle (Nishimura & Forsyth, 1989) is permitted by the data, potentially corresponding to the LAB, although it is not required.A 5.0% velocity increase over 25 km at a depth of ∼140 km provides a good match to the deeper observed positive phase.
It should be noted that many other models could fit the data acceptably.There is no requirement that the trapped melt layers be regularly spaced vertically (e.g., Figure 3b).Lateral variations in the depths, abundance, and thickness of the trapped melt layers, as suggested in geodynamic modeling (Keller et al., 2017) and observations (Ohira et al., 2017), could contribute to variability in individual RFs.The MLD phases are much clearer in the eastern array, which may be partially attributed to lower noise levels, but also could be due to proximity to the Shatsky Rise (Figure 1).The Shatsky Rise formed from the interaction of a hotspot and the Pacific-Farallon spreading center (Nakanishi et al., 1999) beginning just a few Ma after the seafloor beneath the eastern stations was created, so more abundant and volatile-rich melt could have been trapped at that location.In summary, the observed MLD can successfully be explained by the crystallization of partial melts, either with (Figure 3c) or without (Figure 3b) high volatile content.This result is consistent with models (Keller et al., 2017;Ohira et al., 2017) in which volatile-rich partial melts were produced outside the melt focusing region beneath the Geophysical Research Letters 10.1029/2024GL108347 CHEN ET AL. mid-ocean ridge, emplaced in the mantle at the flanks of the ridge, and eventually transported to the western Pacific along with cooling lithosphere (Figure 4).
The observed MLD not only falls within the depth range of nearby wide-angle reflectors (Ohira et al., 2017), but also broadly coincides with SS precursor arrivals from depths of 35-70 km that are found beneath Pacific seafloor (Rychert & Shearer, 2011;Schmerr, 2012;Tharimena et al., 2017).Possible MLD arrivals can also be seen in Sp RFs from 130 Ma Pacific seafloor (Kawakatsu et al., 2009) and in Ps RFs from 75 to 85 Ma in the eastern mid-Atlantic (Hannemann et al., 2017).These findings suggest that oceanic MLDs are not only local features but are possibly widely distributed.In Hannemann et al. (2017), the MLD and the deeper LAB phases are much clearer after a 2-40 s bandpass filter, but become much more ambiguous with a 4-40 s bandpass that may lack enough high-frequency content to isolate the neighboring phases.This comparison may explain why NoMelt OBS data (70 Ma seafloor) detected only a single negative discontinuity at ∼70 km depth from Sp RFs using a 4-35 s filter (Mark et al., 2021), although local active-source data detected a sharp discontinuity at depths of 25-45 km (Lizarralde et al., 2018).Alternatively, the layers detected by the active source data may be too thin to be detected by the RFs.In older oceanic mantle, assuming that MLDs have frozen into similar depths, the thicker lithosphere allows for the isolation of MLD and LAB phases, as observed in this study.Beneath younger seafloor whose age is less than ∼30 Ma, discontinuity detections from SS precursors and RFs are closer in depth (Fischer et al., 2020;Rychert & Harmon, 2018).At these young ages, partial melt emplaced at the flanks of the mid-ocean ridge would result in an MLD signal that overprints the shallow LAB arrival, potentially complicating estimates of LAB depth in RFs.

Partial Melt Within the Asthenosphere
The positive phase we observe at 130-155 km depth, with a peak at 140 km, is consistent with positive velocity gradients that have been previously imaged by Sp RFs beneath the young seafloor in Cascadia (Rychert et al., 2018), beneath Hawaii (Rychert et al., 2013) and in high temperature regions of the upper mantle globally (Hua et al., 2023).This positive velocity gradient is consistent with the base of an asthenospheric layer that contains partial melt.Prior modeling of global observations (Hua et al., 2023) has shown that partial melt above this boundary is necessary to produce a sufficiently vertically localized velocity gradient.We also found that a relatively sharp velocity gradient is required to produce the positive .Schematic model of the lithosphere-asthenosphere system for young-to-old oceanic mantle.Yellow(orange) background: normal oceanic lithosphere(asthenosphere).Red(light-blue) lines: melt(solid) flow from the numerical simulations of Keller et al. (2017).White arrows: spreading direction.Gray patches: partial melt in young oceanic lithosphere.Brown patches: crystallized melts such as dolomite and garnet granulite within the middle of the lithosphere.Zoomed-in circle: grain-scale partial melt at crystal boundaries in the asthenosphere.
phase observed at ∼140 km, but that a 5.0%-6.2%velocity increase over 25 km provides a better fit (Figures 3b and 3c), compared to a velocity step at a single depth (Figure 3a).This result implies an upward increase in partial melt that produces a distributed downward velocity increase within a low-velocity partially molten asthenospheric layer (Figure 4), consistent with high attenuation in the asthenosphere beneath the PLATE stations inferred from P waves (Booth et al., 2014).The melt fraction would be 0.5%-3.3%,given the unknown melt geometry (e.g., Chantel et al., 2016;Hammond & Humphreys, 2000;Kawakatsu et al., 2009;Takei, 2002).These findings are in turn consistent with models for the oceanic asthenosphere where small amounts of partial melt are stable on a widespread basis due to the presence of H 2 O and CO 2 (e.g., Dasgupta, 2018;Gardés et al., 2020;Hirschmann, 2010), the melt is interconnected, and observed localized gradients in seismic velocity correspond to zones of enhanced volatile content (Gardés et al., 2020).The positive velocity gradient we observe at 130-155 km could be the lower boundary of a zone of higher partial melt content, or, if a redox freezing front occurs as shallow as ∼150 km (Stagno & Frost, 2010), the velocity gradient could mark the lower boundary of partial melt in the asthenosphere.
A LAB velocity decrease of 4.5% at 95 km (Figures 3b and 3c) is consistent with the PLATE RF stack, although it is not required.A velocity decrease associated with the oceanic LAB has previously been inferred from Sp RFs in other regions of the Pacific, for example, at 65 km depth with a velocity decrease greater than 4.5% at 70 Ma (Mark et al., 2021), and at 80-85 km depth with a 7%-8% velocity decrease at 130 Ma (Kawakatsu et al., 2009).The LAB at ∼75 km depth at ∼145 Ma inferred from Ps RFs (Olugboji et al., 2016) is smaller than the ∼95 km depth from this study at 155 Ma.However, all of these estimates are more consistent with predictions from CPMs rather than the HSC model (e.g., Rychert et al., 2018).The 4.5% velocity decrease at the LAB suggested in this study could be produced by 0.4%-2.4% melt fractions.Such a boundary could be created by partial melt ponded at the top of the asthenosphere (e.g., Kawakatsu et al., 2009;Rychert et al., 2018;Wang et al., 2020), due to the lower mobility of melt within the cold lithosphere (Sakamaki et al., 2013).

Conclusions
The oceanic lithosphere-asthenosphere system beneath 155 Ma western Pacific seafloor was illuminated by Sp RF using the PLATE OBS data.In the lithosphere, significant negative phases are observed at depths of 33-50 km.These mid-lithospheric depths are similar to those of wide-angle reflectors found in nearby active source experiments, and to the results of other RFs from OBS data and SS precursor studies across the Pacific seafloor, implying that oceanic MLDs may be widespread.The MLDs can be explained by the velocity contrast between background oceanic mantle lithosphere and trapped crystallized carbonated/basaltic melts.The melt was likely emplaced beneath the flanks of spreading centers, eventually crystallizing to form dolomite and/or garnet granulite within old oceanic lithosphere.
A velocity increase of 3%-6% is observed at 130-155 km, and is best explained as the lower boundary of a zone of greater partial melt fraction in the asthenosphere due to higher volatile content.A velocity decrease of ∼3%-4.5% at 95 km is consistent with the data, but is not required.Assuming this boundary represents the LAB, it is more consistent with a CPM for the oceanic lithosphere, versus HSC, and it could signify the upper boundary of asthenosphere that bears a small percentage of partial melt.

Data Availability Statement
Seismic data were accessed through the Incorporated Research Institutions for Seismology (IRIS) Data Management System.This work employed data from the Pacific Lithosphere Anisotropy and Thickness Experiment (PLATE) project under the Z6 network code (No DOI is registered for this network) and all the waveforms are available via the IRIS Data Management Center under the Z6 network code (http://www.fdsn.org/networks/detail/Z6_2009/).The seafloor age data shown in Figure 1 were downloaded from EarthByte (https://www.earthbyte.org/category/resources/data-models/seafloor-age/).The code for synthesizing teleseismic body-wave propagation was downloaded from the Telewavesim webpage (https://paudetseis.github.io/Telewavesim/).The physical properties of garnet granulite, calcite, and magnesite were calculated using the Excel worksheet and macro of Abers and Hacker (2016).The melt and solid flow paths shown in Figure 4 are from the numerical simulations presented in Figure 1a of Keller et al. (2017).

155
Ma western Pacific crust are measured by S-to-p receiver functions using the Pacific Lithosphere Anisotropy and Thickness Experiment oceanbottom seismometers data • A mid-lithospheric discontinuity at ∼40 km depth is observed and can be explained by layers of crystallized partial melt • A partial melt-bearing layer from 95 to 155 km depth is consistent with the receiver function observations Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.Maps showing locations of the Pacific Lithosphere Anisotropy and Thickness Experiment (PLATE) ocean-bottom seismometers and analyzed teleseismic events.Triangles denote PLATE stations; orange show broadly consistent Sp receiver functions (RFs); gray yielded no high-quality RFs; white have pressure data only.Circles mark the epicenters of teleseismic events; orange show identifiable S arrivals for RFs; gray are discarded.Left: Background colors show the age of seafloor.Right: Background colors show bathymetry near the PLATE deployment and the Shatsky rise.Gray lines represent seafloor age isochrons in increments of 5 Ma(Nakanishi et al., 1992).

Figure 2 .
Figure 2. (a) Individual time-domain Sp receiver functions (RFs) from Pacific Lithosphere Anisotropy and Thickness Experiment stations.The 16 traces are sorted by ray parameter and are reversed in time and flipped in amplitude to match the Ps RF convention.Labels list (right) station, event time, back-azimuth, signal-to-noise ratio, and (left) ray parameter.(b) The stack of the 16 traces.Dark and light gray shading on either side of the mean respectively denote ±1 and ±2 standard deviation uncertainties.Positive(negative) phases filled by red(blue) represent velocity increases(decreases) with depth.Background light red(blue) shading highlights the time ranges of significant positive(negative) phases.

Figure 3 .
Figure 3. Modeling of Sp receiver function (RF) stacks for the Pacific Lithosphere Anisotropy and Thickness Experiment (PLATE) region.Left panel: shear velocity structure used in waveform modeling.Right panel: synthetic RF (red line) and stack from PLATE stations in the background for reference.(a) A simple velocity model that matches the data.This model includes sharp velocity contrasts at depths of 44 km ( 3%), and 146 km (3%).(b) A case with crystallized basaltic melts, including three 1-km-thick layers of garnet granulite.A 3% velocity decrease at 100 km depth was included to cancel out a positive phase yielded by water reverberations.Between 130 and 155 km depth, a 6.2% gradual velocity increase provides a better fit to the deeper positive phase.(c) Similar to panel (b), except using a dolomite layer from crystallized carbonated silicate melts.A 4.5% gradual velocity decrease between 95 and 115 km depth is possible at the lithosphere-asthenosphere boundary.At greater depth, a ∼5% gradual velocity contrast is needed to match the deepest observed phase.

Figure 4
Figure 4. Schematic model of the lithosphere-asthenosphere system for young-to-old oceanic mantle.Yellow(orange) background: normal oceanic lithosphere(asthenosphere).Red(light-blue) lines: melt(solid) flow from the numerical simulations ofKeller et al. (2017).White arrows: spreading direction.Gray patches: partial melt in young oceanic lithosphere.Brown patches: crystallized melts such as dolomite and garnet granulite within the middle of the lithosphere.Zoomed-in circle: grain-scale partial melt at crystal boundaries in the asthenosphere.