Constraints on the Upper Mantle Structure Beneath the Pacific From 3‐D Anisotropic Waveform Modeling

Seismic radial anisotropy is a crucial tool to help constrain flow in the Earth's mantle. However, Earth structure beneath the oceans imaged by current 3‐D radially anisotropic mantle models shows large discrepancies. Here, we provide constraints on the radially anisotropic upper mantle structure beneath the Pacific by waveform modeling and subsequent inversion. Specifically, we objectively evaluate three 3‐D tomography mantle models which exhibit varying distributions of radial anisotropy through comparisons of independent real data sets with synthetic seismograms computed with the spectral‐element method. The data require an asymmetry at the East Pacific Rise (EPR) with stronger positive radial anisotropy ξ = VSH2VSV2 = 1.13–1.16 at ∼100 km depth to the west of the EPR than to the east (ξ = 1.11–1.13). This suggests that the anisotropy in this region is due to the lattice‐preferred orientation of anisotropic mantle minerals produced by shear‐driven asthenospheric flow beneath the South Pacific Superswell. Our new radial anisotropy constraints in the Pacific show three distinct positive linear anomalies at ∼100 km depth. These anomalies are possibly related to mantle entrainment at the Nazca‐South America subduction zone, flow at the EPR and from the South Pacific Superswell and shape‐preferred orientation (SPO) of melt beneath Hawaii. Radial anisotropy reduces with lithospheric age to ξ < 1.05 in the west at ∼100 km depth, which possibly reflects a deviation from horizontal flow as the mantle is entrained with subducting slabs, a change in temperature or water content that could alter the anisotropic olivine fabric or the SPO of melt.

Advances in computational power and numerical methods (e.g., Akcelik et al., 2003;Komatitsch & Tromp, 2002;Olsen et al., 1995Olsen et al., , 2003 over the last few decades have made large-scale numerical simulations of the seismic wavefield in 3-D complex media much more feasible than before. This has opened up the possibility of 3-D full waveform tomography (e.g., Chen, Tromp, et al., 2007a), providing higher resolution constraints on 3-D Earth structure. One of the most widely used and accurate forward modeling approaches, the spectral element method (SEM; e.g., Komatitsch & Vilotte, 1998) is being used for such purposes. Fichtner et al. (2009) and Tape et al. (2009) were the first to adopt this approach in regional tomography, in the region of Australasia and Southern California, respectively. In addition, these methods have been applied to image structure in other regions such as the upper mantle beneath the Mid-Atlantic Ridge (e.g., Grevemeyer, 2020), East Asia (e.g., Chen et al., 2015;Chen, Jordan, & Zhao, 2007b), North America (e.g., Zhu et al., 2017), and even globally (e.g., Bozdag et al., 2016;French et al., 2013).
At the global scale, Qin et al. (2008), Qin et al. (2009), and Bozdag and Trampert (2010) used variants of the SEM to test the quality of global tomography models. Moreover, Lekić and Romanowicz (2011a) used a cluster analysis of the global upper mantle radially anisotropic model SEMum (Lekić & Romanowicz, 2011b), which was developed using full-waveform tomography and the SEM, to investigative the age dependence of depth profiles of radial anisotropy. Qin et al. (2008) and Qin et al. (2009) included 3-D anisotropy in their modeling and used a coupled SEM-normal mode approach (CSEM; Capdeville, 2005) to reduce the computational expense of the calculations. Moreover, Lentas et al. (2013) used the SEM along with global tomography models to test the robustness of earthquake source parameters. At the regional scale, Ni et al. (1999) compared real waveforms with 2D ray-based synthetics to constrain the low-velocity anomaly in the lower mantle beneath Africa. Subsequently, Chu et al. (2012) and Chu et al. (2014) used 3-D SEM modeling to place constraints on the geometry of the Juan de Fuca slab and on the layering of the lithosphere beneath the North America craton, respectively. Furthermore, Thorne et al. (2013) used waveform comparisons to evaluate 1-D and 3-D seismic models of the Pacific lower mantle and Parisi et al. (2018) used SEM modeling to understand the effects of isotropic versus anisotropic lowermost mantle structure on waveforms. However, anisotropic full-waveform modeling has not been used yet to constrain the anisotropic structure of the oceanic upper mantle. In this study, we provide constraints on the upper mantle structure of the Pacific by waveform modeling of 3-D radial anisotropy structure. We use the SEM (Komatitsch & Tromp, 2002) to simulate seismic wave propagation for three different 3-D radial anisotropy mantle models that exhibit varying upper mantle distributions of radial anisotropy beneath the Pacific. The synthetic waveforms are compared with independently observed surface waveforms and inverted to obtain a preferred model of radial anisotropy in the uppermost mantle beneath the Pacific. Therefore, we pose the questions: How well do current 3-D anisotropic models fit seismic waveform data not used in their construction recorded beneath the Pacific? What adjustments in radial anisotropy are needed for each model to improve the data fit? Is there an age-dependence to the required radial anisotropy beneath the Pacific?
In the following section, we briefly explain the key features and implementation of the Earth models used here. The independent data set used, waveform comparison criteria and inverse modeling approach are explained in Sections 3-5, respectively. Our results are outlined and discussed in Sections 6 and 7, followed by conclusions.

Earth Models Used
In this study, we use the following global 3-D mantle models in the seismic waveform simulations: (i) S362WMANI (Kustowski et al., 2008); (ii) SGLOBE-rani (Chang et al., 2015) and (iii) SGLOBE-smooth, which was constructed for this explained below. As explained in the previous section, some current 3-D tomography models show linear features in radial anisotropy in the Pacific (e.g., S362WMANI and SGLOBE-rani) while others exhibit a smoother, more distributed signature of positive anisotropy. While S36WMANI is already implemented in the SPECFEM3D_GLOBE package (Komatitsch & Tromp, 2002) used in this study and has been extensively tested by the code's developers and users, we add subroutines to the package to implement the mantle structure of SGLOBE-rani. Due to the challenges in matching model parametrization and the spectral element meshing, we prefer not to implement models that we did not construct since not knowing the full details of the models' construction, notably how the crust is treated (e.g., Ferreira et al., 2010), may bias their implementation in the SPECFEM3D_GLOBE package. Thus, to simulate a smoother anisotropic model we built a new model, SGLOBE-smooth. This model was built using the exact same data set and modeling scheme as SGLOBE-rani (Chang et al., 2015), but using a regularization scheme that included both norm damping (as used in the construction of SGLOBE-rani) and gradient damping. S36WMANI, SGLOBE-rani and SGLOBE-smooth, therefore, reflect the various possible features of radial anisotropy previously reported in the literature for the region (Figure 1 and Figure S1).
S362WMANI uses crustal corrections based on CRUST2.0 (Bassin et al., 2000). Ferreira et al. (2010) and Panning et al. (2010) have shown that different crustal models can strongly affect the retrieval of radial anisotropy in the mantle in tomographic inversions. Moreover, Bozdağ and Trampert (2008) and Lekić et al. (2010) have suggested that imperfect modeling of the crust potentially leads to a significant phase misfit between data and SEM seismograms computed using 3D tomographic models. In order to reduce the effects of the crust on the mapping of radial anisotropy in the mantle, when building SGLOBE-rani, Chang et al. (2015) also inverted for crustal thickness variations with respect to CRUST2.0 by including short-period group-velocity dispersion measurements from Ritzwoller and Levshin (1998). Chang and Ferreira (2017) showed that short-period group velocity data with periods of 20 s or shorter are required in global radially anisotropic inversions to constrain thin oceanic crustal structure. In our analysis, each mantle model is combined with the crustal model that was used in its construction. Hence, we combine the mantle structure of S362WMANI with the crustal velocity structure of CRUST2.0. On the other hand, we combine the mantle structure of SGLOBE-rani and SGLOBE-smooth with the crustal velocity structure of CRUST2.0 and the associated crustal thickness perturbations. The difference in the implemented crustal structures can be seen in Figure S2. The Moho is honored in the spectral-element mesh as a first-order discontinuity.
In our waveform modeling, we use the same P-wave speed and density models that were employed in the construction of the tomography models used here. Hence, S362WMANI is implemented using the following scaling relations for V P : and using the density profile defined by STW105 (Kustowski et al., 2008). On the other hand, SGLOBE-rani and SGLOBE-smooth are implemented using the following scaling relations:

Data
We consider 36 earthquakes that occurred from 2005 to 2018 recorded at over 1,125 different stations from 92 networks (Table 1) to study the Pacific's upper mantle. The vast majority of the events (31 out of 36) are chosen after 2009 so that they are independent, that is, they were not used in the construction of the tomographic models assessed here.
Moreover, events are selected with M w > 5 to ensure a good signal-to-noise-ratio but with M w < 7 to prevent substantial finite-source effects. Shallow earthquakes with a hypocentral depth <50 km are chosen to reduce the excitation of surface wave overtones. For each event, we download 90-min-long three-component broadband waveforms from IRIS recorded within an epicentral distance of 10°-120° to avoid near-source effects, caustics and multiple orbit overlapping phases.

Waveform Comparisons and Azimuthal Anisotropy Correction
We use the SPECFEM3D_GLOBE package (Komatitsch & Tromp, 2002) to simulate the global seismic wavefield for the various tomographic models discussed above. We use 256 spectral elements along each side of a chunk in the cubed sphere of the mesh such that synthetic seismograms are accurate down to T ∼ 17 s. A point source model is used with source parameters from the GCMT catalog.
We seek to quantify the level of agreement between the real surface waveforms and synthetic seismograms. We process the data and corresponding synthetic seismograms in exactly the same way (see previous section). We then isolate the fundamental mode surface waves by windowing the waveforms around the maximum amplitude of the data (with a width of 2.5 times the dominant wave period to encapsulate the phase and avoid the interference of surface wave overtones). We then calculate phase misfits, Δϕ for both Rayleigh and Love waves and for each dominant wave period by cross-correlation between the real and synthetic waveforms whereby positive/negative phase misfits correspond to the synthetic waveforms being faster/ slower than the observed waveforms. In addition, amplitude misfits are calculated using the following equation: which shows the logarithmic ratio between the summed square of the real amplitudes, A real and the summed square of the amplitudes in the synthetics, A synthetic for each sample, i, (1 sps) in the data and synthetic windows. These misfits are computed having shifted the synthetics by Δϕ so that they get aligned with the real data.
To exclude obvious outliers from the analysis, waveforms are accepted if the cross-correlation value between the data and synthetics exceeds 0.7, the absolute phase misfit is smaller than 40 s and if the amplitude misfit is between −3 and 3. Thirty two percent of waveforms are rejected based on these criteria and this leads to a total of 2,307 waveform comparisons used in the analysis. The corresponding source-receiver distributions and great-circle paths are shown in Figure 2.
Upper mantle azimuthal anisotropy has been constrained using fundamental-mode Rayleigh wave data (Debayle & Ricard, 2012;Ekström, 2011;Lebedev & Van Der Hilst, 2008). In radial anisotropy studies, it is often assumed that azimuthal coverage is sufficiently uniform for azimuthal anisotropy to average out. However, this assumption might be inadequate at locations where ray paths exhibit a strong azimuthal preference, such as underneath the Pacific (e.g., Ekström, 2011). Similar to Auer et al. (2014), we correct the data for azimuthal anisotropy by applying ray-theoretical corrections (Boschi & Woodhouse, 2006) based on the global dispersion model GDM52 (Ekström, 2011). We only correct fundamental-mode Rayleigh wave dispersion data since Love-based models of azimuthal anisotropy are typically associated with large uncertainties (e.g., Ekström, 2011

Theoretical Framework
After quantifying the misfits between the data and synthetics, we perform inversions of the phase misfits for isotropic and radially anisotropic V S structure. We parametrize isotropic shear wave speed and shear wave anisotropy as follows: By assuming that structural variations in the Earth are laterally smooth, surface waves can be modeled by solving for so-called "local" eigenfunctions and dispersion relations for a set of 1-D profiles along the ray path (e.g., Dahlen & Tromp, 1998). This approach was previously applied to find global and regional models of the mantle (e.g., Auer et al., 2014;Boschi & Ekström, 2002;Boschi et al., 2009). We write the inverse problem as where Δϕ is the phase misfit measured from the cross-correlation of real and synthetic waveforms for a given 3-D reference model (S362WMANI, SGLOBE-rani, or SGLOBE-smooth). We use phase misfit measurements with and without azimuthal anisotropy corrections to assess their effect on the radial anisotropy retrieved from the inversions. Δ is a step along the ray path, c 0 is a laterally varying phase velocity calculated from the 3-D reference model which is itself directly written from SPECFEM3D_GLOBE, a is the Earth's radius, the kernels V S K and  S K relate the observations to Earth's structure parameters, percentual perturbations in isotropic shear wave speed,  S S V V and perturbations in shear wave anisotropy, δζ S , respectively.
While the sensitivity kernels and the model parameters are given in terms of V S and ζ S , to make comparisons with other studies we convert our model results to the traditional anisotropy parameter   2 2 / SH SV V V . If the real and synthetic waveforms agree perfectly then Δϕ is zero for all paths and the inversion results will show that no perturbation is required from the reference model. If, for example, the Rayleigh/Love waves in the synthetic waveforms are systematically faster relative to the real waveforms then Δϕ is positive and the inversion will indicate positive/negative perturbations in ξ are required to fit the data. Note that the reader can find more details of the theory in Appendix A of the supporting information.

3-D Model Parameterization and Inversion Scheme
In this section, we show the sensitivity kernels and we describe the parameterization, damping and inversion scheme used. The sensitivity kernels for T ∼ 40, 60, and 100 s are shown in Figure S3. Note that the reader can find more details of the inversion in Appendix B of the supporting information.
Lateral perturbations are parameterized on nested shell grids with a node spacing of approximately 650 km. On the other hand, triangle basis functions (linear B-splines) with nodes located at 0. 0, 24.4, 74.4, 129.4, 189.9, 256.4, 329.7, and 410.2 km depth are used to parameterize perturbations in the radial direction. We have two model nodes at the surface and at 24.4 km depth ( Figure S4), so that some crustal perturbations are allowed, which may help to prevent contamination of mantle structure by the crust.
We include exactly the same model norm damping and horizontal gradient norm damping to each model as explained in Appendix B of the supporting information. Finally, the misfit equation is optimized using the least squares algorithm of Paige and Saunders (1982).

Resolution Tests
To illustrate the resolving power of the data set and the robustness of the tomographic models, we present the results from checkerboard tests. We created input models with 5% sinusoidal lateral variations in V S and ξ of size 20° and constant amplitudes in the vertical direction. To check the leakage of isotropic V S structure into anisotropic structure, we created input models that only have V S perturbations and zero ξ perturbations. We performed the converse test to check the leakage of anisotropic structure into V S structure, and we also performed tests where perturbations in both model parameters are present. For all resolution tests, we added Gaussian noise to the synthetic data (see Appendix C in the supporting information).   (Kustowski et al., 2008), SGLOBE-rani (Chang et al., 2015) and SGLOBE-smooth (left to right, respectively) at 100-km depth. (b) Waveform comparisons between the data (gray) and synthetics computed for the 3-D mantle models S362WMANI (blue), SGLOBE-rani (green) and SGLOBEsmooth (red). Waveform comparisons are shown for the vertical (Z), radial (R) and transverse (T) components, respectively, at wave periods T ∼ 40, 60 and 100 s (top to bottom). The gray vertical lines show the surface wave windows considered and the phase misfits are shown in the bottom of each subplot for each model considered.

Illustrative Waveform Examples
T ∼ 60 s) indicate that all 3-D global anisotropic models lead to faster Rayleigh waves than the data and that adjustments in radial anisotropy are required to improve the data fit. Moreover, negative Love wave phase misfits around 10 s (Figure 3; T comp., T ∼ 60 s) show that S362WMANI leads to slower Love waves than the data for this path. Figure 4 shows the distributions of phase misfits for all the paths shown in Figure 2 for each of the 3-D global anisotropic models used here. As seen in the illustrative example above, as the period increases, the phase misfit distributions get narrower and and the median decreases, with the lowest misfits occurring for T ∼ 100 s waves, which have maximum sensitivity around ∼150-km depth (Figure 4, bottom). On the other hand, the largest misfits occur for T ∼ 40 s, which have greater sensitivity to shallower depths and to the crust (Figure 4, top). SGLOBE-smooth leads to the poorest Love wave phase misfits, whereby synthetic Love waves are on average 3-4 s faster than the data (Figure 4; T comp., T ∼ 40 s, 60 s; red). As explained in Section 2, SGLOBE-smooth was built using gradient damping, unlike SGLOBE-rani. This leads to stronger positive anisotropy (V SH > V SV ; Figure 3a) than in SGLOBE-rani and hence to too fast Love waves that do not fit the data well. On the other hand, S362WMANI (Figure 4; T comp., T ∼ 40 s, blue) and SGLOBE-rani (Figure 4; T comp., T ∼ 40 s, green) lead to Love waves that fit the data well, with a median phase misfit of <2.5 s. Regarding Rayleigh waves, S362WMANI fits the data slightly better than SGLOBE-rani and SGLOBE-smooth, with the two latter models showing median misfits of ∼3.5-4 s for T ∼ 40 s and T ∼ 60 s, respectively (Figure 4; Z comp).   Figure 2 for waveforms computed using the Earth models S362WMANI (blue), SGLOBE-rani (green) and SGLOBE-smooth (red) filtered with a dominant wave period of T ∼ 40, 60 and 100 s (top to bottom) for each component (vertical, Z; radial, R; and transverse, T; left to right). A phase misfit of 0 is indicated by a vertical black dashed line for reference. The median and standard deviation can be seen at the top left of each subplot, with the medians also being plotted as vertical colored bars. Positive/negative phase shifts Δϕ mean that the synthetic waveforms are faster/slower than the observations. SGLOBE-rani, which shows stronger, positive radially anisotropic anomalies along and around this path than S362WMANI, synthetic Love waves fit the data within 10 s (Figure 5b, middle). Finally, the model with strong and smooth radial anisotropy beneath most of the Pacific (SGLOBE-smooth) leads to synthetic Love waveforms that are often faster than the data by 10 s near the EPR, South America and for paths between Tonga-Kermadec and North America (Figure 5c, middle). In terms of Rayleigh waves, S362WMANI fits the data well within 10 s except for synthetic Rayleigh waves more than 10 s faster than the data to the west of the EPR (Figure 5a, bottom). A similar pattern of phase misfits can be found for SGLOBE-rani and SGLOBE-smooth, which also lead to synthetic Rayleigh waves more than 10 s faster than the data for paths from Tonga-Kermadec to North America (Figures 5b and 5c, bottom).

Surface Wave Phase Misfits
Note that similar trends can be seen in the Rayleigh and Love phase misfits at T ∼ 40 s, which have main sensitivity around ∼60 km depth ( Figure S5). As the wave period increases to T ∼ 100 s (with main sensitivity KENDALL ET AL.  and Rayleigh (bottom row) wave phase misfits for waveforms filtered with a dominant wave period of T ∼ 60 s, color-coded for synthetics fitting the data within twice the average standard deviation, that is, 7 s (gray), more than 7 s slower than the data (pink-brown), more than 7 s faster than the data (blue). The median and standard deviation for each model are shown at the bottom of each subplot.
around ∼150 km depth), the phase misfits become lower in the whole region ( Figure S6). This is again likely due to the less heterogeneous media being sampled. Figures S7-S9 show Rayleigh wave phase misfits before and after correcting for azimuthal anisotropy at T ∼ 40, 60 and 100 s, respectively. The corrections increase with decreasing wave period as the surface waves are more sensitive to the shallower structure. While the corrections are not insignificant, such as up to 3 s for T ∼ 40 s, the pattern of the corrected phase misfits does not change substantially.

Adjustments in Upper Mantle Radial Anisotropy
Using the procedure described in Sections 4 and 5.2 and the phase misfits corrected for azimuthal anisotropy, we compute the adjustments in radial anisotropy required by the data with respect to the original anisotropy in the various 3-D models considered ( Figure 6). Figure S10 shows the variance reduction versus the model norm for the preferred model (green) for SGLOBE-rani.
Fundamental mode Love waves are mostly sensitive to V SH and therefore synthetic Love waves slower than the data would indicate that the radial anisotropy in the model considered   In agreement with the slow Love wave synthetics of S362WMANI along paths from Tonga-Kermadec to North America shown for example, in Figure 5, we find that the data require an increase in anisotropy of ∼4%-6% in this region (Figure 6a, second row). Fundamental mode Rayleigh waves are mostly sensitive to V SV and therefore synthetic Rayleigh waves faster than the data would indicate that the radial anisotropy is too low in the model considered. Therefore, the fast Rayleigh wave synthetics of S362WMANI west of the EPR shown, for example, in Figure 5 mean that the data require an increase in anisotropy of ∼4%-6% in this region.
The fast Rayleigh wave synthetics of SGLOBE-rani in paths west of the EPR and between Tonga-Kermadec and North America (e.g., Figure 5) suggest that the data require an increase in anisotropy of ∼3%-5% and 2%-4% in these regions, respectively (Figure 6b, second row). Similarly to SGLOBE-rani, the fast Rayleigh wave synthetics of SGLOBE-smooth west of the EPR and between Tonga-Kermadec and North America ( Figure 5) imply that the data require an increase in anisotropy of ∼1%-2% in these regions (Figure 6c, second row). The data also suggest that radial anisotropy in SGLOBE-smooth is too high and that a reduction of ∼4%-6% is required in the south Pacific.
To summarize our findings, the third row in Figure 6 shows the retrieved (i.e., including the adjustments) distribution of perturbations of radial anisotropy for each model at ∼100 km depth. We find that current tomographic models such as S362WMANI and SGLOBE-rani underpredict the magnitude of radial anisotropy in the young Central Pacific and to the east of the EPR. The inversion of the new surface wave phase misfit data presented here reveals that three linear anomalies of positive radial anisotropy along the west coast of South America, west of the EPR and around and south of Hawaii become more pronounced and defined than in the original models. Furthermore, the bottom row in Figure 6 shows the retrieved distribution of absolute radial anisotropy for each model at ∼100 km depth. We find good agreement among these adjusted anisotropic models, which show that the data require ξ to be 1.12-1.15 (dξ/ξ = 5%-7% with respect to PREM) beneath the west coast of South America, 1.13-1.16 (dξ/ξ = 6%-9% with respect to PREM) west of the EPR and 1.12-1.14 beneath the central-Pacific. Radial anisotropy then reduces further to ξ = 1.08-1.1 (dξ/ξ = 2%-4% with respect to PREM) in the northwestern Pacific and to a minimum of ξ < 1.05 in the West Pacific at ∼100 km depth.
We also find good agreement in the retrieved distribution of absolute radial anisotropy at ∼60 km depth ( Figure S11, bottom row). However at ∼150 km depth, the sensitivity kernels for T ∼ 100 s ( Figure S3) become wide and there is a poorer agreement in the magnitude of the retrieved distribution of absolute radial anisotropy ( Figure S12, bottom row).

Azimuthal Anisotropy Contribution
Figures S13-S15 show the contribution that azimuthal anisotropy can have on the retrieved radial anisotropy at ∼60, 100, and 150 km depth, respectively. Including azimuthal anisotropy, corrections reduce anisotropy west of the EPR and increase anisotropy beneath the Central Pacific at 60 and 100 km depth ( Figures S13 and S14), which is similar to the findings of Auer et al. (2014). However, the effect of azimuthal anisotropy corrections on the pattern of retrieved anisotropy is small (max ξ = ±0.01), which is also similar to the findings of Auer et al. (2014) when using a combined data set in their whole mantle inversions. Figure 7 shows the results of the checkerboard tests as described in Section 5.3. We find good recovery of 5% sinusoidal variations in V S (Figure 7, top row) and ξ (Figure 7, middle row) with minimal leakage of V S into ξ and vice versa. When both anomalies are present and inverted for, the larger amount of criss-crossing ray paths near Asia and South America result in higher resolution than in the south and mid-Pacific where there is smearing. Input variations are also reasonably well recovered at shallower depths such as ∼60 km ( Figure S16). As discussed in Section 6.3, the sensitivity kernels for T ∼ 100 s ( Figure S3) become wide and so the resolving power at 150 km depth is weaker ( Figure S17).

Discussion
In this study, we used independent waveform modeling to analyze the robustness of radially anisotropic features in existing tomography models in the upper mantle beneath the Pacific. We presented comparisons of surface waves in real independent waveforms with highly accurate synthetics computed for 3-D radially anisotropic mantle models. Furthermore, we inverted the phase data corrected for azimuthal anisotropy to obtain absolute anisotropy in the uppermost mantle beneath the Pacific. Doing so allowed us to better constrain upper mantle anisotropy structure, which is key for more accurate interpretations in terms of mantle flow.
Our analysis showed that surface wave phase misfits are in the range of about ±15 s for T ∼ 60 s waves (Figures 4 and 5), which require adjustments of up to 6% in the radial anisotropy of the 3-D mantle models considered ( Figure 6, second row). Here we found that the data require strong, positive (V SH > V SV ) lateral variations in radial anisotropy up to ξ ∼ 1.16 at ∼100 km depth (Figure 6, bottom). These findings confirm previous reports since the 1980s of faster SH anomalies in the upper few hundred kilometers beneath the Pacific (e.g., Beghein et al., 2014;Burgos et al., 2014;Cara & Lévêque, 1988;Ekström & Dziewoński, 1998;Gung et al., 2003;Isse et al., 2019;Montagner & Tanimoto, 1991;Nettles & Dziewoński, 2008;Nishimura & Forsyth, 1989). However, we found that current tomographic models such as S362WMANI and SGLOBE-rani underpredict the magnitude of radial anisotropy beneath the east and Central Pacific and to the east of the EPR.
The largest adjustments are required by S362WMANI which exhibits lower anisotropy overall than SGLOBE-rani or SGLOBE-smooth (see e.g., Figure 1). This may be because the latter, more recent models are built using substantially more waveforms and shorter period data than S362WMANI. In particular, at 100 km depth, the region close to the EPR requires the largest adjustments of up to ∼6% in radial anisotropy across the three anisotropic models (Figure 6, second row). This correlates well with synthetic resolution tests by Chang et al. (2015) (Figure 16 in their study), which show a slightly poorer recovery of the amplitude of anisotropic anomalies in the young, eastern Pacific compared with the rest of the Pacific. Therefore the underestimated radially anisotropic anomalies near the EPR are likely the result of poor data coverage. Although the anisotropy resolution of our inversions is limited in the southwest of the EPR (Figure 7) and hence the details of the retrieved structures are not strongly constrained, the need for stronger anisotropy in the region is clearly required by the data. Moreover, we note that the loss in the strength of the anisotropy obtained from our resolution tests is about the same for all regions. Hence, it is appropriate to interpret differences in the amplitudes of the retrieved anisotropy for the various regions considered. To improve the retrieval of radially anisotropic anomalies in the Pacific, particularly west of the EPR and beneath the central Pacific, it is essential to incorporate data from ocean-bottom seismometers (OBSs) such as the Pacific Array (Kawakatsu, 2016) into future data sets. The inversion of the new surface wave phase misfit data presented here reveals more detailed lateral variations in radial anisotropy beneath the Pacific than in previous studies. Specifically, three pronounced linear anomalies of positive radial anisotropy along the west coast of South America, west of the EPR and around and south of Hawaii become more pronounced and defined than in the reference models. These three positive, linear anomalies in radial anisotropy are discussed individually in the following three sections.

Subduction-Related Anisotropy
The data require strong ξ of 1.12-1.15 at ∼100 km depth beneath the west coast of South America. The Andean margin can be subdivided into five main tectonic segments, comprising regions with normal subduction and flat subduction. The majority of fast directions in azimuthal anisotropy studies (e.g., Eakin & Long, 2013;Eakin et al., 2014Eakin et al., , 2015 align trench-parallel, as often found in the forearc of subduction zones (Eakin et al., 2014). This corresponds well with the strong, positive radial anisotropy (ξ = 1.12-1.15) required in this study and could be linked with horizontal flow in the sub-slab mantle. Moreover, a recent surface wave tomography study incorporating a vast amount of OBS data with a focus beneath the Pacific Ocean by Isse et al. (2019) also images strong, positive radial anisotropy with ξ ∼ 1.12 at ∼75-100 km depth beneath the west coast of South America.

Anisotropy Beneath the Young Pacific
To the west of the EPR, a linear region with ξ = 1.13-1.16 at 100 km is required by the data. A positive radially anisotropic anomaly near the EPR is consistent with the study of Gu et al. (2005) and Nettles and Dziewoński (2008) in which an anomaly of dξ/ξ ∼ 6% at ∼100 km is reported. Moreover, the surface wave tomography study by Isse et al. (2019) which incorporated OBS data also images positive radial anisotropy with ξ > 1.1 at ∼75-100 km depth west of the EPR. We propose at least three regions that may contribute to forming this linear, positive anomaly in radial anisotropy along west of the EPR. First, at the north of the EPR we propose that mantle entrained by the Rivera-Cocos subduction zone may result in strong, positive ξ. Second, the EPR is the fastest spreading ridge in the world and strongly deforming horizontal mantle flow can lead to the LPO of anisotropic minerals and therefore strong, positive radial anisotropy. The radial anisotropy could also be attributed to the SPO, for example, from partial melting (e.g., Isse et al., 2019;Tan & Helmberger, 2007;Schmerr, 2012). The presence of melt beneath the ridge axis, particularly to the west, is also supported by Baba et al. (2006).
Despite the Pacific and Nazca plate sharing a boundary in the EPR, several studies have found pronounced asymmetry. For example, the Pacific (west) side is characterized by a higher abundance of seamounts (with a source propagating eastwards at ∼20 cm/yr; Ballmer et al., 2013), lower S-wave velocities (e.g., Dunn & Forsyth, 2003;Forsyth et al., 1998;Toomey et al., 1998), greater shear-wave splitting (e.g., Wolfe & Solomon, 1998), a higher electrical conductivity (e.g., Evans et al., 1999) and slower subsidence (e.g., Cochran, 1986;Morgan & Smith, 1992;Morgan et al., 1995) than the Nazca (east) side. We report an asymmetry in radial anisotropy at the EPR for the first time with ξ = 1.13-1.16 on the west as opposed to ξ = 1.11-1.13 to the east. These results suggest a stronger LPO of anisotropic mantle minerals (e.g., olivine and enstatite) west of the EPR than to the east. This stronger LPO suggests that mantle flow beneath the EPR is not the symmetric corner flow usually assumed in mid-ocean ridge models (e.g., Morgan, 1987). Instead, the stronger LPO could be due to pressure-driven eastern asthenospheric flow from the South Pacific Superswell. This interpretation of the asymmetry across the south EPR (e.g., The MELT Seismic Team, 1998) being the result of vigorous, shear-driven eastward asthenospheric flow opposing plate motion is in line with previous studies (e.g., Ballmer et al., 2013;Conder et al., 2002;Lin et al., 2016;Toomey et al., 2002;Weeraratne et al., 2007). Furthermore, this hypothesis has been supported and explained geophysically and geochemically by some numerical models (e.g., Ballmer et al., 2010Ballmer et al., , 2013Harmon et al., 2011). LPO due to this eastern asthenospheric flow from the South Pacific Superswell could form the southern contribution to the linear anomaly west of the EPR identified in this study.

Anisotropy Beneath the Central Pacific
Radial anisotropy then reduces with lithospheric age before increasing to ξ = 1.12-1.14 beneath the central Pacific, which is consistent with Ekström and Dziewoński (1998) and Nettles and Dziewoński (2008), and the surface wave tomography study incorporating OBS data by Isse et al. (2019). The paths from Tonga-Kermadec to North America tend to be long, so strictly the path-averaged anomalies could be hard to interpret. However, we propose that the positive linear anomaly in radial anisotropy in the central Pacific is due to a combination of (a) anisotropy preserved along the Hawaiian-Emperor chain (e.g., Collins et al., 2012), (b) LPO due to the plume-lithosphere interaction beneath Hawaii (e.g., Auer et al., 2015), a change in olivine fabric as the plume dehydrates (e.g., Karato, 2008) and partially melts or due to SPO from partial melt (e.g., Isse et al., 2019;Rychert et al., 2013;Schmerr, 2012), and (c) poor resolution in the South Pacific.

Age-Dependence of Radial Anisotropy
Radial anisotropy then reduces further to ξ = 1.08-1.1 at ∼100 km depth in the northwestern Pacific, in agreement with Isse et al. (2019), which may reflect the lower levels of deformation beneath the old oceanic lithosphere. In the west, radial anisotropy reduces to ξ < 1.05 at ∼100 km depth. This lower magnitude of radial anisotropy could indicate a deviation from horizontal to vertical flow as flow is entrained with subducting slabs. Alternatively, an age-dependence of radial anisotropy could be linked to a change in rheology with age. Anisotropy near the EPR at 100 km is likely from an asthenospheric source as opposed to a lithospheric source beneath the older ocean. The change in temperature or water content which in part controls the rheology, could in turn alter the anisotropic olivine fabric (e.g., Karato, 2008) or small fractions of melt in the lithosphere may also alter the detected seismic anisotropy (e.g., Tommasi & Vauchez, 2015).

Limitations
By using the advanced SEM for 3-D Earth models, we conduct accurate measurements of phase misfits, which are then interpreted with simpler, computationally efficient inverse modeling. Lekić and Romanowicz (2011a) and French and Romanowicz (2014) also used a similar hybrid approach, which represents a good compromise between accuracy and efficiency. However, our study has some limitations. Moreover, due to the lack of seismic stations in the Pacific, this data set comprises relatively long paths which may average out smaller-scale variations in structure. Future data sets would benefit from including shorter paths to provide more localized constraints on seismic anisotropy such as Lin et al. (2016) and Russell et al. (2018). Ongoing and future seismic deployments in the region should enable further refined analyses in the future. Future work would also include evaluating the retrieved models developed in this study with further SEM simulations. Given that the largest phase misfits are associated with shortest period Love data, in the future if more robust Love-based models of azimuthal anisotropy are available, one should correct the corresponding dispersion curves to constrain the shallower layers more accurately. Although we choose events with Mw < 7 and visually inspect all waveforms, errors can still arise in the measured surface wave misfit from the source phase or cycle skipping, respectively. Finally, source mislocation errors could also contribute to the measured surface wave misfit. For example, considering an event 10 km deeper than that shown in Figure 3 leads to a phase shift of up to ∼8 s to T ∼ 40 s ( Figure S18, top row) and up to 1 s for longer periods such as T ∼ 60 s and T ∼ 100 s ( Figure S18, middle and bottom row, respectively). Future systematic quantifications of phase errors due to source mislocations could be helpful.

Conclusions
In this study, we constrained the radially anisotropic structure in the Pacific upper mantle by using waveform modeling to assess the robustness of anisotropic features beneath the Pacific in 3-D tomographic mantle models exhibiting varying distributions of radial anisotropy. We compared synthetic SEM seismograms with independent surface waveforms recorded in and around the Pacific. The phase misfits for Love and Rayleigh waves were corrected for azimuthal anisotropy and then inverted for the 3-D anisotropy models considered.
We found that the data require an asymmetry at the EPR with stronger positive radial anisotropy ξ = 1.13-1.16 at about 100 km depth to the west of the EPR than to the east (ξ = 1.11-1.13). This asymmetry in radial anisotropy is possibly due to the LPO of intrinsically anisotropic mantle minerals produced by shear-driven asthenospheric flow beneath the South Pacific Superswell. We find that current tomographic models such as S362WMANI and SGLOBE-rani underpredict the magnitude of radial anisotropy in the young-central Pacific and to the east of the EPR. The inversion of the new surface wave phase misfit data presented here reveals that three linear anomalies of positive radial anisotropy along the west coast of South America, west of the EPR and around, and south of Hawaii become more pronounced and defined than in the reference models and are likely related to regional processes. Radial anisotropy then reduces with lithospheric age from ξ = 1.12-1.14 beneath the central Pacific to ξ < 1.05 beneath the oldest seafloor at ∼100 km depth. Lower radial anisotropy could indicate a deviation of the flow direction from horizontal as flow is entrained with subducting slabs, a change in temperature or water content that could alter the olivine fabric or the SPO of melt.

Data Availability Statement
The 3-D radially anisotropic shear wave speed model presented in this study, which we name SPacific-rani, is available online: http://ds.iris.edu/ds/products/emc-spacific-rani/. The data for this research are available through the IRIS Data Management Center: https://www.iris.edu/hq/. Specifically, the authors gratefully acknowledge the availability of broadband seismograms including from the IRIS/USGS (IU), the IRIS/ IDA Network (II); the GEOSCOPE Network (G); the GEOFON Network (GE), and the Pacific21 Network (PS). Source parameters used in the SPECFEM3D_GLOBE simulations were taken from the GCMT catalog (https://www.globalcmt.org/CMTsearch.html).