A Cascadia Slab Model From Receiver Functions

We map the characteristic signature of the subducting Juan de Fuca and Gorda plates along the entire Cascadia forearc from northern Vancouver Island, Canada, to Cape Mendocino in northern California, USA, using teleseismic receiver functions. The subducting oceanic crustal complex, possibly including subcreted material, is parameterized by three horizons capable of generating mode‐converted waves: a negative velocity contrast at the top of a low velocity zone underlain by two horizons representing positive contrasts. The amplitude of the conversions varies likely due to differences in composition and/or fluid content. We analyzed the slab signature for 298 long‐running land seismic stations, estimated the depth of the three interfaces through inverse modeling and fitted regularized spline surfaces through the station control points to construct a margin‐wide, double‐layered slab model. Crystalline terranes that act as the static backstop appear to form the major structural barrier that controls the slab morphology. Where the backstop recedes landward beneath the Olympic Peninsula and Cape Mendocino, the slab subducts sub‐horizontally, while the seaward‐protruding and thickened Siletz terrane beneath central Oregon causes steepening of the slab. A tight bend in slab morphology south of the Olympic Peninsula coincides with the location of recurring large intermediate depth earthquakes. The top‐to‐Moho thickness of the slab generally exceeds the thickness of the oceanic crust by 2–12 km, suggesting thickening of the slab or underplating of slab material to the overriding North American plate.


10.1029/2023GC011088
2 of 22 subduction system transforms into the right-lateral Queen Charlotte and Revere-Dellwood Faults and to the south into the San Andreas Fault (Figure 1).
Farther downdip, the JdF has been identified below the Salish Sea on marine seismic sounding transects through the Juan de Fuca Strait and Georgia Strait.At about 20 km depth, the sharp <2 km thick reflector that marks the top of the slab widens into an up to 10 km wide reflection band, the so-called E-layer (e.g., Clowes et al., 1987a;Nedimović et al., 2003), that extends to depths of at least ∼50 km (Calvert et al., 2006).A similarly thick reflective zone has been identified atop the subducting JdF at 35-40 km depth beneath central Oregon (Keach et al., 1989;Tréhu et al., 1994).It has been argued that the E-layer represents the transition into a wider shear zone that creeps aseismically and hosts episodic tremor and slip (ETS, see e.g., Calvert et al., 2020;Nedimović et al., 2009).
In Cascadia, the subduction zone stratigraphy has also been previously characterized using teleseismic P-wave receiver function data (e.g., Abers et al., 2009;Bostock et al., 2002;Cassidy & Ellis, 1993;Langston & Blum, 1977;Mann et al., 2019;McGary et al., 2014;Nabelek et al., 1993;Nicholson et al., 2005).A recent study employing receiver functions, local tomography and seismic reflection data in southern Vancouver Island suggests that the oceanic crust may reside below the E-layer (Calvert et al., 2020) and that at least part of the E-layer comprises an ultralow S-wave velocity zone (ULVZ), with V P /V S of the order of 2-3 (Audet et al., 2009;Cassidy & Ellis, 1993).In local seismic tomograms, the slab stratigraphy oftentimes appears smeared into a single layer with moderately elevated V P /V S in the order of 1.8-2.0,consistent with basaltic or gabbroic lithologies with some contribution of fluid-filled pores.Interpretation of the oceanic Moho in tomographic models is less ambiguous, where it appears as a strong negative V P /V S gradient to values lower than 1.7 that mark the oceanic mantle below (Guo et al., 2021;Merrill et al., 2020Merrill et al., , 2022;;Savard et al., 2018).
An initial margin-wide map of the top of the JdF was constructed from a mixed data set of earthquake hypocenters, active source seismic profiles, receiver functions and local earthquake tomograms with the aim to model interseismic strain accumulation in the overriding plate (Flück et al., 1997).With increasing data availability over time and a better understanding of subduction processes, the initial model has been updated and extended in space using additional constraints from seafloor magnetic anomalies, deeper seismicity and diffraction of strong earthquake first arrivals (McCrory et al., 2004) and later from relocated earthquake hypocenters and electrical conductivity profiles (Hayes et al., 2018;McCrory et al., 2012).Other slab models are based purely on receiver functions (Audet et al., 2010;Hansen et al., 2012).Despite a broad agreement in recovered slab depths to within ∼10 km, considerable differences exist across these models.These differences are associated with data uncertainties, the fact that the slab models are based on different data types, and with ambiguities in the interpretation of proxies for what constitutes the "slab top" (McCrory et al., 2012).
Here, we construct a margin-wide slab model that honors an oceanic crustal stratigraphy that may consist of up to two layers, including the possibility of subcreted material.Our model is based on the observation that receiver function images of the slab exhibit characteristic successions of positively and negatively polarized conversions that can be explained by interfering forward-and back-scattered seismic wave modes originating at three interfaces.We map these interfaces continuously along dip from the coast to the forearc lowlands (Salish Sea, Willamette Valley) and along strike from Brooks Peninsula on northern Vancouver Island, Canada, to Cape Mendocino, USA (Figure 1).Our results demonstrate how the overall slab morphology is controlled by the Figure 1.The tectonic setting of the Cascadia subduction zone and station distribution employed to determine the slab geometry under the forearc.The convergence of the Juan de Fuca and Gorda Plates relative to stable North America is shown as arrows (DeMets et al., 1994).Terrane boundaries modified after Watt and Brothers (2020).Top inset: Location of the study area on the North American continent.Bottom inset: Earthquake source distribution form 30° to 100° epicentral distance used to compute receiver functions.location of the static backstop.A subduction stratigraphy that is generally thicker than the incoming oceanic crust is testament to complex deformation processes affecting slab morphology along the subduction trajectory.

Data and Methods
A total of 45,601 individual receiver functions recorded at 298 seismic stations distributed across the Cascadia forearc contributed to the slab model.For each station, 100 s recordings symmetric about the P-wave arrival (i.e., 50 s noise and 50 s signal, for convenience) of earthquakes with magnitudes between 5.5 and 8, in the distance range between 30° and 100°, were downloaded (Figure 1).Waveforms with a signal-to-noise ratio smaller than 5 dB on the vertical component or 0 dB on the radial component were excluded.The instrument responses were removed and the seismograms were transformed into the upgoing P-SH-SV modes (Kennett, 1991).The P-component was trimmed to the time window beyond which the envelope fell below 2% of the maximum amplitude and a cosine taper was applied.

Receiver Function Processing
The three-component P-wave spectra were scaled by their signal-to-noise ratio and binned according to their incidence angle in back azimuth bins of 7.5° and horizontal slowness bins of 0.002 s km −1 .Within each bin, radial and transverse receiver functions were computed through frequency domain simultaneous deconvolution (Gurrola et al., 1995), with an optimal damping factor found through generalized cross validation (Bostock, 1998).This approach mitigates the instabilities inherent in spectral division by stacking spectra prior to deconvolution.

Parameter Search
The continental forearc and subducting slab were parameterized as three layers over a mantle half-space, with the subduction stratigraphy bounding interfaces labeled as t (top), c (central), and m (Moho) (Figure 2).Synthetic receiver functions were calculated through ray-theoretical modeling of plane-wave scattering at the model interfaces (Bloch & Audet, 2023;Frederiksen & Bostock, 2000; Figure 2b).The thickness, S-wave velocity (V S ) and P-to S-wave velocity ratio (V P /V S ) of each layer, as well as the common strike and dip of the bottom two layers and the top of the half space (in total 11 parameters), were optimized simultaneously through a simulated annealing global parameter search scheme (Xiang et al., 1997), as implemented in the SciPy package (Virtanen  et al., 2020).In analogy to the annealing process in metallurgy, the scheme samples the misfit function stochastically under a decreasing "temperature" that gradually favors low-misfit parameter combinations.In this way, the algorithm can escape local minima in the misfit function.It has proven efficient in converging toward the global minimum in problems with many independent variables (Kirkpatrick et al., 1983).The misfit was defined as the anti-correlation (1 minus the cross correlation coefficient) between the observed and predicted receiver functions, bandpass filtered between 2 and 20 s period duration.
Initial thickness bounds for the continental crust (Figure 2c) were based on the slab model of Audet et al. (2010) (±10 km).Maximum Layer 1 thickness was constrained by the maximum E-layer thickness of 10 km (Nedimović et al., 2003) and maximum Layer 2 thickness with the thickness of the incoming oceanic crust of 6.5 km (Han et al., 2016).Layer 1 could attain zero-thickness if the E-layer were absent.Because the igneous oceanic crust may be part of the E-layer, Layers 1 and 2 were constrained to have a combined minimum thickness of 6 km.Velocity bounds (Figure 2c) for the continental crust and Layer 2 were based on the 2σ interval of the expected lithologies for continental and oceanic crust, respectively, from the seismic velocity database of Christensen (1996); and for Layer 1 on an analytic poro-elastic model (Bloch et al., 2018) constrained to match the V P /V S observations of the ULVZ (Audet et al., 2009).
To verify convergence toward a global minimum, the global parameter search was initialized with at least three different random number seeds, which affect the distribution from which trial parameter estimates are drawn.The resulting data predictions and models were checked for consistency with neighboring stations, previous tomographic profiles (Guo et al., 2021;Kan et al., 2023;Merrill et al., 2020;Savard et al., 2018), hypocentral locations of low-frequency earthquakes (LFEs) within tremor (Armbruster et al., 2014;Plourde et al., 2015;Royer & Bostock, 2014;Savard et al., 2018Savard et al., , 2020) ) and offshore marine seismic profiles (Suzanne Carbotte, pers. comm;Carbotte et al., 2023).If none of the minimum misfit models of an individual station were consistent with the above constraints, the global search was repeated within narrower bounds around a preferred solution from a neighboring reliable station.Such a model was used only in case it converged toward values away from thickness bounds (Figure 3).For each of the three horizons, a quality and a nominal depth uncertainty were assigned.Quality A denotes a horizon where at least one back-scattered phase in the predicted data correlates with the observed data (Figures 3a and 3b), the predicted data are consistent among neighboring stations and the modeled horizon depth is consistent with the available external constraints.A quality B horizon shows a good phase correlation, but the predicted data are inconsistent with neighboring stations and/or the modeled depth is inconsistent with external constraints.Quality C was assigned to horizons that do not show a convincing correlation between observed and predicted data, usually due to data with low signal-to-noise levels.Stations above the forearc lowlands for which the characteristic slab signature (Figure 2b) is decisively absent and where the onset of eclogitization is expected were marked with a quality X.The nominal depth uncertainty was estimated from the scatter of the local minima in the vicinity of the preferred minimum, as determined in the global search (Figure 3c).

Fitting of Interfaces
In total, 171, 143 and 137 quality A nodes were determined to constrain the t, c and m interfaces, respectively.At the trench, 105 nodes at 3 km below the local bathymetry were inserted to constrain the t and c interfaces, and at 6.5 km deeper to constrain the m interface, representing typical sediment and igneous crustal thicknesses (Han et al., 2016).A spline surface (Sandwell, 1987) was fitted to these nodes to yield margin-wide depth models.The spline coefficients were found using singular value decomposition (Aster et al., 2018;Wessel & Becker, 2008), with the nominal depth uncertainties supplied as weights.The solution was damped by retaining the 116, 117, and 116 largest singular values for the t, c, and m interfaces, respectively, based on the analysis of L-curves and the Akaike information criterion (Figure S1 in Supporting Information S1).

Margin-Scale Slab Morphology
The signature of the subduction stratigraphy can be traced along the forearc from Brooks Peninsula on northern Vancouver Island, across Vancouver Island, the Olympic Peninsula, the Willamette Valley in Washington and Oregon, into the Klamath Mountains and to Cape Mendocino in northern California (Figures S2-S50 in Supporting Information S1).The recovered velocities of the three model layers are consistent for neighboring stations (Figure S51 in Supporting Information S1).Slab morphology suggests a division into four segments: the Klamath, Central, Olympic and Vancouver Island segments (Figure 4).
The Central segment, between 44° and 47°N, reveals the steepest dip, between 10° and 20°, and overall deepest slab, with the t horizon located between 15 and 25 km depth along the coast and dipping to 35-45 km depth before losing expression in advance of the volcanic arc.The Central segment is flanked to the north and south by flatter segments.In the south, the Klamath segment, located between ∼40° and 44°N, displays a more shallowly dipping slab, a contorted t horizon beneath Cape Mendocino and a contorted m horizon along the landward projection of the Blanco Fracture zone.The Olympic segment, located between 47° and 49°N, exhibits a shallow dipping (0-5°) slab beneath the coastal region, and is delimited to the south by a steep downward bend in the t and m horizons near Grays Harbor and by a bend in the slab strike just north of the Juan de Fuca Strait.Along the dip, the slab steepens as it approaches Puget Sound, where it begins to lose expression (Abers et al., 2009).The northernmost Vancouver Island segment is characterized by a moderately dipping slab.Near the northern terminus of subduction, north of Nootka Island, the t and c conversions appear disturbed.In summary, from north to south, the slab (a) dips gently and steepens down dip under Vancouver Island, (b) dips shallowly beneath the Olympic Peninsula, (c) steepens significantly beneath the Oregon Coastal Mountains, (d) subducts in a step-like fashion in front of Klamath Mountains, and (e) becomes contorted in the Cape Mendocino area.A comparison with previous slab models is shown in Figure S52 in Supporting Information S1.

Central Segment
Across the Central segment, the slab has been imaged with seismological methods using data from the CASC'93 experiment that comprised a temporary broadband array of ∼30 stations deployed across the Oregon  forearc (Nabelek et al., 1993).It yielded the first dense receiver function studies targeting the subduction zone structure that clearly revealed subducting oceanic crust (Bostock et al., 2002;Rondenay et al., 2001;Tauzin et al., 2016).The comparison of our model with the teleseismic full-waveform tomogram of Kan et al. (2023) yields a consistent picture of the subduction stratigraphy (Figure 5).As in previous studies, Kan et al. (2023) image the subducting Juan de Fuca plate as a distinctive low-V S zone, which attains velocities as low as 3.3 km s −1 .All three horizons parallel this structure, with t marking the top of the LVZ and c and m marking two steps in the gradual increase toward high V S , characteristic for oceanic mantle of the order of 4.3 km s −1 .This structure has a very clear and characteristic expression in the receiver functions, which weakens near station XZ.A18, beneath the Willamette Valley, as in the tomogram.The entire stratigraphic sequence  (t, c, m) brackets weak slab-related seismicity in the offshore area (Morton et al., 2023).It has a thickness of about 7 km near the coast and thickens arc-ward to about 13 km, with the two layers possessing comparable thickness.

Klamath Segment
Beneath the Mendocino region, the subduction stratigraphy has been imaged as a moderately high-V P /V S zone (1.8-1.9,Guo et al., 2021) complemented by relatively abundant intraslab seismicity defining a tightly confined Wadati-Benioff zone (e.g., Waldhauser, 2009;Wang & Rogers, 1994; Figure 6).The t and m horizons encapsulate the seismically active, moderately high-V P /V S zone, with the m horizon falling in good agreement with the V P / V S = 1.7 contour.Where it projects beneath the Franciscan terrane, the high-V P /V S -zone loses expression and the density of earthquakes diminishes (60 km from coast in Figure 6a).Our slab model here indicates a generally shallower dip that steepens again under the Klamath terrane (100 km from the coast), where it indicates that a low-V P /V S anomaly is located within the subduction stratigraphy.Layer 1 is absent between the coast and the Franciscan terrane and attains a thickness of a few kilometers farther landward.Notably, no seismicity is located within Layer 1.The c horizon defining the base of Layer 1 approximately aligns with the location of LFEs (Plourde et al., 2015).The entire subduction stratigraphy has a fairly uniform thickness of 10 km.The receiver function slab signature is difficult to correlate laterally, presumably due to some combination of variation in overburden and slab properties (Figure 6b).

Olympic Segment
A profile along the dip from the western end of the Olympic Peninsula, across the Juan de Fuca Strait, southern Vancouver Island, and into the Strait of Georgia reveals a flat lying slab beneath the Olympic Peninsula that continues under the Juan de Fuca Strait and gradually steepens under southern Vancouver Island (Figure 7a).The t and m horizons encompass the moderately high-V P /V S zones previously interpreted as the subducting crust in local seismic tomograms (Merrill et al., 2020;Savard et al., 2018).Under the Olympic Peninsula, this zone is seismically active and m agrees well with the V P /V S = 1.7 contour.Beneath southern Vancouver Island, m bounds the top of seismic activity previously interpreted to occur within the subducting mantle (Savard et al., 2018).Layer 1 is absent or very thin beneath the Olympic Peninsula and attains a thickness of about 5 km beneath southern Vancouver Island, where it is aseismic.The c horizon is located 2-3 km above a prominent band of LFE locations (Armbruster et al., 2014;Savard et al., 2018).Tremor hypocenters from Bombardier et al. (2023), see also Kao et al. (2005)) scatter within and above the subduction stratigraphy.The complex overburden structure of the Olympic Peninsula hampers clear identification of c and m; however, correlations of seismic phases along strike and along dip yield a laterally coherent picture.Beneath southern Vancouver Island, the slab reveals a clear and simple receiver function signature that can be traced beneath the Gulf Islands in the Strait of Georgia and loses expression toward the British Columbia Lower Mainland (Figure 7b).

Vancouver Island Segment
The Vancouver Island segment exhibits t and m horizons that bracket NE-dipping regions of elevated V P /V S evident in local seismic tomograms.The m horizon coincides with the V P /V S = 1.7 contour, that also bounds the top of seismicity which has been inferred to reside in the oceanic mantle (Figure 8a and Figures S3-S16 in Supporting Information S1, Merrill et al., 2022;Savard et al., 2018).c can best be seen as a pronounced and distinct horizon in southern and south-central Vancouver Island, where it lies 2-4 km underneath t and decisively above LFE locations (Savard et al., 2018).Toward north-central Vancouver Island, the subduction stratigraphy appears to thicken substantially downdip, from ∼8 km near the coast, to ∼16 km inland.Layer 1 and Layer 2 contribute in equal part to the combined thickness.The c horizon generally follows the LFE locations (Savard et al., 2020).Substantial scatter in the station measurements and difficulties in reconciling phase correlations across closely spaced stations attest to the complex subsurface structures that are also evident in the local seismic tomography and may be related to the subduction of the Nootka Fault Zone as the northern terminus of JdF subduction (Figure 8b, Merrill et al., 2022;Savard et al., 2018).

Interpretation of the Subduction Stratigraphy
The combined thickness of the stratigraphic package comprising t, c, and m horizons exceeds almost everywhere the nominal thickness of the incoming oceanic crust of ∼6.5 km by 2-12 km (Figure 9a).A thickness of ∼7 km is resolved only along the southern Central segment, between ∼43° and 44°N.Model regularization may dampen slab complexity and smooth over interface steps on a ∼20 km scale (e.g., m in Figure 6a), but the excessive thickness of the slab stratigraphy is a robust feature of the model and is almost always underpinned by individual point station measurements.Additional material, other than actively subducting igneous oceanic crust, must therefore make up the subduction stratigraphy.
Layer 2 and the underlying mantle half-space, separated at m, were designed to correspond to igneous oceanic crust and pristine mantle.Where seismic velocities and seismicity images are available, the model appears to have captured this contrast appropriately, so that we confidently interpret m as the oceanic Moho.We cannot exclude the possibility that, where the plate is hydrated, m is biased into the oceanic mantle, lying deeper than the Moho.Signs of mantle hydration may be present under the Cape Mendocino coast and offshore northern Vancouver Island, suggested by a diffuse tomographic Moho, abundant mantle seismicity and the subduction of major fracture zones (Figures 6 and 8, e.g., Chaytor et al., 2004;Merrill et al., 2022;Rohr et al., 2018;Wilson, 1989).Such signatures are, however, not universally present.
The excess thickness is more likely to develop above the plane of active subduction, that is, in Layer 1.Where the thickness of Layer 1 is substantial (i.e., from t to c; Figure 9b), the E-layer (or a reflective zone above the slab) has been detected in reflection seismic surveys (Figure 9b, Clowes et al., 1987a;Keach et al., 1989;Nedimović et al., 2003;Tréhu et al., 1994).Nedimović et al. (2003) suggest that the emergence of the E-layer is related to the occurrence of ETS.The E-layer is typically thicker than Layer 1, which suggests that Layer 1 is part of the E-layer (Calvert et al., 2020).Within the tremor zone, defined by the 0.1 tremor yr −1 km −2 contour (Figures 9b-9d; downloaded from https://pnsn.org;Wech, 2010), the mean and median V P /V S in Layer 1 are 2.49 ± 0.14 (2σ) and 2.44.Outside the tremor zone V P /V S is lower, with a mean value of 2.28 ± 0.14 and median value of 1.95 (Figures 9b,10a,and 10b).A two-sample Kolmogorov-Smirnov test yields a p-value of 5 ⋅ 10 −5 , indicating that the distribution of V P /V S values of Layer 1 from inside and outside the tremor zone are statistically different with >99% confidence.This suggests that the development of Layer 1 as a high-V P / V S ULVZ is related to tremor, in agreement with previous findings (Audet et al., 2009;Song et al., 2009).We interpret t in the tremor zone as the top of this ULVZ.Projecting the tremor epicenters (Wech, 2010) onto the t and c horizons yields tremor depths of 32 ± 10.8 and 38 ± 10.2 km (2σ), respectively (Figures 10c and 10d).Tremor depths are concentrated more tightly when projected to the c horizon, suggesting that tremor occurs closer to the base of Layer 1 (Figure 9c).Inside the tremor zone, where Layer 1 corresponds to the ULVZ, c marks a stark material contrast against the underlying oceanic crust and we interpret c as the base of the ULVZ.
Between the coast and the tremor zone, except between 44° and 45°N, Layer 1 is typically thinner (Figure 9b) and its V P /V S is lower (Figures 9d and 10b), attaining normal values for basaltic material (∼1.8).Layers 1 and 2 still exhibit a combined thickness in excess of the incoming oceanic crust, with Layer 1 displaying properties that are nevertheless similar to oceanic crust.The t horizon is here the top of this excess volume.The c horizon here usually marks a less prominent material contrast than inside the tremor zone.It may seem natural to interpret c as the base of a possible sedimentary blanket above an underlying igneous oceanic crust (e.g., Delph et al., 2018), but we note here that Layer 2 is frequently thicker than oceanic crust, hence the interpretation of c as the base of sediments is possible but not universal.Horizon c may alternatively represent a velocity gradient within a sedimentary layer or the base of altered material belonging to the overriding continental crust.O" mark places where sediment subduction has been detected on marine seismic surveys, "X" where sediment subduction is absent (Han et al., 2016;Tréhu et al., 2012).The thickness of the subduction stratigraphy exceeds the thickness of the igneous oceanic crust.(b) Layer 1 (t-to-c) thickness and tremor zone (Wech, 2010).Downdip thickening of Layer 1 correlates with tremor locations.(c) Depth to c horizon correlates closely with tremor occurrence (Figures 10c and 10d).(d) V P /V S of Layer 1.

Possible Controls on Slab Morphology
The overall slab morphology exhibits a first-order correlation with the location of the static backstops in the Cascadia subduction system (Figure 11, Watt & Brothers, 2020).These are kinematic discontinuities that are related to distinct strength contrasts within the continental crust formed by accreted crystalline terranes.The most important terrane is Siletzia, a basaltic large igneous province that formed offshore as an oceanic plateau, possibly related to magmatism of the Yellowstone hotspot.It can be mapped along coastal Oregon, Washington and British Columbia (Wells et al., 2014).An associated aeromagnetic anomaly indicates that Siletzia is most voluminous in central and northern Oregon (Wells et al., 1998).Reflection seismic together with wide-angle seismic data and geomorphologic markers reveal that the base of Siletzia is up to 35 km deep and possibly extends down to the plate interface (Tréhu et al., 1994).This inference is substantiated by magneto-telluric data that image a voluminous resistive body, interpreted as representing Siletzia, that meets the plate interface in coastal Oregon (Egbert et al., 2022).Kinematically, the thickened Siletz terrane forms a distinct block that rotates clockwise with respect to stable North America (Wells et al., 1998) and displays the lowest interseismic vertical uplift rates along the entire forearc (Mitchell et al., 1994).Where the Siletz terrane recedes far inland on the eastern side of the Olympic Peninsula, giving way to the Olympic complex formed by underthrust marine sediments (e.g., Brandon & Calderwood, 1990), the slab lies shallower and flatter than anywhere else along the entire onshore forearc.Conversely, where the western boundary of the Siletz terrane is located off-shore and Siletzia is thickest, the slab is deepest and has its steepest dip (Figure 11).This suggests that the competence and rigidity of the Siletz block force the descent of the Juan de Fuca slab.It has been suggested that the Kumano pluton influences the subducting Philippine Sea Plate in a similar manner below southwest Japan (Arnulf et al., 2022).
In between the shallowly dipping Olympic and steeply dipping Central segments, a pronounced southward downward bend in the slab is evident along a line extending between Grays Harbor and the southern end of Puget Sound.The bend can be seen in the raw receiver functions, where the timing of the P m S conversion increases, for example, from ∼3 to ∼4 s for rays arriving from NNW relative to those arriving from SSE azimuths at station US.NLWA and again from 4 to 4.5 s just south of that at station UW.WISH (Figure 12).Perhaps significantly, the three largest intermediate depth earthquakes in Cascadia, the 1949 M6.7 Olympia (Nuttli, 1952(Nuttli, ), 1965 M6.7 M6.7 Puget Sound (Langston & Blum, 1977), and 2001M6.8 (Ichinose et al., 2004;Kao et al., 2008) earthquakes occurred near the down-dip continuation of this bend, at depths at or immediately below those projected for the oceanic Moho.Along the Klamath segment to the south (south of 44°N), the slab structure is complex.The Gorda Plate, a relatively young and highly deformed plate (Chaytor et al., 2004;Wilson, 1989), encounters two static backstops, namely the western boundaries of the Franciscan complex and the Klamath terrane (Figure 11, Clarke, 1992;Watt & Brothers, 2020), and is bounded to the south by the Mendocino Fracture Zone.The southern and eastward trending seaward boundary of the thickened Siletz terrane has a reduced impact on slab morphology, resulting in the southward transition to a more gently dipping slab (Figure 11).This geometry is interrupted with the emergence of the Klamath terrane, where the steepest dip of the slab is located near the coast, bending behind the first (seaward) backstop, and unbending beneath the second (landward) backstop.In the Cape Mendocino area, the slab top is contorted in a fashion that yields a flat-lying segment just behind the coast.Because of the generally lower dip in advance of the volcanic arc and its unbending beneath the southern Siletz and Klamath terranes, it appears as if the Gorda plate does not subduct as readily as the Juan de Fuca plate.A possible cause for this behavior is increased buoyancy of the youngest subducting lithosphere (5-6 Ma at the trench, e.g., Wilson, 1993).

Excess Thickness of Subduction Stratigraphy
The nature and origin of the E-layer as a prominent element of the subduction zone stratigraphy that emerges abruptly along the dip trajectory in the vicinity of southern Vancouver Island is a long-standing conundrum in the understanding of the Cascadia subduction zone (e.g., Calvert, 1996;Calvert et al., 2011Calvert et al., , 2020;;Clowes et al., 1987a;Nedimović et al., 2003).Our data show a qualitative correlation between a thick Layer 1 and a thick (>4 km) E-layer where the latter has been imaged (Figure 9b).We also suggest that the reflective zone mapped by COCORP in central Oregon (Figure 9b, Keach et al., 1989) may manifest the presence of a structure with similar origin since it also coincides with a thick Layer 1. Assuming this association holds true along the entire margin, our data suggest that the E-layer is ubiquitous.Its abrupt emergence along dip is likewise reflected in our data: Coastal stations have a tendency to exhibit a thin or absent Layer 1, whereas inland stations generally possess a thick one (Figure 9b), consistent with previous inferences of Layer 1 thickening near the coast line from an amphibious receiver function study (Audet & Schaeffer, 2018).Interestingly, the combined (Layer 1 + Layer 2) thickness of the subduction stratigraphy does not obey the same trend.Places of a thin or absent Layer 1 may have an overall thick subduction stratigraphy coastal Olympic Peninsula and Cape Mendocino) and a significantly thick Layer 1 may correspond to a subduction stratigraphy that does not much exceed the thickness of the incoming oceanic crust (e.g., ∼7 km thickness between 43° and 44°N; Figure 13a).
Sediments entering the subduction system may contribute to the subduction stratigraphy (e.g., Delph et al., 2018), but information about the amount of subducting sediment at the time of writing is scarce.Tréhu et al. (2012) interpret sediments subducting beneath Siletzia on two seismic lines near 45°N (circles on Figure 9a), but not on a third line closer to 44°N (cross on Figure 9a).Within the same latitude interval, the characteristic transition from thickened to normal subduction stratigraphy occurs, suggesting that these subducting sediments make up for the extra thickness (Figure 13b).In contrast, Han et al. (2016) document no sediment subduction at the latitude of the Juan de Fuca Strait, where we image thick (∼11 km) subduction stratigraphy.However, it is possible that sediment subduction occurred at the trench in the latter region at 3 Ma ago, and subsequently ceased.More data are required to conclusively define where sediment subduction contributes to subduction stratigraphy thickness.
We note that Layer 1 emerges at around 30 km depth and gains thickness along the subduction trajectory, and that this thickness is unrelated to the thickness of the subduction stratigraphy updip of this depth (Figures 9a and 9b).This observation suggests that Layer 1 thickens in situ and develops a ULVZ through some depth-activated process.Elevated V P /V S > 2.0 (Figures 9b and 9d) suggest that the medium is fractured and saturated with pressurized fluids (Christensen, 1984), implying that it has lost structural integrity and strength.As a weak zone, the ULVZ is likely to host slip (e.g., Luo & Liu, 2021;Wech & Creager, 2011).LFE hypocenters are located near the base of the ULVZ (Figure 14, Calvert et al., 2020), suggesting that the plate boundary is located near c.Excess thickness may be due to underplating of subducting material, either of sediments atop the oceanic crust (e.g., Delph et al., 2018) or of the upper basaltic crust, which may lose structural integrity through wear (Figure 13c, see also e.g., Calvert et al., 2020;Clowes et al., 1987a).Moderately high seismic velocities (V S > 3.2 and V P /V S < 1.9) indicated by our inverse modeling results for Layer 2 (Figure 2 and Figure S51 in Supporting Information S1) preclude the presence of pervasive fracturing and pressurized fluids (Christensen, 1984).Instead, the presence of slivers of oceanic crust, large enough to not reduce seismic velocities significantly, would be consistent with LFE occurrence inside Layer 2. The subordinate slip represented by the LFEs during ETS episodes is consistent with the process of initiating detachment of the subducting oceanic crust at the LFE horizon (Figure 13c).Slow slip, which makes up the main share of the slip budget at depth (Bostock et al., 2015;Dragert et al., 2001Dragert et al., , 2004;;Kao et al., 2010), may well be located at or above c, that is, at the base or inside the 4-10 km thick ULVZ.
Subcretion and underplating are consistent with earlier inferences made for the onshore Cascadia forearc from a wealth of geophysical data.Calvert et al. (2011) interpret underplating of sediments as taking place south of Puget Sound.Calvert (2004) and Clowes et al. (1987a) inferred that the E-layer constitutes underplated metabasaltic material beneath southern Vancouver Island based on high seismic velocities.The correspondence between these inferred sites of underplating with the thick ULVZ detected here and the widespread distribution of the ULVZ suggest that underplating is occurring through the majority of the entire Cascadia forearc (e.g., Delph et al., 2021).

Conclusion
Receiver functions provide valuable insights into the subduction of the Juan de Fuca and Gorda plates in the Cascadia region.Based on previous studies of receiver-side forward and back-scattered mode conversions, we parameterize subduction stratigraphy in three horizons t, c, and m.Mapping these horizons across the forearc reveals flatter slab segments beneath the Olympic Peninsula and Cape Mendocino, central Oregon exhibits a steeply dipping slab.Below most of Vancouver Island, the slab is marked by modest dips (∼7°-12°).This slab morphology appears to be influenced by the mechanical strength and density of accreted crystalline terranes.A notable Moho step south of the Olympic Peninsula may relate to recurrent, large, intermediate-depth earthquakes beneath Puget Sound.In addition, the presence of a thick topmost layer in the subduction stratigraphy may indicate the widespread occurrence of the E-layer.Previous interpretations suggest that the E-layer represents underplated slab material, including both sediments and metabasalt, implying that underplating occurs through most of the Cascadia forearc.

Figure 2 .
Figure 2. (a) Forearc stratigraphy with the previously identified interfaces.(b) Schematic radial receiver function with the forward and back scattered mode conversions used to constrain the model.Phases may interfere and cancel out in some cases.The absence of specific phase combinations may therefore be meaningful.Upper case letters indicate upgoing rays, lower case letters down-going rays, and subscript the scattering interface.(c) Parameterization of the subsurface model.The possible presence of additional interfaces complicates the phase associations.

Figure 3 .
Figure 3. Global search for subsurface parameters.(a) Receiver function data for station C8.TWBB.(b) Predicted data from the best fitting model with phase labels as in Figure 2b.(c) Local minima encountered in the global search for the 11 subsurface parameters (thickness against V P /V S in the left column and against V S in the right column) using a simulated annealing scheme with the preferred solution marked with a green circle and nominal depth uncertainties with a gray bar.Note the presence of a local minimum.If such minimum proved more consistent with external constraints and neighboring stations, the global search was repeated within bounds around that minimum.

Figure 4 .
Figure 4. Depth to the t, c, and m horizons.Top row: Data points by quality (black frames: A; white frames: B; not used for fitting the interface; and grayed out: C).Stations marked X do not show the respective interface and are interpreted as the location of the eclogitization front.Bottom row: modeled interfaces (Section 2.3) and profile locations (Figures 5-8).

Figure 5 .
Figure 5. (a) Profile A (Figure 4) with slab model and control points superimposed on the V S model of Kan et al. (2023) with seismicity from Morton et al. (2023).Comparison with the V P /V S image is shown in Figure S53 in Supporting Information S1.(b) Receiver function sections of individual stations sorted along the profile, with the receiver functions within each section sorted by the angular distance of the ray back azimuth from the profile azimuth (90°).1.5-20 s bandpass filter was applied.Phase labels correspond as in Figure 2.

Figure 6 .
Figure 6.Red octagon marks the intersection of the profile with the terrane backstop.

Figure 8 .
Figure 8.As Figure 5, but for profile D across Northern Vancouver Island.Tomogram and seismicity from Merrill et al. (2022).Low-frequency earthquake locations from Savard et al. (2020).Comparison with V S is shown in Figure S56 in Supporting Information S1.Receiver functions filtered between 2 and 20 s.

Figure 9 .
Figure 9. Select properties of slab stratigraphy.(a) Combined Layers 1 and 2 (t-to-m) thickness."O" mark places where sediment subduction has been detected on marine seismic surveys, "X" where sediment subduction is absent(Han et al., 2016;Tréhu et al., 2012).The thickness of the subduction stratigraphy exceeds the thickness of the igneous oceanic crust.(b) Layer 1 (t-to-c) thickness and tremor zone(Wech, 2010).Downdip thickening of Layer 1 correlates with tremor locations.(c) Depth to c horizon correlates closely with tremor occurrence (Figures10c and 10d).(d) V P /V S of Layer 1.

Figure 10 .
Figure 10.Properties of Layer 1 in relation to tremor (a, b) V P /V S of Layer 1 at stations (a) inside and (b) outside the tremor zone (0.1 tremor km −2 yr −1 contour; Figure 9) (c, d).Depth distribution of tremor epicenters projected onto the (c) t and (d) c horizons.Numbers at the base indicate the 5%, 50%, and 95% quantiles of the depth distribution.

Figure 11 .
Figure 11.(a) Dip and (b) depth of the t horizon.Static backstop (line with red octagons in panel (a)) and terrane boundaries (thick lines in panel (b)) modified after Watt and Brothers (2020).Shaded area enclosed by white dashed line represents thickened Siletz terrane detected in aeromagnetic data (after Wells et al. (1998)).The location of the terrane backstop correlates with and may exert a first order control on slab morphology.

Figure 12 .
Figure 12.Downwarped Moho from the Olympic Peninsula to Grays Harbor.(a) Map view with Moho depth contours as well as locations and receiver function ray back azimuths of stations shown on the right.Earthquake locations and focal mechanisms from International Seismological Centre (2023).(b) Radial receiver functions sorted by back azimuth.Rays arriving from NNW colored blue and from SSW colored gold.Note the southward down Moho-steps (P m S) at stations coincident with a thickening low velocity zone above at stations NLWA, WISH, WHGC, and RADR.

Figure 13 .
Figure 13.Possible subduction stratigraphies present in the Cascadia subduction zone.(a) Subduction of undisturbed oceanic crust (e.g., Central-South Oregon).(b) Sediment subduction, c may represent the base of the sedimentary later or a horizon within the sediments (e.g., Olympic Peninsula, Northern Oregon).(c) E-layer on top of the subducting crust.Low-frequency earthquake locations may indicate a detachment horizon at or below the base of the ultralow S-wave velocity zone.Low seismic velocities and in situ thickening above suggest ongoing underplating (e.g., Southern Vancouver Island).

Figure 14 .
Figure 14.Histograms of the depth of t, c, and m horizons relative to low-frequency earthquake (LFE) locations for different regions.Bin width is 2 km.LFEs are most closely located to the c horizon.For Vancouver Island, the data indicate that LFEs occur in Layer 2 between c and m.