Seismic evidence for melt-rich lithosphere-asthenosphere boundary beneath young slab at Cascadia

Domain: earth_science
Published: 2024-04-25
Editorial outcome: major_revision

Abstract
The Lithosphere-Asthenosphere Boundary (LAB) beneath oceanic plates is generally imaged as a sharp seismic velocity reduction, suggesting the presence of partial melts. However, the fate of a melt-rich LAB is unclear after these plates descend into the mantle at subduction zones. Recent geophysical studies suggest its persistence with down-going old and cold slabs, but whether or not it is commonly present remains unclear, especially for young and warm slabs such as in the Cascadia subduction zone. Here we provide evidence for its presence at Cascadia in the form of a large (9.8 ± 1.5%) decrease in shear-wave velocity over a very small (<3 km) depth interval. Similarly large and sharp seismic velocity reduction at the bottom of both old and young slabs, as well as along the base of oceanic plates before subduction, possibly represents widespread presence of melts. The melt-rich sub-slab LAB may strongly influence subduction dynamics and viscoelastic earthquake cycles.

Introduction
The Lithosphere-Asthenosphere Boundary (LAB) beneath oceanic plates marks a sharp decrease in seismic velocity. The sharp decrease, loosely but commonly referred to as a negative “seismic discontinuity”, is widely thought to indicate the presence of partial melts 1–5, although there are competing interpretations 6,7. A melt-rich LAB may mechanically decouple the lithospheric plate from the underlying asthenosphere and thus facilitate plate motion 8,9. If it continues into subduction zones, similar decoupling would occur beneath the subducting plate (called the slab) and fundamentally influence subduction dynamics 10 and viscoelastic earthquake cycles 11. Here we refer to the LAB of an oceanic plate before subduction as the plate-LAB and that after subduction as the slab-LAB (Fig. 1). However, it is by no means clear whether the slab-LAB is as common as the plate-LAB. The slab-LAB has been well imaged in the Japan and New Zealand subduction zones 1,8, and in both cases, the age of the subducting plate is very old 12 (~120–130 Ma; Fig. 1a). In contrast, its presence in regions where young plates are subducting remains unclear 3. A magnetotelluric study off the Middle America trench could resolve a melt-rich LAB beneath the young (~23 Ma) Cocos plate but not a slab-LAB in the subduction zone 13. The plate-LAB is seismically detected beneath almost the entire young (<10 Ma) Juan de Fuca plate, but so far, there is no reported evidence for the slab-LAB in the Cascadia subduction zone 14. Whether or not the presence of a melt-rich slab-LAB depends on slab age is a question with important geodynamic implications.

Cascadia is an end-member warm-slab subduction zone where the incoming Juan de Fuca plate is young (<10 Ma) and warm (Fig. 1). Using teleseismic converted waves generated at seismic discontinuities (referred to as receiver functions, RFs), previous studies reported various discontinuities related to the subducting Juan de Fuca slab, including those associated with the fluid-rich subducting oceanic crust, slab-Moho, and anisotropy within the slab crust and/or mantle 15–22. However, no seismic discontinuity associated with slab-LAB has been reported thus far. Seaward of the subduction zone, Rychert et al. detected 14 the plate-LAB throughout the incoming Juan de Fuca plate and considered it strong evidence for the presence of partial melts along the base of the plate. However, whether the inferred melt-rich LAB persists after the plate is subducted remains an unresolved question.

To address this question, we examine P-to-S (Ps conversion) RFs by taking advantage of dense seismic arrays at Cascadia (Fig. 1). One array consists of broadband seismometers with a station spacing of ~5 km and operated in 1993–1994; the other consists of short-period nodal seismometers with a station spacing of ~500 m and operated for one month in 2017 (Supplementary Fig. 1). Data from these arrays have been used in previous studies for investigating the shallow sedimentary basin structures, mid-crustal Conrad discontinuity, continental Moho, slab geometry, fluid-rich slab crust and its dehydration processes, and fossil anisotropy within the slab mantle 15,17,20–26. Below the slab-Moho, some patchy negative RF signals have been detected previously 21,22,24,25, which might be indicative of the presence of a slab-LAB. However, these studies did not consider them to be structural signals because of their intermittent appearance and concerns about artifacts arising from the side-lobe effect of RF deconvolution and multiple reverberations from shallow structures. One study 21 tentatively attributed these negative signals to anisotropy within a very thick slab, yet their findings might be influenced by cross-mode artifacts resulting from the usage of a multiple-mode conversion imaging technique 20. In our present study, by employing a recently developed Bayesian array-based coherent receiver function (CRF) imaging technique 27,28, and conducting multiple frequency analysis and extensive synthetic tests, we are able to detect the slab-LAB from these data in the form of strong and sharp seismic velocity reduction. We consider the imaged slab-LAB to be compelling evidence for the continuation of the melt-rich plate-LAB after subduction.

Results
CRF imaging and comparison with conventional RF analyses
In the distance range (30°–90°) optimal for CRF analyses, the arrays recorded seismic waves generated by distant earthquakes mainly in the southeast, northwest, and southwest directions (Supplementary Fig. 1). Because the Juan de Fuca slab dips to the east, the Ps waves from the southeast are more effective than those from the other directions in constraining slab-parallel seismic discontinuities (Supplementary Figs. 2, 3). We, therefore, mainly use the Ps waves from this direction for our analyses. The results based on waves from the other directions provide corroborative information. We process the data in multiple frequency bands (0.05–0.5 Hz, 0.15–1.0 Hz, 0.25–1.5 Hz, 0.35–2.0 Hz) and employ the CRF technique to conduct a joint analysis of the broadband and short-period data (see Methods; Supplementary Figs. 4–6). The denser spatial coverage of the short-period array allows us to better constrain shallow structures, which aids in constraining deep structures (see Methods).

Compared to previous RF studies 15,22,26, which predominantly focused on relatively low-frequencies (up to 0.3 Hz), we use multiple frequency analyses mentioned above to check the persistence of discontinuities. The frequency dependence of the amplitudes of Ps phases helps to appraise the robustness and sharpness of the discontinuities 29,30. Although increasing the frequency content of the data increases the vertical resolution of subsurface structures, high-frequency analysis is prone to noise and local scattering. To overcome this challenge, we employ the array-based CRF technique, which can leverage the coherency of the wavefield recorded by a dense array to suppress incoherent noise and local scattering (Fig. 2). The CRF method is also less susceptible to interference to the target discontinuity image from nearby strong discontinuities, as shown by the comparison of subsurface images obtained using the CRF and conventional RF methods (Fig. 2 and Supplementary Fig. 7). Another advantage of the CRF is that the reliability of imaged structures is measured using probability distribution, allowing objective assessment of structure identification and interpretation 28.

The CRF images are originally obtained in the time domain. To convert time-to-depth, we employ the seismic velocity model from the tomography study of ref. 31. Although the spatial resolution of this model is much higher than other available regional tomography models 32–34 (Supplementary Fig. 8), its velocity field is still inevitably damped by the smoothing scheme in tomographic inversion and thus distorts the depths of discontinuities in the CRF imaging (Supplementary Figs. 9, 10). As will be detailed in the following section, correction for this distortion is an integral part of the interpretation of our CRF results. The CRF image based on seismic waves from the southeast is shown in Fig. 3. Seismic waves from the northwest are theoretically less optimal for resolving east-dipping discontinuities (Supplementary Fig. 2), but they still result in comparable CRF images for an area immediately to the north (Supplementary Fig. 11), offering further support to the main imaged features shown in Fig. 3. Because of the potential depth distortion discussed above, we refer to the depths shown in Fig. 3 as apparent depths.

The slab-top and slab-Moho discontinuities
Shallower than 60 km, our CRF image reproduces several primary structure boundaries found in previous studies 16,20,21,25,35 and thus demonstrates the efficacy of the CRF technique. At 0–10 km apparent depths, prominent positive discontinuities (Fig. 3a) roughly delineate known sedimentary basins 26,35. The negative intra-crustal discontinuity observed at about 10–20 km depth has also been observed by previous studies 20,25 and may be related to the accumulation of fluids within the continental crust 36,37. In the 20–50 km range, two subparallel east-dipping seismic discontinuities with opposite polarity delineate the slab geometry and the low-velocity subducting crust. By utilizing higher-frequency data in comparison to previous studies 15,24, our proposed slab geometry contains more details; however, it is important to note that the tomographic model used in time-to-depth conversion may introduce slight distortions to the slab geometry, both in our study and in previous ones. The apparent thickness of the low-velocity zone >10 km is larger than what is expected for a normal oceanic crust due to the aforementioned depth distortion caused by velocity smoothing in the tomography model that is used for time-depth conversion. If we follow previous studies in Cascadia to assume a fluid-saturated crust with a very low Vs (about 20–50% Vs reduction) and a high Vp/Vs ratio (about 2–3) 16,17, the thickness is corrected to be ~7–8 km (Supplementary Fig. 10). The CRF amplitudes and our synthetic tests indicate extremely large velocity contrasts across the top (~15–20%) and bottom (~35–55%) of the subducting crust (see details in Methods, Fig. 4), consistent with previous estimates 16,17. The decrease in the amplitudes of the Ps conversions at the slab-top and slab-LAB when the slab depth exceeds ~40 km (Fig. 3a) is likely due to the eclogitization process of the subducting oceanic crust 15,23. Eclogitization is accompanied by dehydration and thus releases aqueous fluid into the overlying forearc mantle wedge, leading to the well-known absence of a detectable continental Moho near the mantle wedge corner 15.

The slab-LAB discontinuity
The most important finding enabled by the newly employed CRF technique is the continuous negative discontinuity at about 60–80 km apparent depths, which we interpret as the slab-LAB (Fig. 3). A direct estimate of the actual depth of the negative discontinuity using observed Ps arrival times shows that it is located about 25 km below the slab-Moho or about 32 km below the slab-top (Supplementary Fig. 10). Patchy negative signals at relevant depths can be seen in some other seismic imaging studies at Cascadia 20,21,24,25, but it was not clear whether they were structural signals or imaging artifacts. Here, we are able to resolve confidently this discontinuity owing to a combination of multiple frequency analyses, the narrow back-azimuth range of earthquakes used in dipping discontinuity imaging, and the advantages of the CRF method. The Bayesian posterior probability distribution of the CRF phases shows a high confidence level for this discontinuity (Fig. 3b), and its presence is persistent across multiple frequency bands (Supplementary Figs. 5, 6).

Our synthetic tests verify that the strong negative signal marking this discontinuity cannot be an imaging artifact due to other structures (Supplementary Figs. 12–16). Basin reverberations occur at shallower depths, and multiples from the slab-top occur at much deeper depths (Supplementary Fig. 12). Multiples from a positive intra-crustal discontinuity (i.e., the Conrad discontinuity) at ~15 km depth in the overriding continental plate would exhibit a pair of positive and negative signals at a depth similar to the observed negative discontinuity (Supplementary Fig. 13). However, such a Conrad discontinuity is not observed in our CRF image (Fig. 3) and the conventional RF images (Supplementary Fig. 5). Although a negative intra-crustal discontinuity is imaged at a depth of 10–15 km in our study, it generates multiples with geometry inconsistent with the slab-LAB discontinuity (Supplementary Figs. 13, 14). More importantly, the RF images constructed with seismic waves from the northwest show a weaker negative signal than that obtained using the waves from the southeast (Supplementary Figs. 3, 15), and the arrival time of the negative phase decreases with the increase of epicentral distance (Supplementary Fig. 16), suggesting that the negative signal is related to the direct Ps phase rather than multiples. We further analyze the RFs as a function of back-azimuth using both synthetic data and real observations, and our results show that this negative signal is better explained by a dipping discontinuity rather than an anisotropy layer (Supplementary Figs. 17–19).

Sharpness and velocity reduction across the slab-LAB
An important attribute of the inferred LAB discontinuity is its sharpness, characterized by the magnitude of the velocity change (dVs) and the depth range (dZ) over which the change occurs (Fig. 4). The amplitude of the Ps phase in RFs is particularly sensitive to dVs, dZ, and the range of dominant frequencies used in the analysis. If dZ is broader than half of the dominant wavelength of the incident P-wave, the amplitude decreases dramatically with increasing frequency 8,29 (Supplementary Fig. 20). The Ps conversions of our inferred slab-LAB consistently exhibit an amplitude of 12±3% of the direct P-arrival across the entire frequency range (Fig. 4). Assuming a typical upper-mantle Vp of ~8 km/s and considering the waveform frequencies used in this study up to 2.0 Hz, the dZ related to this discontinuity should not much exceed ~2 km 29. Further synthetic tests indicate that the slab-LAB corresponds to a sharp discontinuity with dVs of 9.8 ± 1.5% over a depth range of 1.5 ± 1.5 km (Fig. 4 and Supplementary Fig. 21). Such a sharp velocity reduction within 3 km cannot be solely explained by variations in water abundance, temperature, or grain size 1,3,38. A recent study 39 suggests that hydration within the asthenosphere can substantially reduce seismic velocities, offering a potential explanation for the strong and sharp velocity reduction across the LAB. However, the effect of hydration on velocity reduction becomes less pronounced for seismic waveform frequencies considered in this study. Thus, the presence of partial melts stands out as the most likely mechanism for the strong and sharp velocity reduction at the slab-LAB. Similar to previous estimates 1,3,14, a dVs of ~9.8 ± 1.5% across the LAB can be attributed to the presence of ~1–4% melt fraction, depending on the melt geometry 40.

Methods
Data
Seismic data used in this study are from two nearly collocated dense linear arrays in central Oregon, western US, including a broadband (BB) seismic array (with a flat response from ~0.01–50 Hz; https://www.fdsn.org/networks/detail/XZ_1993/) with a station spacing of ~5 km operated from 1993/03/23 to 1994/11/23 and a short-period (SP) nodal seismic array (with a corner frequency of 5 Hz; https://www.fdsn.org/networks/detail/ZO_2017/) with a station interval of ~500 m operated from 2017/06/23 to 2017/08/01 (Supplementary Fig. 1). Compared to the BB array, the SP array is ten times as dense but is limited by a narrower frequency band and a short observation duration. However, after removing the instrument response, frequencies greater than 0.1 Hz are generally well recovered by the SP nodal stations 25,28,50,51. The longer temporal coverage of the BB data and the denser spatial coverage of the SP data complement each other in imaging subsurface structures, as explained in the following section.

Array-based coherent receiver function imaging
We use a recently developed Bayesian array-based coherent receiver function (CRF) imaging technique 27,28. The basic concept of the CRF is similar to that of the conventional single-station receiver function (RF) technique, but the CRF method includes array processing, in which the deconvolution at a single station involves constraints from nearby stations provided that the stations are closely located to detect similar subsurface structures sampled by seismic waves from multiple earthquakes in a narrow back-azimuth range. The basic assumption of the CRF method is that the teleseismic wavefields are coherent and the subsurface structures at nearby stations effectively share similar features. The dominant frequency of teleseismic P-waveforms is around 0.3–0.5 Hz with wavelengths longer than ~10–25 km, which ensures that signals recorded by the ~5 km distributed broadband stations show very similar waveforms, satisfying the coherent assumption required for CRF imaging. In addition, incorporating short-period data with a station spacing of approximately 500 m in CRF imaging improves the resolution of shallow structures, such as sedimentary basins. This, in turn, aids in constraining deeper structures and enables us to conduct higher-frequency CRF analysis. Using CRF, we can suppress incoherent noise and local scattering to retrieve weak signals. Another advantage of the CRF technique is that the deconvolution is performed using a trans-dimensional Markov chain Monte Carlo (McMC) inversion algorithm 52,53. Thus, the obtained RFs are presented in probability formats, enabling us to determine the reliability and uncertainties of the obtained RFs.

We first acquire waveform data from earthquakes with magnitude larger than 5.8 and epicentral distances from 30 to 90 degrees (Supplementary Fig. 1). The recorded seismic waves were generated by distant earthquakes mainly in the northwest, southeast, and southwest directions (Supplementary Fig. 1). In the presence of a dipping structure, RFs depend strongly on the direction of the incoming teleseismic waves 54 (Supplementary Figs. 2, 3). In our study area, the slab dips eastward. To optimize the imaging of the east-dipping slab-parallel seismic discontinuities 55, we mainly use waveforms from the southeast direction (earthquakes colored by orange in Supplementary Fig. 1) in the CRF imaging. The three-component seismograms are rotated into radial, transverse, and vertical components, windowed around the theoretical first P-arrivals, and filtered into four frequency bands (0.05–0.5 Hz, 0.15–1.0 Hz, 0.25–1.5 Hz, and 0.35–2.0 Hz). Compared to the low-frequency analyses (e.g., the 0.05–0.5 Hz used in this study and the up to 0.3 Hz used in previous studies), higher frequencies serve to reduce the wavelength of waveforms, thus increasing the vertical resolution of subsurface structure imaging. However, high-frequency analysis is easily affected by noise and local scattering. Nevertheless, the usage of the CRF method can better suppress incoherent noise and local scattering to improve subsurface structural imaging (Fig. 2), and the integration of multiple frequency data enables us to resolve the robustness and detail structures of the subsurface discontinuities. We then remove low-quality traces which have signal-to-noise ratios (SNR) less than 2 in the vertical component. Here SNR is defined as the ratio of the root mean squared amplitude over an 8-s window after the first P-arrival to that before the first P-arrival. After data preprocessing, 4073 teleseismic waveforms are retained for the BB array, and 2138 waveforms are retained for the SP array. The distribution of the P-to-S conversion points at a depth of 60 km, and the distribution of earthquakes in different back-azimuths are shown in Supplementary Fig. 1.

To image lithosphere-scale structure using CRF, we use a moving subarray of ~15 km with an incremental step for each station. The subarray generally includes 3 BB stations and 25 SP stations (Supplementary Fig. 4a). The chosen subarray size is a compromise between the number of earthquakes recorded at BB and SP stations, as well as the inter-station spacing of the BB and SP stations. We have tested subarray sizes ~5 km larger or smaller but obtained similar results. Within each subarray, the subsurface structure is parameterized as the number of CRF phases, the timing, amplitude, and slowness of each phase, in which the slowness is used to link the subsurface structure beneath different stations within the subarray. We employ the McMC technique to iteratively add, remove, or adjust CRF phases, until the likelihood function (or misfit function) is stabilized 28. Compared to the traditional time-domain iterative deconvolution, in which RF phases are fixed once added, McMC allows for dynamic updates of all phases throughout the process 28. The resultant, more accurate constraint on earlier CRF phases (shallow structures) helps to improve the quality of the later phases (deeper structures). For each subarray, we run McMC Bayesian inversion using 500 CPUs to sample the model space simultaneously and independently. A total of 10,000 iterations are performed on each CPU. The first 5000 iterations are regarded as the “burn-in” period to allow convergence (Supplementary Fig. 4b), and every tenth model in the subsequent 5000 iterations is used to generate the ensemble solutions (Supplementary Fig. 4b). To present the CRF in a similar manner to the conventional RF, we collect all the CRF signals with the posterior probability over a prescribed threshold (e.g., larger than one standard deviation of the mean, that is 68% confidence level) in the density map to construct the 1D trace that contains discrete delta functions with different timings, amplitudes, and slowness (Supplementary Fig. 4b). We then convolve the 1D trace with different Gaussian low-pass filters (Gaussian parameter a = 1.0, 2.0, 3.0, and 4.0) in each frequency band to generate the 1D CRF time series. We stack the 1D CRF waveforms (Supplementary Fig. 4c) from each subarray to produce the 2D stacking images (Supplementary Fig. 4d) via a slow common conversion point (CCP) stacking technique. Furthermore, we extract the posterior probability for the timing of CRF phases at each single-node subarray (Supplementary Fig. 4b), and we then stack the posterior probability distribution of the CRF phases to generate the 2D confidence images. Supplementary Fig. 6 shows the CRF imaging results obtained at multiple frequency bands. We also conduct CRF imaging using BB or SP data alone. Despite very different observation periods and instrument types, the CRF results constructed using the BB and SP data show remarkable agreement (Supplementary Fig. 6).

Time-to-depth conversion using velocity models
To perform the time-to-depth conversion by accounting for the effect of the 3D velocity structure, we use the Vs tomography model from ref. 31 for time-to-depth conversion. The corresponding P-wave velocity and density are obtained through empirical relations 56. However, using this model, as well as other global or regional tomographic models, in time-to-depth conversion results in an unrealistic thick slab-top low-velocity layer (LVL) with thickness >10 km and a relatively thin slab mantle of ~20 km (Supplementary Fig. 9), as explained in the main text. To better constrain the thickness of the slab-top LVL and the slab mantle, we pick arrivals related to the slab-top interface, slab-Moho, and slab-LAB, and conducted synthetic tests to understand the effect of velocity models in time-to-depth conversion (Supplementary Fig. 10). The observed arrival time difference between the slab-top interface and slab-Moho is about 1.3 ± 0.35 s. Following previous studies in Cascadia, assuming a fluid-saturated LVL with a very low Vs (about 20–50% Vs reduction) and a high Vp/Vs ratio (about 2–3) reduces the slab-top LVL thickness to ~7–8 km. Similarly, the observed arrival time difference between the slab-Moho and slab-LAB is 2.3 ± 0.45 s. The thickness of the slab mantle is thickened to ~25 km or larger by considering a 5–10% high-velocity subducting slab mantle and an abnormal Vp/Vs ratio (~1.75).

Synthetic tests
To verify the robustness of our results presented in the main text, we test the effects of shallow structures’ multiple reverberations on slab-LAB imaging in Supplementary Figs. 12–16. The synthetic seismograms are calculated with the Thomson–Haskell matrix method 57,58 using a 1D layered velocity model for each station that is extracted from the 2D velocity model. These synthetics are spatially arranged to mimic the effect of the 2D slab. Given the purpose of testing the effects of shallow multiples, not directly calculating 2D synthetic seismograms reduces the computational cost. Our previous studies have demonstrated the reliability of this approach 27,28. We conduct five types of synthetic tests. Model-1 consists of a 7 km thick low-velocity subducting oceanic crust and a 33 km thick high-velocity subducting oceanic mantle, where the transition from the subducting lithosphere to the underlying asthenosphere is smooth. The slab geometry is obtained from the Slab2.0 model 59 in central Cascadia to mimic the real situations (Supplementary Fig. 12). Model-2 largely follows the structural features in Model-1, but consists of a sharp velocity reduction at the base of the subducting slab (Supplementary Fig. 12). By comparing the synthetics obtained from Model-1 and Model-2, we conclude that the imaged strong negative slab-LAB signal cannot be an imaging artifact due to the multiples from the slab-top and slab-Moho (Supplementary Fig. 12). In Model-3, we test the effects of basin structures on slab-LAB imaging, concluding that the basin reverberations occur at shallower depths than the imaged slab-LAB (Supplementary Fig. 12). In Model-4 and Model-5, we further explored the effects of an intra-crustal discontinuity on slab-LAB imaging (Supplementary Fig. 13). Model-4 contains a positive Conrad discontinuity in the overriding continental plate at 15 km depth, while the Model-5 examines the negative intra-crustal discontinuity (at 10–15 km depth and at 60–130 km distance from coast) imaged in our study (Supplementary Fig. 13). In Model-4, the multiples of the Conrad discontinuity have similar arrivals as the direct Ps conversions of the slab-LAB. However, if the negative phase at around 8 s is the multiples of the Conrad discontinuity, we would expect to observe a similar magnitude positive phase just before the negative phase. In addition, the multiples from the Conrad discontinuity show different amplitude variations with respect to the incident directions (Supplementary Fig. 15), and the arrival time of the negative phase decreases with the increasing of epicentral distance (Supplementary Fig. 16). In Model-5, the multiples of the negative intra-crustal discontinuity exhibit different geometries compared to the imaged slab-LAB.

We further test the effect of an anisotropy layer on RFs using the synthetic data (Supplementary Figs. 17–19). The synthetic RFs are calculated using the Raysum code 54, which is a ray-based method to calculate synthetic seismograms for dipping discontinuities and anisotropic media. Analysis of observed RFs as a function of back-azimuths suggests that a dipping discontinuity provides a more plausible explanation for the negative signal observed at ~8 s.

Multi-frequency receiver function waveform modeling
In this section, we conduct multiple frequency RF analyses to constrain the seismic structural characteristics related to the discontinuities observed in the main text. The amplitude of a P-to-S converted phase in RFs is sensitive to the S-wave velocity contrast across the discontinuity (dVs), the depth range over which the velocity contrast occurs (dZ), and the range of dominant frequencies used in analyses 29 (Supplementary Fig. 20). In this study, we only focus on the discontinuities observed in the distance range of 30–60 km to constrain their seismic characteristics, as they are less affected by the shallow basin structures (Supplementary Fig. 6). We pick the amplitudes related to the slab-top, slab-Moho, and slab-LAB discontinuities shown in the CRF images obtained using the BB data in Supplementary Fig. 6, as the SP data are limited by a narrower frequency band. We then constrain the dVs and dZ through a grid search inversion, in which the misfit is defined as the difference between the observations and synthetics,

misfit = ∑_{i=1}^{4} (obs_i - syn_i)^2 / obs_i^2, (1)

where the i represents the observation or synthetic at a given frequency range.

For a given velocity model, the synthetic RFs are calculated with the Thomson–Haskell matrix method 57,58. We generate a group of synthetic RFs with different bandpass filters similar to that used for real observations. We then pick the amplitudes of RF waveforms associated with each discontinuity and correct the amplitude by assuming the slab geometry from the Slab2.0 model 59. In this region, the slab’s average dip angle and strike direction are 35° and 12°, respectively, resulting in an amplitude correction factor of 1.3 ± 0.1 based on synthetic tests calculated using Raysum 54. Supplementary Fig. 21 shows the relative misfit as a function of dVs and dZ. The relative misfit is defined as:

Relative Misfit = (Misfit_i - Misfit_min) / Misfit_min * 100%, (2)

where the Misfit_i represents the misfit of the ith model and Misfit_min represents the minimum misfit in the grid-search inversion. The velocity structures for the slab-top and slab-LAB are well constrained, as we can fit the multiple frequency RF waveforms simultaneously. The slab-Moho is not well constrained partially due to that the analysis of the Ps phase at low frequencies does not permit to distinguish the slab-Moho from the broader gradient of the high-velocity subducting slab mantle.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript by Wang et al. identifies a coherent negative phase in the receiver functions\nobtained from dense linear array of seismometers deployed at the Cascadia subduction\nzone. The authors attribute the negative phase to the lithosphere-asthenosphere boundary\n(LAB) at the based on the subducting slab. Due to its sharpness, the LAB is interpreted to be\nat the melt-rich condition. The manuscript is overall well-written, and its novelty is clear.\nHowever, the method has some points needed to be further clarified. The requested\nrevisions have the potential to impact the results, and thus, I recommend categorizing these\nchanges as 'major revisions'.\n1. Probability distribution\nThe procedure for calculating the probability distribution needs to be described with greater\nclarity. From my perspective, the process involves two steps: firstly, constructing a 1-D\nprobability distribution for pulse presence by marginalizing the posterior probability across\ntime for each sub-array. Subsequently, this 1-D time series of probabilities is either stacked\nor averaged within each bin on the cross-section using the CCP stacking method. It is\nimportant to note that, during these processes, there is no requirement for 'normalizing' the\nprobability.\nHowever, the manuscript repeatedly mentions normalization at both stages of analysis. The\ninitial reference to normalization is found in the color scale of Fig S4(b), where it states\n\"Normalized PDF (%)\". I am not sure how the normalization is carried out and why it is\nconsidered necessary. It should be emphasized that improper normalization has the\npotential to distort the intended meaning of the probability distribution.\nIn the same figure, there appears to be some ambiguity regarding the connection between\nthe left wiggle (labeled \"PDF of CRF timing\") and the right wiggle (labeled \"CRF (PDF>68%)\").\nSpecifically, the height of the peak at ~5.5 seconds in the left wiggle seems lower than half\nof the peak’s height at ~6.5 seconds. Considering that probabilities should not exceed 100%,\nthe peak at 5.5 seconds should represent a probability lower than 50%. However, in the\nright wiggle, the 5.5-second peak appears as if it carries a significant probability (68%). This\ninconsistency requires clarification.\nFurthermore, the manuscript mentions the use of normalization for imaging at Line 386, but\nits necessity remains unclear to me. It seems that you can obtain the averaged probability in\nthe image domain by simply summing the probabilities and dividing by the number of\npiercing traces at each bin. Therefore, an explanation for why normalization is employed in\nthis context would be beneficial.\n2. Forward simulation methods\nThe manuscript extensively uses forward simulations to demonstrate the robustness of the\nresults but lacks information what specific method is used for the forward simulations,\nexcept for Figure 4, where the Thomson-Haskell matrix method is employed. Such\ninformation is essential to evaluate the validity of the analysis. It is particularly important to\nunderstand to what extent the simulations take into account 3-D heterogeneity and\nmultiple reflections. Some methods based on ray theory can easily overlook the effects of\nmultiple reflection phases.\nRelated to this, it is unclear why Thomson-Haskell’s method was chosen for Figure 4. This\nmethod is applicable to flat layered structures without dipping interfaces. Since you seem to\nuse the method that allows dipping interfaces (e.g., Figure S15-18), it raises the question of\nwhy the same method was not used for the analysis of Figure 4. If there is a specific reason\nfor using the Thomson-Haskell method, it would be helpful to elaborate on how you\ncorrected for any dipping effects.\n3. Iteration numbers\nI have some concerns regarding the relatively small number of iterations, which is set at\n10,000 in your analysis. This number is notably lower than what is typically used in Markov\nChain Monte Carlo (MCMC) analyses. It would be beneficial to justify the adequacy of this\niteration number, for example, by demonstrating the evolution of the likelihood over the\ncourse of iterations.\n4. Other minor comments\n4-1. Line 366: What is the size of the increment for the moving sub-array?\n4-2. Figure 4: I think that not all results from the grid-search is displayed on these Figures\n4(a, b, and c). How did you select the ones to be displayed?\n4-3. Fig. S11: The negative phase interpreted as slab-LAB shows a dip almost zero. Do you\nhave any justification on it?", "reviewer_2": "A review of\nSeismic evidence for a melt-rich Lithosphere-Asthenosphere 1 Boundary along the base of\nyoung slab at Cascadia by Wang et al\nThis paper uses receiver functions to provide images of the LAB beneath the Cascadia\nmargin. This a young subduction system so it is of interest to see if results from here\ncompared with what is found in Japan and New Zealand (NZ) where the subducted plate is\nmuch older.\nThe study presents a compelling case for an image of the LAB. The authors present an\nexhaustive test of different options for the RF signal in terms of potential multiples and\nconvincingly make their case. I think there are some 23 figures in the Supp material !\nThe synthetic modelling is especially good and relevant.\nThe observation that the depth to the LAB is similar to that seen in Japan and NZ is not what\nI would have expected given recent findings showing the LAB increasing its depth with age.\nHowever, we are still very much in the learning stage of the LAB and this is just an added\nnew observation to the mix.\nI was interested to see the authors argue boundary at the top of the LAB is sharp . There is a\nsentence that reads: Our inferred slab-LAB consistently exhibits an amplitude of 12±3% of\nthe direct P-arrival across the entire frequency range (Fig. 4), indicating a sharp discontinuity\nwith a dVs of 9.8±1.5% over a depth range of 1.5±1.5 km (Figs. 4 and S21).\nThe first part of sentence does not make sense ( a slab does not exhibit an amplitude !) and\nin general the whole sentence needs some unpacking. Why not say what the wavelength is\nand then apply the wavelength length rule as given in lines 242 and 243.\nOr maybe use the figure that Stern et al (2015) have in their supplementary material ( Fig.\nS7f) that shows the relationship of Reflection coefficient Vs normalised thickness ( i.e. true\nthickness of transition zone /wavelength). This is an important finding of your paper and\nneeds more explanation.\nHave the authors seen : 128. Herath, P., Stern, T. A., Savage, M. K., Bassett, D., and Henrys,\nS., 2022, Wide-angle seismic reflections reveal a lithosphere-asthenosphere boundary zone\nin the subducting Pacific Plate, New Zealand: Science Advances, v. 8, no. 38, p. eabn5697.\nIn this paper it is argued that the very top of the plate is a thin layer of very high P-wave\nspeed, due to azimuthal anisotropy, an immediately beneath it there is the low speed\nregion. Are your results consistent with this ? or not ?\nEditorial issues:\nLine 42: I would suggest you remove the word “very”. It is meaningless in a science context,\nand the sentence becomes stronger without that word. But your choice.\nLine 55 A note about the origin of the idea of a thin, low-velocity, melt layer at the base of\ntectonics plates facilitating plate motion and allowing plate tectonics to work. In your paper\nyou ascribe this idea to Debayle et al (2020). My sense was that the first group to make the\ndiscovery of this thin layer, and point out that it is this layer that facilitates plates to slide,\nwas the Stern et al (2015) paper. It would be appropriate to credit the first study that made\nthis important link.\nLine 82: another “very” that could be removed.\nLine 269 : I think the word “at” should be “of”.\nIn summary this is an important new finding and should be published. The paper is well\npresented . Figures are exceptional and it is well written.\nMy recommendation that minor edits, as listed above, should be made then it is ready to\ngo.", "reviewer_3": "Review Wang, Chen et al 2024\nThis manuscript presents a new high resolution image of the LAB of subducting plate\nbeneath Cascadia. Interpretation of the phase as the LAB seems sound given that its depth\nagrees with the high resolution tomography in the region. Similarly, much work has been\ndone to demonstrate the robustness of the interpretation (SI).\nThe study is novel in that it is high resolution, requiring a sharp discontinuity at the base of\nthe plate and shows that this is occurring down dip of locations where the LAB has by\nimaged by methods that do not resolve dipping structures as well. This is important to our\nunderstanding of both the LAB and subduction dynamics.\nThe authors interpretation of melt is likely sound. However, it would be nice to see a little\nmore logic to back it up.\nThere has been a lot of talk in the community that a new paper (Liu et al) shows hydration\ncan explain seismic observations and melt is not needed. However, in reality, overall, the\neffects of hydration on seismic waves are mixed depending on study – not just the Liu paper\n(see Cline, Faul, et l), and it remains to be seen which study is correct and why. In addition,\nthe Liu result does not explain much more than a percent or of velocity reduction at shorter\nperiods like those of P-to-S assuming mild hydration (150 ppm), and larger percentages of\nhydration likely cause melt at asthenospheric P/T conditions. That is another reason that\nthis study (using higher freq P-to-S) is important in the debate on the LAB.\nIn terms of whether melt affects the viscosity of the mantle and therefore defines the LAB:\nLab experiments diverge on effect of small amount of melt on viscosity (Hirth & Kohlstedt,\n1995; Mei & Kohlstedt, 2000)(Faul & Jackson, 2007). Theoretical models have argued that\nvery small fractions (< 0.01 %) would have very little effect (Rudge, 2018) or up to 2 orders\nof magnitude of reduction (Holtzman, 2016). However, at higher melt fractions such as\nthose required by this study (>1 %), there is likely a large effect on viscosity. Larger melt\nfractions have been reported by several observationally based papers, but there is still\ndisbelief by many in the community that such percentages could be stable in the mantle\nover time and length scales to be imaged seismically. This is another aspect of the result\nthat is novel – i.e., this is a big deal, something that needs more work. BTW, some 2 phase\nflow models (Katz, Joyce Sim, Turner, etc.) do find melt can pond beneath a permeability\nboundary in large percentages. So, again a hot topic of debate and something that needs to\nbe addressed.\nWhat is the interpretation of the negative intracrustal discontinuity of the upper plate?\nI suggest breaking this into 2 sentences and explaining a little more.\n“The reduction in the amplitude of the Ps conversion from the slab top 188 at ~40 km\ndepth, as well as the absence of continental Moho near the mantle wedge (Fig. 3), is 189\nconsistent with the diminished velocity contrast between the hydrated forearc mantle\nwedge and 190 the subducting or overlying continental crust (Bostock et al., 2002;\nKawakatsu & Watada, 2007).”\nAlso, the lack of a slab Moho sometimes interpreted as eclotigization of basaltic crust at\nleast in the Bostock work.\n“The denser spatial coverage of the short-period array allows us to better constrain shallow\nstructures which aids in constraining deep structures.” What shallow structures? Can you be\nmore specific? Also, how so? I thought that you were using an independent velocity model\nto translate to depth. How does information about shallower structures improve deeper\nstuctures?\n“Recent geophysical studies suggest its persistence with down-going old and cold slabs, but\nit is not known whether it is commonly present, especially for young and warm slabs such as\nin the Cascadia subduction zone.” This is not true. Several authors infer melt beneath the\nslab in Cascadia (e.g Hawley). Please rewrite for accuracy.\n“Nonetheless, 277 neither of these tomographic studies has adequate spatial resolution to\ndetect the melt-rich LAB imaged in our study.” This is vague. Are you sure seismic\ntomography would not detect it, even if blurred by resolution? Please be more specific.\n“We speculate that the melts at both the plate- and slab-LAB have a similar local origin,\nwithout 280 long-distance lateral migration.” This statement and paragraph need more\nlogic. How/why does the observation indicates the melt has not migrated? In general,\nseismology provides a snapshot of the mantle and has difficulty deciphering the past/any\nkind of time differential. Maybe consider results from geodynamic modeling?\nLine 165: Fig S11 does not look comparable to Fig. S3, instead the LAB seems to be flat and\napproaching the Moho toward the right of the figure. Please explain. Is this predicted for\nraypaths from NW?"}, "author_response": "Authors’ Response to the Reviewers’ Comments\nContents:\nResponse to Comments from Reviewer #1 ............................................................. 1\nResponse to Comments from Reviewer #2 ............................................................. 7\nResponse to Comments from Reviewer #3 ........................................................... 10\nNotes:\n Reviewers’ comments are shown in black, and our response are shown in blue.\n The line numbers in this letter correspond to Slab-LAB-maintext-with-changes-tracked.pdf.\n Fig. x - - - - - - - - - - - - - - Figures in the revised main text.\n Fig. Sx - - - - - - - - - - - - - - Figures in the revised supplementary material.\n Fig. Rx - - - - - - - - - - - - - - Figures in the response letter.\nResponse to Comments from Reviewer #1\nGeneral comments:\nThis manuscript by Wang et al. identifies a coherent negative phase in the receiver functions\nobtained from dense linear array of seismometers deployed at the Cascadia subduction zone.\nThe authors attribute the negative phase to the lithosphere-asthenosphere boundary (LAB) at the\nbased on the subducting slab. Due to its sharpness, the LAB is interpreted to be at the melt-rich\ncondition. The manuscript is overall well-written, and its novelty is clear. However, the method\n1\nhas some points needed to be further clarified. The requested revisions have the potential to\nimpact the results, and thus, I recommend categorizing these changes as 'major revisions'.\nResponse: We appreciate the reviewer’s great efforts to review the paper and provide valuable\nsuggestions. Please find detailed point-by-point replies to all the comments below.\nComment 1: Probability distribution\nThe procedure for calculating the probability distribution needs to be described with greater clarity.\nFrom my perspective, the process involves two steps: firstly, constructing a 1-D\nprobability distribution for pulse presence by marginalizing the posterior probability across time for each sub-\narray. Subsequently, this 1-D time series of probabilities is either stacked or averaged within each\nbin on the cross-section using the CCP stacking method. It is important to note that, during these\nprocesses, there is no requirement for 'normalizing' the probability.\nHowever, the manuscript repeatedly mentions normalization at both stages of analysis. The initial\nreference to normalization is found in the color scale of Fig S4(b), where it states \"Normalized\nPDF (%)\". I am not sure how the normalization is carried out and why it is considered necessary.\nIt should be emphasized that improper normalization has the potential to distort the intended\nmeaning of the probability distribution.\nIn the same figure, there appears to be some ambiguity regarding the connection between the\nleft wiggle (labeled \"PDF of CRF timing\") and the right wiggle (labeled \"CRF (PDF>68%)\").\nSpecifically, the height of the peak at ~5.5 seconds in the left wiggle seems lower than half of the\npeak’s height at ~6.5 seconds. Considering that probabilities should not exceed 100%, the peak\nat 5.5 seconds should represent a probability lower than 50%. However, in the right wiggle, the\n5.5-second peak appears as if it carries a significant probability (68%). This inconsistency requires\nclarification.\nFurthermore, the manuscript mentions the use of normalization for imaging at Line 386, but its\nnecessity remains unclear to me. It seems that you can obtain the averaged probability in the\nimage domain by simply summing the probabilities and dividing by the number of piercing traces\nat each bin. Therefore, an explanation for why normalization is employed in this context would be\nbeneficial.\nResponse: Thank you for your comments.\nsolutions in terms of their density distribution as a function of timing and slowness of the Coherent\nReceiver Functions (CRF). In our Markov chain Monte Carlo (McMC) inversion, the ensemble\nsolutions are initially represented by the counts of the accepted models (the CRFs). For\nvisualization, we convert the counts into a density distribution, scaling the maximum value to 1,\nand we use the density plot to represent the Posterior Probability Distribution (PDF) of the CRFs.\nThis approach was chosen to standardize the visual comparison across different subarrays. In\nthe revised version, we have removed the term \"normalization\" and updated the labels to avoid\nany confusion, acknowledging that the use of the density plot to represent PDF is a common\npractice in McMC studies. Please see revised Fig. S4b.\nRegarding the inconsistency between the left and right wiggles in Fig. S4b, we appreciate you\nbringing this to our attention. The right wiggle plot displays the 1D CRF with a 68% confidence\nlevel. To present the CRF in a manner similar to the conventional Receiver Function (RF) which\ncontains discrete delta functions with different timings and amplitudes, we select the peaks with\nthe probability over a certain threshold (e.g., larger than one or two standard deviations of the\nmean) in the 2D density map to construct the 1D CRF. We now realize that using \"normalized\nprobability\" with \"confidence level” could cause confusion, since both are expressed as\npercentage. To avoid confusion, we have revised Fig. S4b and its caption.\nComment 2: Forward simulation methods\nThe manuscript extensively uses forward simulations to demonstrate the robustness of the results\nbut lacks information what specific method is used for the forward simulations, except for Figure\n4, where the Thomson-Haskell matrix method is employed. Such information is essential to\nevaluate the validity of the analysis. It is particularly important to understand to what extent the\nsimulations take into account 3-D heterogeneity and multiple reflections. Some methods based\non ray theory can easily overlook the effects of multiple reflection phases.\nRelated to this, it is unclear why Thomson-Haskell’s method was chosen for Figure 4. This method\nis applicable to flat layered structures without dipping interfaces. Since you seem to use the\nmethod that allows dipping interfaces (e.g., Figure S15-18), it raises the question of why the same\nmethod was not used for the analysis of Figure 4. If there is a specific reason for using the\nThomson-Haskell method, it would be helpful to elaborate on how you corrected for any dipping\neffects.\nFor the synthetic studies shown in Figures S2 and S15-S18, we used the Raysum code\n(Frederiksen & Bostock, 2000), which is a ray-based method to calculate synthetic seismograms\nfor dipping discontinuities and anisotropic media. For the 2D synthetic studies shown in S7 and\nS12-S14, we calculated the synthetic seismograms with the Thomson-Haskell matrix method\nusing a 1D layered velocity model extracted from the 2D velocity model for each station. These\nsynthetics are spatially arranged to mimic the effect of the 2D slab. Given the purpose of testing\nthe effects of shallow multiples, not directly calculating 2D synthetic seismograms reduces\ncomputational cost. Our previous studies have demonstrated the reliability of this approach (Wang\net al., 2021; Zhong & Zhan, 2020). These details have been added in the revised manuscript\n(lines 434-439 and lines 461-464).\nIn the multi-frequency RF waveform modeling section, we use the Thomson-Haskell’s method for\nsynthetic calculation. This method is chosen mainly for its computational efficiency, particularly\nwhen a larger number of layers is required to represent gradual transitions in velocity models.\nAlthough the Raysum code is an alternative with advantages in dealing with dipping layers, its\napplication to a larger number of layers gives rise to extremely long run times and leads to\nsegmentation faults (typically got segmentation fault with more than 8 layers, see Fig. 6 in Bloch\nand Audet (2023)). In using the Thomon-Hakell’s method, we must account for the dipping effect.\nTo address this, we use the slab’s average dip ( ) and strike ( ) from the Slab2.0 model\n(Hayes et al., 2018) and estimate the amplitude correction factor (~ based on synthetic\n12° 356°\ntests (Fig. R1). We have further clarified these points in the revised manuscript (lines 485-487).\n1.3±0.1)\nRF study across different back azimuth. The lower panel shows the variations in Ps amplitude\nratio as a function of back-azimuth for a synthetic model with a dipping LAB interface, in\ncomparison to a refence model with a flat LAB interface.\nComment 3: Iteration numbers\nI have some concerns regarding the relatively small number of iterations, which is set at\n10,000 in your analysis. This number is notably lower than what is typically used in Markov\nChain Monte Carlo (MCMC) analyses. It would be beneficial to justify the adequacy of this\niteration number, for example, by demonstrating the evolution of the likelihood over the course of\niterations.\nResponse: Our Markov Chain Monte Carlo (MCMC) inversion was done using 500 CPUs, with\neach CPU running an independent chain. For each chain, 10,000 iterations were performed.\nConsequently, the total number of the models in the MCMC inversion are five million. Fig. R2a\nshows the evolution of the misfit (normalized to the “best” model) with iterations for the 500 parallel\nchains, and Fig. R2b shows the distribution of misfit for all ensemble solutions. From Fig. R2, we\ncan see that the inversion typically convergences within only 5,000 iterations. We have updated\nFig. S4 in the manuscript to show the evolution of the misfit over iterations.\nrunning in parallel. Each chain contains 10,000 McMC iterations, with the second half of the\niterations used to generate ensemble solutions. The histogram shown in the right panel shows\nthe distribution of relative misfits for all the ensemble solutions.\nComment 4: Other minor comments\n4-1. Line 366: What is the size of the increment for the moving sub-array?\nResponse: We used a moving sub-array of ~15 km with an incremental step for each station.\nSee lines 382-383 in the revised manuscript.\n4-2. Figure 4: I think that not all results from the grid-search is displayed on these Figures 4(a, b,\nand c). How did you select the ones to be displayed?\nResponse: We only displayed those models with relative misfit smaller than 100%. We have\nprovided further clarification in the revised manuscript (lines 262-264 and lines 488-491).\n4-3. Fig. S11: The negative phase interpreted as slab-LAB shows a dip almost zero. Do you have\nany justification on it?\nResponse: The slab-LAB imaged using waveforms from the northwest appears flat, and the\nthickness of the slab seems to vary along the down-dip direction. We have conducted justification\non this feature by considering major factors that may affect it. The Bayesian posterior probability\n6\ndistribution of the CRFs shows a high confidence level for this discontinuity (Fig. S11b),\nsuggesting that noise effects may not be significant on the image. The observed slab-LAB roughly\nresembles the shape of the high-velocity slab shown in the tomographic result (Fig. S11c),\nalthough it manifests at a depth ~10 km deeper at distances of ~0-40 km along the profile. This\ndiscrepancy is unlikely to result solely from the velocity model used in the time-to-depth\nconversion, as recognizing a 10 km decrease in the slab-LAB depth without changing the slab-\nMoho depth would require a slab mantle velocity ~30% lower than that of the current velocity\nmodel, which is unreasonable given the absolute velocities used (Fig. S11c). Further studies are\nrequired to understand the 3D slab effects in RF imaging and to incorporate more accurate\ntomographic models in the time-to-depth migration. If the flat appearance of the imaged slab-LAB\nshown in Fig. S11 reflects a real structural feature, it could suggest the presence of significant\nlateral variations of the slab structure. Indeed, both the strong lateral variation of slab geometry\nshown in the Slab 2.0 model (black dashed lines shown in Fig. S1) and strong variations of the\nplate-LAB depth observed offshore of Cascadia (Fig. S23 in the supplementary material) suggest\na complex slab structure in this region. We have added some words about our analysis on the\nshape of the slab-LAB in the caption of Fig. S11.\nResponse to Comments from Reviewer #2\nGeneral comments:\nThis paper uses receiver functions to provide images of the LAB beneath the Cascadia margin.\nThis a young subduction system so it is of interest to see if results from here compared with what\nis found in Japan and New Zealand (NZ) where the subducted plate is much older. The study\npresents a compelling case for an image of the LAB. The authors present an exhaustive test of\ndifferent options for the RF signal in terms of potential multiples and convincingly make their case.\nI think there are some 23 figures in the Supp material!\nThe synthetic modelling is especially good and relevant.\nThe observation that the depth to the LAB is similar to that seen in Japan and NZ is not what I\nwould have expected given recent findings showing the LAB increasing its depth with age.\nHowever, we are still very much in the learning stage of the LAB and this is just an added new\nobservation to the mix.\nsentence that reads: Our inferred slab-LAB consistently exhibits an amplitude of 12±3% of the\ndirect P-arrival across the entire frequency range (Fig. 4), indicating a sharp discontinuity with a\ndVs of 9.8±15% over a depth range of 1.5±1.5 km (Figs. 4 and S21).\nThe first part of sentence does not make sense (a slab does not exhibit an amplitude !) and in\ngeneral the whole sentence needs some unpacking. Why not say what the wavelength is\nand then apply the wavelength length rule as given in lines 242 and 243.\nOr maybe use the figure that Stern et al (2015) have in their supplementary material (Fig. S7f)\nthat shows the relationship of Reflection coefficient Vs normalised thickness (i.e. true thickness\nof transition zone /wavelength). This is an important finding of your paper and needs more\nexplanation.\nResponse: Thank you very much for your encouraging comments and suggestions. In the revised\nversion, we have revised our wording regarding the “slab-LAB exhibits an amplitude”. In addition,\nwe have incorporated the wavelength rule, as suggested, to estimate the upper bound thickness\nof this discontinuity. Please refer to lines 242-247 in the revised manuscript for these updates. By\nthe way, our results indeed show that the LAB at Cascadia is shallower than those in Japan and\nNew Zealand, consistent with recent findings showing the LAB increasing its depth with age.\nHave the authors seen : 128. Herath, P., Stern, T. A., Savage, M. K., Bassett, D., and Henrys, S.,\n2022, Wide-angle seismic reflections reveal a lithosphere-asthenosphere boundary zone in the\nsubducting Pacific Plate, New Zealand: Science Advances, v. 8, no. 38, p. eabn5697.\nIn this paper it is argued that the very top of the plate is a thin layer of very high P-wave\nspeed, due to azimuthal anisotropy, an immediately beneath it there is the low speed region.\nAre your results consistent with this ? or not ?\nResponse: Thank you for the reference.\nHerath et al. (2022) used active seismic reflection data to investigate the LAB structure in the\nHikurangi subduction zone and identified a series of seismic reflections at depths associated with\nthe LAB. They interpreted these reflections as the presence of a LAB zone consisting of a thin,\nhigh-velocity anisotropic layer atop a low-velocity, melt-rich layer.\nIn our study, we use passive seismic data and the RF technique which can resolve the overall\nvariation of seismic velocity with depth, i.e., a single negative discontinuity corresponding to the\nseismic data, our study does not have the ability to resolve the detailed structure of the LAB zone\nas described by Herath et al.\nSpecific comments:\nEditorial issues:\nLine 42: I would suggest you remove the word “very”. It is meaningless in a science context, and\nthe sentence becomes stronger without that word. But your choice.\nResponse: Here, we feel that the word “very” is a useful part of the scientific discussion because\nit emphasizes that the observed 3 km thickness is unusually small.\nLine 55 A note about the origin of the idea of a thin, low-velocity, melt layer at the base of tectonics\nplates facilitating plate motion and allowing plate tectonics to work. In your paper you ascribe this\nidea to Debayle et al (2020). My sense was that the first group to make the discovery of this thin\nlayer, and point out that it is this layer that facilitates plates to slide, was the Stern et al (2015)\npaper. It would be appropriate to credit the first study that made this important link.\nResponse: Thanks for pointing out this problem. The reference has been added.\nLine 82: another “very” that could be removed.\nResponse: The word “very” has been removed.\nLine 269 : I think the word “at” should be “of”.\nResponse: We have rephrased the sentence (line 272).\nIn summary this is an important new finding and should be published. The paper is well presented .\nFigures are exceptional and it is well written.\nMy recommendation that minor edits, as listed above, should be made then it is ready to go.\nResponse: Thank you very much for your encouraging comments.\nGeneral comments:\nThis manuscript presents a new high resolution image of the LAB of subducting plate beneath\nCascadia. Interpretation of the phase as the LAB seems sound given that its depth agrees with\nthe high resolution tomography in the region. Similarly, much work has been done to demonstrate\nthe robustness of the interpretation (SI).\nThe study is novel in that it is high resolution, requiring a sharp discontinuity at the base of the\nplate and shows that this is occurring down dip of locations where the LAB has by imaged by\nmethods that do not resolve dipping structures as well. This is important to our understanding of\nboth the LAB and subduction dynamics.\nThe authors interpretation of melt is likely sound. However, it would be nice to see a little more\nlogic to back it up.\nThere has been a lot of talk in the community that a new paper (Liu et al) shows hydration\ncan explain seismic observations and melt is not needed. However, in reality, overall, the\neffects of hydration on seismic waves are mixed depending on study – not just the Liu paper\n(see Cline, Faul, et l), and it remains to be seen which study is correct and why. In addition,\nthe Liu result does not explain much more than a percent or of velocity reduction at shorter\nperiods like those of P-to-S assuming mild hydration (150 ppm), and larger percentages of\nhydration likely cause melt at asthenospheric P/T conditions. That is another reason that this study\n(using higher freq P-to-S) is important in the debate on the LAB.\nIn terms of whether melt affects the viscosity of the mantle and therefore defines the LAB:\nLab experiments diverge on effect of small amount of melt on viscosity (Hirth & Kohlstedt, 1995;\nMei & Kohlstedt, 2000)(Faul & Jackson, 2007). Theoretical models have argued that very small\nfractions (< 0.01 %) would have very little effect (Rudge, 2018) or up to 2 orders of magnitude of\nreduction (Holtzman, 2016). However, at higher melt fractions such as those required by this study\n(>1 %), there is likely a large effect on viscosity. Larger melt fractions have been reported by\nseveral observationally based papers, but there is still disbelief by many in the community that\nsuch percentages could be stable in the mantle over time and length scales to be imaged\nseismically. This is another aspect of the result that is novel – ie., this is a big deal, something\nthat needs more work. BTW, some 2 phase flow models (Katz, Joyce Sim, Turner, etc.) do find\n10\nmelt can pond beneath a permeability boundary in large percentages. So, again a hot topic of\ndebate and something that needs to be addressed.\nResponse: Thanks for your insightful comments and constructive suggestions.\nWe agree with your assessment of the findings of Liu et al. (2023). Considering the specific\nfrequencies (up to 2 Hz) utilized in our study, hydration alone may not explain the substantial\nvelocity reductions observed at the slab-LAB. In the revised manuscript, we have added this\nperspective. Please see lines 250-255 in the revised manuscript.\nRegarding the questions about how melt affects the viscosity and the stability of melts over\ngeologically timescales, we agree with the reviewer and have added a few sentences in the\nrevised manuscript. See lines 296-298 and lines 316-320 in the revised manuscript.\nSpecific comments:\nComment 1:\nWhat is the interpretation of the negative intracrustal discontinuity of the upper plate?\nResponse: The negative intracrustal discontinuity may be associated with the accumulation of\ncrustal fluids, resulting from the fluids released by subducting sediments or oceanic crust, or from\nthe heating and dewatering of the mantle wedge, as observed by previous studies in Cascadia\n(Egbert et al., 2022; Wang et al., 2019; Ward et al., 2018) and other subduction zones (Worzewski\net al., 2011). We have added a sentence (lines 169-172) to clarify this.\nComment 2:\nI suggest breaking this into 2 sentences and explaining a little more.\n“The reduction in the amplitude of the Ps conversion from the slab top at ~40 km depth, as well\nas the absence of continental Moho near the mantle wedge (Fig. 3), is consistent with the\ndiminished velocity contrast between the hydrated forearc mantle wedge and the subducting or\noverlying continental crust (Bostock et al., 2002; Kawakatsu & Watada, 2007).” Also, the lack of\na slab Moho sometimes interpreted as eclotigization of basaltic crust at least in the Bostock work.\nResponse: Thank you for your suggestions. We agree that the original sentence does not make\nit clear that the mechanisms for the absence of the slab top at ~40 km depth and the absence of\ncontinental Moho are different. The absence of the slab top at ~40 km depth may primarily be\nattributed to the eclogitization of the subducting oceanic crust, as discussed by Bostock et al.\n11\n(2002)and Rondenay et al. (2008). With the progressive eclogitization of the oceanic crust, the\nslab Moho also becomes invisible as anhydrous eclogite cannot be distinguished from the\nunderlying slab mantle on the basis of seismic wave speeds (Bostock et al., 2002; Chen et al.,\n2005). We have refined our sentences as follows (Lines 186-190): “The decrease in the amplitude\nof the Ps conversions at the slab top and slab-Moho when slab depth exceeds ~40 km (Fig. 3) is\nlikely due to the eclogitization process of the subducting oceanic crust (Bostock et al., 2002;\nRondenay et al., 2008). Eclogitization also releases aqueous fluid into the overlying forearc\nmantle wedge, leading to the well-known absence of a detectable continental Moho near the\nmantle wedge corner.”.\nComment 2:\n“The denser spatial coverage of the short-period array allows us to better constrain shallow\nstructures which aids in constraining deep structures.” What shallow structures? Can you be\nmore specific? Also, how so? I thought that you were using an independent velocity model\nto translate to depth. How does information about shallower structures improve deeper\nstuctures?\nResponse: We appreciate the opportunity to provide further clarification on this matter. It's\nimportant to note that the procedures described does not involve time-to-depth conversion, but\nrather relates to the CRF method used in our study. Specifically, we employ the McMC technique\nto iteratively add, remove, or adjust RF phases, until the likelihood function (or misfit function)\nstabilized. Compared to the traditional time-domain iterative deconvolution, in which RF phases\nare fixed once added, McMC allows for dynamic update of all phases throughout the process.\nThe resultant more accurate constraint on earlier CRF phases (shallow structures) helps to\nimprove the quality of the later phases (deeper structures). The denser spatial coverage of the\nshort-period array significantly improves our ability to constrain the shallow structures, such as\nbasin structures. This will directly contribute to the deeper structures imaging (later CRF phases)\nin the CRF inversion. In the revised version, we have revised our statements to make it clearer\n(line 352 and lines 389-394).\nComment 3:\n“Recent geophysical studies suggest its persistence with down-going old and cold slabs, but it is\nnot known whether it is commonly present, especially for young and warm slabs such as in the\nCascadia (e.g Hawley). Please rewrite for accuracy.\nResponse: Hawley et al. (2016) reported a subslab low-velocity zone (LVZ) through tomography\nstudy along the entire Cascadia forearc and proposed that the LVZ reflects the accumulation of\npartial melts that were present as a thin layer beneath the slab. But in a later study by Bodmer et\nal. (2018), their tomography results show that the LVZs only exist in northern and southern\nCascadia and proposed that they reflect partial melting due to localized mantle upwelling as\nopposed to a thin layer of partial melts below the slab. Thus, whether there is a thin layer of partial\nmelts below the slab in Cascadia remains debated. In light of this, we have revised our statement\n(line 35) to more accurately reflect the current state of knowledge of the sub-slab partial melts\nbelow slab in Cascadia.\nComment 4:\n“Nonetheless, neither of these tomographic studies has adequate spatial resolution to detect the\nmelt-rich LAB imaged in our study.” This is vague. Are you sure seismic tomography would not\ndetect it, even if blurred by resolution? Please be more specific.\nResponse: In our revised manicurist, we emphasize that these tomography models do not have\nenough vertical resolution to detect the sharp slab-LAB seismic discontinuity discussed in our\nstudy. We acknowledge that tomographic imaging may also be sensitive to the LVZs caused by\npartial melts, but its low vertical resolution prevents it from distinguishing whether the seismic\nvelocity sharply or gradually changes with depth. In the revised manuscript, we have revised our\nstatement to make it clearer (lines 287-289).\nComment 5:\n“We speculate that the melts at both the plate- and slab-LAB have a similar local origin, without\nlong-distance lateral migration.” This statement and paragraph need more logic. How/why does\nthe observation indicates the melt has not migrated? In general, seismology provides a snapshot\nof the mantle and has difficulty deciphering the past/any kind of time differential. Maybe consider\nresults from geodynamic modeling?\nResponse: Thanks for the comments. We have refined this sentence to make it clearer:\n“Considering the inference of melt-rich LAB with comparable melt fractions beneath oceanic\n13\nplates both before and after subduction (Rychert et al., 2020; Kawakatsu et al., 2009; Stern et\nal., 2015; this study), we speculate that the melts originate predominantly locally, without long-\ndistance lateral migration from ridges or plumes.”. Please see lines 290-293 in the revised\nmanuscript. We acknowledge this is an important question require further geodynamic studies.\nComment 6:\nLine 165: Fig S11 does not look comparable to Fig. S3, instead the LAB seems to be flat and\napproaching the Moho toward the right of the figure. Please explain. Is this predicted for\nraypaths from NW?\nResponse: We think the reviewer meant Fig. 3 instead of Fig. S3 in this comment. A similar\nquestion has been raised by Reviewer #1. Please see our response to Reviewer #1’s Minor\nComment 4-3 above.\nReference:\nBloch, W., & Audet, P. (2023). PyRaysum: Software for Modeling Ray-theoretical Plane Body-\nwave Propagation in Dipping Anisotropic Media. Seismica, 2(1).\nBodmer, M., Toomey, D. R., Hooft, E. E. E., & Schmandt, B. (2018). Buoyant Asthenosphere\nBeneath Cascadia Influences Megathrust Segmentation. Geophysical Research\nLetters, 45(14), 6954–6962. https://doi.org/10.1029/2018GL078700\nBostock, M. G., Hyndman, R. D., Rondenay, S., & Peacock, S. M. (2002). An inverted\ncontinental Moho and serpentinization of the forearc mantle. Nature, 417(6888), 536–\n538. https://doi.org/10.1038/417536a\nChen, L., Wen, L., & Zheng, T. (2005). A wave equation migration method for receiver\nfunction imaging: 2. Application to the Japan subduction zone. Journal of Geophysical\nResearch: Solid Earth, 110(B11).\nEgbert, G. D., Yang, B., Bedrosian, P. A., Key, K., Livelybrooks, D. W., Schultz, A., et al.\n(2022). Fluid transport and storage in the Cascadia forearc influenced by overriding\nplate lithology. Nature Geoscience, 15(8), 677–682. https://doi.org/10.1038s41561-\n022-00981-8\nFrederiksen, A. W., & Bostock, M. G. (2000). Modelling teleseismic waves in dipping\nanisotropic structures. Geophysical Journal International, 141(2), 401–412.\nhttps://doi.org/10.1046j.1365-246x.2000.00090.x\nHawley, W. B., Allen, R. M., & Richards, M. A. (2016). Tomography reveals buoyant\nasthenosphere accumulating beneath the Juan de Fuca plate. Science, 353(6306),\n1406. https://doi.org/10.1126science.aad8104\nHayes, G. P., Moore, G. L., Portner, D. E., Hearne, M., Flamme, H., Furtney, M., & Smoczyk,\nG. M. (2018). Slab2, a comprehensive subduction zone geometry model. Science,\n362(6410), 58–61.\nreflections reveal a lithosphere-asthenosphere boundary zone in the subducting\nPacific Plate, New Zealand. Science Advances, 8(38), eabn5697.\nhttps://doi.org/10.1126sciadv.abn5697\nLiu, C., Yoshino, T., Yamazaki, D., Tsujino, N., Gomi, H., Sakurai, M., et al. (2023). Effect of\nwater on seismic attenuation of the upper mantle: The origin of the sharp lithosphere–\nasthenosphere boundary. Proceedings of the National Academy of Sciences, 120(32),\ne2221770120.\nRondenay, S., Abers, G. A., & van Keken, P. E. (2008). Seismic imaging of subduction zone\nmetamorphism. Geology, 36(4), 275–278. https://doi.org/10.1130G24112A.1\nWang, X., Zhan, Z., Zhong, M., Persaud, P., & Clayton, R. W. (2021). Urban Basin Structure\nImaging Based on Dense Arrays and Bayesian Array-Based Coherent Receiver\nFunctions. Journal of Geophysical Research: Solid Earth, 126(9), e2021JB022279.\nWang, Y., Lin, F.-C., & Ward, K. M. (2019). Ambient noise tomography across the Cascadia\nsubduction zone using dense linear seismic arrays and double beamforming.\nGeophysical Journal International, 217(3), 1668–1680.\nhttps://doi.org/10.1093gji/ggz109\nWard, K., Lin, F., & Schmandt, B. (2018). High Resolution Receiver Function Imaging Across\nthe Cascadia Subduction Zone Using a Dense Nodal Array. Geophysical Research\nLetters, 45(22), 12–218. ‐\nWorzewski, T., Jegen, M., Kopp, H., Brasse, H., & Taylor Castillo, W. (2011). Magnetotelluric\nimage of the fluid cycle in the Costa Rican subduction zone. Nature Geoscience, 4(2),\n108–111. https://doi.org/10.1038ngeo1041\nZhong, M., & Zhan, Z. (2020). An array-based receiver function deconvolution method:\nmethodology and application. Geophysical Journal International, 222(1), 1–14.\n15"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciate the authors' efforts in enhancing the quality of the manuscript. My\nrequirements have been satisfactorily addressed, and I am confident that the manuscript is\nnow prepared for publication.", "reviewer_3": "Thank you for responding to the requested edits. This is a nice manuscript, of ide interest to\nmany, worthy of publication."}, "author_response": ""}]

Decision letter