Big Mantle Wedge and Intraplate Volcanism in Alaska: Insight From Anisotropic Tomography

We determine high‐resolution tomographic models of isotropic P‐wave velocity (Vp) and tilting‐axis anisotropy of the Alaska subduction zone using a large number of local and teleseismic data recorded at many portable and permanent network stations in and around Alaska. We find a flat high‐Vp slab in the mantle transition zone (410–670 km depths) beneath western Alaska, which is connected with the subducting Pacific slab at 0–410 km depths, suggesting that a big mantle wedge has formed under western Alaska. Our tilting‐axis anisotropy model reveals complex mantle flows in the asthenosphere. Corner flow in the mantle wedge above the subducting Pacific slab and toroidal flow in the big mantle wedge are revealed, which may cause the Cenozoic intraplate volcanoes in western Alaska and the Bering Sea. In central Alaska, the mantle wedge beneath the Denali volcanic gap is characterized by high‐Vp and subhorizontal fast velocity directions normal to the volcanic arc, which may reflect a remnant of the subducted Yakutat slab. In SE Alaska, the shallow subduction of the Wrangell slab is visible above 150 km depth, and hot mantle upwelling through the Wrangell‐Yakutat slab gap may contribute to the Wrangell volcanic field.

In SE Alaska, the Wrangell Volcanic Field (WVF) is located in the Wrangellia Composite Terrane.The WVF is dominated by large andesitic shield volcanoes and calc-alkaline affinity lavas typical of continental volcanic arcs (Richter et al., 1990).However, the WVF also exhibits some unusual characteristics (Martin-Short et al., 2018), such as the presence of adakitic and tholeiitic rocks (Preece & Hart, 2004), and NW-trending progression in activity over time (Richter et al., 1990).In the SE part of the WVF, Mount Churchill is the source of the White River Ash and produced two of the largest explosive eruptions in North America during the last 2,000 years (McGimsey et al., 1992).In addition, the Prindle volcano, a small isolated basaltic cone ∼200 km away from the northeast of WVF, erupted during the Pleistocene (Moll-Stalcup, 1994).
In central Alaska, the Denali Volcanic Gap (DVG), a volcanic paucity area since the Miocene, extends ∼400 km long along the NE strike of the Wadati-Benioff zone (Rondenay et al., 2010).A DVG-parallel tremor zone exists in the Wadati-Benioff zone, which seems associated with the shallow subduction of the Yakutat Terrane (Chuang et al., 2017;Wech, 2016).At the NE end of the DVG, the Buzzard Creek and Jumbo Dome, two typical continental arc volcanoes, erupted in the past 2,000 years and the Pleistocene, respectively (Moll-Stalcup, 1994).
In western Alaska, Cenozoic volcanoes are widespread near the Bering Sea and continental shelf, which are characterized by alkali and tholeiite basalts, basanite, and nephelinite (Moll-Stalcup, 1994).Different from the arc volcanoes, the volcanic rocks in western Alaska commonly contain xenoliths of peridotite, implying an origin from the backarc mantle magmatism (Moll-Stalcup, 1994).The major volcanic fields there include the Imuruk Lake, St. Michael, St. Lawrence, Ingakslugwat, Nunivak Island, and Togiak Basalt (Figure 1a; Table 1).
Seismic instruments across Alaska have been expanded rapidly in recent years (Figure 1c), enabling us to image its crustal and mantle structures in more details and to a greater depth.The 3-D isotropic velocity structure beneath Alaska has been investigated by using various types of seismic data and techniques, including body wave travel times (Burdick et al., 2017;Eberhart-Phillips et al., 2006;Gou et al., 2020;Martin-Short et al., 2016;Qi et al., 2007;Zhao et al., 1995), surface wave dispersions (Feng & Ritzwoller, 2019;Y. Wang & Tape, 2014), full-wave ambient noise (X.Yang & Gao, 2020), receiver functions (Kim et al., 2014;Mann et al., 2022;O'Driscoll & Miller, 2015;Rondenay et al., 2010), as well as seismic refraction and reflection profiles (Christeson et al., 2010;Fuis et al., 2008).The background colors denote elevation whose scale is shown at the bottom-right.The red triangles and red circles denote locations of Cenozoic volcanoes and major volcanic fields, respectively.The blue solid and dashed lines denote major faults and terrane boundaries, respectively (Colpron et al., 2007).The white arrows denote motion directions of the Pacific plate relative to the North American plate (DeMets et al., 1994).AAT, Arctic Alaska Terrane; BC-JD, Buzzard Creek-Jumbo Dome volcanic field; BRF, Border Ranges Fault; CMF, Castle Mountain Fault; DVG, Denali volcanic gap; INF, Iditarod-Nixon Fork Fault; KAF, Kaltag Fault; KOF, Kobuk Fault; KY, Koyukuk Composite Terrane; SEF, Sevier Fault; WCT, Wrangellia Composite Terrane; YCT, Yukon Composite Terrane; YF, Yukon Flats.(b) Distribution of 19,157 local earthquakes used in this study.The dot colors denote focal depths whose scale is shown at the bottom-left.The blue lines denote depth contours of the upper boundary of the subducting Pacific slab (Gou et al., 2019).The green solid and dashed lines outline the inferred Yakutat slab (Eberhart-Phillips et al., 2006) and the Wrangell slab (Gou et al., 2019;Mann et al., 2022), respectively.The yellow circle denotes an inferred aseismic slab edge (Gou et al., 2019).(c) Distribution of 294 permanent network stations (blue squares) and 393 transportable seismic stations (green squares) used in this study.A few recent studies have combined multiple types of seismic data to obtain detailed images of the Alaska subduction zone (Berg et al., 2020;C. Jiang et al., 2018;Martin-Short et al., 2018;Ward & Lin, 2018).However, most of the previous tomographic models are only valid down to ∼200 km depth because the deeper portion is not well sampled by seismic rays from local earthquakes to the seismic stations in Alaska.To investigate the deeper structure under Alaska, earlier studies have used teleseismic data recorded by the Alaska permanent seismic network or the uncompleted USArray Transportable Array in southern and central Alaska (Burdick et al., 2017;Gou et al., 2019;C. Jiang et al., 2018;Martin-Short et al., 2016;Qi et al., 2007).Albeit these previous efforts, the mantle structure and dynamics related to the intraplate volcanism in western Alaska and the Bering Sea are still poorly understood.
Seismic anisotropy refers to the directional dependence of seismic wave velocity.The lattice-preferred orientation of anisotropic minerals, especially olivine, is thought to be the primary mechanism of anisotropy in the upper mantle (e.g., Karato et al., 2008).Meanwhile, the presence of aligned melt and the shape-preferred orientation of isotropic minerals can also produce anisotropy (e.g., Savage, 1999).Therefore, anisotropy can provide valuable insights into mantle flow and the deformation process.In Alaska, two primary types of anisotropy, azimuthal and radial anisotropies, have been studied by making shear wave splitting measurements (Christensen & Abers, 2010;Hanna & Long, 2012;Karlowska et al., 2021;Lynner, 2021;McPherson et al., 2020;Richards et al., 2021;Song & Kawakatsu, 2013;Venereau et al., 2019;Y. Yang et al., 2021), P-wave tomography (Gou et al., 2019;He & Lü, 2021;Tian & Zhao, 2012), and surface wave (ambient noise) tomography (Estève et al., 2021;Feng & Ritzwoller, 2019;C. Liu et al., 2022;Y. Wang & Tape, 2014).The azimuthal or radial anisotropy reflects the directional dependence of seismic velocity on horizontal or vertical plane under an assumption of horizontal or vertical hexagonal symmetry axis (HSA) (J.Wang & Zhao, 2013).However, the two assumptions are not reasonable for subduction zones due to the existence of a dipping slab and related mantle flows (see Zhao et al., 2023 for a detailed review).Recently, a novel tomographic technique has been developed to reveal tilting-axis anisotropy, with the HSA freely orientated in 3-D space (Z.Wang & Zhao, 2021).So far, this type of techniques has successfully revealed the presence of tilting HSA beneath the Japanese region (Z.Wang & Zhao, 2021) and its forearc area (Z.Wang et al., 2022), northern Fennoscandia (Munzarová et al., 2018), Central Mediterranean (Rappisi et al., 2022), and Cascadia subduction zone (Liang et al., 2023).These studies have demonstrated that the tilting-axis anisotropy can reveal more complex 3-D mantle flow by reconciling the contradictory assumptions of azimuthal and radial anisotropies.
In this work, we collect a large number of arrival-time data of local earthquakes and teleseismic events, which are described in Section 2. By enhancing the Tikhonov regularization technique, we improve the 3-D isotropic Vp tomography beneath Alaska, as shown in Section 3. Then we use the new tomographic method (Z.Wang & Zhao, 2021) to determine high-resolution 3-D anisotropic tomography of the Alaska subduction zone down to 800 km depth, as shown in Section 4. Our results shed new light on subduction dynamics and Cenozoic intraplate volcanism in Alaska.

Data
We used P-wave arrival time data recorded by the Alaska Earthquake Center during January 1977 to June 2007, the International Seismological Centre EHB bulletins during January 1970 to April 2021, and the Array Network Facility component of the EarthScope USArray program during January 2015 to April 2021.After merging the three data subsets and removing their duplicated records, we selected the arrival time data of local earthquakes (Figure 1b) according to the following criteria: (a) The study volume is divided into many cubic blocks with a size of 10 km × 10 km × 10 km, and in each block only one earthquake is selected that has the maximum number of P-wave arrivals.(b) Each earthquake was recorded at over 10 stations.(c) All earthquakes are relocated and their uncertainties are less than 5.0 km in hypocentral location and less than 0.5 s in origin time.(d) All P-wave arrivals have an epicentral distance of less than 14° and cut-off residuals of −2.0 and +2.0 s.
We also selected teleseismic events (Figure 2a) according to the following criteria: (a) Each event was recorded at over 7 stations in the study region.(b) All events have an epicentral distance between 30 and 90°.(c) The whole Earth's crust and upper mantle (0-670 km depths) is divided into cubic blocks with a size of 20 km × 20 km × 20 km, and in each block only one event is selected that has the maximum number of P-wave arrivals in our study region.
(d) Because our study region is large (Figure 1), the whole mantle heterogeneity outside the study volume would contaminate the teleseismic relative residuals (Liang, Zhao, Hua, & Xu, 2022;Zhao et al., 2013).Hence, we make the whole-mantle correction (Figures 2b and 2c) to the teleseismic data using the Tohoku2013 global Vp model (Zhao et al., 2013), similar to the teleseismic study of Liang, Zhao, Hua, and Xu (2022) for central and eastern USA.As a result, our final data set contains 1,050,673 P-wave arrivals from 19,157 local earthquakes (Figure 1b) and 702,929 relative residuals of 181,130 teleseismic events (Figure 2).Compared with the data set used by Gou et al. (2019), our data set is certainly better because we have collected more and better data recorded during 2017-2021.

Inversion Strategy
We first use the tomographic method of Zhao et al. (1992Zhao et al. ( , 2012) ) to obtain a 3-D isotropic Vp model beneath Alaska.A P-wave observation equation can be written as: where residual r ij denotes the observed travel time minus theoretical travel time for the ith local event recorded at the jth station; T i , θ i , φ i , and h i denote the origin time, latitude, longitude, and focal depth (i.e., hypocentral parameters) of the ith local earthquake, respectively; Δ denotes a parameter perturbation; s n denotes P-wave slowness at the nth grid node (i.e., structural parameters); ∂T/∂θ, ∂T/∂φ, ∂T/∂h, and ∂T/∂s are partial derivatives of travel time with respect to the hypocentral and structural parameters; ε ij contains higher-order partial derivatives terms and picking error of the observed arrival time.The related tomography parameters are listed in  10, 25, 40, 65, 90, 120, 160, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 800, 900 km 0, 50, 150, 350, 500, 600, 700  In the tomographic inversion, our objective function χ (m) with respect to the model space vector m can be expressed as: where G denotes a sensitivity matrix for the model parameters; d denotes a data vector including travel time residuals of local earthquakes and relative residuals of teleseismic events; λ and L denote the coefficient of smoothing regularization and 3-D finite-difference approximation to the Laplacian operator applied over all model grid nodes (Lees & Crosson, 1989), respectively.
In our previous studies (Liang, Zhao, Hua, & Xu, 2022;Zhao et al., 2012), the smoothing term λ(L • m) is expressed as the difference between the ith target grid node and weighted average of six adjacent grid nodes (northern, southern, western, eastern, top, and bottom) with an index of j: where w ij denotes the distance between the ith and jth grid nodes.The optimal λ value can be obtained by plotting an L-shaped tradeoff curve between the root-mean-square (RMS) of residuals and norm of the 3-D Vp model (Figure S3 in Supporting Information S1).To solve Equation 2, we adopt the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (L-BFGS-B; Morales & Nocedal, 2011) to minimize the objective function.The advantages of the L-BFGS-B algorithm are that it does not require damping regularization and it converges faster with a solution fitting the data slightly better than the LSQR algorithm (Jia & Zhao, 2023;Z. Wang & Zhao, 2019).

Improving the Tikhonov Regularization
Regularization is necessary for solving a geophysical inverse problem (e.g., Menke, 1989), in which a smoothness constraint to the solution is widely applied to stabilize ill-posed inverse problems and avoid overfitting noise or small-scale anomalies in the data.This is typically achieved by minimizing a least-squares objective function to a regularization term that penalizes high-frequency components in the solution, for example, the Tikhonov regularization as shown in Equation 2 (Menke, 1989;Tarantola, 2005;Tikhonov et al., 1995).This approach often leads to a biased and oversmoothed model that lacks sharp features near the boundary of anomalies.Hence, some studies adopted an edge-preserving regularization based on  1 norm such as the total variation regularization (W.Jiang & Zhang, 2018;Lin et al., 2015) and wavelet-based regularization (Loris et al., 2007), but these schemes have not been used in large-scale tomographic studies yet.Alternatively, some studies, on the basis of the Tikhonov regularization (  2 norm), adopted two sets of grids above and below the subducting slab interface and added regularization terms to each set of grids separately (Tian & Zhao, 2012;Zhao et al., 2012), which can also obtain an edge-preserving model near the slab interface.Another problem of regularization is that its parameters need to be set artificially, which may introduce some inevitable errors.The optimal regularization parameter is commonly chosen by plotting an L-shaped trade-off curve between the RMS residual and model norm (e.g., Eberhart-Phillips, 1986;Hansen, 1992;Zhao et al., 1992), though many studies adopt other choosing criteria (e.g., Bayesian framework; Aster et al., 2018;Valentine & Sambridge, 2018).A recent study showed that the solution pattern does not change much for different damping or smoothing parameters if those parameters are in an appropriate range and the tomographic inversion is stable (Gou et al., 2020).Hence, the conventional Tikhonov regularization is still reliable and even plays an irreplaceable role in seismic tomography.
However, in the Alaska subduction zone, the subducting Pacific slab has a sharp upper boundary generating converted seismic waves and other later phases (e.g., Abers, 2005;Ai et al., 2005).The influence of the sharp boundary on the smoothing regularization is poorly understood in tomographic inversions.In addition, there is a big gap in the data volume between different areas (Figure 3) with the outbreaking data recorded at the USArray stations in recent years.Although the regularization can stabilize the minimization process for the heterogeneous ray distribution, accompanied model imbalance in different areas would be more pronounced when the ray density varies in several orders of magnitude (Figure 4a).This is because the sensitivity matrix highly depends on the seismic ray density and paths.The areas with a small ray density are often underfitted (Figure 4a) because their model parameters are less sensitive to the objective function, whereas those areas with a large ray density are often overfitted because their model parameters are more sensitive to the objective function.The regularization holds promise to resolve the model imbalance because it can adjust the sensitivity of model parameters to the objective function.
To overcome these shortcomings of the normal Tikhonov regularization, we modify Equation 3 with three schemes (see Appendix A for details), in which the ray density variation and the sharp slab interface are taken into account in the new regularization (Figure S4 in Supporting Information S1). Figure 5 shows the raw and improved Vp tomographic models.Their differences near the Pacific slab upper boundary are expected, indicating that the amplitudes of the high-Vp slab and low-Vp mantle wedge are underestimated.But some differences also appear in the marginal and deep portions of the study volume where the seismic rays are less dense and do not crisscross well.The low-Vp anomalies in the mantle wedge and subslab mantle have <0.4% differences in Vp perturbation, suggesting that the new regularization has little influence on those areas with denser rays.Compared with the raw Vp model (Figure 5a), the improved model (Figure 5b) shows a clear high-Vp anomaly in the mantle transition zone (MTZ) beneath western Alaska, which demonstrates that the prior information in regularization helps to recover the deep structures in the study region.To investigate the robustness of the tomographic model obtained with the new regularization, we performed a synthetic test (Figure 6).The input model contains two connected high-Vp anomalies to represent the subducting slab in the upper mantle and a flat slab in the MTZ, and two low-Vp anomalies in the mantle wedge and subslab mantle.Then we calculated synthetic travel times for the input model with the same event-station pairs as those in the real data set.Random noise (−0.15 to 0.15 s with a standard derivation of 0.1 s) is added to the synthetic data to simulate the data picking errors.Then we conducted tomographic inversions of the synthetic data set with the normal Tikhonov regularization and the new regularization to obtain output Vp models (Figures 6b-6e).The high-Vp Pacific slab in the upper mantle is clearly visible in these output models.The flat slab in the MTZ can be imaged clearly in the output model obtained with the regularization scheme C, less clearly in the output models with the regularization schemes A and B, but not clearly in the model with the normal regularization (see Appendix for details of the three regularization schemes).Figure 6g shows the difference between the output models obtained with the normal and new regularizations.The mantle wedge and the Pacific slab exhibit systematic negative and positive biases in the Vp amplitude, respectively.Hence, the ray density and sharp velocity boundary should be taken into account in the Tikhonov regularization for the tomographic inversion.

Isotropic Vp Result
Figures 7 and 8 show our final isotropic Vp model obtained by the inversion with the new regularization.In the crust, our isotropic Vp model agrees well with the surface geological features.For example, the Yukon Composite Terrane and the Wrangellia Composite Terrane show high and low Vp anomalies above 30 km depth (Figure 7), respectively.This feature is in line with a previous receiver function study (Miller et al., 2018), which shows an abrupt change in the crustal thickness beneath the Denali Fault.In addition, the Koyukuk Composite Terrane shows a high-Vp crust, whereas the Arctic Alaska Terrane and the Yukon Flats show a low-Vp crust.In the shallow mantle (30-100 km depths), the subducting Pacific slab and the Yakutat Terrane present a high-Vp feature, whereas the subducted Yakutat slab shows a low-Vp anomaly, in line with previous tomographic results (e.g., Eberhart-Phillips et al., 2006).A slab-like high-Vp anomaly appears at 50-100 km depths beneath the Yakutat Terrane and the Wrangellia Composite Terrane, which may represent shallow subduction of the Wrangell slab.Some high-Vp anomalies are also visible at 50-100 km depths beneath the Arctic Alaska Terrane and the Yukon Composite Terrane, which may reflect their lithospheric root.Below 100 km depth, a high-Vp anomaly along the Aleutian Arc is visible, which reflects the subducting Pacific slab.A high-Vp anomaly is also visible at the NE end of the Pacific slab, which may represent its aseismic slab edge (Gou et al., 2019).
In the MTZ (410-660 km depths; Figures 7 and 8), a large high-Vp anomaly appears beneath western Alaska, which may represent a flat slab.We deem that a big mantle wedge (BMW) occurs beneath western Alaska, similar to that in East Asia (Lei & Zhao, 2005;Liang, Zhao, Xu, & Hua, 2022;Zhao et al., 2004).The flat slab is connected with the subducting Pacific slab in its southern and eastern margins, whereas its western margin exceeds our study region (∼168°W longitude) and its northern margin reaches the Kaltag Fault (∼65°N latitude).Our synthetic test (Figure 6) indicates that the flat slab and the BMW beneath western Alaska are robust features.

Tilting-Axis Anisotropy
Because the anisotropic medium with a hexagonal symmetry axis requires five independent elastic moduli, some simplifying assumptions are made by considering the resolution of available seismic data.For example, azimuthal and radial anisotropies assume horizontal and vertical hexagonal symmetry axes (HSAs), respectively, which have been taken by most studies of anisotropic tomography (see Zhao et al., 2023 for a detailed review).However, in the Alaska subduction zone, the dipping subducted slab and 3-D mantle flow around it may cause tilting HSAs, which have been verified by recent studies of anisotropic tomography in Japan (Z.Wang & Zhao, 2021;Z. Wang et al., 2022).The abundant seismic data in Alaska (Figures 1 and 2) enable a tomographic inversion for more complex anisotropies.
Here we take a more general assumption in which the HSA is orientated freely in 3-D space.The tilting-axis anisotropy contains a tilting HSA and an isotropic plane perpendicular to it, so it is also called tilted transverse isotropy.For deformed peridotites or blueschists, their anisotropy is characterized by one fast axis along the lineation and two slow axes perpendicular to it (Nicolas & Christensen, 1987), so we take the fast HSA assumption (Figure 9a).In this case, a fast HSA tomographic model with dipping fast-velocity directions (FVDs) is determined.But in the subduction zone, abundant foliated serpentinite and faults within the slab and mantle wedge exhibit one slow axis perpendicular to the foliation or fault plane (Faccenda et al., 2008;Reynard, 2013).Hence, we also take the slow HSA assumption to obtain a slow HSA model, in which a fast velocity plane (FVP) is perpendicular to the slow HSA (Figure 9b).As a matter of fact, olivine, the most abundant mineral in the upper   , 2006).The blue and red lines denote the slab interface (Gou et al., 2019) and major faults (Colpron et al., 2007), respectively.AAT, Arctic Alaska Terrane; KY, Koyukuk Composite Terrane; NAC, North American Craton; WCT, Wrangellia Composite Terrane; YCT, Yukon Composite Terrane; YF, Yukon Flats. Figure 8. (a-g) Vertical cross-sections of isotropic Vp tomography along seven profiles as shown on the inset map (h).The blue and red patches with black contours denote high and low Vp perturbations, respectively, whose scale is shown at the bottom.In each panel, three black dashed lines denote the Moho, 410 and 660 km discontinuities.The red dots denote seismicity within a 25 km width of each profile.The thick blue line denotes the Pacific slab upper boundary (Gou et al., 2019).The gray patch above each panel denotes the surface topography along each profile.The vertical blue bars and red triangles denote major faults (Colpron et al., 2007) and Cenozoic volcanoes, respectively.BRF, Border Ranges Fault; CMF, Castle Mountain Fault; DF, Denali Fault; DVG, Denali volcanic gap; INF, Iditarod-Nixon Fork Fault; KAF, Kaltag Fault; KOF, Kobuk Fault; TF, Tintina Fault; WCT, Wrangellia Composite Terrane; YCT, Yukon Composite Terrane.(Ismaıl & Mainprice, 1998;Karato et al., 2008), so both assumptions of the slow and fast HSAs are reasonable.In the following figures, we adopt a T-shaped symbol to express the 3-D FVP or FVD projected on the paper plane (Figures 9c and 9d).

Inversion Strategy
We use the tomographic method of Z. Wang and Zhao (2021) to obtain a 3-D Vp tilting-axis anisotropy model beneath Alaska.In a weak anisotropy medium, P-wave slowness s is expressed as: where s 0 denotes isotropic slowness; M denotes a parameter quantifying the anisotropic amplitude; ϕ denotes the angle between the ray propagation direction and the tilting HSA.When M > 0, Vp along the HSA is slowest and the anisotropic medium has a slow HSA normal to the FVP.In contrast, M < 0 represents a fast HSA corresponding to a 3-D FVD.In the Cartesian coordinate (Figure S5 in Supporting Information S1), we adopt three anisotropic parameters (A, B, and C) to express the projection of a titling HSA in the east, north, and depth directions: where q and γ are dip and azimuthal angles of the tilting HSA, respectively.Given that a seismic ray has an azimuthal angle of μ and a dip angle of p, the ray directional vector (d, e, f) can be written as: The included angle ϕ between the ray directional vector and the tilting HSA vector is expressed as: cos = sin sin sin sin + sin cos sin cos + cos cos (7) Thus, P-wave anisotropic amplitude   and slowness s can be expressed as: Most parameters of anisotropic tomography are the same as those of the isotropic tomography as listed in Table 2.
When expressing the observation equation with the first-order Taylor expansion like Equation 1, the linear inverse problem would fall into the trap of locally optimal solution (the isotropic model) because the first-order partial derivatives of travel time with respect to anisotropic parameters are zero according to Equation 8. Hence, we set A = B = C = 0.01 in the starting anisotropic model and adopt the second-order Taylor series to calculate theoretical travel times (Z.Wang et al., 2022).The observation equation for tilting-axis anisotropy tomography can be written as: where m n denotes the nth structural parameter in our model space including isotropic slowness parameter s and anisotropic parameters A, B, and denotes the second-order partial derivatives of travel time with respect to the structural parameters; other symbols are the same as those in Equation 1.The new regularization (scheme C in Appendix) and the L-BFGS-B algorithm are utilized for constructing the objective function and conducting the nonlinear inversion, respectively.The smoothing slope parameter for anisotropy is set to zero according to the L-shaped tradeoff curve (Figure S3 in Supporting Information S1), so the slab interface geometry is taken into account in the regularization.

Recovery of a Synthetic Subduction Zone Model
We assess our tilting-axis anisotropic tomography technique by detecting complex mantle flow patterns induced by a subducted slab.The input model (Figure 10) is a synthetic anisotropy model for a subduction zone proposed by VanderBeek and Faccenda (2021), which was derived from a comprehensive petrological-thermomechanical model extending to 1,000 km depth (for details of its mantle composition, creep mechanism, and fabric development, see Faccenda, 2014).The anisotropic model contains 21 elastic moduli at each grid node, but it has been simplified to a medium with a fast HSA with the lower-order symmetry components neglected (VanderBeek & Faccenda, 2021).In the model space, the subducting slab in the upper mantle and a flat slab in the MTZ induce tilting corner flow and toroidal flow.
To evaluate the performance of our tomographic technique, we employ 16 hypothetical teleseismic events from eight different directions, which are recorded at 770 synthetic, evenly spaced stations (Figure 10).Note that no local or regional earthquakes are included in the synthetic data set (VanderBeek & Faccenda, 2021).The theoretical travel time for the input subduction zone model is calculated by solving the Christoffel equation along the teleseismic ray path.The subslab mantle exhibits trench-parallel and subhorizontal FVDs.However, the amplitudes of isotropic Vp and anisotropy are underestimated to a large extent due to the regularization employed in the tomographic inversion.Although the shallow structure (<∼50 km depth) cannot be accurately resolved due to its small sensitivity to teleseismic data, our tomographic technique is effective in characterizing deep FVD structures (e.g., toroidal and poloidal flows).

3-D Anisotropy Under Alaska
Figures 12 and 13 show our anisotropic model in Alaska under a fast HSA assumption (hereafter we call it the fast HSA model).The isotropic Vp in the fast HSA model is similar to that in the isotropic tomography, whereas fast HSAs (corresponding to 3-D FVDs) can provide more detailed information than azimuthal anisotropy.In the crust (<30 km depth), FVDs in the Koyukuk Composite Terrane, the Wrangellia Composite Terrane, and the Yukon Composite Terrane are subhorizontal (<30° dip angles) and subparallel to surrounding faults including the Tintina Fault, the Denali Fault, and the Iditarod-Nixon Fork Fault.In the shallow mantle (30-100 km depths), the NW portion of the low-Vp Yakutat slab shows NW-SE trending FVDs, whereas the high-Vp Wrangell slab  Beneath the Aleutian arc, the mantle wedge above 100 km depth exhibits trench-normal FVDs with a <30° dip angle (Figures 12 and 13), which may reflect corner flow.In contrast, at 100-300 km depths, the FVDs in western Alaska are horizontal and orientated in the NE-SW direction, which may represent toroidal flow around the subducted Pacific slab.The asthenosphere beneath the Arctic Alaska Terrane has NE-SW trending FVDs with a ∼30° dip angle, which corresponds to the direction of absolute plate motion.The high-Vp anomaly beneath western Alaska at 400-500 km depths shows FVDs with a large dip angle (>30°), which is different from the surrounding low-Vp mantle with subhorizontal FVDs (<20° dip angle).This may reflect the MTZ fabrics due to the lattice-preferred orientation of wadsleyite because horizontal shearing in the MTZ can cause vertical FVDs (e.g., Faccenda, 2014).
We use our fast HSA model above 400 km depth to predict SKS splitting features and compare them with previous SKS splitting measurements.In the prediction, P and S wave anisotropies are assumed to have the same HSA with the lower-order symmetries neglected.The polarized direction and delay time of the SKS splitting are calculated using our multi-layered tilting-axis anisotropy model.Since the SKS splitting represents an integration along a subvertical ray path, the MSAT Matlab toolbox (Walker & Wookey, 2012) is used to synthesize SKS splitting parameters based on a vertical unsplit SV wave with a period of 8 s (Wei et al., 2019).Our predicted SKS splitting shows a circular pattern around the subducting Pacific slab edge (Figure 14a), which aligns well with the mantle flow directions predicted by Jadamec and Billen (2010) (Figure 14b).On average, the predicted SKS splitting shows an angle difference of 22.3° and 32.5° from the observed SKS splitting measurements by Venereau et al. (2019) and McPherson et al. (2020), respectively.At most stations in the continental interior, the angle difference between the predicted and observed FVDs is less than 20° (Figures 14c-14h).However, a large angle difference appears in southern and eastern Alaska.Meanwhile, our predictions tend to underestimate the delay time of SKS splitting probably due to the following reasons.(a) The lack of seismic station in the Pacific Ocean leads to inadequate data coverage, which affects the accuracy of our prediction.(b) The complex mineral composition in the subduction zone may introduce some second-order symmetries that influence the SKS splitting results.(c) The regularization adopted in the tomographic inversion would underestimate the anisotropic amplitude.This comparison verifies the robustness of our fast HSA model and demonstrates that our model reliably captures the anisotropic characteristics in the continental interior.
Figures 15 and 16 show our anisotropic model under a slow HSA assumption (hereafter we call it the slow HSA model).Although the high-Vp subducting Pacific slab in Alaska is less clear than that in the fast HSA model, the low-Vp anomalies in the mantle wedge and subslab mantle are almost the same as those in the fast HSA model.Nevertheless, we can still focus on the FVP features in the mantle wedge and the subducting slab because some minerals with a slow HSA (e.g., serpentine) exist in the subduction zone.At depths of 30-100 km, the mantle wedge is characterized by upright FVPs with their trends oblique to the slab interface, whereas the Pacific slab and the Yakutat slab exhibit trench-parallel and trench-normal FVPs, respectively.Down to 200-300 km depths, the mantle wedge FVPs change to subhorizontal, whereas FVPs of the Pacific slab change to trench-normal with a <30° dip angle.A similar change in azimuthal anisotropy of the Pacific slab at 100-200 km depths is also visible in the previous results (Tian & Zhao, 2012;J. Wang & Zhao, 2013).In the MTZ, the flat slab has tilting FVPs with a NW-SE or E-W trend and a dip angle of ∼40-90°.

Robustness of Anisotropic Tomography Under Alaska
To evaluate the stability of our model with respect to different starting values of anisotropic parameters (A, B, and C), we conducted several tomographic inversions using our real data set.Figure S6 in Supporting Information S1 displays the anisotropic results obtained with different starting values of the anisotropic parameters.When the starting value is >0.006 (corresponding to an anisotropic amplitude >0.01%), the tomographic image shows a consistent FVD pattern.Although the anisotropic amplitude depends on the starting value, the FVD pattern remains the same.Conversely, when the starting value is <0.006 (i.e., the anisotropic amplitude is <0.01%), the solution is trapped in a local optimum, resulting in an isotropic model.Hence, our anisotropic result is stable and consistent when the starting value is >0.006.
To examine the resolution of our slow and fast HSA models, we performed extensive checkerboard resolution tests.The input checkerboard model has a horizontal grid spacing of ∼30 or ∼50 km for the isotropic Vp component, while the spacing is ∼80 or ∼100 km for the anisotropic component.The isotropic Vp perturbation in our input model is 6% and −6%, whereas the anisotropic amplitude is 5% (see Table S1 in Supporting Information S1).When we calculate synthetic travel times for the input checkerboard models, we take into account the real event-station pairs (Figures 1 and 2), 3-D geometry of velocity discontinuities in the crust and mantle, and the data picking errors (−0.15 to 0.15 s Gaussian noise with a standard derivation of 0.1 s).Figures S7-S14 in Supporting Information S1 show the output checkerboard models of isotropic Vp and anisotropy.In the mantle, the resolution of our isotropic Vp model is ∼30 km beneath most areas in the continent and reaches 50 km in some deep areas beneath the continental shelf and ocean.The resolution of our anisotropic models is also high beneath the continent interior but low beneath the ocean and continental shelf due to the limited number of seismic rays there.Near the subducting slab, the isotropic Vp and anisotropy images have resolution of less than 30 and 80 km, respectively.To assess the influence of the data picking errors on the resolution of tomography, we performed three additional checkerboard resolution tests.These tests are the same as the above-mentioned checkerboard tests, but the synthetic travel times contain different levels of Gaussian noise with standard derivations of 0.2, 0.3, and 0.5 s.Figures S15 and  S16 in Supporting Information S1 show the corresponding output models for isotropic Vp and anisotropy, respectively.The resolution of the isotropic Vp component reaches ∼30 km even when the data set has picking errors with a standard deviation of 0.5 s (from −1.5 to 1.5 s).However, the input Vp anomaly amplitude cannot be recovered well.For the anisotropy component, the resolution decreases as the noise amplitude increases.For example, when the noise standard deviation in the synthetic data exceeds 0.3 s, the input anisotropic checkerboard can be recovered only under the dense seismic network.Therefore, a high-quality data set is critical to obtain a high-resolution and reliable anisotropic image.Because our data set has a root-mean-square travel-time residual of 1.28 s for the input checkerboard model, our anisotropic image has ∼80 km resolution when the signal-to-noise ratio of our data set is greater than 4.27.
Dip angle is an important component of the tilting-axis anisotropy, so we performed two checkerboard resolution tests to verify the dip angles of FVDs and FVPs (Figures S17 and S18 in Supporting Information S1).In the two input models, the isotropic Vp and anisotropy grids have the same lateral interval and depth setting as those of the formal inversion.We assigned subhorizontal (10° dip angle) and subvertical (80° dip angle) HSA with a NE trend and 5%  anisotropic amplitude to adjacent nodes of the input anisotropic grid, whereas we did not assign slowness perturbation to the input isotropic Vp grid.The test results show that the dip angles in the input models can be recovered well in most areas beneath the continent, indicating that the dip angle distribution is robust in our tilting-axis anisotropy.
Because our isotropic Vp and tilting-axis anisotropy parameters are coupled in the tomographic inversion (Huang et al., 2015;Zhao et al., 2023), it is necessary to evaluate their trade-off effect.In the first test, the obtained fast HSA model (Figures 12 and 13) is taken as the input model but the anisotropic parameters are set to zero.The isotropic Vp anomalies are normalized to a maximum perturbation amplitude of 2%.The synthetic travel time calculation and tomographic inversion are the same as those of the checkerboard resolution tests.Figure S19 in Supporting Information S1 shows the output models for isotropic Vp and anisotropy.The output anisotropic amplitude is less than 0.5% in most areas, and less than 1.0% in the subduction zone, so weak trade-off occurs between the isotropic Vp and anisotropy.In the second trade-off test, our fast HSA model (Figures 12 and 13) is taken as the input model but the isotropic Vp perturbations are set to zero.The anisotropic amplitudes are normalized to 2% but the FVD pattern is kept unchanged.As a result, the output isotropic Vp anomalies have amplitudes <1% (Figure S20 in Supporting Information S1), indicating that our isotropic Vp image is not an artifact induced by the anisotropic component.

Dip Angle of Tilting-Axis Anisotropy
Figure 17 shows the dip angle distribution of our fast HSA model at different depths, in which most areas exhibit neither horizontal nor vertical FVDs.Above 100 km depth, the rigid lithosphere in Alaska is characterized by subhorizontal FVDs that are distributed in all directions.In most areas, the dip angles of FVDs are smaller than 30°, which is likely a consequence of lithospheric compaction.Hence, in the case of the Alaskan lithosphere, the assumption of azimuthal anisotropy is more reasonable than that of radial anisotropy.In contrast, the mean dip angle of FVDs below 100 km depth is >30° with a predominant azimuth in the NE-SW direction.This may reflect a convective asthenosphere separated by the NE-trending dipping slab that led to trench-parallel flow.At depths of 200-500 km, subhorizontal FVDs appear in western Alaska, whereas eastern Alaska exhibits subvertical FVDs, whose boundary coincides well with the NE end of the subducting Pacific slab (i.e., the aseismic slab edge; Gou et al., 2019).The transition from horizontal flow to vertical flow may be controlled by the geometry of the subducting Pacific slab according to numerical modeling (Jadamec & Billen, 2010;Jadamec et al., 2013).
Figure 18 shows the dip angle distribution of FVPs in our slow HSA model.The mean dip angle is ∼57° at depths of 50-500 km with small variation in depth.Horizontal or vertical FVPs only appear near the slab interface.This result indicates a relatively homogeneous distribution of the FVP dip angle.But the azimuth of FVPs tends to NW-SE below 200 km depth, suggesting that the FVP azimuth is also sensitive to the direction of mantle flow.

Tilting-Axis Anisotropy and Azimuthal Anisotropy
The fast Vp azimuth can be reflected by both azimuthal anisotropy and tilting-axis anisotropy, but their relationship in azimuth has not been investigated.In theory, their azimuths should be consistent.But in practice, the relationship varies with some factors such as the dip angle, anisotropic amplitude and ray coverage.To investigate the azimuthal relationship between the tilting-axis and azimuthal anisotropies, we performed two synthetic tests.
In the first test, our fast HSA model is taken as an input model (Figure 19a) for which synthetic travel times are calculated with the same event-station pairs as those in the real data set.To simulate the data picking errors, we added Gaussian noise (−0.15-0.15s with 0.1 s standard derivation) to our synthetic data set.Then we conducted a tomographic inversion of the synthetic data set for azimuthal anisotropy (J.Wang & Zhao, 2013), in which the ray paths, cut-off residuals, inversion algorithm (L-BFGS-B), and regularization parameters are the same as those in the formal inversion for tilting-axis anisotropy.Figure 19b shows the output azimuthal anisotropy model of this synthetic test.The azimuth of the input tilting-axis anisotropy is similar to that of the output azimuthal anisotropy with a <20° direction difference.But some inconsistency appears at almost all depths, especially in those areas with a lower ray density such as western Alaska and the Pacific Ocean, where the input tilting-axis anisotropy has small amplitudes or large dip angles (Figure 19c).Therefore, azimuthal anisotropy can reflect the direction of nearly horizontal and strong mantle flow with a large anisotropic amplitude, but weak or nearly vertical mantle flow may bring errors to the azimuthal anisotropy model.The second synthetic test is similar to the first test, but our slow HSA model is taken as the input model (Figure 20a).The output 2-D FVDs show a consistent strike with the input FVPs in most areas with a large ray density.The areas near the Bering Sea and the Pacific Ocean (Figure 20b) show a large azimuthal difference between the input FVPs and output 2-D FVDs, probably due to the sparse station distribution and small amplitude projected on the horizontal plane as shown in Figures 9 and 12.In (b), the red and blue background colors denote small and large angular differences between the input 3-D FVD of tilting-axis anisotropy and the output 2-D FVD of AAN, respectively, whose scale is shown at the bottom.The length and orientation of black bars denote the AAN amplitude and its 2-D FVD, respectively.(c) Relationship among the input anisotropic amplitude (horizontal coordinate), input dip angle of FVD (vertical coordinate), and back-azimuthal difference between the input 3-D FVDs of tilting-axis anisotropy and output 2-D FVDs of AAN (color dots).The red and blue dot colors denote small and large direction differences between the input tilting-axis anisotropy and output AAN, whose scale is the same as that in panel (b).The dot size denotes the P-wave ray hit-count at each grid node, whose scale is shown at the lower-right corner of panel (c).
of input anisotropy there.Similar to the first synthetic test, the areas with large anisotropic amplitudes and subvertical FVPs show a high consistency in azimuth between the input FVPs and output 2-D FVDs, whereas those areas with subhorizontal FVPs or small anisotropic amplitudes are likely to exhibit FVP-normal FVDs (Figure 20c).Hence, the tilting-axis anisotropic tomography is much better than azimuthal anisotropy tomography in revealing mantle flow of subduction zones.

Subduction-Induced Mantle Flow
Subduction-induced mantle flows, including toroidal and poloidal flows, have been extensively studied as a main mechanism for explaining seismic anisotropy in subduction zones (e.g., Faccenda, 2014;Faccenda & Capitanio, 2012;Li et al., 2014).To understand the mantle flow induced by the subducted and flat slabs in Alaska, we compare our anisotropy model with previous numerical modeling results.Faccenda and Capitanio (2013) quantified the flow-induced anisotropy in the mantle, specifically related to a subducting plate sinking into the MTZ.Their numerical Model C is useful for understanding seismic anisotropy of the Alaska subduction zone due to the following reasons.(a) The curved subduction zone in their Model C, which spans 2,000 km wide, is comparable to the Aleutian-Alaska arc.(b) The subducting slab in their model, with a thickness of 50-80 km, is analogous to the subducting Pacific slab beneath Alaska.(c) They also simulate a flat slab in the MTZ due to its negative buoyancy, which is identical to the BMW structure we observed in Alaska.
Figures 21a and 21b show the single-layer FVD pattern and predicted SKS splitting (i.e., integrated anisotropy of all layers in the upper mantle) within the Alaska BMW, both of which are derived from our fast HSA model, whereas Figures 21c and 21d display corresponding results from the numerical Model C (Faccenda & Capitanio, 2013).A circular FVD pattern near the slab edge appears in both cases.This consistency suggests that our anisotropic model may reflect the toroidal flow induced by the subducting slab in the upper mantle and the flat slab in the MTZ.
In addition, subduction-induced poloidal flows, including corner flow in the mantle wedge and entrained flow in the subslab mantle, play an important role in mantle dynamics (e.g., Faccenda & Capitanio, 2013).In central Alaska, it is likely that the mantle near the slab edge is predominantly affected by toroidal flow.However, in areas further away from the slab edge, such as beneath southwestern Alaska and the Aleutian arc, corner flow and entrained flow may occur, as shown in the synthetic subduction zone model (Figure 11).In our anisotropic tomography, the mantle wedge and subslab mantle beneath the Aleutian arc exhibit tilting FVDs (Figures 13f  and 13g) and FVPs (Figures 16f and 16g) that are subparallel to the slab interface.This alignment suggests the development of corner and entrained flows away from the slab edge.However, this feature is invisible beneath southeastern Alaska probably due to the termination of the retreat and subduction of the Wrangell slab.

Big Mantle Wedge and Intraplate Volcanism
Earthquakes in the subducting Pacific slab occur down to ∼200 km depth beneath Alaska (Figure 1b), but several previous studies have shown that the Pacific slab has subducted down to a greater depth.Qi et al. (2007)  The present results are generally consistent with the previous findings.The Pacific slab has a small dip angle above ∼100 km depth but a large dip angle at 100-400 km depths with a kink beneath central Alaska (Figures 7 and 8).However, our isotropic Vp image shows a flat high-Vp anomaly in the MTZ that is connected with the subducting Pacific slab in the upper mantle (Figures 7 and 8).A recent study of teleseismic receiver functions has shown that the 410-km discontinuity in western Alaska is deeper than normal (Dahm et al., 2017), which indicates higher-temperature and water-rich upper MTZ there.Hence, we deem that this high-Vp anomaly in the MTZ represents a flat slab connecting with the subducting Pacific plate in the upper mantle, which forms a BMW beneath western Alaska.
The upper mantle structure beneath the Cenozoic volcanic fields in western Alaska is well resolved by our tomographic models (Figures S7-S20 in Supporting Information S1).The flat slab in the MTZ extends to ∼65°N latitude (Figures 7 and 8), so its presence provides an explanation for the low-Vp anomalies in the upper mantle beneath the St. Michael, Ingakslugwat, and Nunivak Island volcanic fields (Table 1).These low-Vp anomalies may reflect hot and wet upwelling flows in the Alaska BMW.It is worth noting that the hydrous upwellings can be driven by processes such as thermochemical convection (Long et al., 2019) or advection-diffusion-reaction (Sun et al., 2019), but their scales are too small to generate significant seismic anisotropy in our fast and slow HSA models.Therefore, the hydrous upwellings beneath the intraplate volcanoes do not contradict the observed subhorizontal FVDs in our anisotropic results (Figures 9 and 10).Instead, as discussed in Section 5.3, our anisotropic results primarily reflect large-scale toroidal flow related to the subducted Pacific slab, so it can provide volatiles to the upwellings and further cause partial melting and volcanic activities.In addition, the low-Vp anomaly beneath the Imuruk Lake volcanic field is less prominent and it exceeds the northern edge of the flat slab in the MTZ.But its upper mantle exhibits FVD and FVP patterns similar to those under the other volcanic fields in western Alaska, suggesting that similar mantle flow occurs beneath the Imuruk Lake volcanic field.
Previous studies of whole-mantle tomography have suggested that the subducted Pacific slab becomes flat in the MTZ beneath East Asia and western Alaska (Figures 22a and 22b) (Zhao et al., 2010(Zhao et al., , 2013)).Regional-scale tomographic models have revealed higher-resolution features of the BMW structure and dynamics (e.g., Lei & Zhao, 2005;Liang, Zhao, Xu, & Hua, 2022;Zhao et al., 2004).The length of the flat slab in the Alaska MTZ (∼800 km) is smaller than that under East Asia (>1,400 km) (Figures 22c and 22d).Notably, both East Asia and western Alaska have Cenozoic intraplate volcanoes, and their volcanic rocks are characterized by alkalic basalts and minor tholeiitic basalts with peridotite xenoliths (e.g., J. Liu et al., 2001;Moll-Stalcup, 1994).Some isotope studies have shown that these intraplate basalts have undergone metasomatism by carbonate melts derived from the deep mantle, which may be derived from subducted and recycled sediments in the MTZ (e.g., Moll-Stalcup, 1994;X.-J.Wang et al., 2017;Xu et al., 2018).These results are important for understanding the formation of intraplate volcanism in western Alaska, suggesting a potential involvement of a volatile-rich mantle source (Figures 22e and 22f).
The isotropic Vp features and presence of Cenozoic intraplate volcanism in the BMW are consistent in Alaska and East Asia, but their anisotropic features have some differences.First, radial anisotropy tomography of East Asia shows a faster vertical Vp than horizontal Vp (Liang, Zhao, Xu, & Hua, 2022), whereas western Alaska exhibits subhorizontal FVDs and FVPs (Figures 13 and 16).In addition, east-west variations in isotropic Vp and azimuthal anisotropy are revealed in the East Asia BMW (Liang, Zhao, Xu, & Hua, 2022;Ma et al., 2019), which are not identified in our Alaskan tomography.These inconsistencies may be attributed to differences in regional tectonics and flow pattern in the BMW.The East Asia BMW is far away from the slab edge, so the poloidal flow may be dominant there.In contrast, western Alaska is close to the slab edge, resulting in toroidal flow in the BMW (see Section 5.3 for details).Because our study region is only part of the Aleutian-Alaska subduction system, the anisotropic structure and mantle flow pattern beneath the Aleutian back-arc areas and the Bering Sea are not well resolved due to the lack of seismic data there.Nevertheless, the MTZ-derived fluids and hot upwelling in the BMW can still cause low-Vp anomalies and the intraplate volcanism in both East Asia and western Alaska, regardless of their specific mantle flow patterns.

The Yakutat Slab and the Denali Volcanic Gap
The Yakutat Terrane, at the NE corner of the Pacific plate, has been accreted to the North American continent since the Cenozoic (Plafker & Berg, 1994).Offshore seismic refraction profiles (Christeson et al., 2010) (Zhao et al., 2013) and (c), (d) regional tomography of the big mantle wedge (BMW) along two profiles in East Asia (Liang, Zhao, Xu, & Hua, 2022)  revealed a dramatic change in the Moho depth from ∼32 km under the Yakutat Terrane to ∼11.5 km under the Pacific Ocean, suggesting that the Yakutat Terrane may originate from an oceanic plateau (Ferris et al., 2003).Its subducted portion, the so-called "Yakutat slab," represents a leading edge of the Pacific plate as it moves down to the mantle beneath the North American continent.Eberhart-Phillips et al. (2006) detected the Yakutat slab, with a thick low-V and high Vp/Vs crust, subducting northwestward to the Buzzard Creek-Jumbo Dome volcanoes (Figures 1b and 23).However, a receiver function study (Mann et al., 2022) revealed another portion of the subducted oceanic plateau (the northern portion of the Yakutat slab) underlying northward to the Wrangell volcanic field, which is bounded by a slab tear with the NW portion of the Yakutat slab inferred by Eberhart-Phillips et al. (2006).A few other studies also inferred that the two portions of the Yakutat slab may be connected at shallow depths (C.Jiang et al., 2018;Wech, 2016), whereas a slab gap or weak zone may occur between them (Gou et al., 2019;X. Yang & Gao, 2020).In the following, the NW and northern portions of the subducted Yakutat Terrane are called the Yakutat slab and the Wrangell slab, respectively, whereas the gap between them is called the Yakutat-Wrangell slab gap hereafter (Figure 23).
The Yakutat slab, despite being older than the Pacific slab, exhibits similar anisotropic features to the Pacific slab.Both slabs display subhorizontal FVDs and subvertical FVPs that tend to the NW-SE direction (Figures 12,13,15,and 16).One possible explanation of this feature is frozen-in anisotropy within the oceanic lithosphere that formed at the mid-ocean ridge.This feature has been observed in the Nazca plate beneath South America (Eakin et al., 2016).Another potential source of anisotropy is the alignment of cracks within the slab due to deformation during its subduc- tion.This alignment can be influenced by slab bending in the forearc areas (Faccenda et al., 2008;Z. Wang et al., 2022).Hence, our anisotropic results support the existence of the Yakutat slab and its deformation during subduction.
The Denali volcanic gap (DVG) in central Alaska, with a width of ∼400 km, is located at the NE end of the Alaska-Aleutian arc.One hypothesis for its absence of volcanic activity is that the subducted Yakutat terrane cools the system and replaces the hot core of mantle wedge (Rondenay et al., 2010).Alternatively, the shallow subduction of the Yakutat oceanic plateau can also explain the DVG because fluids from the subducted oceanic plateau are trapped in its upper crust and such a shallow hydromechanical condition cannot generate melt (Chuang et al., 2017).Furthermore, a combined mechanism is proposed, that is, the low fluid content and low temperature work together due to relative isolation from asthenospheric circulation (Martin-Short et al., 2018).Our isotropic Vp model, in line with previous Vs tomographic models beneath the DVG (Berg et al., 2020;Feng & Ritzwoller, 2019;Martin-Short et al., 2018), shows a high-Vp feature above the subducted Pacific slab (Figure 8d), which supports the above-mentioned mechanisms to some extent.However, our anisotropic model shows that the high-Vp mantle wedge beneath the DVG has NW-SE trending subhorizontal FVDs (Figure 13d) and subvertical FVPs (Figure 16d), which may reflect a remnant of the subducted Yakutat slab (Figure 23).This is because its anisotropy is different from the surrounding mantle wedge with NE-SW FVDs but resembles the subducting Pacific slab above 100 km depth.Therefore, we deem that a remnant of the Yakutat slab underlying the continental crust causes a cool and wedge-shaped fluid-barrier, which further results in the lack of melt beneath the DVG.
The Buzzard Creek and Jumbo Dome volcanoes are located in the northeast of the DVG (Figure 1a and Table 1).Different from the DVG, the two volcanoes have a low-Vp mantle wedge above a high-Vp subducted slab (Figure 8c).This typical subduction structure can explain the presence of the Buzzard Creek and Jumbo Dome volcanoes, which is also supported by a recent 3-D Vs model (X.Yang & Gao, 2020).In addition, they are close to the Yakutat slab edge and the aseismic Pacific slab edge, so toroidal flow (Figure 21) around the slab edges may contribute to the formation of the two volcanoes.

The Wrangell Slab and the Wrangell Volcanic Field
The Wrangell volcanic field (WVF; Figure 1a and Table 1) in SE Alaska has been active since ∼25 Ma ago (Richter et al., 1990) with its last eruption in 1912, but its formation mechanism is still in hot debate.The WVF is located at the Yakutat slab edge (Eberhart-Phillips et al., 2006), so it likely originates from quasi-toroidal flow or partial melting related to the subducted Yakutat slab (Berg et al., 2020;Martin-Short et al., 2018).But many tomographic studies (Feng & Ritzwoller, 2019;Gou et al., 2019;C. Jiang et al., 2018;Y. Wang & Tape, 2014;X. Yang & Gao, 2020) have revealed a slab-like high-V zone (the so-called Wrangell slab) beneath a low-V mantle wedge under the WVF.A low-V accreted sedimentary wedge induced by the subducted Wrangell slab is also visible in some tomographic images (Estève et al., 2021;Ward & Lin, 2018).In addition, the seismicity beneath the WVF outlines a Wadati-Benioff zone that reaches ∼120 km depth (Page et al., 1989;Stephens et al., 1984).Our model also shows the subducted Wrangell slab beneath the Wrangellia Composite Terrane with a flat shape, whose northern edge reaches ∼150 km depth and ∼63°N latitude (Figures 7 and 8b).The WVF is located above the NW corner of the Wrangell slab, indicating that the Wrangell slab subduction contributes to the WVF formation.
Our anisotropic tomography provides new insight into mantle flow beneath SE Alaska.In our fast HSA model, the mantle anisotropy beneath the WVF is characterized by subvertical FVDs with a NW-SE trend above 200 km depth (Figures 12 and 13b), just corresponding to the Yakutat-Wrangell slab gap (Figure 23).This feature may reflect local upwelling induced by the subduction, as predicted by 3-D numerical modeling for Alaska (Jadamec & Billen, 2010;Jadamec et al., 2013).The modeling shows that the Wrangell slab subduction increases the slab-pull driving force and further results in upwelling beneath the WVF in a vertical plane and anticlockwise toroidal flow in a horizontal plane, which explains the WVF formation and other first-order tectonic features.Gou et al. (2019) also found NEE-SWW trending FVDs above 100 km depth beneath the WVF, which is attributed to feeding from the subslab mantle.Therefore, hot mantle upwelling around the Wrangell slab edge (i.e., the Yakutat-Wrangell slab gap) may have contributed to the WVF formation.
The Prindle volcano and Mount Churchill (Figure 1a and Table 1) are located ∼200 km northeast and southeast of Mount Wrangell, respectively.Our isotropic Vp model shows that Mount Churchill is located above the Wrangell slab, whereas the Prindle volcano lies above the northern edge of the Wrangell slab (Figure 8a).Our anisotropic model shows subvertical FVDs with an NWW-SEE trend beneath the two volcanoes, which may reflect hot mantle upwelling triggered by the subducted Wrangell slab (Figure 12a).Hence, the Wrangell slab subduction may have caused not only the WVF but also the Prindle volcano and Mount Churchill in SE Alaska.

Conclusions
We use a large number of seismic data accumulated in the past 50 years to obtain high-resolution tomographic models for isotropic P-wave velocity (Vp) and tilting-axis anisotropy down to 800 km depth beneath Alaska.The major findings of this work are summarized as follows.
(1) In Alaska, the fast velocity directions of P-wave anisotropy are subhorizontal (mean dip angle <30°) in the continental lithosphere but tilting (mean dip angle >30°) in the convective asthenosphere.For a large study region with a highly inhomogeneous ray distribution, strong smoothing (i.e., a large λ value in Equations 2 and 3) is required to stabilize the minimization process, though it often leads to an underfitting model in those less sampled areas.In contrast, a small λ value helps to obtain an edge-preserving model with small-scale anomalies, but it may cause instability in the tomographic inversion (e.g., Zhao et al., 2012).The Alaska subduction zone has a very inhomogeneous distribution of rays, so we adopt a smoothing coefficient λ depending on the ray density N (hitcount, the number of rays around a grid node): where λ ref and N ref denote the smoothing coefficient and ray density, respectively, at a grid node (here N ref = 1,000); λ slope denotes the influence of ray density on the smoothing coefficient, that is, the smoothing slope.Compared with the normal scheme with a constant smoothing value λ ref , we introduce one more smoothing parameter λ slope , which can be determined from an L-shaped tradeoff curve between the RMS travel-time residual and the model norm (Figure S3 in Supporting Information S1).A detailed distribution of the smoothing coefficient in Equation A1 is shown in Figure S4 of the Supporting Information S1.
Figure A1 shows raw and improved Vp tomographic models that are obtained with the constant and hitcount-dependent smoothing regularizations, respectively.As a whole, the two images are similar, but there are some differences in the edge and deep portions where the seismic rays are less dense and do not crisscross well.The low-Vp anomalies in the mantle wedge and the subslab mantle have <0.4% differences in Vp perturbation, indicating that this scheme has little influence on the areas with denser rays.Compared with the raw model (Figure A1a), the improved model (Figure A1b) shows a clearer slab-like high-Vp anomaly in the MTZ beneath western Alaska.

A2. Scheme B: Smoothing Matrix for a Sharp Slab Interface
To retain the sharp feature of the upper boundary of the subducting Pacific slab, we adopt a Vp model with a predefined slab interface (Gou et al., 2019) as prior information for Equation 3. The study volume is divided into two parts by the slab interface, and in each part the normal Tikhonov regularization is applied separately.In this case, the high-Vp slab is isolated from the low-Vp mantle wedge in Equation 3, which is changed to: ( ⋅ ) = Δs Figure A2 shows two models that are obtained with the normal Tikhonov regularization without and with the sharp slab interface.Their difference near the slab interface is expected because the constraints on Vp gradient in six directions are partially canceled out near the slab interface in this scheme.As a result, both the high-Vp slab and low-Vp mantle wedge exhibit higher amplitudes (>3%), suggesting that their perturbations are underestimated by the tomographic inversion with the normal Tikhonov regularization.However, the high-Vp anomaly in the MTZ is unexpectedly enhanced when the sharp slab interface is taken into account in the regularization (Figure A2).Although this scheme changes the smoothing matrix near the slab interface above 400 km depth, the deep structures in the MTZ are also changed probably due to a complex interplay of the structural parameters.Hence, the presence of a sharp edge can also influence other areas in a tomographic image.But the obtained tomographic images depend highly on the robustness of the prior constraint.

A3. Scheme C: Improved Regularization Combining Schemes A and B
In scheme C, we combine schemes A and B to obtain a new 3-D Vp model with the same data set.In the objective function (Equation 2), the smoothing matrix is the same as Equation A2 (Scheme B) but its smoothing weight λ is the same as that in Equation A1 (Scheme A).Our improved Vp model (Figure 5b) not only presents a more balanced distribution in perturbation with hitcount than the raw Vp model (Figure 5a) but also keeps the sharp slab interface.In the formal tomographic inversions (Figures 7 and 8), this scheme C is adopted for the smoothing regularization.

Figure 1 .
Figure1.(a) Tectonic setting of the Alaska region.The background colors denote elevation whose scale is shown at the bottom-right.The red triangles and red circles denote locations of Cenozoic volcanoes and major volcanic fields, respectively.The blue solid and dashed lines denote major faults and terrane boundaries, respectively(Colpron et al., 2007).The white arrows denote motion directions of the Pacific plate relative to the North American plate(DeMets et al., 1994).AAT, Arctic Alaska Terrane; BC-JD, Buzzard Creek-Jumbo Dome volcanic field; BRF, Border Ranges Fault; CMF, Castle Mountain Fault; DVG, Denali volcanic gap; INF, Iditarod-Nixon Fork Fault; KAF, Kaltag Fault; KOF, Kobuk Fault; KY, Koyukuk Composite Terrane; SEF, Sevier Fault; WCT, Wrangellia Composite Terrane; YCT, Yukon Composite Terrane; YF, Yukon Flats.(b) Distribution of 19,157 local earthquakes used in this study.The dot colors denote focal depths whose scale is shown at the bottom-left.The blue lines denote depth contours of the upper boundary of the subducting Pacific slab(Gou et al., 2019).The green solid and dashed lines outline the inferred Yakutat slab(Eberhart-Phillips et al., 2006) and the Wrangell slab(Gou et al., 2019;Mann et al., 2022), respectively.The yellow circle denotes an inferred aseismic slab edge(Gou et al., 2019).(c) Distribution of 294 permanent network stations (blue squares) and 393 transportable seismic stations (green squares) used in this study.

Figure 2 .
Figure 2. (a) Distribution of 25,749 teleseismic events used in this study.The dot colors denote focal depths whose scale is shown below the map.The red box and blue lines show the present study region and tectonic plate boundaries, respectively.(b) Corrected relative travel-time residuals (RTTRs) for the whole-mantle heterogeneity versus raw RTTRs.The colors denote ray density whose scale is shown at the bottom-left.(c) Distribution of raw and corrected RTTRs at each seismic station for all the teleseismic events in (a).The left and right parts of each circle denote the raw and corrected RTTRs, respectively, whose scale is shown at the bottom-right.

Figure 3 .
Figure 3. (a) A vertical cross-section showing the distribution of local rays (red lines) and teleseismic rays (blue lines) within a 25 km width of an E-W profile along 60° north latitude.Black contour lines denote the number of rays per grid node along the profile.The gray patch at the top denotes the surface topography along the profile.Two blue thick lines denote the upper and lower boundaries of the subducting Pacific slab.The black dots denote local seismicity within a 25 km width of the profile.The red triangles denote Cenozoic volcanoes.(b, c) 3-D views of local and teleseismic rays in two cells as shown in (a).

Figure 4 .
Figure 4. (a) Relationship between the Vp perturbation (dVp) and hit-count in our raw Vp model.(b) Histogram of Vp perturbations of our raw Vp model.(c, d) The same as (a), (b) but for our improved Vp model.The colors in (a), (c) denote the density of the Vp perturbation, whose scale is shown at the bottom.See the text for details.

Figure 5 .
Figure 5. Vp images obtained by inversions with (a) the normal regularization and (b) the new regularization.(c) Difference between the images shown in (a) and (b).The upper two rows are map views of the Vp images at 120 and 450 km depths.The lower two rows denote vertical cross-sections along the AA′ and BB′ profiles as shown on the map.In each vertical cross-section, the three dashed lines denote the Moho, 410 and 660 km discontinuities.The black dots denote local seismicity within a 25 km width/depth of each profile/map.The gray patch atop each panel denotes the surface topography.The red triangles and blue lines denote Cenozoic volcanoes and the Pacific slab upper boundary, respectively.MTZ, the mantle transition zone.

Figure 6 .
Figure 6.(a) A vertical cross-section showing the input model of a synthetic resolution test for Vp tomography.(b) Inversion result with the normal Tikhonov regularization.(c-e) Inversion results with new regularization of schemes A, B, and C as described in the Appendix.(f) Difference between the Vp models shown in (b) and (e).The profile location is the same as that in Figure 5 (B-B′).Other symbols are the same as those in Figure 5.

Figure 7 .
Figure 7. Map views of isotropic Vp tomography at different depths as shown at the lower-right corner of each map.The blue and red colors denote high and low Vp perturbations, respectively, whose scale is shown at the bottom.The red triangles denote Cenozoic volcanoes.The green line at 30-100 km depths outlines the inferred Yakutat slab (Eberhart-Philips et al., 2006).The blue and red lines denote the slab interface(Gou et al., 2019) and major faults(Colpron et al., 2007), respectively.AAT, Arctic Alaska Terrane; KY, Koyukuk Composite Terrane; NAC, North American Craton; WCT, Wrangellia Composite Terrane; YCT, Yukon Composite Terrane; YF, Yukon Flats.

Figure 9 .
Figure 9. P-wave anisotropy with (a) fast and (b) slow hexagonal symmetry axes (HSAs).(c) Projection of a fast velocity direction (FVD) on the horizontal and vertical paper planes under a fast HSA assumption.(d) Projection of a fast velocity plane (FVP) on the horizontal and vertical paper planes under a slow HSA assumption.

Figure 10 .
Figure 10.The synthetic subduction zone model from VanderBeek and Faccenda (2021).(a) Map view of isotropic Vp model at 150 km depth.The blue and red colors denote high and low Vp perturbations, respectively, whose scale is shown at the bottom.White triangles denote locations of 770 synthetic stations.(b) Distribution of 16 synthetic teleseismic events.The red rectangle denotes the area with the 770 synthetic stations.(c) Vertical cross-section of the isotropic Vp model along 0° latitude.

Figure 11
Figure 11 shows the output anisotropic tomography obtained with the above-mentioned inversion strategy and the parameters specified in Table 2.The output Vp model exhibits the high-Vp subducting slab in the upper mantle and the flat slab in the MTZ.The output anisotropy also recovers important features of the synthetic subduction zone: (a) The frozen-in anisotropy in the subducting slab is characterized by trench-normal and tilting FVDs with a significant anisotropic amplitude (∼3%).(b) Entrained flow is visible in the mantle wedge near the middle of slab.(c) Toroidal flow is developed around slab edges at 150 km depth but less prominent at 350 km depth.(d)The subslab mantle exhibits trench-parallel and subhorizontal FVDs.However, the amplitudes of isotropic Vp and anisotropy are underestimated to a large extent due to the regularization employed in the tomographic inversion.Although the shallow structure (<∼50 km depth) cannot be accurately resolved due to its small sensitivity to teleseismic data, our tomographic technique is effective in characterizing deep FVD structures (e.g., toroidal and poloidal flows).

Figure 11 .
Figure 11.Map views showing the (a) input and (b) output synthetic subduction zone model at 150 km depth.The blue and red colors denote high and low Vp perturbations, respectively, whose scale is shown at the bottom-right.The black bars denote projections of the fast velocity directions on the paper plane, whose length denotes the anisotropic amplitude with their scales shown at the bottom-left and bottom.(c, d) The same as (a), (b) but at 350 km depth.(e-h) The same (a), (b) but E-W vertical cross-sections along 0° and 4°S latitudes.

Figure 12 .
Figure 12.Map views of Vp anisotropic tomography (the fast HSA model) under an assumption that the fast hexagonal symmetry axis is orientated freely in 3-D space.The layer depth is shown at the lower-right corner of each map.The blue and red colors denote high and low isotropic Vp perturbations, respectively, whose scale is shown at the bottom-left.The T-shaped symbols denote 3-D fast velocity directions (FVDs) projected on a horizontal plane, whose length denotes anisotropic amplitude with its scale shown at the bottom-left.The direction and curvature of a T-shaped symbol denote the back-azimuth and dip angle of a 3-D FVD, respectively, whose scale is shown at the bottom-right.The red triangles denote Cenozoic volcanoes.The green line outlines the inferred Yakutat slab (Eberhart-Philips et al., 2006).The blue and red lines denote the Pacific slab upper boundary (Gou et al., 2019) and major faults, respectively.APM, absolute plate motion; WCT, Wrangellia Composite Terrane; YCT, Yukon Composite Terrane.

Figure 13 .
Figure 13.(a-g) Vertical cross-sections of Vp anisotropic tomography as shown in Figure 12.The profile locations are shown on the inset map (h).The blue and red colors denote high and low isotropic Vp perturbations, respectively, whose scale is shown at the bottom-left.The T-shaped symbols denote 3-D fast velocity directions (FVDs) projected on each vertical cross-section, whose length denotes anisotropic amplitude with its scale shown at the bottom-left.Other symbols are the same as those in Figure 8.

Figure 14 .
Figure 14.(a) Map showing the predicted SKS splitting using our fast HSA model.The orientation and length of blue bars denote the fast direction and delay time of SKS splitting, respectively, whose scale is shown at its lower-right corner.(b) The mantle flow pattern in the Alaska subduction zone from Jadamec and Billen (2010).(c) Comparison of predicted (blue bars) and observed (red bars; Venereau et al., 2019) SKS splitting measurements.The white and yellow colors denote small and large angle differences, respectively, whose scale is shown in (e).Histograms of (e) angle difference and (f) delay time difference between the predicted and observed SKS splitting shown in (c).(d, g, h) The same as (c, e, f) but for the observed SKS splitting from McPherson et al. (2020).

Figure 15 .
Figure 15.The same as Figure 12 but for the slow HSA model under an assumption that the slow hexagonal symmetry axis is orientated freely in 3-D space.The T-shaped symbols denote 3-D fast velocity planes (FVPs) projected on a horizontal plane, whose length denotes anisotropic amplitude with its scale shown at the bottom-left.In each T-shaped symbol, the direction of the longer bar and length of the shorter bar denote the back-azimuth and dip angles of a 3-D FVP, respectively, whose scale is shown at the bottom-right.

Figure 16 .
Figure16.The same as Figure13but for the slow HSA model under an assumption that the slow hexagonal symmetry axis is orientated freely in 3-D space.The T-shaped symbols denote 3-D fast velocity planes (FVPs) projected on the profile, whose length denotes anisotropic amplitude with its scale shown at the bottom-left.

Figure 17 .
Figure 17.(a) Map views, (b) histograms of dip angle distribution, and (c) stereographic projections of tilting fast hexagonal symmetry axes (FHSAs) at different depths.The layer depth is shown at the lower-right corner of each map in (a) and the upper-right corner of each panel in (b).The blue and red colors denote steep and flat FHSAs, respectively, whose dip angle scale is shown at the bottom.The vertical black bar in each panel in (b) denotes the mean dip angle of the FHSAs at each depth.The blue line in (a) denotes location of the slab upper boundary at each depth.

Figure 18 .
Figure 18.The same as Figure 17 but for (a) map views, (b) histograms of dip angle distribution, and (c) stereographic projections of the fast velocity planes (FVPs) in our slow HSA model.

Figure 19 .
Figure 19.(a) Input tilting-axis anisotropy (the fast HSA model) and (b) output azimuthal anisotropy (AAN) of a synthetic test.In (a), the yellow and gray colors denote small and large amplitudes of tilting-axis anisotropy, respectively, whose scale is shown below (a).The T-shaped symbols denote fast velocity directions (FVDs) projected on the horizontal plane as shown in Figures 9 and 12.In (b), the red and blue background colors denote small and large angular differences between the input 3-D FVD of tilting-axis anisotropy and the output 2-D FVD of AAN, respectively, whose scale is shown at the bottom.The length and orientation of black bars denote the AAN amplitude and its 2-D FVD, respectively.(c) Relationship among the input anisotropic amplitude (horizontal coordinate), input dip angle of FVD (vertical coordinate), and back-azimuthal difference between the input 3-D FVDs of tilting-axis anisotropy and output 2-D FVDs of AAN (color dots).The red and blue dot colors denote small and large direction differences between the input tilting-axis anisotropy and output AAN, whose scale is the same as that in panel (b).The dot size denotes the P-wave ray hit-count at each grid node, whose scale is shown at the lower-right corner of panel (c).

Figure 20 .
Figure 20.The same as Figure 19 but for the slow HSA model.In (a), The T-shaped symbols denote fast velocity planes (FVPs) projected on the horizontal plane as shown in Figures 9 and 15.

Figure 21 .
Figure 21.Big mantle wedge anisotropy from (a) our fast HSA model in Alaska and (c) the synthetic subduction zone Model C by Faccenda and Capitanio (2013).The purple areas denote the oceanic plate and the subducting slab.The yellow, red, and brown bars in (a) denote the fast velocity directions (FVDs) at 100, 200, and 300 km depths, respectively.The green and red bars in (c) denote the small and large anisotropic amplitudes, respectively.(b, d) The synthetic SKS splitting results derived from our fast HSA model and Model C by Faccenda and Capitanio (2013).The orientation and length of blue bars in (b) denote the FVD and delay time of SKS splitting, respectively, whose scale is shown at its lower-right corner.
found that the Pacific slab is subducting down to ∼400 km depth and becomes deeper westwards under western Alaska.Martin-Short et al. (2016) obtained a similar result.With the deployment of the USArray portable network in Alaska and the increase in data volume in recent years, the imaging range and resolution have been improved significantly.C. Jiang et al. (2018) conducted a joint inversion of teleseismic S wave travel times and Rayleigh wave dispersions and found that the Pacific slab exhibits a dramatic kink near ∼152°W longitude at ∼120-150 km depths.Gou et al. (2019) obtained high-resolution isotropic Vp and anisotropic tomography showing that the Pacific slab beneath western Alaska reaches 450-500 km depths but does not penetrate into the lower mantle.

Figure 22 .
Figure 22.Vertical cross-sections showing (a), (b) global tomography(Zhao et al., 2013) and (c), (d) regional tomography of the big mantle wedge (BMW) along two profiles in East Asia(Liang, Zhao, Xu, & Hua, 2022)  and Alaska as shown in (e).The blue and red colors denote high and low isotropic Vp perturbations, respectively, whose scale is shown at the bottom-right.Two black dashed lines in (a)-(d) denote the 410 and 660 km discontinuities.The red triangles and yellow dots in (a)-(e) denote intraplate volcanoes and the profile location labels, respectively.(f) A schematic diagram showing that dehydration of the flat slab triggers hot and wet upwelling flows in the BMW feeding the intraplate volcanoes in western Alaska and the Bering Sea.MTZ, the mantle transition zone.
Figure 22.Vertical cross-sections showing (a), (b) global tomography(Zhao et al., 2013) and (c), (d) regional tomography of the big mantle wedge (BMW) along two profiles in East Asia(Liang, Zhao, Xu, & Hua, 2022)  and Alaska as shown in (e).The blue and red colors denote high and low isotropic Vp perturbations, respectively, whose scale is shown at the bottom-right.Two black dashed lines in (a)-(d) denote the 410 and 660 km discontinuities.The red triangles and yellow dots in (a)-(e) denote intraplate volcanoes and the profile location labels, respectively.(f) A schematic diagram showing that dehydration of the flat slab triggers hot and wet upwelling flows in the BMW feeding the intraplate volcanoes in western Alaska and the Bering Sea.MTZ, the mantle transition zone.

Figure 23 .
Figure 23.A 3-D view showing the subducted slab and mantle flow beneath Alaska.The red arrows denote inferred mantle flows including corner flow, toroidal flow, and upwelling.The blue plane with blue contours denotes the subducted slab including the Wrangell slab, Yakutat slab, Pacific slab, and the flat slab in the mantle transition zone (MTZ).The pink dots on the bottom map denote Cenozoic volcanoes including the Wrangell volcanic field (WVF) and the Denali volcanic gap (DVG).Other symbols on the top and bottom maps are the same as those in Figure 1.
(2) In western Alaska, a flat high-Vp anomaly is revealed in the MTZ with its northern boundary at ∼65°N latitude and its western boundary at ∼168°W longitude.A big mantle wedge has formed above the flat slab under western Alaska.Fluids from the slab dehydration and hot upwelling flows in the big mantle wedge have caused the Cenozoic intraplate volcanoes in western Alaska and the Bering Sea.(3)In central Alaska, the Denali volcanic gap is underlain by a high-Vp mantle wedge, whose anisotropy is characterized by subhorizontal and arc-normal fast velocity directions.This result suggests that a remnant of the subducted Yakutat slab becomes a cool and wedge-shaped fluid-barrier, leading to the Denali volcanic gap.(4) In SE Alaska, the subducted Wrangell slab is revealed above 150 km depth.Its mantle wedge exhibits subvertical fast velocity directions, which may reflect hot mantle upwelling induced by the Wrangell slab subduction.The Wrangell volcanic field, Mount Churchill, and Prindle volcano are caused by the hot mantle upwelling around the subducted Wrangell slab.In addition, the Wrangell volcanic field is fed by subslab hot upwelling flow through the Wrangell-Yakutat slab gap.

Figure A1 .
Figure A1.Vp images obtained by inversions with the Tikhonov regularization of (a) constant smoothing coefficient and (b) hitcount-dependent smoothing coefficient.(c) The difference between images shown in (a) and (b).The upper two rows denote map views of Vp images at 120 and 450 km depths.The lower two rows denote vertical cross-sections along the AA′ and BB′ profiles as shown on the map.Other symbols are the same as those in Figure 5.

Figure A2 .
Figure A2.Vp images obtained by inversions with the Tikhonov regularization (a) without and (b) with the slab interface in the smoothing matrix.(c) The difference between images (a) and (b).The upper three rows are map views of Vp tomography at depths of 40, 120 and 300 km.The bottom row shows vertical cross-sections along the AA′ profile as shown on the map.Other symbols are the same as those in Figure 5.

Table 1
Summary of Volcanic Fields in and Around Alaska km