Dual-layered mantle lithosphere beneath southeastern Canadian Cordillera

Domain: earth_science
Published: 2025-11-25
Editorial outcome: major_revision

Abstract
The craton edge beneath southwestern Canada has been characterized as a ‘lithospheric step’ involving the transition from the warm Cordillera to an ancient, mechanically strong North American craton. Using converted shear waves that are sensitive to this boundary zone, this study shows a snapshot of this transition and reveals the existence of westward-dipping interfaces west of the Rocky Mountain Trench (RMT). These interfaces correlate well, both in location and conversion amplitude, with horizontal shear-velocity gradient zones within the craton core. In a span of 400-600 km west of the RMT, the presence of two distinctive horizontal structural gradients (at ~75 km and ~180 km) suggests a dual-lithosphere architecture, where the Precambrian craton core (Laurentia) now underlies the southeastern Cordilleran lithosphere. The morphology of the craton edge suggests influences from uplift and convective erosion over the past 60 million years.

Introduction
Marked by the rugged terrain of the Canadian Rocky Mountains, the mantle lithosphere west of the Cordillera-craton Boundary (CCB) is characterized as thin (~60 km) and uniform, largely formed in response to edge-driven, small-scale mantle convection in a back-arc setting1–3. Based on recent evidence from xenoliths4 and broadband seismic observations1,5, the anomalously thin Cordilleran lithosphere terminates abruptly near the Rocky Mountain Trench6,7 (RMT), where a young upper crustal fault appears to have developed between these two geologically distinct land masses. The proposed CCB is abutted to the east and northeast by an ancient Laurentian lithospheric root, which is 100–200 km thicker6,8–10. The resulting lithospheric step near the CCB underscores sharp changes in subcrustal temperature2,3, morphology6,11, rheology12,13 and composition14,15.

The origin of the present-day craton boundary structure remains debated. Before 2015, it was widely believed that the evolution of the Cordillera was initiated by the break-up and drifting of the Neoproterozoic supercontinent ‘Rodinia’16, followed by episodes of tectonic convergence and landward-dipping subduction. The CCB began to form in the late Devonian16 and continued into the Mesozoic as the Panthalassic Ocean subducted beneath the western margin of North America16,17. For decades, landward subduction was widely regarded as a defining feature of the “accretionary” (non-collisional) model of Cordilleran orogenesis17. This geometry was supported by early electromagnetic and active-source seismic data18, but it was incompatible with a number of subsequent geological and geophysical observations19,20. More recently, results from shear wave tomography6,21 and dynamic calculation22 favored a steeper craton edge and emphasized contributions from edge-driven convection22 or lithosphere detachment21.

The proliferation of regional network data since 2015 (e.g., refs. 23,24) has enabled an appraisal of earlier proposed models of Cordillera-craton interaction (e.g., refs. 16,19,20,25,26). Tomographic imaging based on improved data6,27,28 provides support for an alternative hypothesis centered on a late Cretaceous collision between an exotic Cordilleran ribbon continent and the North American (NA) craton, favoring west-dipping subduction beneath the Cordillera20. The presence of a cryptic collisional suture near the RMT is corroborated by geological observations, e.g., the existence of a carbonatite belt29, and recent geodynamical calculations11,14. Details of westward subduction model (e.g., subduction boundary location and geometry) are crucial for 1) understanding the paleogeographic evolution that culminated in orogenesis; and 2) determining whether orogenesis involved the accretion of exotic lithosphere, simple deformation of existing lithosphere, or some combination of the two.

This study examines potential shallow mantle interfaces beneath the Cordillera-craton transition zone by utilizing S-to-p converted waves resulting from impedance contrasts. We find compelling evidence for an enigmatic suture that is attributable to west-dipping subduction in the upper-mantle. According to our observations, the eastern Cordilleran lithosphere is underlain by the cratonic lithosphere, forming a complex, multi-layered upper-mantle architecture. The geometry of the Moho and the Lithosphere-Asthenosphere Boundary (LAB) across the Foreland margin (the deformation front) of an orogen provides direct constraints on the orogen-autochthon relationship. Commonly, continuity of the Moho and LAB from beneath an orogen out into the autochthon is interpreted as evidence of an autochthonous origin. In contrast, discontinuities affecting either the Moho or the LAB, including dual LABs and stacked Mohos, are commonly associated with 1) suturing of previously disparate lithospheric plates, or 2) significant post-orogenic modification of the lithospheric mantle30,31.

For simplicity, we refer to the lower boundaries of the Cordilleran and NA cratonic lithospheres as the Cordillera Lithosphere-Asthenosphere Boundary (CLAB) and Laurentia Lithosphere-Asthenosphere Boundary (LLAB), respectively.

Results
Our dataset comprises recordings from 107 broadband stations, including 40 from the USArray32 (Fig. 1), which substantially improves the data coverages and volumes examined by earlier regional studies of mantle converted waves (e.g., refs. 1,33). After component rotation (to great-circle path), filtering, and spectral deconvolution, we retained 3968 high-quality SRFs in the frequency range of 0.02–0.2 (or 0.1) Hz and depth-migrated the time axis using PREM34. The predefined discontinuity at 220 km depth, a feature of PREM34, can produce absolute depth errors no greater than 5 km after migrating the S-to-p converted phases (see Supplementary Material “Synthetics”).

The majority of receiver functions after undergoing stacking and bootstrapping analysis (see Supplementary Fig. S6) show robust negative phases in the depth range of 50–350 km, which are typically associated with earlier-reported lithospheric interfaces35–37. Spatially interpolated depths of major negative gradients (phase 1, Fig. 2) approximately follow a Gaussian distribution centered at an average of 207±3 [2σ] km. Approximately 75% of the observations fall within one sigma, whereas outliers near 1.98–1.90 Ga-aged Kiskatinaw and Ksituan domains in northwestern Alberta (block A, Fig. 2) are deeper than 230 km. The azimuthal coverage at our stations is biased toward NE–SW orientations. Despite detecting consistent phases associated with the LLAB and CLAB on individual SRFs (see Supplementary “Resolution”), the limited range of angles prohibits a systematic analysis of the azimuth-dependency of mantle conversions.

From the Buffalo Head Terrane (BHT) to the southern Hearne province, the strongest lithospheric gradient (see Fig. 2a, blocks B, C, and D) progressively deepens by more than 20 km towards the southeast where the cratonic basement is covered by 3 km thick Phanerozoic sedimentary strata38. This conversion occurs at ~180 km depth beneath the center of the BHT, which indicates a locally thinned lithosphere consistent with estimated depths from kimberlite xenolith thermobarometry39,40, electrical resistivity41, and heat flow42 data. The most important outcome of our SRF analysis is the identification of the LLAB beneath eastern Cordillera at ~180 km depth (see Fig. 2a, block E). This deep, sub-Cordilleran LLAB is quantitatively different from the CLAB, which we also identify through a series of west-dipping phases (at ~75 km), and it is not directly associated with the reported ‘lithospheric-step’ (e.g., ref. 14). The most compelling evidence of sub-Cordilleran interfaces and their connections to the adjacent NA craton is regional stacks of SRFs on either side of the RMT. On the craton side, only a single interface at ~220 km, i.e., the LLAB, is detected within the broad depth range of 50–250 km. On the other hand, the upper mantle beneath eastern Cordillera contains distinct interfaces at 180 and 75 km depths. Despite a reduced amplitude and considerable waveform complexities, the CLAB signature is recognizable above ambient noise levels.

Cordillera lithosphere-asthenosphere boundary (region ‘1’)

The SRFs from our study provide unique and reliable constraints on the physical properties of the CLAB (e.g., depth) due to their strong sensitivity to velocity gradient. The average CLAB depth is 75±17 km (see Fig. 2c), after applying timing corrections43 for mantle heterogeneity44, elevation and crustal thickness (CRUST1.045; Supplementary Fig. S3). Three parallel profiles of eastern Cordillera show multiple negative SRF arrivals associated with CLAB west of the RMT (see Fig. 3a, region ‘1’; Fig. 3b–d). Devoid of significant north-south variations (see Fig. 3b–d), a uniform CLAB is consistent with the earlier reported LAB depth (60-85 km) in the western portion of the Cordillera (e.g., refs. 46–48). According to synthetic SRFs calculated from a constant, low P-wave velocity layer (from 60 to 85 km depth; see Supplementary Fig. S9), a thin and hot eastern Cordilleran lithosphere11,28,49–51 is capable of producing a negative conversion with a 20–50% relative conversion amplitude to direct S. Our observed CLAB depth (70–75 km) is within uncertainty of the average depth (65±5 km) of chemical equilibration of Cenozoic to recent mantle-xenolith bearing alkaline lavas of the Canadian Cordillera15.

The simplest explanation for a negative shallow mantle velocity gradient is increased mantle temperature. For decades, a shallow asthenosphere with elevated temperatures has been proposed to isostatically compensate for a 30–48 km shallow, thin and flat Moho (e.g., refs. 2,52,53) and a negative Bouguer gravity anomaly54 beneath the uplifted mountainous terrain west of the RMT (~1500 m high; Fig. 3). Geothermal observations4,11,46,55,56 suggest that the Cordilleran back-arc mantle temperature at the CLAB exceeds 1200°C, nearly 400°C higher than those within the craton to the east at similar depths. Partial melt, small-scale convection and a wet upper asthenosphere (e.g., refs. 7,15,22,55,57) further complicate the thermochemical properties of the Cordilleran upper mantle. These defining characteristics diminish gradually near the eastern edge of the Cordilleran lithosphere, where the CLAB ramps upward towards the Moho of the NA craton near the RMT. Similar mantle conversions have been previously detected by Miller & Eaton33 at 100 km and 180 km depths beneath station SLEB (labelled ME2010, see Fig. 1).

It is worth noting that the long-period nature of S converted waves limits their spatial resolution in mapping the geometry of lithospheric discontinuities. Further complexities could result from corner diffractions, which weaken the SRF’s ability to resolve sharp deflections58, a relatively sparse station density (especially north of 49°N) and back-azimuth coverage. 3-D effects on ray paths may introduce additional uncertainties to our observations (detailed information is provided in the Supplementary material “Resolution”).

Laurentia lithosphere-asthenosphere boundary east of the foothills (region ‘2’)

In a 150 km wide area between the RMT and Cordillera Deformation Front (CDF), which we define as the Cordillera-craton transition zone (Fig. 3a, region 3), the 75-km deep CLAB transitions sharply to the LLAB at a regional average depth of ~207 km. The lithosphere beneath the Archean Hearne Province is notably thick, extending to ~220 km. This estimate is in general agreement with thermal thickness calculations (~200–220 km), defined by the intersection of a geotherm with the mantle adiabat59 (T=1300°C), as well as constraints from garnet xenoliths from depths below 200 km60. These depths are only moderately shallower than those proposed by seismic tomography5,61 (deeper than 250 km) and an earlier magnetotelluric survey62 (~250 km). Differences in data sampling and depth sensitivity37,63 could be partially responsible for the different reported craton depths. The localized LLAB is characterized as a strong amplitude in the SRF data (see Fig. 2c), and its depth is comparable to that of the Precambrian lithosphere64. Regional tomographic models (e.g., ref. 9) indicate shear velocity perturbations in excess of 4% beneath the Laurentian lithosphere, and up to 9% beneath the nearby cratons according to recent studies of SS precursors (e.g., ref. 65). The strong velocity gradients reflect a combination of physical and chemical properties variations within the Laurentian lithosphere, such as dehydration and past partial melting events associated with the cratonic lithosphere. The vertical variations in Mg#, water, and Fe59,66–68, as well as anisotropy variations10,69 could all contribute to a strong mantle gradient near the base of the NA craton. Differences aside, the compelling evidence of a relatively deep and stable LLAB favors a long-lived Archean Hearne lithosphere60,70 that has undergone minimal modifications since the Paleoproterozoic71.

The lithospheric root of the Hearne province shallows to ~200 km beneath the Medicine Hat block (MHB) and ~180 km beneath the Buffalo Head Terrane (BHT). For the former region (i.e., MHB) the Archean-aged lithosphere has undergone significant reworking due to Paleoproterozoic subduction and magmatism62,72,73, which may have resulted in partial removal of the craton root. In the case of BHT, evidence from garnet lherzolite xenolith and magnetotelluric surveys has favored a lithosphere as thin as 180 km39,41, which is in general agreement with our findings (~180 km) from SRFs. Partial melting triggered by hydration from subduction74 or upwelling of low-velocity asthenosphere75 are both capable of causing lithospheric modification and erosion15,76. These processes may have eroded and thinned the lithosphere beneath this region.

Cordillera—craton transition zone (region ‘3’)

An essential feature from the profiles is the west-dipping CLAB near the RMT (see Fig. 3b–d). This morphology is concordant with models of the Mid- to Late Cretaceous Cordilleran orogenesis, which implies that it resulted from interactions between an allochthonous ribbon continent and the NA craton19,20,77. Closure by westward subduction of the ocean basin east of the Cordillera, which continued into the passive west-margin of cratonic NA, and post-collision slab break-off78 were the most significant tectonic events at the former convergent boundary.

To ascertain the orientation of the subduction interface, we linearly fit the depths of the conversions in an area bounded by 150 km west of the RMT and 120 km east of it. The slope of the CLAB is consistent with a westerly dip of 6°, nearly identical to the slope of a gradient zone based on the 0% shear wave velocity contour in model NA1479. This dip angle is significantly lower than has been suggested by other studies using body6 (~70°) and surface28,80 (20–30°) waves at similar positions. We attribute this discrepancy to 1) sparse regional station coverage in the Western Canadian Sedimentary Basin, which hampers the spatial resolution of global surface-wave-based models, and 2) vertical smearing along ray paths8,35,81 from tomographic inversions, especially in the case of body waves. In other words, the boundary between NA and the Cordillera is gentler than a vertical ‘lithospheric step’. Its westward dip contradicts the accretion model, which proposes long-lived1 and landward-dipping (i.e., opposite) subduction beneath the western margin of the Cordillera, during which the Cordilleran lithosphere either underwent thinning or was detached17,25,78.

A related physical description of the Cordillera-craton boundary is ‘suture zone’. Without a consensus, various geophysical and geological studies have placed the eastern edge of the CLAB at the CDF8, RMT6,7, or the Omineca belt4,82. In the averaged SRF profiles (see Fig. 4), conversions associated with the CLAB are relatively diminished at a depth of ~38 km and a position that is ~100 km east of the RMT (near the eastern end of profile). Despite minor variations among different transects (see Fig. 3), e.g., a westward extension of the Cordillera-craton transition zone at ~47° latitude in profile C, the geometry of the CCB relative to the RMT appears stable and uniform across the latitude range from 46°N to 50°N. According to Audet et al.1, this west-dipping interface may extend as far north as northern British Columbia (~60°N), based on their experiments involving P-to-s converted waves, seismic velocities, geothermal modelling and time-dependent thermal simulations.

Western LLAB and dual-layered lithosphere architecture

Judging from the westward continuity of the LLAB (see Fig. 3b–d), the western edge of the intact craton core extends to ~200 km west of the suggested locations from recent tomographic inversions6,28. The upper surface of the intact craton lithosphere is delineated by two semi-parallel, westward dipping (labelled CLAB and phase ‘X’ in Fig. 4) S-to-p conversions with an average separation depth of 50 km. Phase ‘X’ is interpreted as a relic of an oceanic plate that now lies directly above the cratonic lithosphere78. By linearly fitting the depths of phase ‘X’ (see Supplementary material for a robustness estimate) and LLAB (see Fig. 4), which are associated with the top and bottom surfaces of the NA craton, we estimate the western craton margin to be positioned at ~320 km west of the RMT, at a depth of 165 km. This depth is notably shallower than the regional average of ~220 km, a relatively stable value expected for a craton margin that has not undergone major deformations in the past 100 million years50,83,84.

We attribute the progressive westward thinning of the NA craton to erosion, in combination with an uplift caused by a highly convective sub-Cordilleran mantle76,85. Observations of seismic anisotropy86 and geodynamic modelling14,87 suggest that a southwestward mantle wind could exceed 4 cm/yr in the back-arc of the Juan de Fuca subduction zone88,89. In spite of the strength and rigidity of the craton core7, the severity of the convective flow over the past 100 million years is fully capable of shearing the craton margin and elevating the boundary via the pressure differences between the top (low pressure) and bottom (high pressure) of the craton lithosphere14.

A boundary zone dominated by westward subduction resulted in a dual-layered lithosphere architecture in a 300–500 km wide window between the eastern intermontane and the RMT (Fig. 5). This complex geometry includes cratonic Laurentian lithosphere that extends down to ~180 km, directly underlying the Cordilleran lithosphere, which terminates at a depth of ~75 km. The geological evolution that led to the development of the stacked lithospheres includes Neoproterozoic rifting and the subsequent formation of the western Laurentian passive margin; the Devonian to Cretaceous development of a Cordilleran ribbon continent that lay west of cratonic NA, and the mid-Cretaceous closure of the basin separating the ribbon continent from cratonic NA by west-dipping subduction beneath the Cordilleran ribbon continent’s eastern margin. Post-collisional slab break-off at ~95 Ma detached the westward subducting oceanic slab that subsequently foundered into the deep mantle26, leaving behind the westward dipping Laurentian continental lithosphere that had been pulled into and beneath the Cordilleran upper plate by the subducting slab78 and a small remnant of the former oceanic plate (see Fig. 5).

Cratonic lithosphere, when continuous with and attached to subducting oceanic lithosphere, will be dragged into and beneath thinner, arc-bearing upper plate lithosphere. It becomes resistant to subduction due to buoyancy forces, and the resulting duplication of the LAB by lithospheric stacking across a collisional convergent margin is predicted by plate tectonics. The complex lithosphere architecture revealed by our SRF observations is consistent with 1) attempted westward subduction of cratonic Laurentia beneath the Cordillera19,20, and 2) tomographic models that image the Laurentian cratonic lithosphere juxtaposed against and extending beneath the Cordilleran lithosphere1,6. A modern-day parallel is the Tibetan Plateau, where attempted northward subduction of cratonic India beneath Asia gave rise to a lithospheric stack and the development of a complex, multi-layered lithospheric architecture down to 250 km90. Older examples include the Paleoproterozoic stacking of cratonic lithosphere of the western Slave Province beneath the eastern Slave Province to depths greater than 200 km91, and the Archean Sask Craton that was, together with its diamondiferous mantle root, completely overthrust by the Hearne Craton during the Paleoproterozoic Trans-Hudson orogeny92.

Collectively, by taking advantage of SRFs from regional seismic arrays, we are able to reach a few key conclusions pertaining to the Cordillera-craton transition region,

1) The base of the Cordilleran lithosphere directly overlays a wedge-shaped Laurentian lithosphere that bottoms at ~180 km. Together, they form a unique dual-layered lithosphere architecture.
2) The western NA craton boundary dips westward at an approximate angle of 6 degrees. Its present-day morphology may be strongly influenced by Mesozoic subduction, slab break-off and mantle convection.
3) The NA craton edge shows evidence of convective erosion and possible uplift within a window from the RMT to the Intermontane belt.

Our study offers compelling evidence for a multi-layered lithospheric architecture beneath the eastern Cordillera. Granted, many of the inferences from SRFs may be impacted, to varying degrees, by limited azimuthal coverage, depth sensitivity and 3-dimensional nature of the ray paths. They do, however, provide an impetus for a better understanding of the Cordillera-craton boundary evolution, especially when aided by improved constraints from seismic tomography, geological sampling and geodynamic modelling.

Methods
This study utilizes the frequency domain water-level S receiver function technique93,94 to analyze S-to-p converted waves resulting from mantle interfaces (see Fig. S1(a)). Since S-to-p converted waves travel faster than direct S, their waveforms are minimally affected by multiple-reflected S. This is a distinct advantage over P-to-s conversions, which are commonly used in crust and mantle imaging but suffer from contaminations from P wave coda at lithospheric depths69,94,95 (see Fig. S1(c)).

Our S receiver function technique contains two initial steps: 1) signal partitioning that isolates the incident from scattered wave fields, and 2) deconvolution that minimizes source and propagation effects94. The resulting receiver functions are subsequently corrected for epicentral distance (via time shifts) and stacked according to the locations of the conversion points (see below).

Dataset

Our dataset includes 15 permanent stations from the Canadian National Seismograph Network96 (CNSN) spanning 2014 to 2019, 40 temporary stations from USArray operational between 2006 and 201097 (TA), 24 stations from the Regional Alberta Observatory for Earthquake Studies Network24 (RAVEN), and 28 stations from the Canadian Rockies and Alberta Network23 (CRANE) covering the same period. To ensure optimal receiver function results, we limit our earthquake selection to magnitudes equal to or greater than 5.5 and epicentral distances ranging from 55° to 85°. Additionally, only earthquakes with depths shallower than 150 km are included to mitigate the impact of depth phases98,99. All waveform data were obtained from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center based on events listed in the USGS National Earthquake Information Center catalogue. We processed over 20,000 waveforms from more than 2,000 seismic events to produce about 3,900 high-quality SRF traces.

Wavefield isolation

Before deconvolution, we apply a transformation matrix to the recorded wave field to enhance weak converted waves by recasting the coordinate into the L-Q-T system95,100,101 (see Fig. S1(a)). First, we convert local N-E-Z system into R-T-Z system by:

[R T Z]T = [cosγ sinγ 0; -sinγ cosγ 0; 0 0 1] [N E Z]T (1)

where N, E, and Z are the North-South, East-West, and vertical components of the seismogram, respectively. The outputs R, T, and Z are the radial, transverse and vertical component seismograms. Inside the transformation matrix, γ represents the back azimuth of the incident wave. Then, we rotate R–T–Z components into L-Q-T coordinates by:

[L Q T]T = [cos i sin i 0; -sin i cos i 0; 0 0 1] [Z R T]T (2)

where i=arcsin(Pi) is the angle between incident P-wave and the vertical axis. After the rotation, S-to-p converted waves are polarized approximately along the L direction.

Deconvolution

For each event-station pair, we apply a third-order Butterworth filter with corner frequencies of 0.02 Hz and 0.1 Hz, or 0.01 Hz and 0.2 Hz, parameters that are commonly adopted by earlier studies of North America (e.g., refs. 69,98,102). We then deconvolve the L component from the Q component using frequency domain water-level deconvolution103,104. To obtain optimal results, a water-level parameter of 0.01105 was chosen empirically (see Fig. S2) to balance resolution and stability (Supplementary Data 1).

The water-level deconvolution computation follows the equation:

R(ω) = [Z(ω)R*(ω)] / [φ(ω)] G(ω) (3)

Where,

φ(ω) = max(R(ω)R*(ω), c^2) max(Z(ω)Z*(ω), c^2) (4)

G(ω) = ξ exp(-ω^2 / 4α^2) (5)

In these three equations, * denotes the complex conjugate and u(ω) is the seismic signal after transforming to the frequency domain. The constant c is the water level that controls the minimum amplitude allowed in the denominator during deconvolution. In Eq. 5, G(ω) is the conventional Gaussian filter for smoothing, the constant ξ normalizes the filter to unit amplitude, and α controls the width of the Gaussian pulse.

Deconvolved seismograms are then migrated to depth based on the original PREM34. The uncertainty introduced by the 220 km discontinuity in PREM during migration is estimated to be within ±5 km (see Supplementary Figs. S11 and S13). All receiver functions recorded by a given station are then stacked to form a single trace. It is worth noting that only minor differences could be found among different binning procedures (e.g., stacking according to conversion points). The time corrections for lateral heterogeneity43,44 and crust45 (CRUST1.0) are applied based on 1-D ray tracing. The resulting correction value (Fig. S3(a) to (c)) is then applied to all SRFs (see Fig. S3(d)). We visually inspect and reject traces with high pre-signal noise and abnormal reverberations. SRFs have been reversed in their polarity and depth axis to represent a velocity decrease using negative arrivals.

Negative phases after the Moho conversion (~35 km beneath Cordillera and ~40 km beneath craton) are generally caused by reduced velocities (e.g., at the base of the lithosphere). Supplementary Fig. S4 shows an S receiver function example.

Bootstrap resampling. We estimate the uncertainty of SRF by the bootstrap resampling method106. For each station, 20% of quality-controlled SRF traces are randomly selected to create a bootstrapped stack, repeated 100 times to generate multiple stacks. Then we compute the mean and standard deviation of the distribution69,106,107, the latter of which offers an effective estimate of the uncertainty. Well-behaved S receiver functions generally fit within the +2σ (95% confidence level) of the bootstrapped mean. Similarly, errors in depth (d) and amplitude (a) are quantified by their departures from the respective mean values at the base of the lithosphere (see Fig. S1).

Data availability

Seismic data for CNSN network can be requested from Canadian National Data Center (https://www.earthquakescanada.nrcan.gc.ca/stndon/CNSN-RNSC/index-en.php). The original CRANE seismic records can be accessed upon request from the seismic network operator. The IRIS Data Management Center (http://ds.iris.edu/ds/nodes/dmc/) was used to access waveforms, and data from the TA network were made freely available as part of the EarthScope USArray facility. Both data sources are supported by the National Science Foundation under cooperative agreement EAR-1261681. The processed seismic data and S receiver functions generated in this study have been deposited in the corresponding Code Ocean capsule108.

Code availability

The codes of receiver function analysis are available from the Code Ocean capsule108.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript makes an interesting and valuable contribution to ongoing studies of the nature of the craton-Cordillera transition in western North America (specifically across Alberta, Canada, and the northern US). A large number of available seismograph stations is used to track upper-mantle discontinuities using S-wave receiver functions (SRF) and interpret them in the context of the complex interaction between cratonic and Cordilleran lithosphere. The results are of interest to a broad community, especially in terms of their geodynamic significance in terms of reconstructing the nature and evolution of the boundary at depth.\n\nThe SRF dataset is robust and carefully processed, with attention to characterising uncertainties in the waveforms (via a bootstrapping technique) and assuring that the final waveforms are of the highest possible quality. The methodology for calculating the SRF and making comparisons with appropriate synthetic data is valid and clearly explained.\n\nThere are a few issues that I flagged while reading the manuscript for which I have scientific questions that need to be more clearly addressed in order to fully assess the results and judge the full validity of the interpretations/conclusions, as follows:\n(1) The authors migrate the SRFs from the time domain to the depth domain using the PREM (Dziewonski & Anderson, 1981) global reference model. The original version of this model includes a significant discontinuity in shear-wave velocity at 220 km depth, i.e., fairly close in depth to the inferred Laurentian lithosphere-asthenosphere boundary (LAB). Many seismic studies using PREM in some way (e.g., as a starting model for inversions) use a smoothed version in which the discontinuity is removed, to avoid any potential artefacts. Was a smoothed PREM model used in this study, or the original version including its discontinuity?\n(2) The Cordilleran LAB is identified as a slightly dipping interface by the changes in its depth for stations along the east-west profile. In the original single-station SRFs, are there any back-azimuthal variations that could also shed some light on the dipping structure? Or energy on the Q-component? In the supplemental material we see an example of SRFs sorted by epicentral distance; how would such SRFs look when sorted by back-azimuth, for a station where the dipping CLAB is visible in the stacked data?\n(3) The S20RTS global tomographic model (Ritsema et al., 1999) is used to apply timing corrections to the SRF. What is the motivation for using such an old and low-resolution model rather than either a more modern global model (e.g., S40RTS of Ritsema et al., 2011 or some of the newer ones) or continental-scale model? This is not clear from reading the manuscript.\n(4) Phase ‘X’ is identified in the figures and mentioned in the last paragraph of the discussion. From the images in Fig 3 (and equivalents in supplementary material), the phase is very difficult to see, and assessing its robustness (e.g., compared with other peaks and troughs in the section) is therefore also difficult. From the figures, it looks as if the phase ‘X’ peaks are no higher in amplitude than some other peaks in the SRF that are not attributed to a significant feature. An extra figure (e.g., the SRF section like Figs 3, S9 and S10 but without the tomography background), and extra justificative text, in the supplementary material, would be valuable.\n\nRelated to the questions posed above, it may be instructive to test the effects of (i) the original versus a smoothed PREM model, if appropriate (depending on the version used in the current migrations), and (ii) the effects of using a different mantle tomographic model to carry out the time corrections on the SRFs.\n\nThe citations are appropriate and comprehensive in general, though some minor corrections are necessary, and there are a few superfluous references (see below for details).\n\nIn general, the manuscript is clearly written and well presented. Minor corrections (e.g., typographic or grammatical errors, or text that needs a little extra clarity) are listed below.\n\nOverall, I recommend moderate revisions for this manuscript, with particular attention to the scientific questions detailed above.\n\nSpecific (line-by-line) comments\nLines 8-9 (Abstract) – How relevant to this abstract is the reference to an upper-crustal fault, given that the focus of the paper is the LAB and deep lithospheric architecture? A statement like this would be fine in the introduction, but it risks giving readers the wrong impression on what the paper is really about. If the authors decide that it is important to include, the word “that” must be removed to correct the grammar.\nLine 12 (Abstract) – The term ‘underplate’ is usually used in the context of magmatic processes (e.g., mafic underplating of the crust due to plumes, rifting, etc.). “underlies” would be more appropriate here.\nLine 17 – The term ‘flat’ is fine for a boundary like the Moho or LAB, but does not work for the lithosphere as a whole; perhaps “uniform” here?\nLine 23 – Not sure how useful the reference to Bedle & van der Lee is, given the more recent references to continental-scale models that follow it.\nLine 24 – “H. Yuan” -> “Yuan”.\nLine 41 – “Gu et al. 2010” -> “Gu et al. 2011” (see remark in References section below).\nLine 54 – “lithosphere boundary” -> “lithosphere-asthenosphere boundary” (?)\nLines 68,290 – See scientific questions above regarding the use of PREM for migration.\nLine 69 – “stcacking” -> “stacking”\nLines 71-73 – The Artemieva, Priestley (“2017” -> “2018) and Rychert references are directly relevant here, as they specifically talk about the LAB, but the others are global tomography models that do not discuss (and often do not directly resolve) interfaces, so their inclusion is misleading. Also, the Masters citation does not correspond to any existing publication (see remark in References section below).\nLine 82 – “thin lithosphere” -> “relatively thin” or “relatively thinner”, since a 180-km lithosphere is still significantly thicker than the global average.\nLine 87 – See scientific questions above regarding dipping phases.\nLine 91 – “These” -> “This”.\nLines 98,294 – See scientific questions above regarding the use of the Ritsema et al. (1999) model for time corrections.\nLines 102, 115 – Citation for Hansen should be 2015 (see remark in References section below).\nLine 104 – “p-wave” -> “P-wave”.\nLine 106 – “relative strength to” -> “strength relative to”.\nLine 107 – Needs some further comment; how does one define ‘marginally’ in this context? Is an 8-13km difference negligible or significant, and why?\nLines 110-111 – A negative velocity gradient is explained here by a reduced mantle temperature, but velocity and temperature have an inverse relationship (higher temperature = lower velocity). Is this a mistake, or an aspect that I am missing?\nLine 123-124 - This statement is a bit confusing. From the current wording, it gives the impression that at 200 km depth there is a 3% Vs difference between Hearne & Medicine Hat, but this is not evident in Fig 3. Is this a reference to a velocity difference between Hearne & MH, or between the Cordillera and Hearne/MH region?\nLine 130 – “consistent with” requires some clarification. It is consistent in that the elastic thickness for cratonic lithosphere is significantly greater than that for the Cordillera, but the absolute values (220 km for seismic, ~100 km for Te) are very different, and this should be made clearer to avoid confusion.\nLine 130 – the Artemieva & Mooney paper cited here does not deal with elastic thickness, but thermal thickness.\nLine 131 – “xenolith” -> “xenoliths”.\nLines 142-143 – Indeed there is evidence for processes that could lead to lithospheric erosion, but the difference in lithospheric thickness quoted here between Hearne and MHB is only 20 km (220 vs 200), which could be argued to be ‘marginal’ and perhaps within uncertainties of seismic methods. How significant is this decrease really?\nLine 155 – “Cordillera” -> “the Cordillera”.\nLines 160, 658 – The tomographic model NA14 refers directly to the one described by Yuan et al. (2014), and I am not sure that there is any direct link to the Yuan & Romanowicz (2010) publication.\nLine 170 – “DeCELLES” -> “DeCelles”.\nLine 184 – Perhaps worth specifying “earlier body-wave tomographic inversions”? Figures 3, 4 and S9 actually give the impression of high-velocity lithosphere extending quite a long way west beneath the Cordillera in the surface-wave/waveform tomographic models.\nLines 186, 215-221 – See scientific questions above regarding phase ‘X’.\nLine 204 – See comment for Line 12 regarding the term ‘underplating’. Perhaps “underthrusting of the cratonic lithosphere beneath the Cordilleran lithosphere” would be a more appropriate wording here.\nLine 225 – “of Cordillera” -> “of the Cordillera”.\nLines 234-235 – “beneath eastern Cordillera” -> “beneath the eastern Cordillera”.\nLines 236-237 – regarding ‘profound implications’; the new observation is certainly important, but the current wording is rather vague. Is it possible to expand on this?\nLines 250-251 – there are official citations/references in existence for CNSN and USArray; these should be included.\nLine 253 – “2010” -> “2011” (see remark in References section below).\nLines 274-275 – There are two different filter passbands mentioned. What decides which of the two options is used? Is the application of a given band applied systematically by station, or individually on a seismogram-by-seismogram basis?\nLine 289 – What value of alpha was used in this study, and to what maximum wavelet frequency does it correspond?\nReferences:\nLine 324 – What is “932” in the Aki & Richards reference? This reference also needs publisher information.\nLine 325 – Journal name is missing.\nLine 344 – Spurious word “Received”; also, this reference should be placed with the “K”s.\nLine 358 – Missing volume/number.\nLine 394 – Remove superfluous “(December 2022)”.\nLines 402-404 – Duretz et al reference is not cited in the text.\nLine 426 – “mantlelithosphere” -> “mantle lithosphere”.\nLine 433-435 – Year should be 2011; also the journal and doi details are incorrect. Should be: \"Seismological Research Letters, 82(4), 575-588, doi:10.1785/gssrl.82.4.575\".\nLine 436-438 – This appears to be at least 2 different references mashed together.\nTitle belongs to Hansen, Dueker & Schmandt 2015 (EPSL); doi belongs to Yu et al 2022 (JGR). Neither are in CJES and, in 2010, CJES was only up to volume 47.\nLine 445 – “cretaceous-tertiary” -> “Cretaceous-Tertiary” (these are proper names).\nLine 449 – “love wave” -> “Love wave”\nLine 496 – Specify that this is an MSc thesis.\nLines 498-499 – Remove “Linked references are available on JSTOR for this article”. Journal name is missing. Google Scholar search finds a paper from 1996 (Phil Trans Roy Soc London) with this title. Guy Masters has no 2019 publications.\nLine 500 – Needs capital letters: “The Upper Devonian Aley carbonatite”.\nLine 528 – Incomplete reference; is this a book, journal, other publication? Remove “February 2002”.\nLine 530 – Year should be 2018. Needs book title, editor names & publisher.\nLines 559-560 – Journal and doi details are incorrect. Should be: \"Journal of Geophysical Research: Solid Earth, 125(8), e2020JB019570; doi:10.1029/2020JB019570\".\nLines 561-562 – What is “Development, January, 1-182”? This reference is a PhD thesis and therefore needs correction.\nLines 568-569 – Reference should be: \"The Regional Alberta Observatory for Earthquake Studies Network (RAVEN), CSEG Recorder, 40(8), 34-37\"\nLine 571 – “alberta” -> “Alberta”\nFigures:\nFig 1, legend – “Archen” -> “Archean”; “accreated” -> “accreted”.\nFig 1 caption, lines 626-629 – Abbreviations should go with the explanation for panel (a).\nFig 2 – Do any of the SRFs shown in panel (b) contain the ‘X’ phase? If so, it would be useful to indicate it. If not, it would be useful to add a figure (in supplementary information?) showing the identification of this phase, and justifying the choice of the particular peaks. The waveforms are not easy to see against the dark blue of the high-Vs regions of the tomographic model, so the phase identification is not really clear from the figures.\nFig 3 caption, line 663 – The ‘X’ phase is explained as a “major velocity decrease in each profile”. However, in profiles (b) and (c), the X-line seems to correspond to an *increase* in Vs according to the surface-wave model, and in profile (d) it's pretty much in the middle of the high-Vs region.\nFig 4 – Top of main panel: “NA raton” -> “NA craton”.\nFig 4 – The cyan circles used to show the LAB picks get lost against the light blue of the background velocity model. Suggest changing to a different colour that shows up better.\nSupplementary material:\nFigs S2-S4: Specific station names (EDM and BRU) are given in reference to SRF waveforms. For the benefit of readers unfamiliar with Canadian seismographs, please indicate these two stations in one of the maps in Fig S3.\nFig S3a,c – As stated by the authors, most corrections are the same sign (negative for crust/topography, positive for the combined effects) for the reasons noted. There are a few stations with opposite corrections that seem spurious given their proximity to ‘normal’ stations (e.g., at ~54N, 119W in (a), 54N, 117W in (c)). What could cause these discrepancies over such a short lateral distance?\nFig S3 caption, line 28 – Why use S20RTS for the mantle corrections, when it is a 25-year-old tomography model that has been superseded many times over by higher-resolution models?\nFig S4 – There is quite a lot of variation, which does not seem systematic, in the arrival time of the LAB phase in the individual RFs when sorted by epicentral distance and depth-converted. For this type of plot, it would be more instructive to see the SRFs plotted as a function of time instead of depth, and to add a second plot where the y-axis is back-azimuth instead of epicentral distance.\nFig S6 caption – It would be useful to add that the synthetics are shown for vertical and radial components (at least, I assume that is what the “Z” and “R” denote).\nFig S7 caption – For extra clarity, state that the synthetic SRFs are computed from the model shown in Fig S6.\nFig S9 caption, line 121 – The velocity model SL2013NA should be cited as Schaeffer & Lebedev 2014. The 2013 publication was the presentation of the global model SL2013.\nFig S9 caption – For image (b), please explain the magenta horizontal lines in the depth section. For the ±1% velocity contours, specify “The solid white lines”.\nLine 191 – Please provide the citation (and corresponding reference) for “the NA14 velocity model”.\nLine 192 – “2013” -> “2014” (see comment about Fig S9 caption).\nSupplementary references:\nLine 211 – What is “932” in the Aki & Richards reference? This reference also needs publisher information.\nLine 232 – Missing journal and volume information; month is not required.\nLines 240-241 – Replace with Schaeffer & Lebedev (2014, EPSL) reference.", "reviewer_2": "This manuscript presents a dual-lithosphere structure, consisting of the Cordillera LAB and Laurentia LAB, beneath the eastern Cordillera. It also identifies a westward-dipping Cordillera LAB, as opposed to a \"lithospheric step,\" across the transition from the Cordillera to the North American Craton in southwestern Canada using Sp receiver function imaging. The manuscript is well-organized, and the observations are clearly articulated. However, my primary concern pertains to the novelty of the reported findings, as similar features have been documented in these tectonic regions in prior studies, albeit with different datasets. Detailed comments are provided below.\n\nPrevious studies have reported comparable structures in the same tectonic regions. For instance, Audet et al. (2019) identified a westward-dipping upper mantle discontinuity using Ps receiver functions (see their Fig. 2b). Additionally, Miller and Eaton (2010), employing the same Sp receiver function technique as the authors, documented a dual-discontinuity structure (refer to their Fig. 4, station SLEB). While these studies are cited, the authors should include a more in-depth comparison and discussion of the similarities in results. This would help readers appreciate how the current work confirms and extends previously observed structures. Including these prior observations in Fig. 1 would also provide readers with a clearer contextual background.\n\nFurthermore, the final model proposed in Fig. 5 of the manuscript closely mirrors the models presented in Audet et al. (2019) (their Fig. 5). The authors should discuss the similarities and differences between these models and highlight any new constraints their study offers. Additionally, as noted by Audet et al. (2019), seismic imaging alone may not fully resolve these models. It would be beneficial for the authors to address this limitation in their discussion. Moreover, could the authors kindly provide a more detailed discussion on why a dual lithosphere matters for the broader scientific questions regarding the structure of the lithosphere and its origins? A review of current hypotheses on this topic might help place their work within the context of ongoing debates in the field.\n\nOther minor comments:\nFig. 2: I don’t see any red rectangles in a (location of stations used in b).\nFig. 3b-d: the legend labels seem to be off.\nSp-RF resolution: Lekic and Fischer (2017) provided some useful information on the relationship between station spacing/frequency content and the ability to resolve dipping structures using Sp-RFs (e.g., their Fig. 5). The authors may want to consider mentioning this to prove the robustness of their observation.\n\nReferences:\nAudet, P., Currie, C. A., Schaeffer, A. J., & Hill, A. M. (2019). Seismic evidence for lithospheric thinning and heat in the northern Canadian cordillera. Geophysical Research Letters, 46(8), 4249–4257.\nLekić, V., & Fischer, K. M. (2017). Interpreting spatially stacked Sp receiver functions. Geophysical Journal International, 210(2), 874–886.\nMiller, M. S., & Eaton, D. W. (2010). Formation of cratonic mantle keels by arc accretion: Evidence from S receiver functions. Geophysical Research Letters. https://doi.org/10.1029/2010GL044366", "reviewer_3": "The submitted manuscript provides an in-depth analysis of the geometry of the lithosphere-asthenosphere boundary beneath the transition from the North American cordillera to the Laurentian craton using Sp receiver function analysis. This analysis is generally considered to be the gold standard of LAB imaging. The authors use their results to argue for a dual lithosphere architecture where the craton underplates the Cordilleran lithosphere.\n\nMy overall assessment of the manuscript is as follows:\n• I think that the manuscript and the results are interesting, and I enjoyed learning about the work of the author(s). If this were to be published, I would certainly read it again and reference it in my own work.\n• The manuscript itself is well written, and the figures are generally clear and easily understood even though there is a lot of information contained within them.\n• Although I am not particularly familiar with the debate surrounding the tectonic evolution of the Cordillera in North America, I thought that the authors did a good job of discussing other constraints (other seismological methods, other geophysical methods, geochemical, etc.) within their paper. I am not going to comment on the particulars of their conclusions as they relate to the history of the Cordillera as I imagine others are better versed in the particulars.\n• The methodology as presented is straightforward and without significant controversy. I use Sp receiver function analysis frequently and found no issues in the (pre)processing, deconvolution, bootstrapping or migration.\n• I do have two major concerns related to the results themselves. My first concern is what the results imply about the structure of the lithosphere-asthenosphere boundary beneath the cratons and how the results of this study compares to similar studies from the same region as well as other cratons globally. My second major concern is the lack of discussion regarding resolution and ability to draw conclusions about structure as a result. I provide more detailed discussion of these two concerns below.\n\nConcern #1: LAB beneath cratons\nThe base of the cratonic lithosphere is commonly thought to mark a gradual transition from lithosphere to asthenosphere, the result thermal cooling. The resulting velocity gradient is thought to be small, with only a few percent velocity change occurring over 50-100 kilometers (e.g., Fischer et al., 2010, Annual Reviews). As a result, it is thought that receiver functions, including Sp receiver functions, are relatively insensitive to this boundary (Mancinelli et al., 2018; GRL). This makes it difficult for receiver functions to image the boundary and many studies (e.g., Hopper and Fischer, 2018, G-cubed) show a dearth of converted phase energy at LAB depths beneath the cratons. Therefore, I was quite surprised to see how substantial in amplitude the LAB phase associated with the craton is within this study and had hoped that the authors would comment on this observation. In fact, in their summed stacks for the two regions, the cratonic LAB appears to be nearly equal in amplitude to the Moho phase, suggesting comparably sharp velocity gradients.\n\nThe authors do provide some cases of velocity modeling in the supplement but never directly relate this to a modeled velocity gradient for the cratonic LAB. At a minimum this needs to be explored and discussed. I ask that the authors consider why they observe large amplitude LAB phases beneath the craton, in contrast to what is typically observed globally. I also ask that if they do model velocities, use this to put constraints on what this might imply regarding the physical properties of the boundary.\n\nRelated to this point, the study of Hopper and Fischer (2018, G-cubed) overlaps with the study area of this manuscript. Specifically, their cross-section A-A’ in Figure 1 and Figure 3. In their paper, there is little negative phase energy associated with the LAB within the craton east of the Cordillera. Instead, they image a shallower boundary (referred to as a mid-lithospheric discontinuity). A discussion comparing these two sets of results, reasons for their discrepancies, and any ways in which the results can be reconciled (and what that may mean tectonically) is necessary to make this manuscript strong.\n\nConcern #2: Sp receiver function sensitivity, resolution and finite frequency effects\nIt is clear that there is lateral complexity in the structure of the lithosphere beneath the study area. The authors use this complexity to deduce something about the tectonic history of the region. However, there is little discussion of the limitations of the methodology or how issues such as sampling off of the ray path, or the frequency limits the resolution of the analysis, may affect their interpretation. For example, in the first paragraph of section 3 the authors describe the dip of the CLAB is a westerly 6°, which they argue is nearly identical to the slope of the 0% shear velocity contour of model NA14. I think that the authors need to be more careful when discussing this point. Similarly, the discussion of the location of the suture zone in the same section should acknowledge limitations. As I suggestion, I would encourage the authors to look at a paper that focuses on the 2D modeling of Sp receiver functions (Lekic and Fischer, 2017 GJI). In their paper they demonstrate that things such as finite frequency effects in the presence of dipping or lateral complex structures, can complicate the interpretation of Sp receiver function results. Given the relatively coarse station spacing and binning by station of this study, it is possible that other structures off of the assumed ray path may be influencing the interpreted structure.\n\nMinor comments and suggestions:\nAt the start of the second paragraph of the section titled “S receiver function imaging”, “stacking” is spelled incorrectly. In the same sentence the authors refer to boostrapping under a “three-sigma criterion” but every other instance of uncertainties in the paper refer to two sigma. Clarification on what is meant would be appreciated.\nUnder the section on “Data Availability” the authors should refer the EarthScope DMC and not the IRIS DMC. More specifically, they should use the wording provided by EarthScope as it relates to the SAGE award (https://www.earthscope.org/how-to-cite/).\nBe more careful when referring to amplitude, strength, gradient, etc. throughout the manuscript. Assume not all readers are well versed in what “strength” might mean when referring to amplitude.\nA comment on the availability of code and data – not sure if Nature Communications has different standard than other journals but “available upon request” is a less commonly accept format than in the past. I understand that data is sometimes under an embargo, but the deadline for when it should be made available, and at which data center, should be made clear.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We are grateful to the four anonymous reviewers and the Editor, Laura Frahm, for their constructive comments and professional handling of our manuscript. We took extra time and care in our response letter and revisions to address all comments. This document contains our point-by-point responses to the questions and comments raised by reviewers. Please also refer to the revised manuscript, which includes marks highlighting detailed modifications. The original comments and questions are shown in red, our replies are in black, and the quotes (texts and figure captions) from the revised manuscript are shown in blue.\nEditor (Remarks to Author):\nThank you again for submitting your manuscript \"Dual-layered Mantle Lithosphere Beneath Southeastern Canadian Cordillera\" to Nature Communications. We have now received reports from 4 reviewers and, after careful consideration, we have decided to invite a major revision of the manuscript. You will see that the reviewers are interested in your study but also raise the need for improvement. In particular, they mention the need for a better comparison with other studies while highlighting the advances of this study (ref#2), a discussion about the observed large amplitude LAB phases (ref#3) as well as the limitations of the methodology (ref#3) and in general more information about the used models (ref#1).\nReply to Editor\nWe sincerely thank the Editor for the professional handling of our manuscript. We appreciate the insightful questions and comments raised by both the Editor and the reviewers, and we have made a concerted effort to address them thoroughly in the revised manuscript. In this revision, we have:\n Added further comparisons with previous studies in the study area.\n Examined the LAB phase amplitude by using synthetic S receiver functions (SRFs).\n Examined errors arising from the migration model’s structural assumptions.\n Addressed the methodological limitations of the S receiver function approach."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have provided a detailed and comprehensive response-to-reviewers document, comprising explanations, descriptions of modifications, justifications of approach, revised text/figures used in the new manuscript, and figures created for the response document. I appreciate this thorough approach and the care taken with the responses, including several new tests of the seismic data and models.\n\nThe new version of the manuscript has been extensively reworked, with new information that helps readers assess the robustness and novelty of the study. In my opinion, it is almost ready for publication; I have just a few very minor wording changes to suggest; see below.\n\nRequested modifications (line numbers refer to tracked-changes manuscript version):\nLines 84-86: The additional information about the effect of the PREM 220-km discontinuity is useful, but it needs reference to the Supplementary material so that readers can see that the stated error value has been determined by testing (this is done in the Methods section, lines 381-382, but should be added here too).\nLine 91: Given the vertical resolution of the SRF, the value of “±2.75 km” for the uncertainty seems unrealistically precise; suggest quoting it to the nearest kilometre (i.e., “±3”).\nLine 101: “relatively thin” would read better as “relatively thinner” for extra clarity.\nLines 103, 121, 129: “eastern Cordillera” -> “the eastern Cordillera”.\nLine 134: “cordilleran” -> “Cordilleran”.\nLine 137: “Canadian Cordilleran” -> “Canadian Cordillera”.\nLine 138: Repetition; suggest changing to “…explanation for a negative shallow mantle velocity gradient”.\nLines 180, 183: “Laurentia […] lithosphere” -> “Laurentian […] lithosphere”.\nLine 181: “SS precursor” -> “SS precursors”.\nLine 186: “thick and stable” -> “deep and stable” since the text refers to the depth of the boundary rather than the thickness of the lithosphere above it.\nLine 202: “survey” -> “surveys”.\nLine 210: “imply” -> “implies that”.\nLine 246: With the new text quoting a depth of 165 km, the “(< 170 km)” in the next sentence is no longer needed.\nLine 286: Since tomographic models cannot directly resolve faults, I suggest removing the word from “fault juxtaposition”.\nLine 395: “Then compute” -> “Then we compute”.\nLine 411: Is there a suitable citation/reference for Code Ocean? I have not heard of this repository before, and other readers may also not be aware of it, so a reference or a URL is needed here.\n\nSupplementary material:\nLines 258, 322: “the PREM” -> “PREM” or “the PREM model”.\nLine 288: “that our SRF” -> “than our SRF”.\nLine 308: “ocean-ward” -> “oceanward”.", "reviewer_2": "This paper has been significantly revised and we think it is ready for publication.", "reviewer_3": "This is the second review that I am providing for this manuscript. As I noted in my original review, the manuscript is well-written, thoughtful and presents interesting results. Looking through the comments to the reviewers, as well as the revised manuscript, it is clear that the authors put significant thought and effort into responding to the comments. While I may not wholly agree with their interpretation in some instances, I think that any disagreement is worth follow-up study at another time and should not impede the publication of the manuscript. I fully support the quick publication of this work and thank the authors for their careful consideration of the reviewers comments.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Subject: Nature Communications manuscript NCOMMS-24-75865\n“Dual-layered Mantle Lithosphere Beneath Southeastern Canadian Cordillera”\nPoint-by-point response to the reviewer’s comments\nWe sincerely thank the four anonymous reviewers and the Editor, Laura Frahm, for their careful review and professional handling of our manuscript. We appreciate all the kind words from the reviewers, and are glad that our efforts in the last round of reviews were appreciated by all four reviewers. Since Reviewers #2 and #3 are satisfied with our manuscript and did not raise additional questions/suggestions, this document will mainly focus on the additional comments from Reviewer #1 and #4.\nPlease also refer to the revised manuscript, which includes marks highlighting detailed modifications. The original comments are shown in red, our replies are in black, and the quotes (texts and figure captions) from the revised manuscript are shown in blue. All the line numbers in our response will refer to manuscript with markups."}]

Decision letter
Editor (Remarks to Author):
Thank you again for submitting your manuscript "Dual-layered Mantle Lithosphere Beneath Southeastern Canadian Cordillera" to Nature Communications. We have now received reports from 4 reviewers and, after careful consideration, we have decided to invite a major revision of the manuscript. You will see that the reviewers are interested in your study but also raise the need for improvement. In particular, they mention the need for a better comparison with other studies while highlighting the advances of this study (ref#2), a discussion about the observed large amplitude LAB phases (ref#3) as well as the limitations of the methodology (ref#3) and in general more information about the used models (ref#1).