Subducted Lithosphere Under South America From Multifrequency P Wave Tomography

We analyze mantle structure under South America in the DETOX‐P1 seismic tomography model, a global‐scale, multifrequency inversion of teleseismic P waves. DETOX‐P1 inverts the most extensive data set of broadband, waveform‐based traveltime measurements to date, complemented by analyst‐picked traveltimes from the ISC‐EHB catalog. The mantle under South America is sampled by ∼665,000 cross‐correlation traveltimes measured on 529 South American broadband stations and on 5,389 stations elsewhere. By their locations, depths, and geometries, we distinguish four high‐velocity provinces under South America, interpreted as subducted lithosphere (“slabs”). The deepest (∼1,800–1,200 km depth) and shallowest (<600 km) slab provinces are observed beneath the Andean Cordillera near the continent’s northwest coast. At intermediate depths (1,200–900 km, 900–600 km), two slab provinces are observed farther east, under Brazil, Bolivia and Venezuela, with links to the Caribbean. We interpret the slabs relative to South America’s paleo‐position over time, exploring the hypothesis that slabs sank essentially vertically after widening by viscous deformation in the mantle transition zone. The shallowest slab province carries the geometric imprint of the continental margin and represents ocean‐beneath‐continent subduction during Cenozoic times. The deepest, farthest west slab complex formed under intra‐oceanic trenches during late Jurassic and Cretaceous times, far west of South America’s paleo‐position adjoined to Africa. The two intermediate slab complexes record the Cretaceous transition from westward intra‐oceanic subduction to eastward subduction beneath South America. This geophysical inference matches geologic records of the transition from Jura‐Cretaceous, extensional “intra‐arc” basins to basin inversion and onset of the modern Andean arc ∼85 Ma.

that cap the Andean mountain belt, and are a product of the subduction of oceanic plates. The modern Andes (Figure 1) span the west coast of South America, stretching from the Caribbean Plate to the northern Scotia Plate (Ramos & Folguera, 2009). Studies of the area have shaped geologic thinking so that the term "Andean-style subduction" often refers to ocean-beneath-continent subduction as a general process.
Seismic illumination of the mantle beneath South America is variable, a consequence of very uneven distribution of seismic sources and receivers across the region. The active western margin produces a coast-parallel belt of abundant, strong seismicity beneath the Andes. Instrumentation efforts have focused on this region, often motivated by its serious seismic hazards. The largely aseismic interior and the eastern coastal parts of the continent are much more sparsely instrumented, in particular the vast and remote Amazon and Orinoco basins (Figure 1).
Most tomography studies of South America have produced regional, upper-mantle models beneath various parts of the Andes (Bianchi et al., 2013;Heintz et al., 2005;Marot et al., 2014;Scire et al., 2016;Ward et al., 2013). Along strike they have resolved the actively subducting Nazca slab in the upper mantle, which plunges eastward beneath the Andes at normal and flat dip angles (Marot et al., 2014;Porter et al., 2012;Scire et al., 2016;Wagner et al., 2005); and seismically slow anomalies above and below the Peruvian and Chilean flat slabs (Portner et al., 2017;Scire et al., 2016). Tomography has also illuminated the upper mantle around the South American-Caribbean plate boundary (Bezada et al., 2010;Miller et al., 2009) and beneath the Atlantic/South America boundary region (Celli et al., 2020).
Only a few tomography studies in South America have been conducted at the continental scale, notably the surface-wave tomography of Heintz et al. (2005) and the body-wave tomographies of Grand (1994), Ren et al. (2007), and Portner et al. (2020). Much of what is known about lower-mantle structure under South America has been learned from global-scale tomography models (e.g., Amaru, 2007;Bijwaard et al., 1998;Li et al., 2008;Lu et al., 2019;Montelli et al., 2006;Obayashi et al., 2013;Simmons et al., 2012). These studies have traced the eastward-dipping, margin-parallel slab curtain of the Nazca plate (and ancestral Farallon plate) down to ∼900 km depth, with significant fragmentation along strike.
geometries. The Northern Andes contain sequences of oceanic terranes that accreted from Upper Jurassic to Eocene times (e.g., Braz et al., 2018;Kennan & Pindell, 2009;Spikings et al., 2015). This contrasts with the Central Andes, which lack evidence for Jurassic and younger terrane accretions (e.g., Sempere et al., 2008), instead containing sedimentary and volcanic formations interpreted as reflecting deposition in a convergent-margin setting since early Paleozoic times (Pepper et al., 2016).
The model DETOX-P1 (Hosseini et al., 2020) that is described, resolution-tested, and analyzed here is a global, teleseismic P wave tomography model, constructed with particular attention to South American regional data. We were motivated to better resolve and understand the discontinuous and scattered geometries of lower-mantle slabs that had become apparent in global tomographies.
This study represents the first regional-scale and "tomotectonic" interrogation of this recent global inversion, against a backdrop of only a few comparable, global P-velocity models published over the past decade (Lu et al., 2019;Obayashi et al., 2013;Tesoniero et al., 2015), during which the number of seismic broadband stations has grown enormously. Section 2 gives an overview of the DETOX-P1 model, its strengths, and limitations. Section 3 presents the P wave traveltime measurements made, with a focus on South America, and Section 4 summarizes the tomography method. Section 5 offers a descriptive survey of slab geometries under South America. In Section 6, we attempt to reconcile these observations with plate kinematic models and geological records from mid-Jurassic (∼180 Ma) to Oligocene times (∼30 Ma), exploring the null hypothesis that subducting lithosphere essentially sank in place, that is, displaced laterally by no more than a few hundred kilometers since entering its paleo-trench (Steinberger et al., 2012;Sigloch & Mihalynuk, 2013van der Meer et al., 2010), such that the locations and strikes of slabs faithfully record their paleo-trench geometry. This is the simplest, physically reasonable slab sinking hypothesis, but the implications of such uniform mantle behavior and the powerful predictive capabilities that arise from it motivate us to test whether it can be rejected.

Teleseismic Global Tomography Model DETOX-P1
Seismic tomography, and more specifically the inversion of teleseismic body waves, has long been the preferred tool for inferring mantle structure below lithospheric depths in contexts where high spatial resolution is of the essence. For mantle regions well sampled by body waves, P wave inversions usually offer higher resolution thanks to their shorter wavelengths and much more abundant observations. Global tomography models are very large data gathering and processing efforts. Therefore new issues or fundamental revisions of models are relatively infrequent. Recent global P wave tomographies published before DETOX-P1 (Hosseini et al., 2020) were GAP-P4 by Obayashi et al. (2013), SPani-P (Tesoniero et al., 2015), TX2019slab-P by Lu et al. (2019) and LLNL-G3D-JPS (Simmons et al., 2019). Between updates or new issues, global models typically see extensive use for tectonic and geodynamic interpretation, either with or without re-inversion of the data in order to test for hypothesized mantle features ("resolution tests").
DETOX-P1 was produced as the first and least comprehensive of three global P-velocity models (Hosseini et al., 2020). DETOX-P2 made a major methodological advance by adding core-diffracted P waves ("Pdiff") to the "conventional" teleseismic P wave observations of DETOX-P1; on top of those, DETOX-P3 included PP observations at (two times) teleseismic distances. The purpose was to fundamentally improve the sampling of the mantle's lowermost third, where P waves and their interaction with the Earth's core had been poorly modeled previously, casting doubts on any deep, structural P wave results. Hosseini et al. (2020) accordingly focused on DETOX-P1 as a starting point against which the enhanced deep resolution of DE-TOX-P2 and DETOX-P3 compared favorably.
Yet DETOX-P1 is usually the superior model in the upper half of the mantle, in that slabs of subducted lithosphere are more crisply resolved, less blobby, and less affected by smearing artifacts. The reason is that the teleseismic-only observations of DETOX-P1 are more straightforward to regularize, compared to the heterogeneous P/Pdiff/PP data combinations in DETOX-P2 and-P3, where all parts of the mantle are not easily brought into focus simultaneously. Roughly speaking, DETOX-P1 usually seems clearer above ∼1,500-1,800 km depth, the most relevant depth range for this study. Below 2,000 km depth, DETOX-P2/3 is clearly better because DETOX-P1, robustly damped by design, fades below these depths.
MOHAMMADZAHERI ET AL. Hence for the investigation of upper and mid-mantle structure, we recommend relying on DETOX-P1 first and foremost, with a supporting role for DETOX-P2/-P3 where slabs need to be followed into the deepest mantle. Here we demonstrate the first such application, supported by extensive resolution testing.
Below 1,500 km depth, the mantle under South America, the focus of our study, contains very little seismically fast material (interpreted as subducted slabs), in marked contrast to the surrounding ocean basins and to Central/North America (Hosseini et al., 2020). Hence, we focus on DETOX-P1 for the technical part of this paper and for the bulk of its interpretations. In the Discussion, we draw on two slabs resolved by DETOX-P3 under the South/Equatorial Atlantic to extend our tomotectonic analysis based on DETOX-P1 to deeper mantle depths.

Traveltime Data for Tomography
Our teleseismic P wave tomography maximizes resolution under the sparsely instrumented South American continent by combining modern broadband data in a waveform-based, multifrequency inversion with the much longer instrumental record of high-frequency, picked traveltimes. Like regional studies, we put special emphasis on South American seismometer networks. Like global tomographies, we attempted to include all suitable broadband stations located within teleseismic distances of South America, usually when recording South American earthquakes. This permits imaging to greater depths than typical regional studies, to approximately 1,800 km. Figure 2 shows the teleseismic sources and receivers included in the DETOX-P1 data set, together with the number of traveltime measurements obtained at each receiver. Figures 2a and 2b show the earthquake locations used for finite-frequency waveform measurements and for ISC-EHB picks (Weston et al., 2018), respectively. Figures 2c and 2e show global stations and South American, respectively, that yielded finite-frequency waveform measurements, as described in Section 2.2. Figures 2d and 2f represent analyst-picked P-arrival time data, a subset of the new ISC-EHB catalog (Weston et al., 2018) selected to be as complementary as possible to the waveform measurements (details in Hosseini et al., 2020).

Traveltime Picks From the International Seismological Centre
In South America, both earthquakes and stations are clustered along the active, western margin, and are sparse across large swaths of the interior and eastern passive margin. Arrival time picks by human analysts for regional network data extend back three decades longer than broadband stations, and they include data from short-period seismometers and from broadband networks not set up to deliver waveforms to international data centers. Hence these picked traveltimes, curated and bundled by the ISC, can plug significant spatial data gaps in broadband coverage-compare Figures 2e and 2f. It is technically straightforward to invert traveltime picks, modeled by ray theory, jointly with finite-frequency measurements (Montelli et al., 2004;Sigloch et al., 2008).
From the new ISC-EHB catalog (Weston et al., 2018), we adopted ∼128,000 P-arrival picks measured by 1,075 stations in South and Central America (Figure 2f), and another ∼812,000 P wave picks from 4,549 stations at teleseismic distances outside the study area ( Figure 2d). The events occurred between 1967 and 2014, source magnitudes exceeded 5.6, and only arrival time readings of high precision (<0.1 s error) were accepted. The ISC stations partly coincide with stations used for broadband measurements (compare Figures 2c-2f).

Finite-Frequency Traveltime Measurements
For DETOX-P1, we measured the largest data set of teleseismic, multifrequency P wave traveltime anomalies to date (Hosseini et al., 2020). Measuring a multifrequency traveltime anomaly involves cross-correlation of an observed waveform with its synthetic counterpart in multiple frequency passbands. such earthquakes recorded by regional and global networks between 1990 and 2016. The usability of an earthquake for tomography depended on the deconvolution of a viable STF (see below). For all networks that were retrievable through the "web services" protocol of the International Federation of Digital Seismograph Networks (FDSN), the broadband seismograms were downloaded from international data centers, bandpass-filtered between 0.01 and 3.5 Hz, and corrected for instrument responses using the open-source, Python-based tool suite obspyDMT (Hosseini & Sigloch, 2017).
Within South America, our broadband data set contains all permanent and temporary networks available through the IRIS DMC, through European data centers including GEOFON Data Centre (1993) and GEO-SCOPE (1982), and through regional centers such as the Universidade de Sao Paulo Seismological Center. Table S1 summarizes the South American networks used, comprising 529 stations in South and Central America, between longitudes 111°W-35°W and latitudes 20°N-55°S.
Where studies with regional station coverage make traveltime measurements by cross-correlation, they usually correlate seismograms with spatially adjacent seismograms (Portner et al., 2020;Scire, 2015;Scire et al., 2016). This straightforward, computationally inexpensive approach is not possible for the global spread of our aggregated "network" because waveforms for a given earthquake vary strongly across the Earth's surface as a function of the radiation pattern. Specifically, over 75% of events occurred at depths shallower than 50 km. For shallow earthquakes, the first-arriving P pulse is followed within seconds by its depth phases pP and sP, which have large amplitudes of their own, and different radiation patterns. Traveltime measurements from data-to-data correlations would be biased by these systematic, unmodelled shifts of the three radiation patterns, especially in areas of large interstation distances, which are common across South America. Hence, we explicitly model the depth phases and correlate each observed seismogram with its predicted counterpart, consisting of a synthetic Green's function convolved with an estimate of the broadband STF. This measurement method is adapted from Sigloch and Nolet (2006) and Hosseini and Sigloch (2015).
Broadband Green's functions were obtained using the Instaseis software (van Driel et al., 2015), a look-up database of synthetic seismograms for regional and global distances, which were computed by the full-waveform, axisymmetric spectral element method AxiSEM (Nissen-Meyer et al., 2014), up to a dominant period of 2 s. Our velocity reference model is IASP91 (Kennett & Engdahl, 1991), complemented by the Q model of PREM (Dziewonski & Anderson, 1981).
For the moderate to large events we use (M w between 5.8 and ∼7.8), earthquake ruptures and hence STFs are several seconds in duration. Such STFs blur together the tightly spaced P, pP, and sP pulses in the Green's functions they convolve. Hence, a well-correlated synthetic seismogram needs to incorporate an accurate estimate of the event's STF. Our iterative, linearized estimation (deconvolution) of the broadband STF also yields estimates for source depth and for the moment tensor (Sigloch & Nolet, 2006;Stähler & Sigloch, 2014).
Mohammadzaheri (2019) estimated broadband STFs, depths, and focal mechanisms for ∼2,000 earthquakes of magnitude M w > 5.8, recorded between 2008 and 2016. They were deconvolved from seismograms of the Global Seismographic Network, which has a relatively even spatial spread. This data set was complemented by STFs estimated by Sigloch (2011) for earthquakes recorded between 1990-2008, processed by the same method. In total, 3,472 STFs were acceptable (epicenters shown in Figure 2a).
For each source-receiver pair, the broadband Green's function was convolved with the event's broadband STF to produce a synthetic seismogram. P wave traveltime anomalies were estimated by cross-correlating observed waveforms with synthetic waveforms in eight frequency passbands, using a Python-based tool suite by Hosseini and Sigloch (2015). The center periods of the Gabor (log-Gaussian) passband filters ranged between 30 and 2.7 s, spaced by factors of 2 (i.e., filter of constant fractional bandwidth).
MOHAMMADZAHERI ET AL.  The finite-frequency traveltime anomaly is defined as the time shift between observed and modeled (bandpassed) seismograms for which the cross-correlation function attains its maximum value xcorr (Dahlen et al., 2000), a value that can in principle range between 1.0 and −1.0. We accept traveltime measurement for tomography if xcorr > 0.8. Figure 3a displays histograms of the cross-correlation coefficients obtained for each frequency band. The number of successful measurements decreases with the center period and drops precipitously for periods shorter than 2.7 s. This high-frequency limit seems to be mainly due to systematic noise (source-side and receiver-side reverberations) in the observed waveforms. Measurement errors on the traveltime anomalies were assigned semi-empirically based on the value of xcorr (Sigloch, 2011).
The global, teleseismic, finite-frequency data set thus produced contains ∼5 million acceptable traveltime delays in the eight overlapping frequency passbands. Only ∼2.5 million traveltime measurements from four nonoverlapping frequency bands (center periods 30.0, 15.0, 7.5, and 3.8 s) were retained to manage computational cost during inversion. These traveltime anomalies were corrected for station topography, crustal structure, and ellipticity before inversion. Their source and station locations are shown in Figures 2a and 2c, respectively. This teleseismic P wave data set was also used by Hosseini et al. (2020) for joint inversion with Pdiff and PP waves (producing models DETOX-P2 and DETOX-P3).
The set of successful measurements within South America (Figure 2e) consists of 61,605 traveltime delays from 529 broadband stations in the four frequency bands. 603,465 additional traveltime anomalies were obtained from receivers at teleseismic distances recording 684 earthquakes in and around South America (Figure 2c).

Linear Inversion Under the Finite-Frequency Approximation
In finite-frequency tomography, we solve a system of linear equations with data index i running over the range of source-receiver-passband combinations. dT i is a traveltime anomaly and dV P /V P is the quantity to solve for-the deviation of P wave velocity V P from spherically symmetric reference earth model IASP91 (Kennett & Engdahl, 1991) Laske et al. (2001). The temporal spread of these raw data is due not just to 3-D velocity anomalies in the mantle but also, and indeed mostly, to systematic mismatches between source depth and origin time at this processing stage. Compared to initial seismicity catalog values, source depth has already been re-estimated (as part of the STF deconvolution), whereas origin time gets re-adjusted only later during the tomography stage.
mantle. The Fréchet, finite-frequency or "banana-donut" kernel   Κ i x quantifies the sensitivity of the observation dT i to mantle structures dV P /V P . We compute kernels using the method of Dahlen et al. (2000) and software by Tian et al. (2007).
Discretizing Equation 1 onto an irregular grid of tetrahedra with linear interpolation between vertices, we obtain the linear system where index j runs over grid vertices, Vj is the P-velocity anomaly at vertex j and A ij results from the quadrature of K i onto vertex j (Montelli et al., 2004). The tetrahedral mesh is adaptive (Persson & Strang, 2004), with vertices spaced more closely in densely sampled areas, as discussed in Hosseini (2016) and Hosseini et al. (2020). The facet length of the tetrahedra ranges between ∼80 km in the upper mantle under highly sampled regions and ∼400 km near the core in the lower mantle (see Figure 5 in Hosseini et al., 2020). Figure S3 shows map views of the volumes of Voronoi cells encasing the tetrahedra' vertices j. The volume enclosing each vertex depends on the edge length of its closest tetrahedra. A Voronoi cell is closely similar to the volume that contributes to the sensitivity kernel values A ij . Figure S4 displays, for all vertices j located beneath South America, the kernel   sum i ij A , which is a rough proxy for the goodness of sampling and thus resolution of a location j by the seismic wave paths available to this study. Comparison of Figures S3 and S4 shows how cell volumes are smaller where the seismic wave paths sample more densely (the proxy for sampling density is the column sums of the sensitivity kernel matrix). The purpose of this adaptive densification of the inversion grid is to honor all spatial detail delivered by the data in well-sampled regions while stabilizing the ill-posed part of the problem arising from poorly sampled areas.
The tomography simultaneously inverts for the V j and for four scalars "event corrections" per earthquake (dlat, dlon, ddepth, origin time), expressed as deviations from initial catalog values and subject to empirical prior uncertainties as follows:   lat     lon depth 10 km and   0 2 T s. For multifrequency measurements, the effective prior uncertainty in source depth is smaller than 10 km (because much care was taken to estimate depth during the STF inversion step) and ∼10 km for ISC measurements. We do not implement station corrections but rather correct for crustal structure a priori using model CRUST2.0 by Laske et al. (2001). On average, all earthquakes moved ∼3 km after inversion in the DETOX-P1 model.
The partly overdetermined and partly underdetermined linear system is solved in a least-squares sense by adding a weighted combination of norm damping and second-derivative smoothing (detailed in Hosseini et al., 2020). The full linear system becomes   Figure S5 displays the trade-off curve between data fit and regularization ("L-curve") at 1,000 km depth, showing how solutions change with changing regularization parameters. The preferred solution presented in Sections 4 and 5 features   2 1.152. Figure S5 shows it to be in the regularization regime where the solution has robustly stabilized. This is the solution shown by Hosseini et al. (2020) as the "DETOX-P1" model and contrasted to models DETOX-P2 and-P3, which were regularized to almost identical chi-square misfit values and additionally inverted Pdiff and PP wave traveltimes.

Resolution Tests
We performed synthetic resolution tests using regularly spaced Gaussian balls that fill the depth of the mantle column (Figures 4 and 5). The maximum velocity anomaly  3 % j V , occurs at a ball's center and decreases with distance according to a bell curve. The "radius," chosen as 200 km throughout the mantle column, is the distance where the velocity anomaly has decreased to a value of  max 1 , where e is the base of the natural logarithm. Before inverting the synthetic traveltime data, we added a realistic amount of white Gaussian noise, with a standard deviation of 0.35 s.
Two patterns have been used in the "checker-ball" resolution tests shown in Figures  , we conclude that the DETOX-P1 model offers its best resolution down to ∼1,800 km. For structures below 1,500 km depth, we compare to DETOX-P3, which specializes on resolving the lower and lowermost mantle (Figures S12-S27). Figure S11 presents a checker-ball test for DETOX-P3 from 300 to 2,700 km depth. Figure 4 shows that in the upper mantle, anomalies of 200 km are well resolved beneath the western half of the continent but gradually less so toward the east and south. Coverage and resolution improve from the transition zone downward, where the P waves travel increasingly horizontally with more abundant crossings. In the lower mantle, resolution and coverage remain very good down to 1,500-1,700 km under the continent ( Figures 5 and S4).
Elongated smearing artifacts are observed beneath the noninstrumented Atlantic, parallel to the dominant ray directions, which run toward the northeast and southeast, away from Andean seismicity. Similar smearing effects are present beneath the Pacific basin between 1,400 and 1,500 km depth. We adjust for these oceanward artifacts in our interpretations.

Tomography Results: Mantle Provinces Under South America
This section presents the results of our tomographic inversion, the DETOX-P1 model. Figure 6 shows map sections between 300 and 800 km depth, Figure 7 between 900 and 1,600 km. Figure  In the map sections of Figures 6 and 7, fast anomalies tend to be narrowly elongated, consistent with the geometries of (paleo-)trenches. Along strike, fast anomalies are often clearly segmented on lateral scales of 1,000-3,000 km, consistent with the segmentation of present-day trench networks into arc segments. Individual fast segments are still linearly elongate, or sometimes arcuate (beneath the Lesser Antilles subduction zone; Figures 6a-6c). Fast anomalies in the shallow upper mantle show a tight spatial correlation with active subduction zones (Lesser Antilles, Central America, Northern Andes, and Central Andes), consistent with their interpretation as subducted lithosphere ("slabs"). Integration of linear anomalies across depths shows slabs to have a tabular form in 3-D, consistent with subduction of lithosphere along a trench, over a period of geologic time. By extrapolation, we also interpret deeper fast anomalies as subducted lithosphere, including those now disconnected from the surface.
We consider a fast anomaly as a robustly interpretable slab of subducted lithosphere if and where its amplitude exceeds dVp/Vp > +0.25%, relative to reference model IASP91. This relatively weak threshold takes into account the observation that few fast anomalies exceed 1% amplitude while staying spatially connected, especially not below the transition zone. This is true for our own and most other global models, see model comparison Figures 9-12 and Figures S12-S27. For conceptual simplicity and clarity, we choose a single slab interpretation threshold of dVp/Vp > 0.25% for all depths, which is about the strongest dVp/Vp level at which fast anomalies generally tend to present as linearly elongated and interconnected belts at any mantle depth. Such geometric characteristics match those of present-day trench networks at the surface and of their associated, recently subducted slabs in the uppermost mantle. The conceptually simplest interpretative approach entails a baseline expectation that lower-mantle slabs are just older, sunken versions of such recent, shallow slabs. At larger velocity thresholds, fast anomalies appear spatially disconnected and blobby, no longer lending themselves to interpreting paleo-trenches. Figure S5, a visualization of candidate tomography solutions along a regularization trade-off curve, further details this rationale and its robustness.
The slabs visible in Figures 6 and 7 fall into four main provinces according to their clustering in depth, latitude, and strike. Sections 5.1-5.4 are dedicated to these four slab provinces, progressing from shallow to deep, or from the more recently subducted lithosphere to the older plates. We introduce nomenclature for the arc-scale segments of the massive, eastward-dipping slab complex under South America, since not all of the segments have been described separately. Unlike most previous authors, we do not name slabs based on their presumed plate origins, such as Nazca, Farallon, Cocos, because future, more complete paleo plate models are likely to involve additional plates. Rather we name slab units after surface features that have the same locations, strikes and extents, presumably due to a shared origin. We newly name three slabs beneath the Atlantic, which have either not been described or not named previously. Table S2 maps our regional slab units to the global slab nomenclature of van der Meer et al. (2018), who in turn discuss the nomenclature of earlier workers.

Slabs in the Upper Mantle (300-600 km Depth)
The mantle at 300 km depth (  Unlike beneath Peru, the Chilean flat slab segment is imaged by tomography, albeit with moderate smearing westward and downward. Note that each cross-section is visually subdivided by five evenly spaced green markers. Since the absolute lengths of the profiles vary, the absolute spacing between two green markers varies across sections.
interpret slabs above 300 km because the densities of wave-path crossings are unsatisfactory in most areas, due to a lack of wave-path coverage.  Top row: DETOX-P1 (this study); regional P wave model SAMS_P_2019 by Portner et al. (2020); LLNL_G3Dv3, a global body wave model (Simmons et al., 2012). Bottom row shows three global P wave models: MITP08 (Li et al., 2008), PRI-P05 (Montelli et al., 2006), and UU-P07 (Amaru, 2007 , the CA slab shifts eastward and the CoMa slab shifts south-eastward, that is, both slabs dip beneath the curved continental margin. At 300-400 km, a marked gap ("P") is present between ∼9°S and 1°N, that is, between the Ca and CoMa slab segments (Figures 6b and 6c). The gap disappears between 500 and 600 km depth (Figures 6c and 6d). At the surface, this Peruvian slab gap coincides with a paucity of seismic stations, raising doubts whether the apparent lack of slab might be an imaging artifact caused by norm damping. However, both our resolution ( Figure 4) and kernel densities ( Figure S4 and Figure 6) are good. Moreover, in this margin segment (between ∼6°S and ∼4°N), the ISC-EHB catalog (Weston et al., 2018) lists earthquakes from the surface to 200-250 km depth and again below ∼600 km, but none between 300 and 600 km (Figures 8b and 8c).
To further increase confidence that the Peruvian slab gap is a real feature, we implemented hypothesis-driven resolution tests ( Figure S7). The input model features a slab paralleling the Andean margin that is vertical and 100 km wide between 0 and 300 km depth, and eastward-dipping, 150 km wide between 300 and MOHAMMADZAHERI ET AL.  Figure S7) cut through the input model across the location of Peruvian gap. Sections a′ and b′ show that the continuous input anomaly is faithfully recovered under Ecuador and northern Peru. Although stations are lacking in northern Peru, wave paths propagating from regional seismicity (above 250 km depth) still sample the slab gap region, as Figure S4 shows. Hence all indications are that the gap is real, that is, that no slab is present at 300-600 km under Ecuador and northern Peru.
Below 600 km depth, slab beneath the Andean margin is continuous and strikes NNW along its entire length (Figures 6e and 6f  right angle. Hence, the slab gap seems to have accommodated trench divergence over the subduction times corresponding to 600-300 km depth. This inference is spatiotemporally consistent with an Eocene change in Farallon plate motion, from north-eastward to eastward, as recorded by magnetic isochrons. Hence, the southern flank of the Farallon plate may have veered toward the more easterly motion of the future Nazca plate well before 23 Ma (Lonsdale, 2005; and references therein). In contrast, Taboada et al. (2000)   Other tomography models at 400 km depth (Figure 9) essentially agree on the shallow slab arrangements described. Models PRI-P05 (Montelli et al., 2006) and LLNL-G3Dv3 (Simmons et al., 2012) show the Peruvian slab gap less clearly than other models, and they hint at southward slab continuation beneath the Atlantic, offshore Argentina. Similarly, the regional tomography by Portner et al. (2020) shows less of a slab gap than a thinning of fast anomalies. Additional tomography models (nine in total) are compared in Figures S12-S14, at 400-600 km depth.

Slabs From ∼600 to ∼900 km Depth
The dominant structure at 600-900 km is a continuous, NNW-striking slab inboard of the Andean margin ( Figures 6 and 7), which we term Trans-Andean slab (TA). This slab is the downward continuation of the CoMa and CA slabs, and with the Peruvian slab gap filled in below ∼600 km depth-we rename this merged structure in order to emphasize its latitudinal continuity (Trans-Andean). The slab extends ∼3,800 km between latitudes 5°N and 28°S and it dips eastward between 600 and 900 km depth. The simple linear geometry of the TA slab bears little resemblance to the present-day curvature of the continental margin, again in contrast to the shallower slab segments of Section 5.1. Most evidently absent is the imprint of the oceanward-concave Bolivian orocline, which is clearly expressed in CA slab above 500 km depth ( Figure 6, roughly between 10°S and 22°S).
Under Central America and the Caribbean, the CAM and ANT slabs are well expressed between 600 and 900 km; the Antilles slab is displaced to the south compared to shallower depths. No slab is imaged south of 30°S beneath Patagonia and the South Atlantic, with the exception of the South Sandwich arc. Imaging resolution south of 30°S is significantly better and broader in this depth range than in the upper mantle ( Figure 4 vs. Figure 5); the main weakness is SE-ward smearing. A targeted resolution test (Figures S8 and S9) shows that shape recovery for a slab below 600 km depth is achieved south of 30°S, although with diminished amplitude.
The model comparison at 800 km depth ( Figure 10) shows somewhat larger variance than in the upper mantle. The single, long Andean slab is evident in all global models, consistently terminating close to 30°S, beneath Argentina. None of the models shows well-defined fast anomalies farther south, except near the Scotia Sea (Sc slab of van der Meer et al., 2018) and the South Sandwich arc. The strike of the Andean slab varies perceptibly between models, especially along its northern half, where some models show it diffuse or trending north, instead of NNW. Presumably this reflects resolution limitations, specifically the smearing together of the NNW-trending Andean slab and the N-trending Antilles slab. Both slabs are well resolved individually by our study. The Portner et al. (2020) model shows more scattered small-scale structure across the continent than the other models. Since a physical interpretation in terms of widely scattered subduction products seems unlikely, this model's overall appearance may reflect less conservative regularization (smoothing and damping) of the inverse problem.

Fast Anomalies From ∼ 900 to ∼1,200 km Depth
Below 900 km, imaging resolution is good under the continent (except easternmost Brazil) and also under most of the adjacent ocean areas (Figures 5 and S6). The slab assemblage becomes markedly more fragmented and complex (Figures 7 and 8a). A first-order observation in the depth range of 900-1,200 km is that all slab geometries change. The Antilles (ANT) slab shrinks with depth, and its southern part merges with TA slab into a longitudinally widened slab body around the equator, presumably smearing together contributions from spatially close subduction zones. With depth, the southern tip of the TA slab fades: at 900 km depth, TA extends to ∼28°S, at 1,200 km it is limited to ∼18°S.
The deep part of the TA slab does not dip down to the east, unlike shallower TA; instead this main Andean slab is essentially vertical from 900 km to its lower end at ∼1,200 km. We return to this subtle but arguably important observation in Section 6.2. The farthest eastward extent of TA, which thus occurs at 900-1,200 km and at the southernmost TA segment, is ∼57°W.
Below 1,000 km, a different, westerly slab appears beneath the coast of Colombia, here named the Western Cordillera slab (WC, Figures 7c-7f). WC slab strikes north or NNW, and links up with the deep CAM.

10.1029/2020JB020704
20 of 35 In the model comparison at 1,200 km depth (Figure 11), all models pick up on these multi-slab complexities. The degree to which slabs are resolved separately and with elongated (trench-like) geometries is variable. No model shows TA slab continuing to south of ∼30°S, and in most models, it terminates by ∼20°S at these depths. The models also agree that TA is localized under west-central Brazil and does not extend under eastern Brazil. The mantle under Argentina and the Argentine basin remains free of well-defined slabs in all models. Some models feature a blue haze in the Argentine Basin, which forms a well-defined high-velocity domain at greater depth.
In light of plate reconstructions and widely held interpretations of persistent Andean arc magmatism, the area east of ∼57°W might be expected to hold slab as well (as discussed in Section 6.2). Hence, we need to be confident that TA slab does not in reality extend farther east than imaged, at the 900-1,200 km or deeper levels. The resolution tests in Figures 4 and 5 indicate that, for the relevant region of (eastern) Brazil, recovery amplitudes are somewhat weakened, but that shape recovery is excellent at 900-1,200 km, in fact better than at any other depth in DETOX-P1. A targeted slab resolution test bears this out ( Figures S8 and S9). It is implemented as horizontal slabs based on the broad longitudinal spread of slab material at 1,000-1,200 km depth (interpreted as "stagnating" in the lower mantle by Portner et al. (2020); we will offer the alternative interpretation of three concurrent trenches (WC, TA, and ANT) during a fundamental plate reorganization that presaged the modern Andean arc). The test results ( Figures S8 and S9) show that if TA slab extended farther east, this should be easily resolved. Resolving a dipping slab would be equally feasible since resolution remains strong to ∼1,800 km depth (Figures S12-S27).

Fast Anomalies Below ∼1,200 km Depth
Below 1,200 km, the combined Andean-Antilles slab (TA-ANT) rapidly fades, while the westerly WC slab persists beneath the Colombian coast and offshore. In Figure 7, the easterly, NNW-striking TA slab still extends south to ∼18°S at 1,200 km, but by 1,600 km depth it has retreated to around the equator as a very weak anomaly. Figure 5 demonstrates excellent resolution, implying that ∼1,600 km is the true lower limit of the Trans-Andean-margin slab. This result is confirmed by model DETOX-P3 in Figure S24 (Hosseini et al., 2020), which images more deeply into the mantle by including core-diffracted P-diff waves.
CAM slab remains clearly expressed to 1,600 km depth in DETOX-P1; WC slab remains robust to ∼2,200 km. In DETOX-P3, slab in the CAM/WC region continues to the core-mantle boundary in DETOX-P2 (Hosseini et al., 2020).
The model comparison at 1,500 km ( Figure 12) shows slabs to be confined to latitudes north of ∼20°S (WC/ CAM), except for an isolated, seismically fast blob under Central Brazil (centered on ∼25°S and ∼60°W) in PRI-P05 and UU-P07. The extended model comparison down to 2,200 km ( Figures S24-S27) shows how DETOX-P3, with its deep Pdiff wave constraints, evolves very differently from all other (teleseismic-only) models. Toward 2,200 km, all but DETOX-P3 become fainter under the action of regularization/norm damping and show less mutual agreement (compare WC slab, among the most robust features). DETOX-P3 features strong anomalies with elongated, slab-like strikes beneath WC, the southern Atlantic and eastern seaboard of Brazil. The direct comparison between DETOX-P1 and DETOX-P3 is particularly informative because they share everything, including the same robust regularization, except that DETOX-P3 adds constraints from the deepest-diving Pdiff (and PP) waves. Gutscher et al. (1999Gutscher et al. ( , 2000 related the seismicity gap under Peruvian flat slab to subducted topographic anomalies, i.e., Nazca ridge and Inca Plateau. In this scenario, the shape and the length of subducted topographic anomalies are similar to their conjugates on the preserved Pacific plate, having formed during symmetric seafloor spreading at the East Pacific Rise. Slab flattening and extinction of the overlying Andean volcanic arc followed the subduction of the thick plateau crust, with a compositional buoyancy that delayed sinking into the mantle. Although geodynamic modeling has shown that the subduction of a buoyant anomaly (e.g., an oceanic plateau, aseismic ridge or seamount chain) can cause a narrow flat slab, MOHAMMADZAHERI ET AL.

10.1029/2020JB020704
21 of 35 an additional mechanism was proposed to be required to explain flat slab segments under Peru and Central Chile (van Hunen et al., 2002). Geodynamic modeling by these authors shows that the thermal structure of the subducting plate also impacts its buoyancy, and critically, its strength and ability to bend.
The Inca Plateau had been proposed to have subducted ∼10-12 Ma (Gutscher et al., 1999). However, our depth-time conversion suggests that the slab gap appears under Ecuador and northern Peru in the Eocene, implying a different scenario. We propose that the slab gap reflects a tectonic reorganization, as follows. The Ecuadorian margin accreted an oceanic plateau (Piñon plateau) which plugged the subduction zone, forcing it to jump outboard to the western edge of the plateau during Eocene times (Kerr et al., 2002), or during Late Paleocene (∼58 Ma) (Jaillard et el., 2009). Its buoyant lithosphere may have clogged the trench instead of sinking like the seafloor to its north and south, during the times that correspond to slab now located between 300 and 600 km depth. Newer slab above 200 km in Figure 8b, as outlined (only) by seismicity, would then represent lithosphere subducted since the times that a rejuvenated trench stepped outboard (west), leaving the Piñon plateau accreted to the Ecuadorian margin, and a hole in the Andean slab reflecting the extent and time span of its nonsinking.

Synopsis of Slab Configurations Over Time: Two Fundamentally Different Regimes?
Figures 13a and 13b show all seismically fast anomalies of dVp/Vp > 0.25% located below 400 km and below 600 km, respectively. Tabular domains greater than ∼5° long are presumed slabs. The differences between panels (a) and (b) indicate where material subducted during the time period corresponding to 400-600 km slab depth; slab present in panel (b) but not in panel (c) subducted during the time period corresponding to 600-800 km slab depth, and so on. If slab sank at a uniform rate of, say, 10 mm/a, then 600-800 km deep slab would have subducted between 60 and 80 Ma. We do not assert any such depth-to-time mapping (or sinking rate) a priori. Rather we will estimate it from the observations, independently in various locations across the slab assemblage.
Slabs in blue and green shades (i.e., down to ∼1,000 km, Figures 13a-13c) resemble present-day trench geometries and subduction polarities at the surface. The Andean CA/TA slabs are elongated parallel to the western coast and have a concave-to-the-west curvature, similar but more open than the Bolivian orocline. The slab curvature increases to shallower depths (aqua to dark blue levels), consistent with a ∼45° clockwise rotation of the orocline's southern limb since Late Cretaceous (Maffione et al., 2009), and 29° since Eocene (Arriagada et al., 2008). The CA/TA slab dips eastward down (dark blue through green color levels). The Antilles slab (ANT) is arcuate like its overlying Lesser Antilles arc and dips westward. The shallow CAM matches the strike of the Cocos trench and dips eastward. Hence, the blue-green assemblage in Figure 13c shows the spatiotemporal continuation of present-day deposition styles of ocean lithosphere in the Andean, Cocos, and Antilles trenches, back to 1000-km depth-equivalent time.
At deeper, earth-toned levels in Figures 14a-14d, the slab assemblage is completely reconfigured. Below 1,200 km (Figure 14a), the Andean CA-to-TA slab is essentially gone. Instead, the dominant feature is a massive vertical slab wall extending north-northwest from the deep CAM and WC slabs under western Colombia. It crosses the Caribbean and tracks up the eastern seaboard of the United States, to Nova Scotia. Along its strike, this slab wall terminates upward at the salmon-pink level (800-1,000 km depth, best appreciated in Figures 13a-13c). Among the most voluminous and visible features in the entire mantle, this super-slab under the Americas has been resolved to varying detail by probably every global tomography in the past 20+ years, in addition to pertinent regional studies (e.g., Grand, 1994;Ren et al., 2007;Sigloch, 2011). Almost all authors have interpreted this slab as eastward subduction of the Farallon/Kula plate under the North American continent and Central American blocks, until Sigloch andMihalynuk (2013, 2017) reinterpreted the lower, vertical parts as intra-oceanic, westward subduction of Mezcalera Ocean (part of the North American plate) under the Insular/Guerrero microcontinent. North America, pulled westward by those trenches, ultimately collided with their arc terranes, which forced a flip of subduction polarity to eastward beneath the reconfigured continental margin, i.e., the establishment of "Andean-style" Farallon subduction beneath North America. At South American latitudes, the super-slab's southernmost segment, WC slab, has MOHAMMADZAHERI ET AL.

Problem Statement: Attributing Complexity to the Surface Versus the Mantle?
The challenge here is to reconcile the earth-toned slab assemblage (below 800-1,000 km) in Figures 13d  and 14 with Andean subduction history ( Figure 15). The widely held perception is that the Andean margin hosted a continental arc since the Jurassic (e.g., Ramos, 2018) or latest Triassic (209 Ma; Spikings et al., 2015), although widely acknowledged to be in extension until Late Cretaceous times (e.g., Ramos, 2018). Figure Torsvik et al., 2019), locating South America 4,000 km to the east of the present coastline by 140 Ma. At the same time, slab that reaches deeper than the east-dipping TA slab is found only in the opposite direction, far to the west. Hence, the continent and the presumed areas of subduction seemingly diverge into the distant past. The longitudinal translation of the slab system, from farthest east to farthest west, occurs quite abruptly across the depth range of 1,000-1,200 km, suggesting a rapid and fundamental reorganization of the plate system in the depth-equivalent time period, from an earlier system that was closely tied to Central and North American events, to the modern, Andean-style trench system. Resolution tests ( Figure 5 and Figures S8-S9) show that the Andean-margin TA slab does not reach farther east than ∼57°W because a hypothetical eastward continuation would be readily resolved. Yet quantitative plate reconstructions in absolute hotspot reference frames place the Andean continental paleo-margin significantly farther east than ∼57°W for absolute times exceeding ∼75 Ma (Seton et al., 2012) (Matthews et al., 2016), as seen in Figure 15. When South America broke away from Africa ∼120-140 Ma, its western margin lay >1,500 km east of ∼57°W and the TA slab. The intervening area is slab-free at 1,000-1,600 km in DETOX-P1. Which slab, therefore, corresponds to presumed marginal subduction at times before 75 or 85 Ma? This is especially puzzling given that the TA slab (east of 57°W, above 900 km) conforms closely to the expected shape of slab deposited beneath the continental margin. If a trench existed along the Andean margin prior to the times during which the blue-aqua slab levels in Figures 13 and 15 were deposited, then mantle advection behavior must have changed massively, for unknown reasons. Yet slab geometries below 900 km (Figures 13d-16a) appear equally structured (elongate, narrow, vertical) as slabs above 900 km, and plausible, especially as they represent seamless extensions of the North/ Central American configuration.
Plate reconstructions (Matthews et al., 2016;Seton et al., 2012) first position the Andean margin above the easternmost south TA slab between 75 and 85 Ma. If lithosphere had been deposited beneath the Andean margin before 75-85 Ma, then that lithosphere could not have sunk nearly vertically, unlike more recent Andean-margin slabs-unless the reconstructions had South America's position wrong (too easterly) by >1,500 km. This would imply that the Indo-Atlantic hotspots were not nearly fixed before ∼80 Ma, for which there is little evidence (Butterworth et al., 2014;O'Neill et al., 2005). This postulation would discard the most reliable workhorse of paleo-reconstruction. In summary, neither highly variable slab sinking styles nor highly mobile hotspots (only) before ∼80 Ma are attractive hypotheses due to their complexity and the cascading issues they would cause for paleo-reconstruction globally.

Vertical Sinking Hypothesis for Viscously Thickened Slabs
As always, one should attempt to reject the simplest hypothesis first. From a mantle physics perspective, the simplest hypothesis is that slabs below 900 km reflect trench geometries equally faithfully as slabs above 900 km. The simplest hypothesis that can in turn be made about the deposition of the aqua to dark blue slab levels (above 800 km and below ∼300 km) is that they sank vertically down and are still located directly beneath the absolute coordinates where they originally entered the mantle (within observational uncertainties of a few hundred kilometers laterally, as will become clear). Slabs that did not displace laterally relative to the lower mantle also did not displace laterally relative to the hotspots (to the extent that the hotspots were fixed). Under this hypothesis, a mantle column that is slab-free should not have experienced a trench passing above. A stationary trench would produce a dip-less, vertical slab beneath it. For a dipping slab, the vertical sinking hypothesis implies that its dip was produced by the absolute lateral motion of its paleo-trench-westward in the case of the Andean margin and Central America, eastward in the case of the Antilles arc (Matthews et al., 2016;Seton et al., 2012). Specifically, if i denotes the slab's dip angle measured relative to the surface, then tan(i) = v sink / v trench , where v sink is the slab sinking velocity (vertically down) and v trench the trench migration velocity relative to the lower mantle (oceanward rollback). Note that the trench convergence rate v conv , the velocity at which the incoming plate converges on the trench, does not feature. This is taking a viscous view of the subduction process: if we applied this approximation starting from the surface, it would imply that the moment the subducting plate enters the mantle, it turns viscous and drips vertically down beneath the trench at the rate v sink , while the trench moves horizontally away from this drip at the rate v trench . The trench convergence rate v conv of a plate is typically 3-10 times faster than v trench , where typical values would be 50 mm/a versus 10 mm/a (Goes et al., 2011;Sdrolias & Müller, 2006). In this example, a cross-sectional 500-km long swatch of lithosphere would enter the mantle every 10 m.y. and drip vertically down in the 100-km wide column of accommodation space created by the retreating trench during this period. In reality, the transition from MOHAMMADZAHERI ET AL.  Figure 13). The quantitative plate reconstruction of Matthews et al. (2016) is superimposed at 30 Ma (dark blue coastlines), 70 Ma (aqua green), and 120 Ma (dark red). The line color associated with a reconstruction time is chosen to match the color of slab that subducted at this time, assuming the slab sinking rate was 10 mm/year. Hence slab that subducted at 120 Ma would have sunk to 1,200 km depth (dark red). The Andean margin lines up well with slab at the blue and aqua green levels but not at the red level since all slabs are located much farther west than the Andean margin at 120 Ma. Figure S28   Present-day continents are shaded gray. White arrows denote two ray-shaped smearing artifacts that should be ignored but make no difference to the interpretation. The basin (shaded cyan) occupies the space between the 140-Ma margin and the closest slab(s) to the west, i.e., AgB and TA. (b) Cartoon showing the subduction sequence inferred along the red cross-sectional line in (a). From ∼140 to 0 Ma, SAM moves continually westward in a lower-mantle reference frame. From 140 to ∼85 Ma, TA trench sits stationary above deepest TA slab (salmon-colored arrow). SAM approaches from the east while TA trench consumes the basin lithosphere into the deepest TA slab. Following collision/accretion by SAM ∼85 Ma, subduction into TA trench flips to eastward beneath the SAM margin-the birth of the modern Andean arc. All slab geometries are explained by vertical sinking (after viscous thickening). Slab dip (or not) is the consequence of trench motion, not of lateral slab displacement within the mantle. rigid to viscous will be gradual, but it remains the case that each vertical column ultimately needs to accommodate N-times more lithosphere than was originally overlying it at the surface, where N = v conv /v sink . Ribe et al. (2007) and Gibert et al. (2012) modeled how softening and thickening can occur by slab folding in the transition zone, against the high-viscosity backstop of the lower mantle. The tightness and stacking of the folds adjust to accommodate the thickening ratio N. This process is cartooned in Figure 16b for the Andean margin slab, followed by vertical sinking once the folds have been laid down (Čížková et al., 2012).
Vertical sinking will not apply in the upper few hundred kilometers, where the slab is still rigid and where its seismicity shows no signs of slab thickening. However, an inspection of the combined seismicity and tomography of the Andean slab in Figure 8 shows the slab fundamentally changing between 300 and 600 km, consistent with the hypothesized behavior. In every panel (except for b and c, in the Peruvian slab gap), we observe a shallower slab dip in the uppermost mantle, clearly delineated by seismicity, and a much steeper dip deeper down, as indicated by the imaged slab and in most panels also by deep seismicity. In Figures 8d and 8e, abundant seismicity shows the transition from shallow to steep dip particularly clearly in a seismicity kink occurring around 200 km depth. For the vertical viscous sinking hypothesis, this seismicity kink likely signals the point from which the slab starts to deform viscously.
In panel 8a, the slab is imaged through all depths and appears much fainter (light blue) and thinner above 300 km depth than below. In the remaining panels, slab in the uppermost mantle is essentially not resolved, whereas from the transition zone down, it appears as a robust fast anomaly, hundreds of kilometers wide.
Clearly, the slab is widening several-fold as its passes through the transition zone. This cannot represent an imaging artifact where a non-thickened slab became laterally smeared at deeper depths. In principle, near-surface structure is easier to resolve than deeper anomalies. In the Andes, the near-trench regions also benefit from the region's best source and receiver coverage (Figure 2). Despite these advantages, the shallow slab is not resolved whereas the deep slab is. In addition, the fast dVp/Vp anomaly in Figure 8a has a higher amplitude in the deep, thickened part than in the thinner near-surface part-this difference must be due to slab thickening, since it cannot be due to slab cooling with depth (rather the contrary). Under the densely instrumented United States, slab thickening is confirmed more directly because the dense USArray deployment resolves a shallow slab of the expected lithospheric thickness (<100 km), but shows it widening several-fold from ∼400 km down (Sigloch et al., 2008, their Figure 1).
The more voluminous such viscous slab pile-ups become, the more nearly vertically they would tend to descend, as massive Stokes sinkers under the pull of gravity (Morgan, 1965). The slabs under the Americas are among the most massive in the entire mantle. Hence, it seems reasonable to proceed with this vertical sinking hypothesis to explore what exactly it implies for the subduction history of South America, and whether observations exist that would permit to reject this simple hypothesis.
The east-dipping CA-TA slab above ∼900 km was almost certainly deposited beneath the Andean margin. We can estimate the vertical sinking rate from its dip. Figure 15 shows the dark-blue, 30-Ma margin contour aligning with the dark-blue slab contour at ∼300-km depth, and the aqua-colored, 70-Ma margin contour roughly aligning with the aqua-colored TA slab contour at ∼800-km depth. An optimized fit is achieved between the 75 Ma contour and the 800-km slab. This makes the eastward dip of the Andean slab consistent with vertical sinking at v sink =(800 − 300)km/(75−30)m.y. = 11.1 mm/a. Note that this estimate requires slab geometries and a plate reconstruction in a lower-mantle frame, but no geological constraints such as ages of arc volcanism. There is an obvious reading uncertainty; a plate reconstruction uncertainty; a spatial (imaging) uncertainty due to blurriness where the slab edge is located; and an uncertainty about the exact shape of the paleo-margin (due to subduction erosion, shortening, extension, and coastwise terrane translation). The latter uncertainty is often the largest; for the similar case of the North American Cordillera, an error propagation analysis by Sigloch and Mihalynuk (2013) yielded a timing uncertainty of ∼7 m.y., and a sinking rate uncertainty of 1-2 mm/a (relative error of ∼10%-20%) for times since 80 Ma. A slab sinking rate of 11 ± 2 mm/a for CA/TA slab would be the same (within error) as values across the lower-mantle slab assemblage under North America, which ranged from 9 ± 1 mm/a to 12 ± 2 mm/a (Sigloch & Mihalynuk, 2013). Van der Meer et al. (2010) used the alternative method of associating slab tops and bottoms around the globe with arc volcanic dates, and obtained slab sinking rates of 12 ± 3 mm/a for the lower mantle (although subsequently van der Meer et al., 2018 proposed larger sinking variability or uncertainty).

Slab Constraints on a Jura-Cretaceous "Intra-Arc" Basin Offshore the Andean Margin
Under the vertical sinking hypothesis, eastward termination of TA slab at 57°W implies that the trench was never located farther east in absolute coordinates. Hence no trench should have lined the Andean margin for the times when hotspot frame reconstructions place the margin (shelf break) east of 57°W, which is before ∼80 Ma according to Seton et al. (2012) or before ∼75 Ma according to Matthews et al. (2016), cf. Figure 15 (±7 Ma, as discussed above). We denote this time of first longitudinal overlap between TA slab and the reconstructed SAM margin as T 0 . We note the good temporal correlation of T 0 with the onset of vigorous arc magmatism in southern Peru (Toquepala arc in the region of Arequipa) ∼86-82 Ma constructed on older South American basement (Ramos, 2018;Spikings et al., 2015), marking the onset of the modern (Central) Andean arc. Before ∼85 Ma, Cretaceous magmatic arc rocks are only found in accreted terranes of the northern Andes, which are "suspect," separated by major fault zones and ophiolitic slivers from older Triassic-Jurassic arc rocks flanking cratonic South America. Thus, land geology is consistent with the separation of the Early to mid-Cretaceous arc from stable South America by an oceanic basin of uncertain width.
For times before T 0 the reconstruction constrains the basin's width as the space between 57°W and the South American margin. This space widens back in time as the South Atlantic is reconstructed to its zero-width state within Gondwana. By ∼120 Ma, in Figure 15, the space between the TA slab and the dark red coastline is 1,500 km wide at its narrowest, around the latitudes of northern Peru (∼1,800 km to central Peru; Figure 16b). This width rises well above the longitudinal reconstruction uncertainties for SAM, which increase as absolute positioning constraints from the Indo-Atlantic hotspots peter out between 100 and 120 Ma (O'Neill et al., 2005;Seton et al., 2012).
The geometries of two additional slabs, located beneath today's Atlantic, independently support the basin width implied by the hotspot frame reconstruction. These slabs are shown in Figure 16a, in a slab stack rendering of the DETOX-P3 model, since they are too deep to be confidently resolved by DETOX-P1. The Eastern Brazil/Equatorial Atlantic slab (EBEA) strikes northeast at ∼1,600 km to 2,200-2,400 km depth. The Argentine Basin slab (AgB) strikes southeast at ∼1,300 and ∼2,000 km depth. Both slabs are flat-topped and have no obvious dip. They can be further inspected and compared in Figures S21-and S27, and resolution tests of Figure S11. Figure 16a reconstructs South America and Gondwana to 140 Ma (Matthews et al., 2016), with continental positions almost unchanged to 120 Ma ( Figure 15). The lithosphere of the inferred, >1,500-km wide offshore basin is shaded blue.
The EBEA slab strikes parallel and moderately inboard of the 140-Ma paleo-margin of northwestern SAM-a candidate slab for a continental arc. This is reassuring in light of a first-order geologic constraint: before the modern Andean arc, confidently interpreted arc magmas last intruded the continental margin strata of SAM in Late Jurassic times, in Colombia (e.g., Ibaqué batholith, 183-149 Ma according to Spikings et al., 2015). Well-dated, calc-alkaline arc magmatism commenced in Colombia in Late Triassic (e.g., ∼209 Ma; Spikings et al., 2015). Arc termination at ∼149 Ma and the upward termination of EBEA slab at ∼1,600 km implies a sinking rate of 10.7 mm/a for EBEA. A lower slab termination at 2,200-2,400 km (variable along the strike of EBEA) and arc onset ∼209 Ma imply sinking rates of 10.5-11.5 mm/a. Within the error limits of ∼2 mm/a, these rates are indistinguishable from those obtained for the TA slab, the latter without reference to arc volcanism, painting a self-consistent picture.

Basin Closure by Intra-Oceanic Subduction and Establishing the Modern Andean Arc
No slab is observed beneath or aligned with more southerly segments of the SAM margin in Figure 16. The AgB slab strikes north-northwest, ∼60° to the paleo-Chilean margin, and terminates flush against it, implying intra-oceanic subduction. The relative depths of the slabs suggest that, at roughly uniform sinking rates, EBEA slab (>1,600-km deep) would have gone extinct before AgB slab (>1,300 km), and that TA slab (extending to ∼1,200-km depth in the south, ∼1500 km in the north) started depositing around the time that AgB slab shut down. Honoring these geometric constraints, the blue-shaded lithosphere must have subducted westward into the AgB and TA slabs, accommodating South America's westward motion as the South Atlantic opened. This basin lithosphere could not have subducted eastward beneath the Andean margin because no candidate slabs are imaged beneath the basin-the area across which a marginal Andean trench would have swept. Westward slab pull originating at the AgB and TA trenches would have favored an extensional regime on the continental margin. Indeed, the Peruvian margin was shaped by ex-MOHAMMADZAHERI ET AL.
10.1029/2020JB020704 28 of 35 tension and rifting during Jurassic and Early Cretaceous times, and it is difficult to observationall separate purely rift-related magmatism from arc-related. No arc was built between 145 and 110 Ma, and not all magmatic rocks of the Jurassic or Late Triassic age are clearly subduction-related (Spikings et al., 2015). A sinistral transpressional margin may have prevailed in the Jurassic, forming pull-apart basins; for example, the Ilo batholith in Peru (173-152 Ma) was emplaced into a slumping and extending margin (Boekhout et al., 2012). According to Ramos (2018) "the magmatic arc products were juvenile poorly evolved in intra-arc basins up to 9 km (6 miles) thick." Extension-related batholiths in Equador are dated 185-145 Ma (Spikings et al., 2015. For the Central Andes, Ramos (2018) summarizes the literature as indicating that subduction "developed under an extensional regime from Early Jurassic to Early Cretaceous times." Absence in Figure 16 of margin-aligned slab under the 140-Ma margin of Peru, Ecuador, and Chile margin, as well as west of it beneath the blue-shaded basin, suggests that subduction contributions to the Jurassic magmatic margin were minimal or absent. In Chile, Cretaceous volcanic rocks are dominantly bimodal and primitive (Muñoz et al., 2018) with isotopic geochemistry very different from "Andean type" magmatism and fitting to a model mixing curve with Pacific MORB and Jurassic plutonic end-members (Morata & Aguirre, 2003).
Time T 0 ∼75-85 Ma, when the reconstructed margin reaches southern TA slab, is the geophysically predicted time of basin closure, trench override, subduction flip, and conversion from intra-oceanic to continental arc ( Figure 16b). The first-order geologic observable of the Late Cretaceous is the Ayabacas giant submarine collapse ∼90 Ma (Ramos, 2018). This marks a change to collisional orogenesis in southern Peru, 5-10 m.y. later than in northern to central Peru (Ramos, 2018), and ∼185 km of shortening (Anderson et al., 2017). These observations fit the geophysically predicted scenario of the juvenile, westward-dipping AgB and TA arcs shedding their products into the blue-colored basin, while debris scraped off its subducting lithosphere piled up a subduction complex along the AgB and TA trenches, which was inverted as the basin approached closure. The AgB/TA arcs would have accreted or subcreted upon closure, and arc polarity was forced to flip eastward beneath the continental margin, to enable the continued westward drift of South America. This is marked by the onset of the modern Andean continental arc ∼85 Ma (85-45 Ma Toquepala arc in South Peru and North Chile; Ramos, 2018;Spikings et al., 2015).
Our interpretation of the TA slab topology, showing a transition at ∼900 km depth, from vertical to sloping slab deposition, is consistent with South America having collided with the intra-oceanic TA arc starting around 90 Ma. The timing of this collision is also consistent with both sedimentary and structural evidence in the Central Andes showing a transition from post-rift subsidence (during Late Jurassic-Early Cretaceous) to Late Cretaceous shortening. A north-to-south decrease in shortening (20°S-30°S, e.g., Horton, 2018) coincides spatially with the southern extent of the proposed collider, the TA intra-oceanic arc. Based on the slab extents, the collision of the TA arc would have been limited to north of ∼35°S, providing an answer to the long-standing question of why the Central Andes followed a distinctly different evolution path from that of the Southern Andes. DeCelles et al. (2015) suggested that the majority of Triassic and Jurassic arc rocks in the Western Cordillera of Chile could have originated west of South America, and were subsequently accreted. In the case of the TA arc, this accretion began ∼90 Ma, and relicts of the TA arc have since been translated northward, laying bare the cratonic margin at the latitude of southern Peru. Classical Andean-type continental arc magmatism is a product of eastward subduction beneath the continental margin since TA collision.
In the mantle, the polarity switch of lithospheric deposition into the TA trench is recorded by the changing character of the TA slab: below 800-1,000 km depth, the TA slab is vertical (evident by comparing Figures 14 and 15: the ochre and red levels are represented in Figures 14a and 14b but remain hidden behind the salmon layer of the slab stack in Figure 15). A vertical slab is indicative of a stationary trench (necessarily intra-oceanic, given that the SAM margin was moving). Above 900 km, the TA slab dips eastward, as expected from vertical sinking under a trench dragged westward by the continent (Figure 16). Hence the observation of a deep vertical part of TA slab suggests an intra-oceanic origin for this slab, and geometric constraints (slab free areas west of TA) require such an origin in the vertical sinking scenario. Associating TA slab inflection at 900 km with the Ayabacas collapse (∼90 Ma) or with modern arc onset (∼85 Ma) implies an averaged sinking rate of 10.0 mm/a to 10.6 mm/a for the deep, vertical part of TA slab. MOHAMMADZAHERI ET AL.

Summary: Complexity Resided at the Surface, Not in the Mantle
In summary, the hypothesis of vertical slab sinking and hotspot fixity infers a >1,500 km wide oceanic basin between the Andean margin and the volcanic arc in Jurassic and Early Cretaceous times. The basin's geophysically predicted evolution is in agreement with first-order geologic constraints, which represent completely independent observations. Hence, we perceive no need or motivation to invoke more complex hypotheses that involve non-vertical slab sinking (more than several hundred kilometers relative to the lower mantle) or reference frames that do not honor hotspot fixity (van der Meer et al., 2010). Across the lower-mantle slab assemblage, we consistently found sinking rates of 10-11 mm/a, with and without invoking knowledge of the ages of associated arcs.
Geometric constraints require that the "intra-arc" basin closed by westward (oceanward) subduction. Hence, the closed "Peru Basin" was not a backarc basin (formed east of an eastward-dipping subduction zone), an inference that differs from earlier interpretations (e.g., Braz et al., 2018;Kennan & Pindell, 2009). Establishing the polarity of intra-oceanic arcs from land-geologic constraints is notoriously difficult, especially if a polarity flip occurred, with younger continental arc overprinting the older intra-oceanic arc and suture. We acknowledge the remaining ambiguities in the geologic record, which have led to diverging interpretations. Specifically, a confident demonstration of arc volcanism on the Andean margin in Early Cretaceous times (as opposed to remelting of arc crust during extension) could seriously challenge our hypothesis.
A closely analogous controversy exists in the North American Cordillera on whether a dozen Jura-Cretaceous "collapsed basins", located inboard of accreted terranes, represent former backarc basins of eastward Farallon subduction (Monger, 2014;Pavlis et al., 2019), or a mature, westward subducted ocean basin, the Mezcalera Ocean (Sigloch & Mihalynuk, 2013. Labeled in Figure 16, this inferred Mezcalera Ocean represents the northward continuation of the Andean intra-arc basin constrained here. Hence further investigation of geologic constraints in the Andes will have implications for the understanding of mantle and of tectonic processes across the American Cordilleras. The South American situation is more favorable than for North America because the Peru basin sutured during the more recent times when hotspots still yield good paleo-positioning constraints. Arc terranes that accreted to the northern Andean margin in Cretaceous times must be associated with the deep, intra-oceanic Western Cordillera slab (WC), and with the deepest parts of the (northern) TA slab. Further work is needed and ongoing toward reconciling these shifting slab geometries below ∼1,000 km with the complex sequence of terrane accretions in the western Cordillera of Colombia and Ecuador (e.g., Braz et al., 2018;Kennan & Pindell, 2009;Kerr et al., 2002;Reynaud et el., 1999;Spikings et al., 2015).
The hypothesis that a hemisphere-spanning, intra-oceanic archipelago pulled in the Americas by westward subduction during Jura-Cretaceous times was put forward by Moores (1970) a half a century ago, motivated by the quest to explain ophiolite distributions in the southwest United States. Moores (1998) extended the evidence to Central and South America, an effort revived by Hildebrand and Whalen (2014). The root cause of such a vast and mantle-anchored arc complex has been attributed to a meridional girdle of degree-2 mantle downflow that persisted for the last 200 m.y. or more (Spencer et al., 2019).

Conclusion
The DETOX-P1 teleseismic P wave tomography has brought the elongated geometries of subducted slabs under South America into sharper focus, especially in the lower mantle down to ∼1,800 km depth. Mantle structures above ∼900 km depth correlate well with the modern (Early Cenozoic to present) geology of the Andes. Confirming previous studies, we resolve an eastward-dipping slab that conforms to the curvature of the Andean margin, i.e., has carried a faithful record of its paleo-trench geometry to the depths of the mantle. By contrast, slab geometries below 1,000 km are decorrelated from the shape of the continent and from its reconstructed paleo-position.
We explored the hypothesis that these deeper slabs reflect paleo-trench configurations equally faithfully as shallower slabs. More specifically, we asked whether existing observational evidence that might permit to MOHAMMADZAHERI ET AL.

10.1029/2020JB020704
30 of 35 reject the simple hypothesis that slabs sank vertically in the lower mantle, after deforming and widening viscously in the transition zone. This hypothesis predicts the geometry and temporal evolution of an "intra-arc" basin (Peru Basin) offshore the Jura-Cretaceous Andean margin. The Peru Basin must have closed by westward, intra-oceanic subduction, followed by a forced polarity flip to eastward, modern-style Andean subduction ∼85 Ma. Much evidence for the inferred basin has been reported in the geologic literature but its interpretations have been overprinted by expectations of keeping a trench close to the Andean margin at all times. The vertical slab sinking hypothesis is not compatible with this interpretation, which should therefore be scrutinized. The vertical slab sinking hypothesis does appear to accommodate all first-order constraints from the geologic record. Hence, we currently perceive no reason for rejecting it in favor of hypotheses that postulate more complex slab sinking behaviors or non-stationary hotspots. Our reinterpretation of subduction geometries bears close resemblance to those of the North American Cordillera.

Data Availability Statement
All the data used in this study are freely retrievable at all international data centers that implement the "web service" protocols of the International Federation of Digital Seismograph Networks. Specifically, we used the obspyDMT software of Hosseini and Sigloch (2017) to retrieve from the data centers listed at this URL: https://docs.obspy.org/packages/obspy.clients.fdsn.html. See Supporting Information S2 for details of the networks located in South and Central America. Our tomography model DETOX-P1, and the models compared to in this paper, are freely available through the "SubMachine" tomography web portal at: https://www.earth.ox.ac.uk/∼smachine/cgi/index.php?page=tomo_depth (where the user can select models, depths, and other parameters). The DETOX models can also be downloaded from the open-access repository Zenodo (https://doi.org/10.5281/zenodo.3993276).