Seismic evidence for subduction-induced mantle flows underneath Middle America

Laboratory experiments and geodynamic simulations demonstrate that poloidal- and toroidal-mode mantle flows develop around subduction zones. Here, we use a new 3-D azimuthal anisotropy model constructed by full waveform inversion, to infer deep subduction-induced mantle flows underneath Middle America. At depths shallower than 150 km, poloidal-mode flow is perpendicular to the trajectory of the Middle American Trench. From 300 to 450 km depth, return flows surround the edges of the Rivera and Atlantic slabs, while escape flows are inferred through slab windows beneath Panama and central Mexico. Furthermore, at 700 km depth, the study region is dominated by the Farallon anomaly, with fast axes perpendicular to its strike, suggesting the development of lattice-preferred orientations by substantial stress. These observations provide depth-dependent seismic anisotropy for future mantle flow simulations, and call for further investigations about the deformation mechanisms and elasticity of minerals in the transition zone and uppermost lower mantle.

M apping mantle flows induced by descending oceanic lithospheres has profound implications for understanding the dynamic, thermal and chemical evolution of subduction systems 1,2 . Two end-member flow circulation patterns have been proposed based on geodynamical, geochemical and seismological observations. The first end-member suggests that the downdip movement of slabs entrains mantle flows to great depths, causing 2-D corner flows perpendicular to the trench in the mantle wedge and beneath the slabs, also known as the poloidal-mode flow 3 . The second end-member suggests that, to preserve the mass displaced by trench migration and slab rollback, 3-D return flows develop in the sub-slab region and surrounding the lateral edge of sinking slabs, also known as the toroidal-mode flow. The toroidal-mode flow typically happens for convergent plate margins with fast trench migration speeds 4,5 , and is important for controlling trench curvature 6 . To date, the behaviors of these two end-members have been extensively investigated using laboratory experiments 7,8 and mantle flow simulations 9,10 .
Seismic anisotropy is a valuable tool to detect mantle flows in subduction zones 11,12 . Olivine and orthopyroxene, two dominant upper mantle minerals, are intrinsically anisotropic in terms of seismic wavespeeds 13,14 . The strain-induced lattice-preferred orientation (LPO) of olivine leads to the direction and polarization dependence of seismic wavespeeds 15 , which can be detected using shear wave splitting and surface wave tomography. Previous shear wave splitting studies document a trench-parallel pattern beneath the fore-arc and the subducted slab 4,16 , as well as a trench-perpendicular pattern in the back-arc mantle 3 . However, shear wave splitting measurements represent the path-integrated effects of anisotropy over the entire raypath, leading to poor vertical resolution and controversial interpretations 17 . Surface wave tomography, on the other hand, enables us to constrain depth-dependent azimuthal anisotropy by taking advantage of the dispersive characteristics of Rayleigh and Love waves [18][19][20][21] . However, its lateral resolution is typically limited to several hundreds of kilometers due to the use of long-period signals, making it difficult to resolve small-scale features, such as the toroidal-mode mantle flows in the vicinity of sinking slabs. The different lateral and depth resolutions of shear wave splitting and surface wave tomography also make it challenging to reconcile these two independent measurements 20,22 .
In this study, we characterize detailed, depth-dependent seismic anisotropy for the Middle American and Caribbean subduction systems using full waveform inversion [23][24][25] . The Middle American Subduction Zone (MASZ) is a 2700 km long, active convergent margin, stretching from Mexico to Costa Rica (Fig. 1), which involves the subduction of the young Rivera and Cocos Plates underneath the North American and Caribbean Plates. Both the Rivera and Cocos Plates are the remnants of the ancient Farallon Plate, which broke into several fragments in the Oligocene 26 . Currently, the 10 Myr oceanic crust of the Rivera Plate is subducting underneath the Jalisco Block of Mexico at a speed of 30 mm/yr. Along the strike of the Middle American Trench (MAT), the convergence rate of the Cocos Plate relative to the overriding Caribbean Plate increases from 50 mm/yr in southern Mexico to 90 mm/yr offshore Nicaragua and Costa Rica 27 .
As forming in the Cretaceous, the Caribbean Plate has undergone complicated interactions with the North American, South American and Farallon Plates 28 . It is now bounded to the north by a 1200 km east-west orientated, left-lateral Septentrional-Oriente Fault System, and to the south by another >1000 km long, right-lateral fault, the San Sebastian-El Pilar Fault System. To the west, it is confined by the eastward descending Cocos Plate, and to the east by the westward sinking Atlantic Oceanic lithosphere to form the Lesser Antilles Volcanic Arc (LAVA). The transition from the subduction of the Atlantic slab to transform motion along the San Sebastian-El Pilar Fault System has been suggested as an example of a Slab-Transform Edge Propagator (STEP) fault 29 . Two tectonic reconstructions have been proposed to explain the origin of the Caribbean Plate: the Pacific 30 and intra-Americas origin 31 , placing this plate to the west and east of the Farallon Plate in the Cretaceous, respectively. The arcuate-shaped Antilles Subduction Zone is still an enigma in Middle America, which includes the Greater Antilles in the north underneath Hispaniola and Puerto Rico, and the Lesser Antilles in the southeast (Fig. 1). The shallow portion of the subducting slab can be directly inferred from the Wadati-Benioff zone 32 . In the lower mantle, the slabs underneath the Caribbean Plate merge with the Cocos Plate into the prominent fast Farallon anomaly, which has been detected in global tomography models 33,34 .
The Middle American and Caribbean Subduction Zones are ideal locations to investigate the dynamics of subduction systems because of strong variations in the slab morphology 35 , fast trench migration speeds (20-30 mm/yr) 5 , the concave-shaped LAVA, as well as heterogeneous geochemical observations in the Trans-Mexico Volcanic Belt (TMVB) 36 and Central American Volcanic Arc (CAVA) 37 . All these features indicate complex mantle flows and subduction processes underneath Middle America, which are explored by using a new azimuthal anisotropy model US 32 in this study. 32 . With the framework of full waveform inversion, we combine three-component shortperiod (15- Table 1 gives the numbers of stations in major contributing arrays (with stations > 30) located in Middle America and used for constructing model US 32 . The current model is the result of 32 preconditioned conjugate gradient iterations. The first 22 iterations are utilized to constrain radially anisotropic model parameters, including horizontally propagating, vertically (β v ) and horizontally (β h ) polarized shear wavespeeds 38 . The last ten iterations are used to further constrain radially (L and N) and azimuthally anisotropic (G c and G s ) model parameters. These four model parameters are simultaneously updated using the Fréchet derivatives calculated by the adjoint state method, details for constructing the model can be found in Zhu et al. 39 and methods. Figure 2 demonstrates the capability of full waveform inversion to match three-component observed (black) and predicted (red) seismograms for a shallow earthquake that occurred in 2012 underneath Costa Rica, both observations and predictions are bandpass filtered from 25 to 100 s. High-quality data from the high-density USArray enables us to capture complex, triplicated body waves passing through the upper mantle discontinues, and offers us opportunities to harness these signals for delineating detailed structures near these phase transition boundaries. In Fig. 2, three-component body waves, as well as Rayleigh and Love waves are matched very well, suggesting that the mantle structures underneath Middle America are well constrained in the current model US 32 . More waveform fittings can be found in Supplementary Note 1, and comparisons with previous tomographic studies can be found in Supplementary Note 6.

3-D azimuthal anisotropy model US
Resolution analysis. Analyzing resolution or uncertainty of full waveform inversion is important and challenging, considering computational costs for performing iterative inversion [40][41][42][43] . Here, we use the approximated diagonal Hessian and point-spread functions (PSFs) to analyze illumination, resolution, and more importantly tradeoffs among different model parameters for several key features discussed in the following sections. Figure 3 shows the horizontal distribution of the approximated diagonal Hessian underneath Middle America at depths ranging from 100 to 600 km, which is a good indicator for the ray density coverage 44 . Owing to the configuration of stations and earthquakes ( Supplementary Fig. 2), the illumination for North America, the MAT and the Lesser Antilles is good at all depth ranges. In contrast, the illumination for the central Caribbean is not as good as other regions, especially at greater depths. Figure 4 presents one PSF test around the Rivera slab offshore Mexico, a location with interesting features for discussing toroidal-mode mantle flows. In this calculation, a Gaussian anomaly with a half width of 120 km is added to G s at a depth of 350 km, and three other model parameters (G c , L and N) are held unchanged. Results for this PSF test suggest that there are good constraints for the fast axis orientations at this location, which depend on the ratio between G c and G s . For instance, there is limited leakage from G s to other three model parameters, and the resulting PSF for G s perturbation is well focused. More PSF tests can be found in Supplementary Note 2. These PSFs give us confidence about the distributions of wavespeeds and anisotropic fabrics at these locations.     Subduction-induced mantle flows. The narrow Rivera slab is imaged dipping~40°northeast underneath the Jalisco Block in Mexico (profile A-a in Fig. 6). It appears that the Rivera slab is detached from the surface, despite the occurrence of Mw > 5 earthquakes at depths shallower than 150 km, which is consistent with global tomography results 48 . Its mantle wedge is characterized by slow wavespeed anomalies, suggesting the ascent of hot asthenosphere from 400 km towards the TMVB, which might explain uplift and low Bouguer gravity anomalies (~−100 mgal) in this region (Supplementary Fig. 1). The subducted Rivera and Cocos slabs are separated by a slab window, filled with pronounced slow anomalies (<−2%) from 300 km down to the transition zone around 12-18°N (Fig. 5), in agreement with previous body and surface wave tomographic results 49,50 . Furthermore, the fast axes are revealed following the  . This pattern reflects the development of return flows surrounding the narrow (~200 km), fast rollback (~31 mm/yr) Rivera slab, as well as escape flows moving across slab windows, owing to large pressure gradients created by slab fragmentations as demonstrated in previous physical and numerical experiments 6,7 . The return flows and slab detachment might contribute to the observed mixed subduction-related calcalkaline rocks with intraplate-like lavas and adakites along the TMVB 36 , as well as the unusual oblique angle (16°) between the MAT and TMVB. In profile B-b (Fig. 6) and 3-D visualization (Fig. 7), the Cocos slab plunges~60°down to the transition zone underneath Nicaragua and Costa Rica. Its northern boundary is truncated sharply along the extension of the Septentrional-Oriente Fault System (Fig. 5). Flat subduction of the Cocos Plate might exist underneath central Mexico ( Supplementary Fig. 14), however, it is still challenging to detect thin, horizontal contrasts at shallow depths using seismic tomography because of its limited vertical resolution. The observation of different dip angles of the Rivera and Cocos slabs along the strike of the MAT is consistent with previous studies from seismicity distribution 35 , local tomography 49 and receiver functions 51 . The steeper dip of the Cocos slab can be attributed to its older age (23 Myr) and therefore greater density. Similar to the Rivera slab, the Cocos slab is also detached from the surface, with the mantle wedge filled with strong slow seismic anomalies connecting to the CAVA and Cayman Trough. Moreover, the trench-parallel anisotropic pattern is observed in the sub-slab region of the Cocos slab at depths ranging from 300 to 450 km (Fig. 5). There is a sharp transition in the mantle wedge from the trench-parallel pattern at 50 km to the trenchperpendicular pattern at greater depths. These complex, depthdependent anisotropic fabrics suggest the development of toroidal-mode flows in the sub-slab region, in response to fast slab rollback (25 mm/yr) of the MAT 5 . Toroidal-mode flows are forced to tilt toward shallower depths in the mantle wedge due to fast slab rollback, and interfere with the corner flows entrained by the downdip motion of the Cocos slab, in agreement with 3-D geodynamic flow simulations 9 . Furthermore, the fast axes in model US 32 run across the Panama slab window 28 between the Cocos and Nazca Plates at 300-400 km depth. This observation supports the hypothesis about the subduction of the Cocos-Nazca Ridge and Galápagos hotspot track in this region 28 , and offering an explanation about the occurrence of Neogene alkaline and arctholeiitic volcanism in Panama 52 .

Shallow structures.
Underneath the LAVA, the Atlantic lithosphere sinks toward the west with~50°dip down to the transition zone (Figs. 6C-c and 7). In contrast to the Rivera and Cocos slabs, the Wadati-Benioff zone beneath the LAVA reflects subduction from the surface down to 200 km, with a slow anomaly imaged in the mantle wedge underneath the LAVA. In Fig. 5, there are two fast bodies at 300-400 km depth, one underneath Puerto Rico and the other underneath the southern LAVA, which are separated by a slab window around 16-18°N. Previous P wave tomography studies report a similar feature, which was attributed to the subducted North American-South American plate boundary 28,53 . Similar to the anisotropic fabrics observed around the Rivera slab, the fast directions surround the lateral edge of the descending Atlantic lithosphere, and move across the slab window (Figs. 5 and 7). This pattern suggests the development of return flows warping around the narrow, rapidly advancing (17 mm/yr) Atlantic slab, which pushes materials within the mantle wedge and induces flows around its edge. The current model fails to delineate east-west oriented anisotropic fabrics along the San Sebastian-El Pilar Fault System as demonstrated in shear wave splitting 54 , suggesting limited resolution due to the lack of stations at this location ( Supplementary Fig. 2).
In profile D-d (Fig. 6), the Caribbean slab dips southeast underneath the Southern Caribbean Deformation Zone and Maracaibo Block. Its dip increases from 30°at depths shallower  Fig. 9 The histograms of angle differences between the fast axis orientations in model US 32 and the plate motion model NUVEL-1A 58 . a-f Results at depths ranging from 70 to 400 km. Blue and red histograms are results for the hotspot and no-net rotation reference frames, respectively. Numbers in the brackets to the right corner of each histogram denote mean angle differences between these measurements. than 200 km to around 90°at greater depths, which might be affected by the presence of the thick (~250 km), craton-like overriding South American Plate. To the northwest, the Cayman Trough is underlain by a triangular slow anomaly (<−1.5%) down to 300 km, reflecting strong deformation owing to the relative motion and pull-apart between the North American and Caribbean Plates along the Septentrional Fault System since the Eocene 28 . By pushing mantle materials in the sub-slab region, the fast northwestward rollback of the Caribbean slab (22 mm/yr) 55 might also contribute to the observed northeast-southwest oriented flow field underneath the Caribbean Plate at depths ranging from 300 to 500 km (Fig. 5).
At 600 km, the fast wavespeed anomalies discussed above merge into a continuous, arcuate-shaped feature surrounding the Gulf of Mexico and Caribbean, which then turns to a northsouth oriented fast Farallon anomaly at 700 km 34 , stretching from the southern U.S. to South America (Fig. 5). Strong slow anomalies are observed to the west and east of the Farallon slab, which might reflect water entrained down to the transition zone by subduction of the Farallon Plate since the Mesozoic. In comparison to the upper mantle and transition zone, the fast orientations turn perpendicular to the strike of the Farallon anomaly, reflecting poloidal-mode mantle flows in the uppermost lower mantle. This flow field might be driven by strains produced when the Farallon slab penetrated the 660-km discontinuity to enter the lower mantle. This process results in large stresses due to a significant viscosity jump across the 660-km discontinuity 56 , and might produce detectable straininduced seismic anisotropy.

Discussion
Comparisons with plate motion data. Debayle et al. 57 concluded that surface wave azimuthal anisotropy models correlate well with global plate motion models 58,59 for regions with a relatively simple tectonic history, such as oceanic plates. However, for regions with a complicated tectonic history, such as continents, their correlation is much poorer. Figure 8 presents the angle differences between the fast axis orientations in model US 32 for the study region with the global plate motion model NUVEL-1A 58 at depths ranging from 70 to 400 km. Both hotspot (Fig. 8) and no-net rotation ( Supplementary Fig. 17) reference frames are used for these comparisons. Figure 9 shows the histograms for these angle differences at different depths. Overall, the current seismic model fits plate motion results better for the hotspot reference than the no-net rotation reference at all depths. For instance, from 150 to 300 km, the mean angle differences for the hotspot reference are <30°, whereas these numbers are >35°for the no-net rotation reference. It appears that the fast axis orientations in the current model fit plate motions with the hotspot reference frame better for the North Atlantic and Caribbean. In contrast for the Cocos Plate, the no-net rotation frame fits seismic measurements better. Major differences come from regions with complicated fast axis orientations, especially in the vicinity of subducting slabs, which are not included in large-scale plate motion results.
Implications for mantle flow simulations and deformation mechanisms in the deep Earth. LPO of anisotropic mineral aggregates deformed by dislocation creep is the major cause of seismic anisotropy. Experiments demonstrate that the major minerals within the mantle transition zone and the uppermost lower mantle, such as wadsleyite and bridgmanite, might also be intrinsically anisotropic 60,61 . Diffusion creep or superplasticity is considered as the dominant deformation mechanism in the deep mantle 62 , and this prevents the development of LPO. Therefore, the existence and nature of seismic anisotropy within the transition zone and lower mantle are still highly debated. A recent study reported strong radial anisotropy (up to 2%) in the vicinity of subducting slabs down to 1000-1200 km 63 , which raises questions about the deformation mechanisms in the deep Earth. In addition, a number of previous surface wave tomography 64,65 and shear wave splitting 66,67 studies speculated about whether or not azimuthal anisotropy exists in the deep mantle, especially around subducted slabs. Anisotropic signatures (with magnitude around 1%) observed at greater depths in this study provide new seismic evidence for the development of complicated mantle flows induced by subduction systems with fast trench migration speeds and strong variations in slab morphology, such as tearing and detachment. All these seismic observations suggest that dislocation creep might have a more important role in the deep Earth than we thought before. In addition, all interpretations in this article assume that we have A-, C-and E-type LPO for olivine aggregates, other types of LPO, such as type-B fabric affected by water distributions, may change current interpretations 14,68 , especially for the mantle wedge regions.
Furthermore, depth-dependent azimuthal anisotropy observed in this study offers more constraints for investigating the physics of subduction, such as using mantle flow simulations 9,10 and physical experiments 7,8 . Most recent mantle flow simulations allow us to construct geodynamic models to fit different geophysical and geological observations, such as shear wave splitting 54,69 and plate motions 58 . However, shear wave splitting represents the integrated anisotropic effects over the entire depth range 17 . Therefore, it is a non-unique solution to match splitting measurements through LPO calculations 70 and flow simulations. The depth-dependent azimuthal anisotropy patterns delineated in this study provide future opportunities to better constrain 3-D deep mantle flows and investigate physical processes through geodynamic modeling.
Misfit function. The details of constructing model US 32 can be found in Zhu et al. 39 Here, we summarize the key components of this study. The spectral-element method 71,72 is used to solve the anisotropic/anelastic wave equation, and compute threecomponent synthetic seismograms. Each forward simulation takes 45 minutes using 144 CPU cores, therefore, in total, it took 1.87 million CPU hours to construct the current model. Frequency-dependent phase measurements for three-component short-period (15-40 s) body waves and long-period (25-100 s) surface waves are combined to simultaneously constrain deep and shallow structures 73 . No crustal correction 74,75 is required in the inversion owing to the simultaneous update for the crustal and mantle structures. A multi-scale inversion strategy 76 , i.e., starting with long-period signals for long-wavelength updates and gradually incorporating short-period signals to improve resolution, is utilized to mitigate the cycle skipping problem 77 . The total misfit function includes six categories: P-SV body waves on vertical and radial components, and SH body waves on transverse component; Rayleigh surface waves on vertical and radial components, and Love surface waves on transverse component. A multi-taper technique 78 is employed to calculate frequency-dependent phase discrepancies between observed and predicted seismograms within automatically selected windows from FLEXWIN 79 . The total misfit is defined as follows where N c denotes the total number of categories, here N c = 6 for three-component body and surface wave contributions. w c is the weighting factor used to balance the magnitudes of each individual category. w m represents the weighting factor used for the multi-taper measurements. Functions Δτ m (ω) and σ ϕ m ðωÞ represent the frequency-dependent phase anomalies and their associated uncertainties, respectively.
Model parameterization. There are 13 independent model parameters in a typical tomographic study for simultaneously constraining radial and azimuthal anisotropy 18,80 . However, due to the non-uniqueness of geophysical inverse problems and limited sensitivity of seismic data, not all of them can be well resolved in real applications. Previous theoretical studies have suggested that Rayleigh waves are mostly sensitive to 2θ azimuthal dependent coefficients G c and G s 81,82 , whereas Love waves have weak sensitivity to E c and E s . Here subscripts c and s stand for the cosine and sine terms, respectively, and θ denotes the local azimuth. In this study, we combine three-component surface and body waves to simultaneously constrain four model parameters: two radially anisotropic parameters L and N, and two azimuthally anisotropic parameters G c and G s 39 .
Thus, the relative perturbation of the total misfit, δJ, can be expressed as the following volume integral These four model parameters are related with the fourth-order elastic tensor C ij as Misfit gradients. The misfit gradients with respect to the four model parameters L, N, G c and G s can be calculated as follows 83,84 K L ¼ K c 44 þ K c 55 ; K N ¼ K c 66 À 2K c 12 ; Here the primary gradient K cij can be computed via the adjoint state method 85 as where ε j and ε y i are the elements of the strain tensors for the forward and adjoint wavefields, respectively.
A preconditioned conjugate gradient method 86 is used to iteratively update the above four model parameters. The final isotropic shear wavespeed V s and radially anisotropic parameter ξ can be calculated from L and N as With the distributions of G c and G s , the peak-to-peak anisotropic strength G 0 and fast axis direction ϕ are computed via Data availability All the continuous seismic data are collected from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (http://www.iris.edu/ds/nodes/dmc/). The digital model for US 32 can be downloaded from the attached data file.

Code availability
The open source spectral-element software package SPECFEM3D_GLOBE and the seismic measurement software package FLEXWIN used in this study are freely available for download via the Computational Infrastructure for Geodynamics (CIG; geodynamics. org). This is UTD Geoscience contribution number 1357.