The deforming Nazca slab in the mantle transition zone and lower mantle: Constraints from teleseismic tomography on the deeply subducted slab between 6°S and 32°S

The Nazca–South America subduction zone is a long-lived plate boundary system that shows significant along-strike variations in subduction angle. Subduction of the Nazca plate under the central Andes of South America occurs at an angle of ∼30° down to depths of ∼300 km along much of the margin. There are two regions of flat or shallow subduction, one occurring under Peru, north of 14°S, and another located under central Chile and Argentina, at ∼30°S (Cahill and Isacks, 1992; Hayes et al., 2012). North of 14°S, the Nazca slab subducts at a shallow angle for several hundred kilometers inland (Cahill and Isacks, 1992; Hayes et al., 2012; Bishop et al., 2013) before resuming subduction at a very steep angle (∼70°) down to the bottom of the mantle transition zone (MTZ) (Scire et al., 2016). We combined data from 14 separate networks deployed in the central Andes at various times from 1994 to 2013 to image the subducting Nazca slab down to 895 km depth between 6°S and 32°S using finite-frequency teleseismic P-wave tomography (Fig. 1). Previous studies by the authors (Scire et al., 2015, 2016) used subsets of this data set to examine variations in upper mantle structure down to 660 km depth. By combining these data sets, we are able to increase the depth to which we can resolve velocity anomalies in the mantle due to the increased aperture of the array at the surface. In this paper we will focus on the geometry of the Nazca slab in the MTZ and lower mantle, providing new insights into how the subducting Nazca slab interacts with the high-viscosity lower mantle. The tomographic images presented here offer new constraints on the morphology of the subducted Nazca slab on a regional scale previously unmatched by most tomography studies. The along-strike variations in the Nazca–South America subduction zone and the regional scale of our tomographic images allow us to look at the relationship between slab angle and slab deformation in the MTZ and lower mantle in a single system. Our work offers new insights into slab-MTZ interactions and how along-strike variations in subduction zones influence the deformation that subducted slabs undergo in the MTZ and lower mantle.

Variations in the width of the Nazca slab in the MTZ have been observed using teleseismic tomography, indicating that some deformation of the Nazca slab is occurring in the MTZ, possibly related to resistance to continued subduction through the 660 km discontinuity into the lower mantle (Scire et al., 2015, 2016). Increases in slab thickness in the MTZ have also been observed in global tomography studies and have been interpreted as evidence for temporary stagnation of some subducting slabs in the MTZ, resulting in thickening and deformation of the slab in the MTZ before subduction continues into the lower mantle (e.g., Bijwaard et al., 1998; Li et al., 2008). The results from these studies are consistent with receiver-function studies (Liu et al., 2003) and modeling studies (Quinteros and Sobolev, 2013) that indicate a correlation between an abnormally thick MTZ and resistance to subduction of the Nazca plate into the lower mantle.

The imaged morphology of subducting slabs in the MTZ varies globally, with some slabs penetrating into the lower mantle with only minor broadening in the MTZ while other slabs flatten out and stagnate in the MTZ over broad regions (e.g., Bijwaard et al., 1998; Li et al., 2008; Fukao et al., 2009). In the Nazca–South America subduction zone, global tomographic studies have observed a diffuse fast anomaly in the lower mantle that corresponds to the presumed location of the subducting Nazca slab, indicating that while deformation of the Nazca slab may occur in the MTZ, the slab continues to subduct into the lower mantle and does not experience long-term stagnation in the MTZ (e.g., Bijwaard et al., 1998; Li et al., 2008; Zhao, 2004; Fukao et al., 2001, 2009; Zhao et al., 2013). Similarly, Engdahl et al. (1995) observed an amorphous fast anomaly in the lower mantle associated with the subducting Nazca slab using teleseismic and regional P phases to image the Nazca–South America subduction zone between 5°S and 25°S down to ∼1400 km depth. Using teleseismic travel time tomography, Rocha et al. (2011) imaged the Brazilian mantle on the eastern side of South America between 38°W and 58°W. At the very westernmost edge of their study area (57°–58°W), they imaged a fast anomaly from 700 km to 1300 km depth which they interpreted as being the Nazca slab in the lower mantle, although its geometry it difficult to resolve due to vertical smearing on the edge of their model space. South of 30°S, Pesicek et al. (2012) used a nested global-regional tomography model to image the Nazca subduction zone down to 1500 km depth. They imaged a detached segment of Nazca slab below ∼200 km depth to the south of 38°S, which they interpreted as being a result of a tear corresponding to the location of the Mocha fracture zone. North of the hypothesized tear, Pesicek et al. (2012) imaged the Nazca slab in the transition zone, but the slab passes out of the eastern edge of their model space before reaching the lower mantle.

The data set used in this study consists of 384 short-period and broadband seismometers deployed in the central Andes as part of 14 separate networks between 1994 and 2013 (Fig. 1). The networks include 13 temporary deployments with varying numbers of stations. The data set also includes data from permanent stations of the Plate Boundary Observatory network that were deployed before 2013. Two permanent stations from the Global Seismograph Network (stations LVC and NNA) and one permanent station from the Global Telemetered Seismograph Network (station LPAZ) are also included. Arrival times were picked for direct P phases for 402 earthquakes with magnitudes >5.0 between 30° and 90° away from the study region (Fig. 2A). Additional arrivals were picked for PKIKP phases for 144 earthquakes with a similar magnitude limit between 155° and 180° away from the study region. As in previous studies by the authors (Scire et al., 2015, 2016), arrivals were picked using the multi-channel cross correlation technique described by VanDecar and Crosson (1990) and modified by Pavlis and Vernon (2010). A total of 27,435 direct P arrivals were picked for the broadband stations in three frequency bands (0.04–0.16 Hz, 0.1–0.4 Hz, and 0.2–0.8 Hz) and for the short-period stations in one frequency band (0.5–1.5 Hz) due to the limited frequency content of data from the short-period stations. An additional 10,055 PKIKP arrivals were picked in the same three frequency bands for the broadband stations. Distribution of the picks in the three frequency bands is dominated by the 0.2–0.8 Hz band, with 44% of all picks in the 0.2–0.8 Hz band, 24% in the 0.1–0.4 Hz band, and 21% in the 0.04–0.16 Hz band. The 0.5–1.5 Hz frequency band contains the fewest picks, 11% of the total number of picks, due to the limited duration of the deployments of most of the short-period stations and the lack of PKIKP picks for these stations. The azimuthal distribution of incoming rays is uneven, with the majority of events occurring in either the Middle America subduction zone to the northwest or the South Sandwich subduction zone to the southeast (Fig. 2B).

Travel time residuals were calculated relative to the IASP91 reference model (Kennett and Engdahl, 1991) and then demeaned for each event to calculate the relative travel time residuals. Travel times were corrected for crustal variations using Moho depth estimates from receiver functions from Bishop et al. (2013) and Ryan et al. (2016) and crustal thickness estimates from Tassara and Echaurren (2012) in order to account for crustal heterogeneity (Schmandt and Humphreys, 2010). Crustal corrections were made using a homogenous layered velocity model corresponding to the average velocity of the central Andean crust (Zandt et al., 1996; Swenson et al., 2000; Dorbath and Masson, 2000). Crustal corrections for stations in the forearc were calculated using a faster velocity model based on local tomography studies (Graeber and Asch, 1999; Koulakov et al., 2006; Schurr et al., 2006) and ambient noise tomography that covers our study region (Ward et al., 2013).

This study uses a finite-frequency teleseismic tomography algorithm discussed in detail by Schmandt and Humphreys (2010) and used in previous studies by the authors (Scire et al., 2015, 2016). The finite-frequency approximation defines the sampled area around the geometrical ray path by the first Fresnel zone, whose width is dependent on both the frequency and the distance along the ray path (Hung et al., 2000; Dahlen et al., 2000). The sampling in each model layer is determined by the width of the Fresnel zone, which increases with distance from the source or receiver, and the differential sensitivity within the Fresnel zone, resulting in theoretical “banana-doughnut” sensitivity kernels (Dahlen et al., 2000).

The model space is parameterized into a series of nodes in a nonuniform grid that dilates with increasing depth and distance from the center of the model. In the shallowest layers of the model, nodes are spaced 35 km apart in the center of the model where sampling is densest, increasing to 56 km apart at the edges of the model. Vertical node spacing increases with depth, from 35 km at 60–95 km to 60 km at >715 km depth. The horizontal node spacing in the shallowest layers varies from 35 km to 56 km and also increases with depth such that the horizontal node spacing varies from 52 km to 83 km at a depth of 940 km. The uppermost model layer (60 km) and the lowermost model layer (940 km) are removed from any interpretation because these layers absorb effects from anomalies outside of the model space that are not accounted for in the crustal corrections or the global mantle velocity model respectively. In addition, station and event terms are incorporated in order to compensate for perturbations outside of the model space. Because a homogenous, layered velocity model is used to calculate the crustal corrections, the station terms absorb local perturbations in the velocity structure as well as any errors in the a priori crustal thickness. The event terms account for variations in the mean velocity structure under the subset of stations that record a specific event.

The LSQR algorithm of Paige and Saunders (1982) is used to invert the frequency-dependent relative travel time residuals to calculate velocity perturbations within the modeled volume by obtaining the minimum length model that satisfactorily explains the observed data. The inverse problem is regularized using norm and gradient damping as well as model smoothing to account for uncertainties in the location of the ray paths, details of which are discussed in Schmandt and Humphreys (2010). A tradeoff analysis was performed between the variance reduction and the Euclidean model norm to choose the overall damping and smoothing weights (Supplemental Material, Fig. S1¹). The chosen damping (D5) and smoothing (S4) parameters result in a variance reduction of 82.3%.

The sampling of the model space is determined by calculating normalized hit quality maps in which each node is assigned a hit quality value between 0 and 1 which is a function of both the number of rays sampling that node as well as the azimuthal distribution of those rays (Fig. 3; Biryol et al., 2011; Schmandt and Humphreys, 2010). This calculation relies on the idea that better sampling of a node is achieved with the intersection of rays from multiple azimuths. A node that has been sampled by multiple rays from all four of the geographical quadrants is assigned a hit quality of 1 while a node that is not sampled at all has a hit quality of 0. Comparisons of the hit quality maps for the combined data set used in the current study with the maps for the subsets of data used previously by the authors (Scire et al., 2015, 2016) show that combining the two smaller data sets has resulted in a substantial increase in hit quality, particularly at MTZ depths (Fig. 3). This indicates that our combined data should be better able to resolve velocity anomalies in the MTZ as well as in the lower mantle due to the increased sampling at these depths.

Synthetic resolution tests were performed to determine the ability of our model to resolve structures in the MTZ and lower mantle. A synthetic “checkerboard” input was created with alternating fast and slow velocity anomalies with amplitudes of +5% and –5% respectively defined in eight-node cubes (Fig. 4). The size and location of the synthetic input anomalies change with depth due to the dilation of the node spacing with depth and distance from the center of the model. In order to assess the extent of any vertical or lateral smearing, two zero-anomaly nodes separate the positive and negative synthetic “checkers”. The output from the checkerboard tests shows that lateral resolution is good in the center of the model space, with increased lateral smearing and decreased amplitude recovery observed toward the edges of the model space (Fig. 4). On the eastern edge of the model space in the deepest layers where we expect the Nazca slab to be located in the MTZ and lower mantle, amplitude recovery is decreased and lateral smearing is increased from the center of the model space, but ∼30% of the input anomaly is still being recovered. Lateral smearing at the eastern edge of the model space in the lower mantle depth slices is dominantly to the southeast due to large number of rays coming from the South Sandwich subduction zone (Fig. 2B). Neutral layers (505 km and 660 km) show evidence of vertical smearing with low-amplitude anomalies being erroneously resolved into neutral background layers.

Additional synthetic tests were performed to assess the ability of our inversion to resolve the subducting Nazca slab in the MTZ and lower mantle. The Slab1.0 global subduction zone model (Hayes et al., 2012) was used to create a synthetic slab anomaly that is ∼70–100 km thick with amplitudes of +5% to input into the model in order to test the ability of the inversion to resolve the Nazca slab (Fig. 5). Because the Slab1.0 model ends between 600 km and 660 km depth due to limited information about the structure of the Nazca slab in the lower mantle, the synthetic Nazca slab anomaly was taken from the 660 km depth slice and projected into the lower mantle depth slices (715–940 km). The recovered slab anomaly from 300 km to 900 km depth is continuous with amplitudes of ∼60% of the input amplitude. A slight decrease in amplitude recovery is noted with increasing depth, but the continuity of the recovered slab anomaly remains at all depths. Limited broadening of the recovered slab anomaly is observed, indicating that variations in the width of the slab anomaly in our results probably represent real variations in the thickness of the Nazca slab. As in the checkerboard tests (Fig. 4), some smearing of the slab anomaly is observed at the eastern edge of the model space in the lower mantle depth slices (e.g., 895 km) although the amplitude of the smeared slab anomaly decreases rapidly with increasing distance from the synthetic slab anomaly (Fig. 5).

We present our results in horizontal depth slices (Figs. 6 and 7) and vertical cross sections (Fig. 8). Because detailed discussions of the shallow mantle have previously been published for subsets of this data set (Scire et al., 2015, 2016), we focus our discussion on the morphology and deformation of the Nazca slab in the mantle transition zone (MTZ) and the lower mantle.

The Nazca slab is observed as a trench-parallel fast anomaly (∼1%–4%) in the MTZ that generally follows the Slab1.0 contours (Hayes et al., 2012; Fig. 6). The subducting slab in this region is mostly aseismic in the depth interval of 300–500 km, so these are the first detailed images of the slab in this depth range. South of ∼26°S there is a parallel slab-like anomaly beneath the Chile-Argentina border. Although this feature is within the “well-resolved” area of our study, it lies mostly south of the southern edge of our array and will be better resolved in a separate ongoing study incorporating stations located further south. We refrain from interpreting this feature until the study with better resolution of this feature is complete. The fast slab anomaly continues into the lower mantle as a fairly continuous anomaly at 715 km and 775 km depth (Fig. 7). Resolution of the Nazca slab is complicated by increased lateral smearing along the eastern edge of the model space as noted above and shown in Figures 4 and 5. The fast Nazca slab is still observable in the lower mantle, but lateral smearing to the southeast must be taken into account when interpreting the geometry of the slab at 835 km and 895 km depth.

Our observations of the Nazca slab in the lower mantle (Figs. 7 and 8) are consistent with global tomographic studies, many of which have imaged a diffuse fast anomaly in the lower mantle associated with the subducting Nazca slab (e.g., Bijwaard et al., 1998; Li et al., 2008; Zhao, 2004; Fukao et al., 2001, 2009; Zhao et al., 2013). However, we are able to offer higher-resolution images of the slab and along-strike variations in the slab geometry in the lower mantle. We find that the slab anomaly generally broadens in the lower mantle. While the slab anomaly is fairly continuous at 715 km and 775 km depth, segmentation of the slab anomaly is observed below 800 km depth (Fig. 7). This segmentation occurs primarily south of the bend in the slab anomaly associated with the Bolivian orocline. While some of the eastward smearing of the slab anomaly in the southern part of our study region is likely due to the ray path distribution as discussed previously, the along-strike segmentation of the slab does not appear to be an effect of the decreasing resolution with depth because similar segmentation is not observed in the synthetic slab recovery tests (Fig. 5). These recovery tests also showed that we have the resolution to distinguish between a slab that terminates at the base of the MTZ and one that continues into the lower mantle. North of the orocline, the slab anomaly is more coherent and exhibits higher amplitudes than in the southern part of the study region. In addition, the slab penetrates nearly vertically into the lower mantle north of the orocline, compared to a low angle of penetration south of the orocline (Fig. 8).

A large high-amplitude slow anomaly (–3% to –4%) is observed in the MTZ to the east of (above) the slab anomaly north of 18°S (anomaly A, Fig. 6). This anomaly appears to continue vertically into the upper mantle where we have good resolution east of the slab (anomaly A, Fig. 8, cross sections D1 and D2). North of the orocline, multiple smaller but vertically elongated low-velocity bodies are located west of (below) the slab in the upper mantle and were discussed in more detail in Scire et al. (2016). South of the orocline the most prominent low-velocity anomalies in the MTZ also occur west of (below) the slab. In particular, sub-slab low-velocity anomalies are observed immediately adjacent to the slab between ∼20°S and 26°S (Figs. 6 and 8). Previous work by the authors (Scire et al., 2015) also imaged these sub-slab anomalies and interpreted them as either local thermal anomalies or regions of hydrated material in the MTZ. We note that the slab recovery tests (Fig. 5) show a small-amplitude low-velocity “halo” artifact around the high-velocity slab, but the amplitude of this artifact is several times smaller than that of the slab anomaly whereas the low-velocity features described above are about the same amplitude but opposite sign as the slab anomaly.

In order to better illustrate the changing morphology of the subducted Nazca plate within the MTZ, we isolated just the high-velocity anomaly associated with the slab (Fig. 9). We also plot the slab width relative to the average width at 410 km at depths of 550 km and 660 km along strike from north to south in Figure 10. The slab width was measured by plotting the +1% velocity anomaly contour on the 410 km, 555 km, and 660 km depth slices. Measurements of the width of the slab anomaly defined by this contour were then taken perpendicular to the strike of the slab at every ∼200 km starting at 8°S and moving south along strike. We took the average of the slab width at 410 km depth and divided the measurements at 555 km and 660 km depth by this average to plot the slab width in the MTZ relative to the slab in the upper mantle. Along-strike changes in slab width are observed, particularly to the north of 16°S where the slab anomaly is as much ∼2× broader with a steeper dip than to the south (Fig. 10). As we observed in our earlier study (Scire et al., 2016), the slab anomaly becomes much wider with depth in the northern part of the study region. Near 14°S the slab anomaly width at 410 km is ∼200 km, and at 605 km the anomaly is ∼400 km wide. Slab recovery tests for this region (Fig. 5) show that we can recover ∼100-km-wide slab anomalies at these depths. The slab anomaly bends in the MTZ between 15°S and 17°S, paralleling the change of the strike of the trench at the surface associated with the Bolivian orocline. Between ∼19°S and 26°S, the slab anomaly morphology changes from a relatively thin and segmented appearance in the upper 50 km of the MTZ to a more continuous and thicker appearance at deeper depths. South of 26°S, the slab anomaly appears relatively similar at different depths (Fig. 9). The strongly segmented appearance of the slab at 660 km depth may be due to diminishing resolution at this depth, but we note that our slab recovery test (Fig. 5) does not show this effect.

The thickening of the Nazca slab in the MTZ that we observe in the northern part of our study region (Fig. 6) is similar to observations made in global tomography studies which interpreted the observed widening of the slab in the MTZ as a result of temporary stagnation as the slab resists subduction into the lower mantle (e.g., Bijwaard et al., 1998; Li et al., 2008; Zhao, 2004; Fukao et al., 2001, 2009; Zhao et al., 2013). Because we do not see any evidence for continued stagnation to the east, we conclude that the stagnation is only temporary, unlike for the Pacific slab under East Asia which remains stagnant in the MTZ forming a laterally extensive body for several hundred kilometers (e.g., Zhao, 2004; Fukao et al., 2001, 2009; Li et al., 2008; Zhao et al., 2013; Fukao and Obayashi, 2013). Several studies have suggested that folding or buckling of the slab occurs in the MTZ to accommodate a decreasing sinking velocity of the slab in the lower mantle due to the viscosity increase across the 660 km discontinuity (e.g., Ribe et al., 2007; Běhounková and Čížková, 2008; Lee and King, 2011; Garel et al., 2014). This mechanism is consistent with the observed thickening of the slab anomaly in the MTZ toward the northern end of our study region (Figs. 6 and 8). Little to no thickening of the slab anomaly is observed in the MTZ to the south of 16°S despite similar resolution across this region, indicating that deformation of the slab in the MTZ varies along strike (Figs. 9 and 10). This type of variation in the subducted slab deformation could be due to the interplay of the slab temperature, slab strength, viscosity of the slab, and the viscosity increase at the 660 km discontinuity (Stegman et al., 2010; Garel et al., 2014).

It seems unlikely that the viscosity contrast at the 660 km discontinuity changes over the relatively short distance between the differing slab morphologies that we observe. The morphology of the slab anomaly in the north falls into the vertical folding category of Garel et al. (2014) or folding mode of Stegman et al. (2010). Folded slab formation is enhanced by weak slabs (based on the ratio of the viscosity of the slab to the viscosity of the upper mantle) with little trench retreat, based on modeling by Stegman et al. (2010). Hence, one possibility is that the slab to the north is weaker than the slab south of 16°S. The age of the subducting slab is not very different between these two regions (Müller et al., 2008), but near the Nazca Ridge subduction in the north, the tomography shows a strong low-velocity anomaly beneath the Nazca Ridge (Figs. 8A and 8B) that may have added heat to the base of the slab resulting in a warmer and weaker slab. Another possibility is that in the northern segment the trench is not retreating as the slab sinks into the mantle, but rather the horizontal part of the slab is nearly overriding the deep portion of the slab that is anchored in the transition zone, producing the near-vertical dip of the slab to the north of the Bolivian orocline. The near-vertical penetration of the slab into a more viscous lower mantle is a geometry more conducive to slab thickening than a more gradual angle of slab penetration in the lower mantle according to Garel et al. (2014). The subducted slab further south has a more continuous dip that is consistent with trench retreat. Further modeling is needed to understand the along-strike changes in the Nazca slab subduction.

Focal mechanisms for events in the MTZ (Ekström et al., 2012) indicate that the Nazca slab is generally in a state of downdip compression at depths between 410 and 660 km (Figs. 6 and 9), consistent with a viscosity increase across the 660 km discontinuity that results in increased resistance to subduction into the lower mantle (e.g., Isacks and Molnar, 1971; Vassiliou, 1984). While the causes of deep seismicity are uncertain, several models have been proposed to explain the occurrence of seismicity in the MTZ, including transformational faulting associated with phase changes of olivine in cold subducting slabs in the MTZ (e.g., Wiens et al., 1993; Kirby et al., 1996; Guest et al., 2003, 2004) and thermal shear instabilities (e.g., Ogawa, 1987; Kanamori et al., 1998; Karato et al., 2001). The fault planes of the focal mechanisms change orientation along strike, following the strike of the imaged slab as it bends in the MTZ. Seismicity in the transition zone in our study region is mostly clustered into two distinct regions, with a cluster of events occurring north of 12°S where the hypocenters are concentrated within the area of broadening of the slab anomaly in the MTZ, and a second cluster occurring south of 26°S where the hypocenters are aligned near the top edge of the narrower slab anomaly. Work by Myhill (2013) indicates that such bands of denser deep seismicity could be associated with the localization of strain in the hinge zones of folds in the slab in the MTZ. If, as previously discussed, this zone of thickening of the slab in the MTZ is related to buckling or folding of the slab in response to increased viscosity across the upper-lower mantle boundary, then the distribution of the associated deep seismicity could similarly be affected by the deformation of the slab in the MTZ. The cluster of seismicity associated with the broader slab anomaly in the north could be explained by the slab thickening and folding. However, the seismicity south of 26°S does not correspond to a similar thickening of the slab anomaly, and its cause is uncertain.

An isolated region of very deep seismicity (>625 km depth) in the MTZ corresponds to the location of the 9 June 1994 M_w 8.3 deep Bolivia earthquake which occurred at ∼640 km depth at ∼14°S (International Seismological Centre, 2013; Fig. 6). This event was the largest deep-focus earthquake recorded until the more recent 2013 M_w 8.3 Sea of Okhotsk earthquake. Slip occurred along a subhorizontal plane (e.g., Kikuchi and Kanamori, 1994; Beck et al., 1995) perpendicular to the observed strike of the slab anomaly in our tomograms. Prior to the deep Bolivia earthquake, few deep-focus earthquakes were observed in this region. The focal mechanism of the 1994 Bolivia earthquake deviates from the more consistent pattern of downdip compression observed in the clusters of deep seismicity to the north and south (Fig. 6). Detailed analyses of the fault rupture pattern of this earthquake have indicated that the 1994 Bolivia earthquake was probably caused by a combination of factors (e.g., Ihmlé, 1998; Zhan et al., 2014). Zhan et al. (2014) concluded that the initial rupture occurred in the cold core of the slab due to transformational faulting as metastable olivine underwent a phase change (e.g., Kirby et al., 1996). A large sub-event during the first phase of rupture triggered shear melting, allowing the rupture to continue to propagate away from the core of the slab into the warmer slab material (Ogawa, 1987). The alignment of the east-west nodal plane with the strike of the slab anomaly in our results (Fig. 6) agrees with the idea that most of the slip was perpendicular to the strike of the slab as the rupture propagated away from the cold core of the slab into the warmer slab material, where the melting point is more easily reached and slip increased due to positive feedbacks during shear melting.

Our models demonstrate that the Nazca slab shows along-strike variations in morphology in the lower mantle (Figs. 7 and 8). The high-amplitude, continuous slab anomaly north of 16°S corresponds to the steeply dipping (∼70°) segment of the Nazca slab inboard of the Peruvian flat slab as well as the thickened slab in the MTZ (Scire et al., 2016). This corresponds to the region where Engdahl et al. (1995) observed ponding of high-velocity material immediately below the MTZ. Where the dip of the Nazca slab is more gradual when entering the MTZ (Fig. 8), less widening of the slab anomaly is observed. The slab anomaly, particularly where we observe less widening of it, also appears to begin flattening to the east after it enters the lower mantle (Fig. 8), although this could be influenced by lateral smearing at the edges of our model. Because several other global and regional tomography studies (e.g., Bijwaard et al., 1998; Fukao et al., 2001; Zhao, 2004; Li et al., 2008, Engdahl et al., 1995) imaged the diffuse fast anomaly in the lower mantle associated with the Nazca slab broadening to the east, however, it is possible that this apparent flattening is not solely a result of lateral smearing. Běhounková and Čížková (2008) concluded that buckling or folding of the slab related to the viscosity increase between the upper and lower mantle would explain the imaging of the slab as diffuse fast anomalies in the lower mantle rather than the narrower slab anomalies which are observed in the upper mantle. Rocha et al. (2011) imaged a fast anomaly that they interpreted as the subducting slab at 900 km depth at the edge of their study region under Brazil at 57°–58°W, which also indicates that at least some of the eastward flattening that we observe is genuine. The along-strike segmentation of the slab in the lower mantle south of the Bolivian orocline also appears to correspond to the more gradual dip angle of the slab as it enters the MTZ. In the north where the incoming dip angle is high and we observed extensive deformation in the MTZ, the slab anomaly remains coherent into the lower mantle, unlike further south. This potentially indicates that the slab dip angle and subsequent deformation in the MTZ are influencing the strength and cohesion of the slab in the lower mantle.

Hydration of the mantle transition zone by the transportation of hydrous phases to depth by subducting slabs has been suggested by a number of studies (e.g., Bercovici and Karato, 2003; Smyth and Jacobsen, 2006; Ohtani et al., 2004; Magni et al., 2014). Direct evidence for hydration of MTZ minerals has been observed recently in a sample of hydrated ringwoodite discovered as an inclusion in a diamond (Pearson et al., 2014). Outer rise normal faulting is thought to allow for hydration of the oceanic mantle lithosphere (e.g., Peacock, 2001; Garth and Rietbrock, 2014), although the degree of mantle lithospheric hydration is uncertain. The high-amplitude low-velocity anomaly observed in our tomograms in the MTZ to the north of 18°S (anomaly A, Fig. 6) could be associated with dehydration of the oceanic mantle lithosphere in the MTZ or just below the 660 km discontinuity in the lower mantle, leading to a region of locally elevated hydration in the MTZ. Alternatively, this low-velocity anomaly could possibly reflect dehydration melting related to the downward flow of hydrous minerals due to the penetration of the Nazca slab into the lower mantle. Schmandt et al. (2014) theorized that downward flow into the lower mantle would result in dehydration melting due to the low water storage capacity of lower mantle minerals (e.g., Bolfan-Casanova et al., 2000). Such a hydrous melt is likely to be slightly less dense than the top of the lower mantle and would therefore rise upwards back into the MTZ (Sakamaki et al., 2006), so the presence of this melt in the lower mantle is likely short lived. We suggest that this low-velocity anomaly could reflect the presence of a hydrous melt that has returned to the MTZ driven by positive buoyancy after dehydration melting has occurred in the lower mantle. In addition, dehydration melting has been suggested to occur when hydrous wadsleyite moves upwards across the 410 km discontinuity into the olivine stability field in the upper mantle (Bercovici and Karato, 2003). Both Schmerr and Garnero (2007) and Contenti et al. (2012) have observed a disrupted 410 km discontinuity under the central Andes, which both have attributed to compositional heterogeneity in the upper part of the MTZ in this area. Schmerr and Garnero (2007) suggested that their observation of a deepened 410 km discontinuity could be due to the presence of a hydrated “lens” of wadsleyite, while Contenti et al. (2012) suggested that the decreased reflectivity of the 410 km discontinuity could be due to either the effects of hydration of the zone directly above MTZ or the presence of a layer of neutrally buoyant melt. While some studies have suggested that this hydrous melt would be neutrally buoyant at the 410 km discontinuity (e.g., Sakamaki et al., 2006; Bercovici and Karato, 2003), Hirschmann (2006) argued that the density contrast between the hydrous melt and solid peridotites is such that the hydrous partial melts could potentially percolate up rather than pooling at 410 km depth. Continued upward flow of hydrous melts from the MTZ could explain the apparent continuation of the low-velocity anomaly into the upper mantle (anomaly A, Fig. 8, cross sections D1 and D2), although the effect of this hydrous melt on the upper mantle and potentially even the cratonic lithosphere under which it is rising is unknown.

By combining data from 14 separate networks deployed in the central Andes, we are able to image the Nazca–South America subduction zone to depths of ∼900 km between 6°S and 32°S. Our images of the Nazca slab in the MTZ and lower mantle provide new insights into the deformation of the slab and its interactions with the boundary between the upper and lower mantle. We observe significant thickening of the slab in the MTZ in the northern part of our study region (north of 16°S), which is due to temporary stagnation and possibly folding or buckling of the slab in the MTZ in response to the decreasing of the sinking velocity of the slab due to increasing viscosity in the lower mantle. In contrast, limited thickening of the slab is observed to the south where the slab appears to subduct into the lower mantle with little to no deformation in the MTZ. While the Nazca slab penetrates into the lower mantle through the entirety of our study area, significant along-strike variation in the amount of deformation both within and beneath the mantle transition zone is observed. The character of the slab anomaly in the MTZ and lower mantle changes from a broad, higher-amplitude, cohesive anomaly in the north to a more segmented and weaker slab anomaly in the south. This variation appears to at least partially correspond to changes in the dip of the slab where it enters the MTZ. We observe a low-velocity zone above the slab in the MTZ to the north of 18°S, which we attribute to the effects of hydration of the MTZ. Our images of the subducting Nazca slab provide new constraints on slab morphologies in the mid-mantle, which could help guide interpretations of modeling studies exploring the dynamic interactions between the slab and the MTZ. We have documented significant along-strike variability in the morphology and deformation of the Nazca slab in the MTZ and lower mantle, indicating that along-strike variations in the plate boundary have a significant effect on the deformation and fate of the slab at depth. This has implications not only for the long-lived South America subduction system which has had significant variations in subduction style both along strike and throughout its history, but also for other subduction systems that show variations in the plate boundary both spatially and temporally.

Support for this research was provided by the U.S. National Science Foundation (NSF) as part of the CAUGHT (NSF award EAR-0907880), PULSE (NSF award EAR-0943991), and MIMOSA (NSF award EAR-1415914) projects. ML acknowledges support from NSF award EAR-0943962. Additional support for this research was provided by ExxonMobil as part of the COSA project. The seismic instruments for the CAUGHT and PULSE arrays were provided by the Incorporated Research Institutions for Seismology (IRIS) through the PASSCAL Instrument Center at New Mexico Tech, Yale University, and the University of North Carolina, Chapel Hill. Data collected are or will be available through the IRIS Data Management Center. Additional data were also obtained from the IRIS Data Management Center. The facilities of the IRIS Consortium are supported by the NSF under Cooperative Agreement EAR-1063471, the NSF Office of Polar Programs, and the U.S. Department of Energy National Nuclear Security Administration. The authors would like to acknowledge the GEOFON Program at GFZ Potsdam as a source of much of the waveform data. The Global Seismographic Network (GSN) is a cooperative scientific facility operated jointly by IRIS, the U.S. Geological Survey, and the NSF. We thank Rob Clayton and Paul Davis for providing PeruSE data. The authors thank Brandon Schmandt at the University of New Mexico and Gene Humphreys at the University of Oregon as the original authors of the teleseismic tomography code used in this study. We would also like to thank two anonymous reviewers for their comments, which helped greatly improve this manuscript. Chevron-Texaco and Conoco-Phillips provided additional support for AS.