Seismic Imaging Beneath Cascadia Shows Shallow Mantle Flow Patterns Guide Lower Mantle Upwellings

The mantle transition zone (MTZ) plays an important role in modulating material transport between the upper mantle and the lower mantle. Constraining this transport is essential for understanding geochemical reservoirs, hydration cycles, and the evolution of the Earth. Slabs and hotspots are assumed to be the dominant locations of transport. However, the degree of material transport in other areas is debated. We apply P‐to‐S receiver functions to an amphibious data set from Cascadia to image the MTZ discontinuities beneath mid‐ocean ridges, a hotspot, and a subduction zone. We find a MTZ thinned by 10 ± 6 km beneath the ridges and by 8 ± 4 km beneath the base of the slab, closely resembling the 660 discontinuity topography. Depressions on the 410 discontinuity are smaller, 5 ± 2 km on average, focused in the north and the south and accompanied by supra‐410 discontinuity melt phases. The depressions occur away from locations of uplifted 660 discontinuity, but near slow seismic velocity anomalies imaged in the upper mantle. This suggests lower mantle upwellings occur beneath ridges and beneath the base of slabs but stall in the transition zone, with upper mantle convection determining upward material transport from the transition zone. Therefore, upper mantle dynamics play a larger role in determining transfer than typically assumed.


Introduction
The mantle transition zone (MTZ) is the region between the upper mantle and lower mantle bounded by two seismic discontinuities at approximately 410 (the 410) and 660 km (the 660) depth (Dziewonski & Anderson, 1981).These are typically interpreted as the pressure-induced transformation of olivine grains into denser crystal structures, as predicted by laboratory experiments: α olivine to β-spinel (wadsleyite) at ∼410 km depth and ringwoodite to bridgmanite and magnesiowüstite at ∼660 km depth (Ringwood, 1975).The MTZ is arguably one of the most important parts of the Earth for understanding its evolution and behavior because: (a) any material moving between the upper and lower mantle has to pass through the MTZ, (b) the MTZ has the capacity for storing huge amounts of water, potentially several oceans worth (Hirschmann & Dasgupta, 2009), and (c) the top of the MTZ (and/or the shallow lower mantle) may be a region of widespread deep melts (Bercovici & Karato, 2003;Schmandt et al., 2014), which may act as geochemical reservoirs and/or alter mantle viscosity.
It is typically assumed that material transfer between the upper and lower mantle occurs where hotspots upwell vertically and where subducting plates descend vertically or sub-vertically.Seismic tomography models also find slow velocities interpreted as ascending material beneath hotspots and fast velocities associated with descending material beneath subduction zones (Nolet et al., 2007).Seismic imaging of MTZ discontinuities generally Abstract The mantle transition zone (MTZ) plays an important role in modulating material transport between the upper mantle and the lower mantle.Constraining this transport is essential for understanding geochemical reservoirs, hydration cycles, and the evolution of the Earth.Slabs and hotspots are assumed to be the dominant locations of transport.However, the degree of material transport in other areas is debated.We apply P-to-S receiver functions to an amphibious data set from Cascadia to image the MTZ discontinuities beneath mid-ocean ridges, a hotspot, and a subduction zone.We find a MTZ thinned by 10 ± 6 km beneath the ridges and by 8 ± 4 km beneath the base of the slab, closely resembling the 660 discontinuity topography.Depressions on the 410 discontinuity are smaller, 5 ± 2 km on average, focused in the north and the south and accompanied by supra-410 discontinuity melt phases.The depressions occur away from locations of uplifted 660 discontinuity, but near slow seismic velocity anomalies imaged in the upper mantle.This suggests lower mantle upwellings occur beneath ridges and beneath the base of slabs but stall in the transition zone, with upper mantle convection determining upward material transport from the transition zone.Therefore, upper mantle dynamics play a larger role in determining transfer than typically assumed.

Plain Language Summary
The mantle transition zone (MTZ), the region between the upper mantle and lower mantle, plays an important role in material transport during mantle convection.Subduction zones and hotspots are regarded as the dominant locations where material moves across the MTZ, but it is debatable in other areas.We apply P-to-S receiver functions which reflect the inner seismic discontinuities of the Earth to a seismic data set from Cascadia.We find the MTZ is thinned beneath the ridges and beneath of the base of slabs, closely resembling the 660 discontinuity topography, likely caused by upwellings from the lower mantle.Depressions on the 410 discontinuity, inferred to be upwellings to the upper mantle, are smaller and occur away from locations of upwelling indicated by the uplifted 660.However, they are near slow seismic velocity anomalies imaged in the upper mantle.Lower mantle upwellings occur beneath ridges and beneath the base of slabs but instead of ascending vertically as typical diagrams demonstrate they may stall in the transition zone and be guided by the upper mantle convection.The results suggest upwellings may take place in regions other than hotspots, and upper mantle dynamics play a larger role in determining transfer than typically assumed.
Several lines of geochemical evidence suggest greater complexity than large, long-lived full mantle convection, or in other words, a well-mixed mantle.These include isotopic differences between ridges and hotspots and bulk compositional discrepancies between the predictions and observations of Earth's interior layering and atmosphere (Hofmann, 1997).Possible explanations include pockets of enrichment, chemical piles, layering of stagnant slabs (Ballmer et al., 2015;Jellinek & Manga, 2004;Kellogg et al., 1999;Marquardt & Miyagi, 2015; J. P. Morgan & Morgan, 1999), or stable high-viscosity lower-mantle convective cells (Ballmer et al., 2017).Alternatively, the mantle may pervasively rise across the MTZ, but typically be compositionally filtered during the process (Bercovici & Karato, 2003).Although more intricate convection and alternative locations of material transfer across the MTZ may be implied by these conceptual models, tight constraints have remained elusive.
Seismic observations have also been used to argue for greater complexity than the simple vertical ascent of material beneath hotspots.Deflected upwellings or tree-like upwelling structures have been inferred from global seismic tomography images (French & Romanowicz, 2015;Tsekhmistrenko et al., 2021).Imaged slow velocity regions or MTZ thinning have been interpreted as material upwelling from the lower mantle in alternative locations away from hotspots (Agius et al., 2021;Suetsugu et al., 2006).However, it has also been suggested that some of the observations may also be related to nearby plume material (Portner et al., 2017;Zhao et al., 2013) or be artifacts of tomographic inversion (Bezada et al., 2016).
Several studies have found seismic phases interpreted as conversions from a melt layer just above the MTZ.This is thought to be the result of hydrated material that upwells through the MTZ, but melts when it enters the lower water solubility conditions shallower than the MTZ.However, the phases have been used to link the presence of a melt layer to a variety of different tectonic environments (Hier-Majumder et al., 2021;Revenaugh & Sipkin, 1994;Song & Kim, 2011;Vinnik & Farra, 2007) and also to suggest the presence/absence of the melt layer is sporadic without necessarily having a tectonic correlation (Tauzin et al., 2010;Wei & Shearer, 2017).
High-resolution in situ imaging of MTZ discontinuities provides independent and tight constraints on material transfer and melt layering.Yet, this imaging has proven challenging given that most seismic stations are located on land, that is, limited to ∼30% of Earth's surface, and typically away from many tectonic regions of interest.
The Cascadia Subduction Zone lies in the northeastern Pacific where the Juan de Fuca and Gorda Plates subduct beneath the North American Plate (Figure 1).The Cascadia Initiative experiment deployed an amphibious seismic array from the intermediate spreading Juan de Fuca and Gorda Ridges, including the Cobb Hotspot, to the Cascadia Arc and Backarc region (Toomey et al., 2014) offering an excellent opportunity to investigate material transfer at a range of tectonic environments.Body wave tomography images high-velocity slabs descending to MTZ depths and slow velocities in the upper mantle beneath the base of the slabs (Bodmer et al., 2018).However, connections between these shallow mantle dynamics, the MTZ, and the lower mantle have yet to be investigated.Here we image the MTZ discontinuities beneath the Cascadia region using P-to-S (Ps) receiver functions (RF)s from an amphibious data set and compare the variations of these discontinuities and MTZ thickness with the seismic velocity anomalies.

Receiver Functions
In this study, we used an amphibious data set recorded by 202 stations including 136 ocean bottom seismometers (OBS)s and 66 onshore instruments.The data were collected by the Cascadia Initiative (OBSIP, 2011), USArray Transportable Array (Array, 2003), the Pacific Northwest Seismic Network (Washington, 1963), Berkeley Digital Seismic Network (Center, 2014), and the United States National Seismic Network ((ASL)/USGS, 1990).We used earthquakes of magnitudes greater than 5.5 Mw with epicentral distances between 35° and 80° that happened during the deployment period of the Cascadia Initiative experiment (Figure 1).These initial parameter cutoffs resulted in seismograms from 29,683 event-station pairs.4 and 6.Black circles along the cross sections correspond to a spacing of 100 km.The inset map shows the locations of earthquakes (red circles) used in this study.Juan de Fuca Ridge (JdFR), Blanco Fracture Zone (BFZ), Gorda Ridge (GR), and Mendocino Fracture Zone (MFZ) are as labeled.
The seismograms located on land were rotated into the P-and S-wave components using a transformation matrix (Cerveny, 2005) for the free surface: where U R and U Z are the radial and vertical components, p is the ray parameter, α and β are the P and S wave velocities, and η x is given by  √  −2 −  2 ( =  ) .Stations located on the seafloor were rotated using the transformation matrix for the solid Earth-ocean interface: where η x is given by , ρ is the density, and subscripts 0 and 1 correspond to α, β, and ρ for the water layer and seafloor layer, respectively.
Parameters used for rotation of the onshore stations were as follows: P-wave velocity = 5.5 km/s, S-wave velocity = 3.2 km/s, and density = 2,900 kg/m 3 .For OBS data, we used the sediment properties calculated from Ps sediment-crust conversion delay times for each individual station using empirical relationships (Rychert et al., 2018).The OBS data were not corrected for tilt or compliance since these corrections are applied to the frequency band of 0.005-0.03Hz (Bell et al., 2015), that is, below that of this study (0.05-0.2 Hz), and therefore they do not influence the Ps receiver function imaging of MTZ discontinuities (Agius et al., 2021).We show examples of coherent RFs with adjacent piercing points at 410 km from on-and off-shore data in Figure S1 in Supporting Information S1.
We applied a band-pass filter between 0.05 and 0.2 Hz to the waveforms and then calculated the signal-to-noise ratio (SNR).The SNR of each signal was determined by comparing the largest amplitude in a 20-s-long window after the predicted P-wave arrival time to the standard deviation of the amplitudes in a window beginning 60 s prior to the predicted P-wave arrival.For land stations, waveforms with SNRs smaller than 2 were discarded.We did not discard OBS data using the SNR approach.Instead, we performed an initial visual inspection of the OBS waveforms and performed an initial elimination.We can confirm by visual inspection that the largest amplitudes also correspond to the P-waves and the initial SNR cutoff is a valid approach.The remaining waveforms, corresponding to 21,356 event-station pairs, were deconvolved using the extended-time multi-taper frequency domain deconvolution technique (Agius et al., 2017(Agius et al., , 2021;;Chichester et al., 2018;Helffrich, 2006;Lavayssière et al., 2018;Possee et al., 2021;Rychert and Harmon, 2018;Rychert, Harmon, & Armitage, 2018;Rychert et al., 2013Rychert et al., , 2014Rychert et al., , 2021Rychert et al., , 2020Rychert et al., , 2018aRychert et al., , 2018b;;Shibutani et al., 2008).
Each RF was migrated to depth with respect to sea level, correcting for station elevation.The mantle migration model assumes velocities from the IASP91 model (Kennett & Engdahl, 1991).For crustal velocities, a modified version of the CRUST1.0model was used (Laske et al., 2013).For land stations, the sediment and crystalline crust layers from CRUST 1.0 were used.For the OBS stations, we used the sediment thicknesses and velocities derived from the Ps sediment-crust conversion delay times at each station (Rychert et al., 2018) and crystalline crust values from CRUST1.0.
We inspected each waveform and corresponding RF visually to pick reliable waveforms, eliminating unstable deconvolutions characterized by pure ringing.The P410s and P660s conversions arrived with coda and in some cases strong noise.In some cases, only one of the conversions was clear, while the other was suppressed by noise or ringing.Establishing two separate data sets for the P410s and the P660s can maximize the potential use of as many RFs as possible (Agius et al., 2017).For RFs in each data set, we required Ps phases that clearly stand out from the Moho and either the 410 or 660 discontinuities, within 40 km of the theoretical.We considered the converted phases from 410 to 660 in comparison to those at all adjacent piercing points within ∼100 km of the conversion, after migration to the depth domain, and eliminated outlier RFs with conversions at apparent depths different (>∼20 km) than the neighboring waveforms.Thereafter, we obtained 1919 RFs in total.The P410s data set consisted of 952 waveforms with 205 from OBS data, whereas the P660s data set was comprised of 1,261 waveforms with 392 from OBS data.
We found the amplitudes of RFs from OBS data were systematically larger than those from land data (Figure 2).This has been observed in previous amphibious experiments, and it is expected for typical sediment/ocean interfaces (Kim et al., 2021).We divided the RF amplitudes by the ratios of the median amplitudes of the offshore and onshore P410s and P660s data, respectively, as done in previous work, in order to meld the data sets (Kim et al., 2021).Our testing with synthetic seismograms indicate this is a reasonable way to deal with the two data sets (Supporting Information S1).Synthetic testing shows that ocean bottom receiver functions are predicted to show a different character than land receiver functions at shallow depths because of the very different shallow layering (Figure S2 in Supporting Information S1).However, this does not affect the resolution of discontinuities at MTZ depths.In addition, reverberations in ocean receiver functions from shallow discontinuities are not expected to strongly influence MTZ discontinuity resolution (Text S1 and Figure S2 in Supporting Information S1).

Error Bars, Depth Migration, and Tests
Selected RFs were migrated to depth and back-projected along the theoretical raypath and stacked onto a 3-D grid that has a lateral spacing of 1° by 1° and a 1 km depth spacing.The grid smoothing was determined by the Fresnel zone where the radius is , and λ is the wavelength and d is the depth, centered in the grid.Since we used separate data sets for P410s and P660s, two 3-D grids were generated first and then merged into one by using a linear weighing between 410-to 660-km depth of the grids.A weighted average is applied to the 410 and 660 data sets.The weights for the P410s data set grid decrease linearly from 1 to 0 from 410 to 660 km depth, while weights for the P660s data set grid increase linearly from 0 to 1 over the same depth range.There are sufficient waveforms (>5) to resolve discontinuities in all bins plotted (Figure S4 in Supporting Information S1).There are no gaps within the model, therefore sampling does not affect our interpretation.The region of overlap of resolved 410 and 660 is >10° longitude and >20° of latitude.
The reported error bars in the abstract are the standard errors of the means of each bin across subregions.In the main text, the error bars of the 410 and 660 km discontinuity depths represent the standard error of the mean of the depths of the peaks of the individual waveforms in each bin.These vary across the study area from 1 to 4 and 1 to 7 km for the 410 and the 660, respectively (Figure S5 in Supporting Information S1).The error bars on transition zone thickness reported for a given region in our study area correspond to the errors propagated from the 410 and the 660 discontinuity depths.
We performed tests using different migration models to demonstrate the robustness of our observations (Figure 3).In all tests, we left the sediment and crustal corrections intact, given that these shallow layers are well-constrained by independent observations.We report maximum differences compared to depths assuming the IASP91 model for the mantle, noting that the observed phases and depth anomaly patterns presented and discussed in the following sections remained intact.
We first tested the effect of using the 1-D mantle velocities from the PREM model (Dziewonski & Anderson, 1981) instead of the IASP91 model.The average differences in the depths of the 410 and the 660 between these two models are 1.69 and 3.88 km, respectively, and 3.39 km for the difference of the MTZ thickness.
We also tested a series of 3-D models.We tested the effect of using the regional P-wave velocity anomalies from model CAS2018_P (Bodmer et al., 2018) (extrapolating in the western and southern edges to cover our study area), relating the anomalies to P-velocity assuming an IASP91 background model and calculating S-wave velocities assuming the Vp/Vs values from IASP91.The results have a good agreement with migration using the 1-D mantle model.Ninety-five percent of pixel-to-pixel differences of the depths of the 410 and the 660 and the thickness of the MTZ were within 5, 11, and 10 km respectively.
We tested the global model PRI (Montelli et al., 2006), which has both Pand S-wave velocities.The average differences in the depths of the 410 and the 660 between these two models are 2.50 and 10.24 km, respectively, and 7.08 km for the difference of the MTZ thickness.
We also tested the global S-wave velocity model SEMum2 (French et al., 2013) assuming P-wave velocities calculated using the Vp/Vs ratio from the IASP91 model.The average differences in the depths of the 410 and the 660 between these two models are 4.40 and 9.61 km, respectively, and 5.58 km for the difference of the MTZ thickness.
Changes in depth resulting from migration model changes are sometimes larger than the reported formal error bars.However, this is a systematic error or bias.The trends are robust regardless of the migration model, and therefore do not affect our interpretation or our understanding of the Earth.
Since the overall observed patterns remain intact, for example, the depressions of the 410 in the north and south and the uplifted 660 beneath the ridges and beneath the base of the slabs, these tests show that the main interpretations of this manuscript are robust.We prefer the IASP91 migration model given that it avoids inherent additional uncertainties related to resolution in tomographic models.In addition, it avoids uncertainty related to choosing a Vp/Vs ratio when using local models that are only for Vp.

Temperature Estimates
We estimated temperature anomalies for the MTZ using the relationships of temperature with discontinuity depths applied in previous work (Agius et al., 2021).These include a +2.9 MPa/K slope for the 410 (Bina & Helffrich, 1994) and a −2.5 MPa/K slope for the 660 (Ye et al., 2014) based on experimental relationships of pressure and mineral phase transitions.The average 5 km depression of the 410 and the 15 km elevation of the 660 correspond to 60 and 190 K thermal anomalies, respectively, using the assumed Clapeyron slopes.
A wide range of Clapeyron slopes has been reported, from 1.5 to 2.9 MPa/K for the phase transition at the 410 (Akaogi et al., 1989;Bina & Helffrich, 1994) and from −4.0 to −2.0 MPa/K at the 660 (Bina & Helffrich, 1994;Ito et al., 1990).The thermal anomalies associated with the average depression and uplift reported above would be 60-115 and 138-276 K, respectively, for these slopes.However, this does not affect our interpretation of upwellings and downwellings based on the topography of the discontinuities and body wave tomography (Bodmer et al., 2018).

Results
The main phases that we image are positive phases, or velocity increases caused by the 410 and 660 discontinuities, and these are imaged across the study region.We also image negative phases at ∼ 370 km depth, strongest in discrete locations beneath the Mendocino Fracture Zone and the north-central Juan de Fuca Plate (Figures 4a  and 4c).Time-domain deconvolution tests suggest that the strongest supra-410 phases are not sidelobe artifacts, but rather robust features (Text S2 and Figure S3 in Supporting Information S1).We image strong positive phases beneath the 660 at ∼ 750 km depth east of 122°W and north of 47°N (Figure 4c).Similar phases were detected in this region by a previous RF study and interpreted as a region of high basalt content (Frazer & Park, 2021).We do not interpret these phases further, nor do we interpret the negative phases between the 660 and the positive phases at ∼750 km depth, given the proximity of those depths and the potential for the interference of sidelobes.Strong negative phases beneath the 410 are imaged but at the very southeastern corner of our study area (Figure 4a).These could potentially be related to the edge of the slab; although, dipping structures can cause complex receiver function imaging, making definitive interpretation challenging (Lekić & Fischer, 2017).Therefore, we do not interpret further since the phases are not the topic of this work.
The average thickness of the mantle transition of our study area is 245 ± 1 km in agreement with the global average observations of 242 ± 2 km from Ps RFs (Lawrence & Shearer, 2006), 242 ± 20 km from SS precursors (Gu & Dziewonski, 2002;Lawrence & Shearer, 2008), and also the theoretical 250 km thickness (Kennett & Engdhal, 1991), with a total range between 231 and 258 ± 6 km.East of the Cascadia Volcanic Arc, the MTZ is thinned by 5-10 km in comparison to the theoretical.In the middle of our research area, roughly along the coast, the MTZ is of normal thickness with zones of slight thickening, up to 10 km, in the north beneath south Vancouver Island to the Olympic Peninsula and in the south beneath the Cascadia Arc in northern California.The thickening is mainly because of the depression of the corresponding 660.The greatest amount of thinning (15-20 km) is located beneath the ridges and transform fault system in the west from the junction of the Gorda Ridge and the Mendocino Fracture Zone to the inner corner of the intersection of the Juan de Fuca Ridge and the Blanco Transform Fault.We do not observe focused thinning, uplifted 660 or depressed 410 directly beneath the Cobb hotspot.
The 660 topography is uplifted by up to 15 ± 4 km in the west and 10-20 ± 3 km in the east.The 410 is normal to slightly uplifted (4 ± 2 km) in the west of our research area and uplifted by 10-15 ± 2 km in the east.The 410 is depressed in the central parts of our study area, most notably by up to 13 ± 3 km in the northernmost section in a region below the trench and the ridge and by up to 11 ± 3 km in the south beneath the trench east of the Gorda Ridge (Figure 5).The depression of the 410 is 5 ± 2 km on average.The thickened MTZ (by 2-12 km) in the center of our study region, primarily related to a depressed 660, corresponds to subtle P-wave velocity anomalies (0.5%) (Figure 5b).Since this is >600 km east of the slab inferred from body wave tomography, it would need to be an ancient torn piece of Farallon slabs (Bodmer et al., 2018).This is not a major feature of our model and we have no strong interpretation of it.

Discussion
The most prominent anomaly in our result, the thinned MTZ beneath the western portion of our study area by 10-20 km agrees with thinning of 16 km reported in this region by a previous global SS precursors study (Huang et al., 2019).Finer scale variability in our MTZ discontinuity topography is not likely resolved by the global study.The MTZ discontinuities have been imaged beneath western North America by multiple RF studies with resolutions that are similar to ours.Overall, the magnitudes and patterns of the discontinuity topographies and MTZ thickness are variable among these studies, likely owing to differences in data and/or approach, such as data selection (manual vs. automatic) and/or migration assumptions (Eagar et al., 2010;Gao & Liu, 2014;Schmandt et al., 2012;Zhou, 2018).However, we do find some general agreement with previous studies.For instance, two studies found MTZ thinning of 5-10 km in the eastern part of our study region in agreement with the 5-20 thinning observed in our result (Eagar et al., 2010;Schmandt et al., 2012).

Abnormal Discontinuity Topography and Estimated Temperature Anomalies
The 410 and 660 depths may be affected by temperature and/or composition (Ghosh et al., 2013;Hirose, 2002;Ito & Katsura, 1989;Katsura & Ito, 1989;Smyth & Frost, 2002).However, we do not find particularly strong evidence for compositional effects in general, for example, uplifted 410 or depressed 660 accompanied by slow seismic velocities, which would correspond to a water interpretation, or very depressed or broad 660, which would correspond to a basalt interpretation.Instead, the variations of the 410 and the 660 depths are in good agreement with the seismic tomography.The locations of depressed the 410 correspond to low-velocity anomalies.The uplifted and depressed 660 locations correspond to low-velocity anomalies and high-velocity anomalies, respectively (Figures 5a and 5b), which is more consistent with thermal anomalies caused by upwellings or downwellings.
There is a small region on the eastern side of our study area that is characterized by an uplifted 410 and slow velocities in the transition zone (Figures 6b and 6c), that is, indicative of hydration.This could be a possible explanation.However, the thickening could also be related to the nearby slab.If the uplifted 410 in the eastern section of the study area is caused by the slab it would correspond to an anomaly of 121-181 K (Figures 5a  and 6).We do not find a corresponding depressed 660 in the east, likely because the slab is located outside our study region at those depths, as suggested by the body wave tomography anomalies (Figure 6) (Bodmer et al., 2018).Given the ambiguity of the slab location, we do not have a strong interpretation of the observed uplifted 410.
Very large depressions on the 660 are thought to be caused by very high temperatures or basalt contents, for example, where garnet would dominate over olivine at the 660 at 200-300 K above the global average mantle temperature (Hirose, 2002).This has been observed for instance beneath Iceland (Jenkins et al., 2016).However, we do not observe this, nor do we necessarily expect very large temperatures or basalt contents.
Therefore, we proceed assuming the temperature-sensitive olivine phase changes are related to the Clapeyron slopes (Bina & Helffrich, 1994), and that compositional effects on the remainder of the study area besides the eastern portion, if any, are minimal.
The thickened MTZ (by 2-12 km) in the center of our study region, primarily related to a depressed 660, corresponds to subtle P-wave velocity anomalies (0.5%) (Figure 5b).Since this is >600 km east of the slab inferred from body wave tomography, it would need to be an ancient torn piece of Farallon slabs (Bodmer et al., 2018).This is not a major feature of our model and we have no strong interpretation of it.
The lack of strong MTZ thinning beneath the Axial seamount, the volcanically active surface realization of the Cobb Hotspot, could be explained by a variety of possibilities.A conduit that is less than ∼100 km in diameter, the size of one of our bins, is below our resolution.The anomaly could also be muted in magnitude.Gravity arguments have been used to suggest it has an excess temperature of 30°-40°C (Hooft & Detrick, 1995), which would correspond to a 3-5 km transition zone thickness anomaly, that is, near the size of our error bars and at the edge of resolution.Alternatively, the hotspot may take a non-vertical path through the upper mantle, crossing the 410 at the location of our anomaly 100 km to the northeast.
The thermal anomalies predicted for the discontinuity topographies are consistent with a previous result from the Mid-Atlantic Ridge, more than that inferred at the East Pacific Rise, and less than inferences from hotspots.
The more muted 410 thermal anomalies associated with the depressions in the north and south of the study region (60-97 K) similarly agree with 410 observations at the Mid-Atlantic Ridge, which suggested a 60 K anomaly (Agius et al., 2021), and the predictions are also smaller than estimates based on RF imaging of the 410 beneath hotspots including, for example, 100-165 K in Hawaii (Agius et al., 2017) and 200 K in Iceland (Jenkins et al., 2016).The result is different than the inferred lack of a thermal anomaly beneath the East Pacific Rise (Shen et al., 1998), which suggests either that faster-spreading centers are not underlain by deep upwellings, or that the upwellings were at the edge of the resolution of the previous work (Agius et al., 2021).

Dynamics
Several global tomography models also support deep upwellings in this region, finding slow velocities at 660-1,000 km depth (Figure S5 in Supporting Information S1).The variations of the 660 topography contribute more to the MTZ thickness differences as they are correlated at 0.73 with the MTZ thickness (Figure 7).The thinned MTZ and uplifted 660 in our study area suggest that material upwells from the lower mantle beneath the intermediate spreading ridges in Cascadia and beneath the base of the Cascadia Slabs.This expands the interpretation that upwelling occurs beneath slow-spreading environments (Agius et al., 2021) to intermediate-spreading environments.It suggests that sub-slab upwellings inferred from slow velocity anomalies need not necessarily be artifacts of tomography resolution (Bezada et al., 2016;Bodmer et al., 2020).While many geodynamic models predict downwelling sub-slab material, some also suggest that buoyant upwelling material may be possible beneath the base of a subducting slab, which could potentially uplift the 660 (Yang & Faccenda, 2020).
Surprisingly, the 410 depressions do not correspond to the locations of the uplifted 660, and instead, occur roughly between the uplifted 660 regions.This observation is consistent with a non-vertical ascent of material through the MTZ, which is much different from the typical paradigm which includes the simple vertical ascent of material through the MTZ (W.J. Morgan, 1971).However, the observation agrees with observations from the Mid-Atlantic Ridge where the peak 410 topography was offset from that of the 660 (Agius et al., 2021).In our study region, the strongest 410 depressions occur in the north and the south, suggesting a more 3-D flow.In the south, the 410 depression and the eastern 660 uplift occur just beneath the base of the slab as defined by the 0.5% fast contour from body wave tomography (Bodmer et al., 2018) (Figure 5).In the northern sections, the 660 is uplifted beneath the base of the slab but the 410 depression occurs roughly midway between the ridge and the slab, that is, ∼500 km west of the base of the slab.One possibility is that the upwellings from the lower mantle in our study region stagnate within the MTZ before reaching 410.The stagnation could occur because the upwellings dehydrate during ascent through the MTZ, as predicted by geodynamics (Long et al., 2019).
Instead, the locations of thermal upwellings suggested by our 410 depressions are coincident with the locations of slow seismic velocity anomalies from the body wave tomography (Figure 7).The slow velocities were interpreted as buoyant upper mantle flow at <100-350 km depth.In the north, the inferred flow originates from the west, bending eastward toward the slab as it nears the surface (Figure 6, light gray arrows) (Bodmer et al., 2018), which is also generally consistent with west-east shear wave splitting directions reported east of the Cobb Hotspot (Martin-Short et al., 2015).In the south, the flow occurs from east to west up the base of the slab (Figure 6, light gray arrows) (Bodmer et al., 2018), which could also be consistent with complex shear wave splitting directions reported from the eastern side of the Blanco Fracture Zone.Our 410 depressions are coincident with deeper realizations of the slow anomalies associated with the interpreted flow, suggesting that the flow originates from within the MTZ (Figure 6, arrows pointing upwards at 410 km depth).
There are a few possible causes for the inferred reversed flow direction in the north in comparison to the south.The difference may be related to additional buoyancy from the Cobb Hotspot (Bodmer et al., 2018).Neither our MTZ discontinuity depths nor the body wave tomography detects a strong hotspot anomaly.However, we cannot preclude the possibility that the hotspot causes the variation in dynamics with latitude.An alternative is that the difference in flow direction occurs because the southern edge of our study area is located near the hypothesized southern edge of the Gorda slab, as inferred from tomography (Bodmer et al., 2018).Toroidal flow around the slab edge could result in more complex mantle flow in the region (Eakin et al., 2010).The intermittently imaged negative discontinuities at 365-385 km depth may be caused by conversions from the top of a supra-410 melt layer.In the north, the negative phase is strongest directly above the inferred location of upwelling based on 410 topography.In the south, the negative phase exists just to the west of the location of inferred upwelling.Supra-410 melt is predicted to occur when hydrous material upwells from the MTZ into the less-soluble overlying mantle (Bercovici & Karato, 2003), and discontinuities related to such a layer have been imaged by a variety of studies (Frazer & Park, 2021;Wei & Shearer, 2017).The observations are consistent with this model and suggest that upper mantle flow may enter the MTZ and entrain hydrous material to shallower depths.A supra-410 discontinuity is not imaged above the inferred upwelling based on the 410 topography at middle latitudes (Figure 4b).Therefore, MTZ hydration may be variable and/ or supra-410 melt may be transient in time and/or space possibly owing to melt migration (Hier-Majumder et al., 2021).Either could explain the wide variety of hypothesized preferential tectonic locations and/or sporadic imaging of supra-410 melt layers (Bercovici & Karato, 2003;Wei & Shearer, 2017).We note that our testing indicates the supra-410 phases persist regardless of the deconvolution method and migration models used (Text S2, Figure S3 and S7 in Supporting Information S1).However, the melt-layer discontinuities are not necessarily required to support the conclusions of this paper, although they are intriguing and aligned with our interpretation.

Conclusions
We use extended-time multi-taper Ps RFs to image the MTZ discontinuities beneath the Cascadia region using an amphibious data set.We interpret phases that are above the formal error bars of the stack, away from the far edges of the model, and at depths deeper than those contaminated by crustal reverberation artifacts.Although migration model assumptions can impact the absolute depths of the imaged discontinuities, the existence of the MTZ phases and the trends of the MTZ discontinuity depths and MTZ thicknesses that we interpret are robust regardless of migration model assumptions.
The MTZ is thinned beneath the Juan de Fuca and Gorda Ridges and beneath the Cascadia Slabs and corresponds to the locations of the uplifted 660.The most depressed 410 lies beneath the northern Juan de Fuca Plate and beneath the Gorda Plate, offset from where the 660 is largely uplifted.The depths of the 410 and the 660 beneath the Cascadia are not anticorrelated as in classical diagrams for regions where ascending/ descending happens.However, the variations of the 410 and the 660 topography are in good agreement with thermal predictions based on previous seismic tomography observations.We do not find strong evidence for compositional effects in the central portions of our study region and the locations where our interpretations are focused.
Overall, this suggests that deep upwellings from the lower mantle may occur beneath tectonic locations other than hotspots, including beneath ridges and beneath slabs.However, the upwellings are sluggish and may stagnate in the MTZ, possibly owing to dehydration, instead of continuing to ascend vertically.Instead, upper mantle convection can interact with the shallow MTZ, entraining hydrated MTZ material and transporting it into the mantle above.This provides a new transport mechanism for the redistribution of hydration into the mantle above the MTZ.It suggests that in regions distant from a major hotspot like Hawaii and Iceland, that is, the majority of the Earth, upper mantle dynamics may play a larger role in dictating material transfer from the MTZ into the upper mantle.

Data Availability Statement
The methods used are standard and widely used (Helffrich, 2006)

Figure 1 .
Figure 1.Map of the study area.The background shows bathymetry/topography.White inverted triangles represent seismometer locations offshore and white squares represent seismometer locations onshore.Red volcano symbols indicate the locations of major volcanos of the Cascadia Arc and the Axial Seamount (the active volcano associated with the Cobb Hotspot), which is on the Juan de Fuca Ridge.Thick red lines show plate boundaries, and the serrated line shows the trench.Three thick black lines represent the locations of the cross sections in Figures4 and 6.Black circles along the cross sections correspond to a spacing of 100 km.The inset map shows the locations of earthquakes (red circles) used in this study.Juan de Fuca Ridge (JdFR), Blanco Fracture Zone (BFZ), Gorda Ridge (GR), and Mendocino Fracture Zone (MFZ) are as labeled.

Figure 2 .
Figure 2. Piercing points of RFs at 410 and 660 km depth and corresponding P410s P660s RF amplitudes.(a) Amplitudes of P410s versus longitude.Blue circles show piercing point longitudes at 410 km of RFs from offshore data.Orange squares show piercing point longitudes at 410 km of RFs from onshore data.The red thick line shows the median of P410s RF amplitudes from offshore data, and the black thick line shows the median from onshore data.Numbers on the right side of the lines indicate the corresponding median values.(b) Similar to (a) but for P610s.(c) Spatial distribution of piercing points (blue circles for ocean data and orange squares for land data) at 410 km depth.Red lines show plate boundaries.The black circle indicates the bin used in Figure S2 in Supporting Information S1.(d) Same as (c) but for the P660s data set.

Figure 4 .
Figure 4. Vertical cross-sections.(a) Top: bathymetry/topography along 40.5°N from 131° to 115°W.Bottom: Vertical cross-section from 3-D migrated RFs along 40.5°N with location shown in Figure 1.Dashed black lines indicate 410 and 660 km depths.Black circles along the cross sections correspond to a spacing of 100 km.Black arrows indicate supra-410 phases.Small black arrows indicate significantly depressed 410 and uplifted 660.(b, c) the same as (a) but along 45°N (b) and 47.5°N (c), respectively.The plate boundaries are labeled as follows: Gorda Ridge (GR), Mendocino Fracture Zone (MFZ), and Juan de Fuca Ridge (JdFR).

Figure 5 .
Figure 5. Maps of the 410 topography and the 660 topography and the thickness of mantle transition zone (MTZ).(a) Colors show the depths of the 410.Thick black lines represent plate boundaries, and the serrated line shows the trench.Red triangles indicate the locations of major volcanos of the Cascadia Arc and the Axial Seamount.Red contours show −0.5% P-wave velocity anomalies and blue contours show +0.5% P-wave velocity anomalies from body wave tomography (Bodmer et al., 2018).Gray circles indicate the locations of bins used for time-domain tests (Figure S3 in Supporting Information S1).(b) same as (a), but for the 660.(c) similar to (a), but for the MTZ thickness.Thin black lines show contours of MTZ thickness.

Figure 6 .
Figure 6.Vertical cross-sections compared to anomalies from tomography and inferred dynamics.(a) Top: bathymetry/topography along 40.5°N from 131° to 115°W.Bottom: Vertical cross-section from 3-D migrated RFs (gray wiggles) along 40.5°N.Background colors are P-wave velocity anomalies from CAS2018_P(Bodmer et al., 2018).Semi-transparent gray and black arrows with shades indicate inferred pathways of upwellings and downwellings (slabs).Semi-transparent gray arrow shows suggested upper mantle upwellings(Bodmer et al., 2018).(b, c) Are the same as (a) but along 45°N (b) and 47.5°N (c), respectively.Plate boundaries are labeled as follows: Gorda Ridge (GR), Mendocino Fracture Zone (MFZ), and Juan de Fuca Ridge (JdFR).Black circles at 800 km depth are plotted every 100 km along the cross-sections and correspond to those plotted in Figure1.

Figure 7 .
Figure 7. Maps of the 410 and the 660 differential topography, the differential thickness of mantle transition zone (MTZ), and average P-wave velocity anomalies and their correlations.(a) The 410 topography difference compared to 410 km.(b) Average P-wave anomalies (dVp) from 100 to 410 km (Bodmer et al., 2018).(c) Regions of 410 depressions and slow average Vp from 100 to 410 km depth (orange) and regions of uplifted 410 and fast average Vp from 100 to 410 km depth (green).(d) The 660 topography difference compared to 660 km.(e) MTZ thickness difference compared to 250 km.(f) Regions of MTZ thinning and uplifted 660 (orange) and regions of MTZ thickening and depressed 660 (green).