Receiver function mapping of mantle transition zone discontinuities beneath Alaska using scaled 3-D velocity corrections

The mantle transition zone is the region between the globally observed major seismic velocity discontinuities around depths of 410 and 660 km and is important for determining the style of convection and mixing between the upper and the lower mantle. In this study, P -to- S converted waves, or receiver functions, are used to study these discontinuities beneath the Alaskan subduction zone, where the Paciﬁc Plate subducts underneath the North American Plate. Previous tomographic models do not agree on the depth extent of the subducting slab, therefore improved imaging of the Earth structure underneath Alaska is required. We use 27 800 high quality radial receiver functions to make common conversion point stacks. Upper mantle velocity anomalies are accounted for by two recently published regional tomographic S -wave velocity models. Using these two tomographic models, we show that the discontinuity depths within our CCP stacks are highly dependent on the choice of velocity model, between which velocity anomaly magnitudes vary greatly. We design a quantitative test to show whether the anomalies in the velocity models are too strong or too weak, leading to over- or undercorrected discontinuity depths. We also show how this test can be used to rescale the 3-D velocity corrections in order to improve the discontinuity topography maps. After applying the appropriate corrections, we ﬁnd a localized thicker mantle transition zone and an uplifted 410 discontinuity, which show that the slab has clearly penetrated into the mantle transition zone. Little topography is seen on the 660 discontinuity, indicating that the slab has probably not reached the lower mantle. In the southwest, P410s arrivals have very small amplitudes or no signiﬁcant arrival at all. This could be caused by water or basalt in the subducting slab, reducing the strength at the 410, or by topography on the 410 discontinuity, preventing coherent stacking. In the southeast of Alaska, a thinner mantle transition zone is observed. This area corresponds to the location of a slab window, and thinning of the mantle transition zone may be caused by hot mantle upwellings.


I N T RO D U C T I O N
The mantle transition zone delineates the upper from the lower mantle, and its characteristic phase transitions both affect and reflect the style of mantle convection. Subducting plates leave an imprint on the boundaries of the mantle transition zone, the seismic discontinuities at approximately 410 and 660 km depth. By studying these we can constrain the depth extent of subducting slabs. In Alaska, tomographic models do not agree on the depth extent of the subducting Alaskan slab. Here, we use receiver functions, which are waves that convert at discontinuities, to image the mantle transition zone to better constrain the Earth structure underneath Alaska.

The Alaskan subduction zone
Alaska, a U.S. state, is located on the northwestern edge of the North American Plate and along the northern part of the Pacific subduction zone. This subduction zone accommodates convergence between the northernmost part of the Pacific Plate and the Alaskan part of the North American plate at the Alaskan-Aleutian megathrust (see Fig. 1). The Pacific Plate moves 5.2 cm yr -1 to the NNW with respect to the North American Plate (DeMets & Dixon 1999). In the east, the Pacific Plate terminates against the Fairweather-Queen Charlotte transform system. At the eastern part of the subduction zone lies the Yakutat terrane, a region of unusually thick oceanic crust (Gulick et al. 2007), moving with the Pacific oceanic lithosphere and being  (Bird 2003). Red dots denote active volcanoes. The orange lines show the subducted Pacific plate at depth with contour intervals of 20 km as imaged in the Slab2 model by Hayes et al. (2018) and the subducted Yakutat crust as imaged by Eberhart-Phillips et al. (2006) is shown by the thin dashed line. accreted to the Alaskan margin (Fuis et al. 2008). Collision and underthrusting of the Yakutat terrane is believed to have caused some unusual features, such as a region with absence of seismicity known as the Denali volcanic gap and the Wrangell volcanism (Eberhart-Phillips et al. 2006). This Wrangell volcanic field lies in the north of the Yakutat block, and is a region with a decrease in seismicity below 50 km depth. The existence of a Wrangell slab has been debated in several studies (e.g. Jadamec & Billen 2012;Bauer et al. 2014;Jiang et al. 2018;Martin-Short et al. 2018;Gou et al. 2019), but is thought to only be present in the upper mantle and is not expected to have a signature on the 410 km discontinuity.
Seismic tomographic results are also used to debate the extent of the Pacific slab into and across the mantle transition zone (MTZ, 400-700 km depth) further west. Qi et al. (2007) found fast velocities associated with the slab reaching depths of around 400-450 km, but also disconnected fast velocities lower in the MTZ, which they propose to be remnants of an older plate, the Kula Plate. The models of Martin-Short et al. (2018) indicate that the Pacific slab continues to depths below 400 km, but do not have good resolution below that depth. In the upper mantle Vs model constructed by Jiang et al. (2018) the subducting plate is most prominent at depths until 390 km depth, but weaker fast anomalies appear in the MTZ. In the P-wave velocity model for North America by Burdick et al. (2017) the end of the slab seems to penetrate into the lower mantle, but resolution is lost at greater depth. The P-wave velocity model by Gou et al. (2019) shows the slab down to 450-500 km. Comparing both S and P velocities, 'Atlas of the Underworld' (van der Meer et al. 2017) identify the end of the slab within the MTZ.
Plate tectonic reconstructions have long shown that convergence has occurred between the North American Plate and the Pacific and previously Kula and Farallon plates for >100 Myr (Engebretson et al. 1985;Müller et al. 2016;Seton et al. 2012;Madsen et al. 2006). However, the depth extent of the slabs observed in tomographic models is much less than the several thousand kilometres of lithosphere that is expected to be subducted beneath the Aleutian trench based on these plate reconstruction models. This may suggest that the subduction configuration in the past may have been slightly more complex than generally perceived.
The disagreement across different tomographic models on the depth extent of the subducting slab could be due to reducing resolution of the models across the MTZ and the potential smearing along body wave ray paths. Therefore we need additional methods to image the MTZ and the possible extent and existence of slabs in the Alaskan region. In this study, receiver function analysis is used to study the topography of the global mantle discontinuities, to gain more insight in the Earth structure underneath Alaska and the depth extent of the subducting slab.

Mantle discontinuities
Since tomographic models do not provide a consistent answer to the slab structure underneath Alaska, we study the MTZ in this region. The MTZ divides the Earth's upper and lower mantle and is characterized by abrupt changes in seismic velocities. The globally observed major seismic velocity discontinuities occur around depths of 410 and 660 km (hereafter referred to as the '410' and '660') and mark the top and the bottom of the MTZ. These discontinuities have been interpreted as polymorphic phase changes in the Mg 2 SiO 4system, the olivine system. The 410 is associated with a phase transition from olivine (α-(Mg,Fe) 2 SiO 4 ) to wadsleyite (β-(Mg,Fe) 2 SiO 4 ) in an exothermic reaction (Katsura & Ito 1989). The 660 is associated with a phase transition from ringwoodite (γ -(Mg,Fe) 2 SiO 4 ) as it disassociates into perovskite ((Mg,Fe,Al)(Si,Al)O 3 ) and magnesiowüstite ((Mg,Fe)O) in an endothermic reaction (Ringwood et al. 1975;Ito & Takahashi 1989). The non-olivine or garnet-pyroxene system also predicts phase transitions around the depth of the 660, and it depends on temperature which transition dominates (Hirose 2002;Yu et al. 2011).
The phase transitions do not occur at the exact same depth everywhere, but they vary depending on temperature, composition (Cobden et al. 2008;Xu et al. 2008;Hernández et al. 2015) and water content (Smyth & Frost 2002;Thio et al. 2016;Buchen et al. 2018). The behaviour of a phase transition is described by its clapeyron slope, defining the temperature-pressure range in which the transition occurs. For the olivine system, the transition at the 410 has a positive clapeyron slope, while the transition at the 660 has a negative clapeyron slope. Due to the opposite signs of the clapeyron slopes describing the olivine phase transitions around 410 and 660 km depth, the reactions at the top and the base of the MTZ behave opposite to temperature anomalies. This means that in colder regions an uplifted 410 and depressed 660 are expected, which will thicken the MTZ. Conversely, in hotter regions a depressed 410 and uplifted 660 are expected, thus creating a thinner MTZ.
Care should be taken in interpreting the transition zone thickness as temperature anomalies alone, because other factors may also control the reactions. The 660 discontinuity is thought to be especially variable, since not only the transition from ringwoodite to perovskite and magnesiowüstite occurs at this depth, but phase transitions in the non-olivine components play a role as well.
In relatively low temperature regions, the transition in the olivine system is thought to be dominant, but ringwoodite dissociates to ilmenite first before transitioning to perovskite around the depth of the 660 (Stixrude & Lithgow-Bertelloni 2011). Garnet instead of olivine dominates the phase transition at relatively high temperatures (Stixrude & Lithgow-Bertelloni 2011). While temperature estimates of where ilmenite or garnet dominates the 660 vary widely, the phase diagram in Stixrude & Lithgow-Bertelloni (2011) suggests the ilmenite transition dominates the 660 at less than ∼1500 K while the garnet transition dominates the 660 at temperatures larger than ∼2100 K. Geodynamical estimated temperatures for plumes or slabs in the MTZ also vary. For example, Zhang et al. (2013) find temperatures of 1100 K at transition zone depth for a geodynamical slab model, and Maguire et al. (2018a) model a plume excess temperature of δT=400 K which relates to temperatures of ∼2200 K at 660 km.
These phase transitions can cause a complex and broad 660 discontinuity with sometimes multiple signals (Deuss et al. 2006;Andrews & Deuss 2008). The presence of water can further complicate the phase transitions in the MTZ. In the olivine system, wadsleyite and ringwoodite can store significant amounts of water, which alter the kinetics of the phase transformations (Ohtani et al. 2004).
Receiver functions contain P-to-S converted phases from discontinuities. Many studies have used teleseismic receiver functions to image an anomalous MTZ. Some examples are Maguire et al. (2018b), Tauzin et al. (2017), Jenkins et al. (2016), Owens et al. (2000), Li et al. (2000), Andrews & Deuss (2008) and many more. In this study receiver functions are used to map the topography on the major seismic mantle discontinuities. Receiver functions are sensitive to the local Earth structure underneath a seismic station by representing the response to the incoming compressional (P) wave, and can therefore be used to study discontinuities at different depths.
To use receiver functions to map the topography on the discontinuities, one needs to correct for anomalous crust and upper mantle structure. The exact velocity structure is usually not known, and accurately correcting for all velocity heterogeneities is challenging. Differential times between the arrivals from the 410 and 660 are less sensitive to variations in upper mantle structure, and many studies therefore interpreted the MTZ thickness instead of the discontinuity depths. Tauzin & Ricard (2014) propose a method to seperate the competing effects of velocity structure and temperature on the discontinuities, and use this to estimate the clapeyron slopes of phase changes in the transition zone. Also, some receiver function studies test the robustness of the obtained discontinuity depths using synthetic velocity models (e.g. Jenkins et al. 2016).
Several receiver function studies have been conducted in Alaska. For example, Ferris et al. (2003), Veenstra et al. (2006) and Miller et al. (2018) looked at the subducting plate at shallow depth and crustal thickness and structure, and Ai et al. (2005) studied the crust and upper mantle discontinuities in a small area in central Alaska. Using receiver functions to study the transition zone discontinuities underneath a large part of Alaska has only previously been done by Dahm et al. (2017). They find a broad zone with a thicker than normal MTZ and shallower than normal 410 and conclude that the subducting Alaskan slab has reached the MTZ where it is probably broken into fragments. They used the velocity model by Martin-Short et al. (2016) for the time-to-depth conversion. Much more data has become available compared to their data set which runs to late 2016. This increase is mainly due to the ongoing deployment of the USArray Transportable Array (TA), which has also greatly increased coverage in the more remote corners of Alaska. Besides the improved data coverage, our study also compares using two recently published S-wave velocity models for upper mantle corrections, which have also benefited from the improved data coverage. By comparing these two regional velocity models, we do show the sensitivity of the MTZ discontinuity topography to the velocity models, and particularly the degree of damping in the inversion, used for the mantle corrections.

P-to-S receiver functions
When an incoming seismic compressional (P) wave meets a boundary or discontinuity the energy of the wave is partly converted to shear (SV) wave energy. This produces waves referred to as Pds phases, where the d denotes the depth of the conversion. The P410s and P660s, the waves converted at the 410 and 660 discontinuity, arrive 40 to 70 s after the direct P-arrival on the radial component of a seismogram. Their delay times can be used to determine the Downloaded from https://academic.oup.com/gji/article-abstract/219/2/1432/5543222 by Utrecht University user on 19 November 2019 depths of the discontinuities assuming a velocity model between the discontinuity and the seismic station (Langston 1979). To highlight the converted phases, we stack the signals. To do this, we first deconvolve the source component. Assuming that the P-wave arrival on the vertical component seismogram represents the source, we deconvolve the vertical component from the radial component. The time series resulting from this is called a receiver function and is sensitive to sharp structures underneath a station.
Creating the receiver functions can be done with various deconvolution techniques. We choose to use the iterative time domain deconvolution described by Ligorría & Ammon (1999). In this method the receiver function is constructed by adding Gaussian pulses defined as G(ω) = e −π 2 ω 2 a 2 with a = 0.4 rad s -1 . This filter is also used to pre-filter the two components. The resulting Gaussian pulses in the receiver function in the time domain have a full-width at half maximum of ∼4 s. Using higher frequencies will result in lower amplitude conversions, for example shown by Jenkins et al. (2017), while using lower frequencies will broaden the pulses further, potentially capturing less detail in the receiver functions.
The timing and amplitude of iteratively added Gaussian pulses to the receiver function is determined with the following steps: The current receiver function is convolved with the vertical component to produce the predicted radial component. This predicted radial component is then subtracted from the observed radial component to produce the unmodelled signal in the data. The new Gaussian pulse is added where the unmodelled signal is largest, determined by cross-correlating the unmodelled signal with the vertical. The Gaussian pulse is scaled by the amplitude of the maximum of the cross-correlated signal. These steps are iterated until adding more pulses shows insignificant improvement or 200 pulses have been added.

Data acquisition and processing
Data are collected from the Incorporated Research Institutions for Seismology (IRIS) database. ObsPy (Krischer et al. 2015) is used to download and process the data. Stations are selected with latitudes of 56-72 • N and longitudes between 130 • and 170 • W. The events used occurred between January 2000 and October 2018 with magnitudes (M w ) between 5.5 and 8.3. Only events within an epicentral distance range of 30 • to 90 • from the stations are selected. Earthquakes with an epicentral distance smaller than 30 • do not have a ray path with a turning point deep enough to image the MTZ discontinuities, while for events with an epicentral distance larger than 90 • the P wave is affected by the core-mantle boundary. The final data set before quality control consists of 405 348 events to station pairs from 477 stations across Alaska. Figs 2(a) and (b) show the distribution of events, stations and pierce points (the location at depth where the P-to-S conversion occurs).
The trend and mean in the data are removed and the north-south and east-west traces are rotated into radial and transverse component seismograms with 10 samples per second. Data is prefiltered between 0.01 and 0.2 Hz. Seismograms are then cut from 25 s before to 150 s after the predicted direct P arrival and these traces are used to create the receiver functions.
An automated quality check on the created receiver functions is done using the following criteria: (i) The seismograms should have a signal-to-noise ratio of 2.5, that is the energy of the P-wave arrival should be larger than 2.5 times the average energy in the time window used for the receiver functions.
(ii) The receiver function convolved with the vertical component seismogram should at least reproduce 60 per cent of the radial component.
(iii) The direct P arrival in the receiver function needs to be within 1 s from zero, that is the main energy on the radial and vertical should be correlated.
(iv) There should be no peaks before the main P wave which are more than 30 per cent of the amplitude of the P wave and no peaks after the main P wave which are more than 70 per cent of the amplitude of the P wave.
(v) There should be signals after the main P arrival in the receiver function which are at least 4 per cent of the amplitude of the P wave.
Data for all stations are checked visually to ensure the automated quality check did not leave receiver functions with anomalous characteristics in the data set. After the quality check the data set contains 27 800 receiver functions.

Time-to-depth conversion
Next, a time-to-depth conversion is performed on every trace in order to interpret the receiver functions in terms of discontinuity depths. We have used three velocity models to perform this conversion. First, the 1-D radial velocity model PREM (Dziewonski & Anderson 1981) is used. The seismic travel time calculator Taup (Crotwell et al. 1999) as implemented in ObsPy (Krischer et al. 2015) calculates the predicted arrival times for the direct P wave and P-to-S conversions. The velocity model calculates the predicted arrival time for a range of Pds phases using PREM and interpolates for all conversion depths, taking into account epicentral distance and event depth.
Because the Alaska region contains 3-D upper mantle velocity variations, we need to account for these in the time-to-depth correction. Therefore, two recently published regional 3-D velocity models are used: (ii) The model by Jiang et al. (2018) is an upper mantle Rayleigh and S-wave tomography model, constraining absolute velocities in the top 240 km and relative velocities beneath, and is hereafter called model JSWLW.
Where the models give relative velocities, these are converted to absolute velocities using IASP91 (Kennett et al. 1995). The P-wave model for the relative velocity models is in both cases obtained by scaling the shear wave velocity perturbation by assuming that dln V S /dln V P increases from 2 at the surface to 3 at the core-mantle boundary, following Ritsema et al. (2011) Where absolute V S values are given in the uppermost 200 km in MABPM and the uppermost 240 km in JSWLW we apply the same dln V S /dln V p scaling as for greater depths, using the average 1-D V S in each model as the reference velocity. The V P reference velocity is obtained by rescaling the average 1-D V S by the V P /V S ratio in IASP91.
The MABPM model includes its own crustal model. For model JSWLW crustal model CRUST1.0 is added (Laske et al. 2013). These two models and PREM are each used to trace back the P and converted S phases within the 2-D plane between the station and earthquake. We use the pierce points and delays for the converted phases in our stacking procedure. Fig. 3 shows a slice through both 3-D velocity models at 390 km depth. Note that JSWLW covers a wider area than MABPM. The velocity anomaly patterns of these models are similar, with at the same location the most prominent feature being the fast velocity anomaly in central Alaska, which corresponds to the subducting Pacific slab. However, the anomalies in model JSWLW are approximately twice as strong as the ones in model MABPM. This has large implications for the time-to-depth correction, as we will show in Section 3. In Section 4, we explore rescaling of the time-to-depth corrections for both models.

Common conversion point stacking
To improve the signal to noise ratio the time-to-depth converted traces are used to make common conversion point (CCP) stacks (Dueker & Sheehan 1997;Kosarev et al. 1999). CCP stacks are made with receiver functions that sample the same region at depth while recorded at different stations and originating from different sources. We discretize the area underneath Alaska with a grid spacing of 0.5 • in latitude and longitude and 2 km in depth, for depths between 60 and 1300 km. Every receiver functions is then traced back from the recording station to the source, using one of the velocity models described above. The amplitudes in the receiver functions at each depth are stacked into nearby grid points using a weighting factor based on the distance to the piercepoints. The weighting factor follows Lekic et al. (2011) and  and depends on the Fresnel zone half width (FHZW), defined as where λ(z) is the wavelength at depth z. The FZHW is used to normalize the distance between the conversion point at a given depth and the grid point in the stack: The weight γ (¯ ) is then given by: In this way some smoothing is applied beyond the width of the Fresnel zone, that is at the edge of the Fresnel zone the weighting factor is 0.25. The smoothing increases as a function of depth, as the Fresnel zone widens with depth. Additionally, we smooth generously by using the FZHW for a 10-s shear wave, which is a dominant period in shear waves and at the half-width of our Gaussian filter. The sum of the weights at 410 and 660 km depth can be seen in Fig. 4 for the CCP stack using PREM (Dziewonski & Anderson 1981). The coverage is wider for 660 km depth, because the ray piercing points at 660 km depth have a wider distribution around the stations and the Fresnel zones are wider. The coverage is densest beneath the southern part of Alaska, due to the large number of seismic stations in that region. The summed weight maps for the 3-D models look very similar.
Errors in the receiver function stacks are calculated using a running standard deviation at each grid point, which is then used to calculate the standard error. We consider arrivals as significant when they have a stacked amplitude larger than twice the corresponding standard error and therefore non-zero with 95 per cent confidence.

R E S U LT S
The receiver function data set is used to make CCP stacks to study the topography on the transition zone discontinuities underneath Alaska using PREM, MABPM and JSWLW. Here we will present our results on the MTZ thickness and the topography of the 410 and the 660, for the different velocity models used.

Topography of the 410 and 660
Fig. 5 shows the transition zone thickness for the CCP stacks corrected with PREM, MABPM and JSWLW respectively. The MTZ thickness is calculated by subtracting the depth of the 410 from the depth of the 660, represented by the peak amplitudes around those depths in the stacks. The thickness of the MTZ varies considerably (Fig. 5). The most prominent and robust feature is the thicker transition zone beneath central Alaska. Around that region the MTZ thickness seems to be close to the globally considered average thickness of 250 km, with thinner parts in the southeast and southwest. Fig. 6(a) shows the depths of the peak amplitudes between 382 and 442 representing the 410 discontinuity. Fig. 6(b) shows the same for peak amplitudes between 639 and 699, representing the 660 discontinuity. This scaling is chosen in such a way that the mean depth is located in the middle of the colour scheme. The mean depths are 413 and 671 km for MABPM and 411 and 667 km for JSWLW. Since the strongest velocity variations are present in the upper mantle, the transition zone thickness is less dependent on the velocity model used for the mantle corrections than the absolute  discontinuity depths of the 410 and 660. Indeed, the transition zone thickness maps corrected with the different models (Fig. 5) are more in agreement with each other than the discontinuity depths maps (Figs 6a and b). The thickness variations seen here are thus robust features.
Next, we compare the discontinuity topography maps, starting with the 410 topography using MABPM with that using PREM (Fig. 6a). The most prominent feature in both images is the significantly uplifted 410 in central Alaska. The location of these uplifted P410s arrivals corresponds well with the location of the slab at 410 km depth in the 3-D tomographic models MABPM and JSWLW (see Fig. 3). Since MABPM accounts for the fast velocities inside the slab, the time-to-depth correction shifts the arrivals deeper than with PREM where there is no correction for the fast velocity anomaly. Correcting with this velocity model takes out some of the extreme values present in the PREM results. Due to the positive clapeyron slope of the transition from olivine to wadsleyite an elevated 410 is expected at the location of a cold subducting slab. As this anomaly largely correlates with fast shear velocities in the tomographic models, the uplifted 410 in the stack using the MABPM model is most likely a signature of the slab at that depth. The 410 further has a deep part in the south east (around 60 • N and 135 • W) and a shallow part in the north east (around 69 • N and 135 • W). The 660 discontinuity corrected by MABPM (Fig. 6b), shows less topography than the 410. The P660s arrivals appear slightly shallower beneath central Alaska, and slightly deeper surrounding Downloaded from https://academic.oup.com/gji/article-abstract/219/2/1432/5543222 by Utrecht University user on 19 November 2019 this area, with an uplifted area in the north east and a deeper part in the south east. Again, if we compare it with the 660 map corrected with PREM, some extreme topography has been subdued resulting in less variation on the corrected 660 topography.
As opposed to MABPM, the JSWLW covers our entire study area. Additionally, the velocity anomalies are stronger in this model, therefore we also observe larger corrections with respect to the topography in the PREM stack, resulting in less topography on the 410 and 660. The 410 does not show a clear uplifted part around the slab location (around 64 • N and 152 • W), although this area is surrounded by a ring of uplifted 410 correlated with slow seismic velocities in the JSWLW model (see Fig. 3). The topography shows weaker variations in depth overall and does not have large consistent structures. Contrarily the 660 shows a prominent deep part in the south around -152 • W, which is not seen when correcting with PREM or MABPM. Around the edges of our study area the JSWLW maps generally show more average depths for the 410 and 660 compared to the other two maps. Shifts are particularly strong in the southeast. The slow velocities in the JSWLW model here, cause what appears to be a deep 410 and deep 660 in the other two maps, to appear as a slightly deep 410 and shallow 660.
Apparently, the resulting depth of the 410 and 660 discontinuities in the CCP stacked receiver functions is very dependent on the velocity model used for the mantle corrections. Implications of using these different models and interpretation of the imaged features will be discussed in the Section 4.

Cross sections
We also present two cross sections to show the Pds stacked amplitudes and waveforms (Fig. 7). Relative velocities in MABPM and JSWLW are plotted in the background to locate the large fast velocity anomaly corresponding to the subducting Alaskan slab. Both cross sections intersect the location of the subducting slab in the velocity model. In general, P410s is usually characterized by a simple single peak in the seismic observations, while P660s often has a broader and more complicated waveform (Deuss et al. 2006;Andrews & Deuss 2008). Indeed, in the cross sections it can be seen that P660s has a broader and more asymmetric peak than P410s, although in the frequency band used here, there is no observable splitting of the P660s.
We start with the cross sections for the stacks corrected with MABPM. On cross section AB, it can be seen that the 410 discontinuity is uplifted around the slab location (between angular distance 4 • and 6 • ), and deeper in the rest of the cross section. Looking at the 660, there is less topography on this discontinuity than on the 410. In the south (near B), both the 410 and the 660 appear deeper than in the rest of the profile. In cross section CD, located to the south of AB, there is again not much topography observed on the 660 discontinuity. In addition, in a part of this cross section amplitudes of the CCP stacked traces decrease significantly around 410 km depth.
The cross sections corrected with JSWLW show differences with the cross sections described above. Around the slab location in cross section AB, we do not see an uplifted 410, but instead the 660 appears deeper there. In the south, the 410 is somewhat deeper, while the 660 is shallower than in the rest of the profile. On cross section CD the 660 also seems to be deeper around the slab location, as opposed to little variation seen on the 660 when using the MABPM model. The 410 arrivals decrease in amplitude in a part of the cross section, just as we observe for the MABPM results. So, although the depth of the discontinuities shifts when using a different velocity model, the waveforms stay the same.
The CCP stacks show hints of additional arrivals in the upper mantle, within the MTZ, and at the top of the lower mantle, some of which can be seen in the cross sections. We have analysed the incidence angles or slownesses of these arrivals, but the results show they are more likely multiples resulting from shallow complexities in the upper mantle than robust arrivals from depth.

Dependency on velocity models
As shown in Section 3, the depth of the discontinuities is very dependent on the velocity model used for the mantle corrections. The strong differences between the two recent models applied here are due to the large amplitude differences in the tomographic models. In this Section we will show through quantitative analysis that one model is undercorrected while the other is overcorrected. We also show how tomographic velocity corrections can be rescaled to improve upper mantle wave speed corrections within CCP stacks using this test.
Looking at the two 3-D regional models in Fig. 3, the pattern of fast and slow anomalies is similar, but the anomalies at a given depth are much stronger for JSWLW than for MABPM. For example, at 390 km depth, the largest anomaly in the slab in MABPM is 2.3 per cent fast while the largest anomaly in the slab in JSWLW is 5.1 per cent fast, more than twice as large. At 250 km depth, the largest anomaly in JSWLW has a value of 11.7 per cent, while this is only 3.6 per cent in MABPM. Both models have relative traveltime anomalies at deeper depths (below 200 km for MABPM and 240 km for JSWLW). For such models the mean 1-D absolute reference velocities is unknown, and our assumed background model might cause artefacts (Bastow 2012). For example, what we use in our corrections as slow anomalies, that is the anomalies around the fast slab in Fig. 3, might only be slow with respect to the regional mean and not the global 1-D reference model we use here. Further, the ill-posed nature of tomographic inversion problems leads to nonunique solutions. For example, different methods (e.g. using a full waveform inversion scheme or not), choices of model parametrization, the amount of damping and smoothing and the amount of iterations can all lead to different resulting tomographic models, with anomalies of widely varying amplitudes. De Wit et al. (2012) studied the non-uniqueness of tomographic models by examining the null-space of the forward operator and found a broad range of acceptable solutions compatible with the data. Quantitative assessment of the posterior model covariance and synthetic sensitivity analyses could help assess the difference and accuracy in observed amplitudes in the models (Rawlinson & Spakman 2016).
Here we present a different way to assess the amplitudes of velocity anomalies in tomographic models through a number of correlation tests. Some receiver function studies use the correlation between the 410 and 660 depths to test whether there are unaccounted for velocity perturbations after mantle corrections (e.g. Dueker & Sheehan 1998;Chevrot et al. 1999;Gao & Liu 2014;Dahm et al. 2017). We add to this calculations of the correlation between the discontinuity depths and the correction applied to the discontinuity depths for the different velocity models. Fig. 8 shows the various correlation plots for discontinuity topography in the CCP stacks for the velocity models PREM (Fig. 8a),  and JSWLW (Figs 8e-f). For this analysis we  (Fig. 3). Our CCP grid spacing every 0.5 • in latitude and longitude direction results in unequal sampling. We therefore interpolate the observed topographies to a grid of 1336 equally spaced points (914 for MABPM) for use in these statistical tests. We used the Pearson correlation coefficient, or Pearson's r, as a measure of correlation. Here we describe the three different correlation tests shown, what we expect for each of these correlations with the help of the cartoon in Fig. 8(h), and what we observe in the case of each model. (Figs 8a, b and e): If the topography on the 410 and 660 is purely due to temperature variations in a olivine-dominated system, we expect anti-correlation between the topography on the 410 and 660 (solid lines in Fig. 8h). In the case that the 410 and 660 are affected by different compositions and temperatures, that is anomalies do not extend across the MTZ, there could be an overall lack of correlation. For these reasons we do not know if after 3-D corrections the 410 and 660 should be anticorrelated, or decorrelated. Additionally, a garnet dominated transition zone at higher temperatures results in local deepening of the 660 and therefore could introduce a degree of positive correlation with the 410 topography.  Lack of corrections for upper mantle structure has the same topographic effect on the 410 and 660 and causes correlation between the two (dashed lines in Fig. 8h). Any incorrect shifts by the velocity model across the upper mantle could increase the correlation between the 410 and 660, while incorrect shifts from velocity anomalies within the MTZ would decorrelate the 410 and 660. As expected, PREM results show the highest correlation between the depths of the 410 and 660 discontinuity, with a Pearson's r of 0.79 (Fig. 8a). Corrections using MABPM lower the Pearson's r to 0.53 (Fig. 8b), but as correlation remains significant, the topography is likely undercorrected. The results corrected with JSWLW have a very low Pearson's r of only 0.08 (Fig. 8e), potentially implying suitable corrections, but this remains inconclusive since the expected correlation given perfect corrections is unknown.

(i) Correlation between the 410 and 660
(ii) Correlation between the 410 and corrections on the 410 (Figs 8c and f): Assuming perfect 3-D velocity corrections, we could expect three possible cases here. If there is no correlation between the temperatures in the MTZ and shallow mantle, we expect no correlation between the resulting topography and the corrections. If there is some correlation between the temperatures in the MTZ and shallower structures, that is in the case of continuous hot upwellings or cold downwellings, the positive Clapeyron slope of the 410 phase transition leads to negative correlation with the corrections. In other words, a deep 410 in a hot area was shifted to shallower depths, while a shallow 410 in a cold area was shifted to deeper depths (Fig. 8h). If there is some anticorrelation between the MTZ and upper mantle, that is where a slow upper mantle sits above a cold slab in the MTZ, the corrections would be opposite in sign, and a positive correlation is expected.
(iii) Correlation between the 660 and corrections on the 660 (Figs 8d and g): Due to the opposite sign of the Clapeyron slope of the 660 phase transition compared to that for the 410, we expect opposite behaviour in the correlation coefficients. Where there is correlation between temperatures in the MTZ and shallow mantle, we expect positive correlation between the topography of the 660 and the corrections thereof. Where the 660 is deep, corrections have pushed it deeper, and where the 660 is shallow, the corrections have pulled it shallower. When there is anti-correlation between the MTZ and upper mantle, a negative correlation is expected.
In the case that there is a degree of anti-correlation between the temperatures in the MTZ and shallow mantle, we would see positive correlation between the depth of the 410 and corrections thereof, while at the same time negative correlation between the depth of the 660 and corrections on the 660. Since we do not see this in either one of the models when we apply the statistical test, we argue that in this area we do not have such an anticorrelation. We therefore argue that we either have zero, or a positive correlation between temperatures in the MTZ and the shallow mantle.
Looking at the correlation between the 410 and the correction on the 410 for the two models (Figs 8c and f), the MABPM results show a degree of negative correlation, while those for JSWLW show positive correlation. If we take the thought experiment that there is a degree of positive correlation between the MTZ and upper mantle, we would expect to see a negative correlation here. The positive correlation for the JSWLW then indicates that the corrections on the 410 are either too strong and thus placing the discontinuity too deep (dotted line in Fig. 8h) or are introducing incorrect topography through incorrect anomalies.
For the correlation between the 660 and the correction on the 660 (Figs 8d and g), the Pearson's r for the two models again has an opposite sign. In the thought experiment of positive correlation between the MTZ and upper mantle, the unexpected negative Pearson's r for MABPM means this case is undercorrected (dotted line in Fig. 8h). The Pearson's r value of 0.82 for JSWLW is also suspiciously high, and combined with the positive correlation for the 410 in Fig. 8(f) could be explained by topography introduced by over-correcting of the model. The strong positive correlation could also, combined with the decorrelation between the 410 and 660 in Fig. 8(e), suggest strong corrections across the MTZ, introducing topography on the 660 only.
In summary, the negative correlation coefficient between the 660 and corrections on the 660 for the MABPM results suggests undercorrection, while the positive correlation between the 410 and the corrections on the 410 for the JSWLW results suggests overcorrection. This is consistent with the stronger amplitudes of the velocity variations seen in the JSWLW over MABPM (see Fig. 3).
We investigate the implications on the correlation coefficients by rescaling the topographic corrections by a range of factors. In order to do so, we take the difference between the topography of the 3-D model and PREM, multiply that with a specific factor and add it to the PREM topography again. A major assumption in this approach is that the corrections everywhere are over-or underestimated by the same degree, while in reality resolution and amplitude recovery will vary laterally and with depth in a tomographic model. Fig. 9 shows behaviour of the three correlation coefficients exploring a rescaling factor for the MABPM model between 1.0 and 2.4, and for the JSWLW model between 0.5 and 1.2. The correlations vary with scaling factor in ways we expect. Undercorrected models lead to negative correlations between both the depths of the discontinuities and the corrections on the depths, while overcorrected models lead to positive correlations for both. It should be noted that the correlation coefficient for the 660 increases faster than that for the 410, which is due to there being less topographic variations on the 660 compared to the 410 causing any topography introduced by overcorrections to lead to strong correlations.
As expected, the increase in correction rescaling factor generally leads to a decrease in the correlation between the 410 and 660, but if corrections become too strong these can increase again as seen for the JSWLW case, where the minimum correlation between the 410 and 660 lies around a factor of 0.85-0.9.
We show the range of Goldilocks values for the rescaling factor (grey shaded area in Fig. 9), where the correlation between the depth of the 410 and the correction on the 410 depth is zero or negative and the correlation between the depth of the 660 and the correction on the 660 is zero or positive. We suggest the optimal values are where these correlations are roughly of the same magnitude at opposite signs, these correlations result from a degree of correlation between anomalies in the MTZ and shallower structure affecting both the 410 and 660. This gives factors of 1.3 and 0.7 for MABPM and JSWLW, respectively. For these factors the correlation coefficient between the 410 and 660 is 0.39 for MABPM and 0.17 for JSWLW. As we do not expect any positive correlation between the 410 and 660 due to anomalously hot, garnet-dominated, MTZ beneath Alaska, we attribute this residual correlation due to imprecision of the tomographic models in general, and our simplified rescaling in particular. While the rescaling values found are largely driven by the difference in amplitude of the velocity anomalies in the two models, the exact values found also depends on our assumptions when scaling to Vp velocities.
This exploration additionally shows that their are no choices of rescaling that suggest decorrelation or a degree of anti-correlation between the corrections and the topography seen on the MTZ, that is their are no values where both the 410 and 660 are decorrelated with the corrections, or values where the 410 has positive correlations with the corrections while at the same time the 660 has negative ones. We could however have chosen rescaling factors, close to our current chosen factors, where corrections are only correlated with the 410 topography and decorrelated with the 660 topography.
The resulting topographies using these scaling factors are shown in Fig. 10. Compared with Fig. 6, the topographies corrected with the different models are, as expected, more similar. The correlation between the two models changes from 0.66 in the old scaling to 0.86 in the new scaling for the depth of the 410, and from -0.16 to 0.42 for the 660. Interpretation of the final topographies is described in the next section. With the tests above we provide a quantitative assessment of the amplitudes of velocity anomalies in tomographic models based on an independent data set. It should be noted that this assessment is not directly sensitive to the amplitude of the anomalies as such, but to the vertically integrated anomalies, and can therefore have less resolution on the degree of smearing that occurs in a tomographic model. The rescaling approach provides a new avenue to use the tomographic models for a better time-todepth correction, but needs to be tested at different regional and global scales. In our current application we use a single rescaling factor, ignoring the fact that amplitude recovery can vary laterally and with depth for a tomographic model.

Interpretation of the imaged features
The two features that have little dependence on the velocity model used for the corrections are the MTZ thickness and the waveforms of the converted phases. Indeed, the MTZ thickness only shows minor changes after rescaling, and is therefore not shown again in Fig. 10. After rescaling we are able to interpret the absolute discontinuity depths with more confidence than without the rescaling before.
As shown in Section 3 the MTZ thickness for all three models increases in the area corresponding to the location of the subducting slab. A thicker transition zone implies either only an uplift of the 410, or an uplifted 410 and depressed 660. In Fig. 10(a) it can be seen that both rescaled models show an uplifted 410 discontinuity to depths of 390-395 km around the location of the slab. Since the transition from olivine to wadsleyite has a positive clapeyron slope, this uplifted 410 is interpreted as the interaction of the cold slab with the phase transition with an estimated temperature temperature anomaly of 180-290 K, assuming clapeyron slopes between +4 MPa K -1 (Katsura et al. 2004) and +2.5 MPa K -1 (Katsura & Ito 1989).
We do not observe a deeper 660 in the same region as the uplifted 410. Some studies have shown that the clapeyron slope for the 660 might be more gentle than that of the 410 (e.g. Katsura et al. 2003;Fei et al. 2004), leading to less topography on the former discontinuity. However, potential signatures of slabs on the 660 have been found in various places (e.g. Cao & Levander 2010;Suetsugu et al. 2010;Gao & Liu 2014;Lee et al. 2014), with depression of the 660 as much as 35 km. Particularly strong signatures on the 660 have also been found in places of ponding slabs (e.g. . Overall in the area studied here, the 660 discontinuity shows little topography (Fig. 10b) with variations up to ∼10 km, compared the ∼25 km on the 410. The P660s in the stacks shows consistent amplitudes and waveforms, with no observations of splitting. These observations imply that the slab has not reached the 660.
The results on both discontinuities suggest that the Pacific slab, together with the Yakutat block in the eastern part of the subduction zone, is indeed subducting there and that the slab is penetrating Downloaded from https://academic.oup.com/gji/article-abstract/219/2/1432/5543222 by Utrecht University user on 19 November 2019  Results are plotted if the sum of weights for the CCP stack is more than 40 stacked receiver functions and if the arrivals are significant. Note that MABPM results are only shown where the 3-D model is constrained (see Fig. 3). The colour scale is consistent with Fig. 6. through the 410 discontinuity and ending somewhere in the MTZ. This is in agreement with the depth extent of the slab concluded by Dahm et al. (2017), who used the velocity model by Martin-Short et al. (2016), of which the model MABPM is an update. It is also in agreement with the depth extent of the slab in the velocity models that we use, and with the 'Atlas of the Underworld' (van der Meer et al. 2017), and the P-wave velocity model by Gou et al. (2019). These observations raise the question where the thousands of km of subducted lithosphere in common plate reconstruction models (Madsen et al. 2006;Seton et al. 2012;Müller et al. 2016) currently resides in the mantle. Part of the convergence may be taken up by intra-oceanic subduction, as suggested by Domeier et al. (2017) and Vaes et al. (2019). The slab associated with such a subduction zone may correspond to anomalies identified in previous studies (e.g. Sigloch 2011;Müller et al. 2016;Domeier et al. 2017) such as the North Pacific Slab in the 'Atlas of the Underworld' (van der Meer et al. 2017). Our results are consistent with the previously published tomographic images that suggest that the Alaskan slab does not reach the lower mantle and may thus inspire future tectonic reconstructions that reconcile this with large-scale plate convergence.
However, the signal of uplifted 410 extends further to the east (around 65 • N and 142 • W) than we might expect from tomography or present-day surface subduction. Also, looking at the location of the subducting slab in the velocity models (Fig. 3) one would expect the uplifted signal on the 410 to continue west of the Aleutian Island arc. This signal can be observed in JSWLW, but is relatively weak in MABPM.
Further south, beneath the Kenai Peninsula, the area around Katmai National Park and Kodiak Island, shows a different and more confusing topography on the 410. This area lies beneath the subducting slab. A large part of this shows a deep 410, juxtaposed against a shallow 410 beneath the Kodiak Island. Overall the P410s here has weak amplitudes and both maps of MABPM and JSWLW having large regions where no significant arrivals can be picked (shown by the grey areas). The entire area has reasonable data coverage ( Fig. 4) and data quality, as exemplified by the quality of the maps for the 660 in this area. In general though, the observed P410s are weaker in amplitude and have narrower arrivals compared to the 660, making this phase more sensitive to incoherent stacking either due to strong topography in the area, or due to incorrect 3-D velocity corrections of the complex tectonic region above this area. The entire velocity structure may contribute to amplitude variations, for example velocity gradients near a discontinuity or multiples that contribute energy near the delay times of the Pds phases (Andrews & Deuss 2008). It is therefore difficult to judge if the topography and the amplitudes in this area are a true representation of a complex discontinuity or an artefact.
In cross section CD in Fig. 7, there is a region where the uplifted P410s has very weak amplitudes for both velocity models. Again the weak amplitudes could be an artefact of stacking across an area of strong topography. If it is a true signature of the discontinuity, a possible explanation for a weak or absent 410 is the presence of water. The existence of water expands the stability field of wadsleyite to a lower pressure. This causes a broader seismic discontinuity at 410 km depth (van der Meijde et al. 2003;Litasov et al. 2006), and thus smaller amplitudes in the arrivals on the receiver functions. Presence of basalt can also reduce the signal at the 410, because it weakens the phase changes in the olivine system. The model by Xu et al. (2008) showed that an increasing basalt fraction leads to a smaller velocity jump around 410 km depth. Further, the presence of ferric Fe in the mantle may play a role in broadening the 410 (Frost & Dolejš 2007).
A thinner than normal MTZ is observed in the southeast of the area studied, where thicknesses of ∼230 km are observed. The location corresponds to a low velocity anomaly in the 3-D models of Martin-Short et al. (2018), Jiang et al. (2018), Qi et al. (2007) and Burdick et al. (2017). In the geometric plate-tectonic model of Madsen et al. (2006) this location corresponds to the location of a slab window. Qi et al. (2007) suggest this low velocity anomaly is associated with upwelling mantle in the slab window. Both the velocities in the model, and the thinner MTZ are consistent with the interpretation of a hot mantle upwelling with olivine dominated phase transitions. If we take 20 km thinning of the MTZ and look at an example where the 410 is 10 km deeper and the 660 is 10 km shallower, this corresponds to a temperature increase between 94 and 150 K on the 410 (assuming clapeyron slopes between +4 MPa K -1 (Katsura et al. 2004) and +2.5 MPa K -1 (Katsura & Ito 1989)) and between 118 and 160 K on the 660 [assuming clapeyron slopes between -3.4 MPa K -1 (Hernández et al. 2015) and -2.5 MPa K -1 (Ye et al. 2014)].
In the northeast, the maps using PREM show an extremely shallow 410 and 660. This area lies at the edge of the North American craton. Cratons consist of highly stable continental crust, which are chemically depleted areas with a lower than average temperature (Yuan & Romanowicz 2010) and have a higher than average velocity. This area lies outside the region where the velocity model MABPM has good resolution. The JSWLW model shows fast velocities in this region, which, as expected, shifts the 410 and 660 deeper, resulting in average depths in the original stack (Fig. 6) and close to average depth after rescaling (Fig. 10). In all models the MTZ seems slightly thicker than average, around 260 km thick. In the global model UU-P07 (Amaru 2007;Hall & Spakman 2015) the fast anomaly below northwestern North America is interpreted as the Yukon slab, visible within the upper part of the lower mantle. Thickening of the MTZ here might be due to the presence of the Yukon anomaly.

C O N C L U S I O N S
Using 27 800 high quality receiver functions to map the MTZ beneath Alaska, we demonstrate that receiver function common conversion point stacks are highly dependent on the velocity model used for the mantle corrections by using two different recent models with differing velocity anomaly amplitudes. We show that correlations between the depth of the discontinuities and the topographic corrections from the 3-D velocity model can be used to test the degree of under-or overcorrection. We find one model undercorrected and the other model overcorrected, and rescale the topographic corrections to obtain comparable topographic maps. The final maps show the effects of the slab on the 410, confirming that it has penetrated into the MTZ. As the 660 appears unaffected, the slab must end within the MTZ. In some parts, the P410s arrivals have a very small amplitude, which can be caused by compositional variation in the mantle, or by too much topography on the discontinuity. A thinner than normal transition zone is observed in the southeast, which may be caused by hot mantle upwellings associated with a slab window.

A C K N O W L E D G E M E N T S
The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata and derived products used in this study. The largest seismic networks contributing data to this research were the USArray Transportable Array from the EarthScope USArray project and the Alaska Regional Network operated by the Alaska Earthquake Center. We thank R. Martin-Short and C. Jiang for providing the 3-D velocity models for Alaska, and we thank Alistair Boyce and two anonymous reviewers for helpful comments. AvS and AD were supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 681535 -ATUNE) and a Vici award number 016.160.310/526 from the Netherlands organization for scientific research (NWO). SC was supported by the European Research Council (ERC) under Downloaded from https://academic.oup.com/gji/article-abstract/219/2/1432/5543222 by Utrecht University user on 19 November 2019