Non-steady-state subduction and trench-parallel flow induced by overriding plate structure

Earth and Planetary Science Letters 401 (2014) 227–235

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

Non-steady-state subduction and trench-parallel flow induced by overriding plate structure Juan Rodríguez-González a,∗ , Magali I. Billen b , Ana M. Negredo a a b

Department of Geophysics and Meteorology, Universidad Complutense de Madrid, Instituto de Geociencias (UCM, CSIC), Faculty of Physics, 28040 Madrid, Spain Department of Earth and Planetary Sciences, University of California–Davis, Davis, CA 95616, USA

a r t i c l e

i n f o

Article history: Received 23 January 2014 Received in revised form 10 May 2014 Accepted 8 June 2014 Available online xxxx Editor: Y. Ricard Keywords: subduction dynamics slab dip trench-parallel flow seismic anisotropy thermo-mechanical numerical modeling Middle and South America

a b s t r a c t The direction of plate tectonic motion and the direction of mantle flow, as inferred from observations of seismic anisotropy measurements, show a good global correlation far from subduction zones. However, this correlation is poor near subduction zones, where below the slab seismic anisotropy is aligned parallel to the trench and above the slab has a complex pattern, which has not been fully explained. Here we present time-dependent three-dimensional (3D) fully-dynamic simulations of subduction to study the effect of overriding plate structure on the evolution of slab geometry and induced mantle flow. We find that along-strike variation in thermal thickness of the overriding plate causes increased hydrodynamic suction and shallower slab dip beneath the colder portion of the overriding plate; the variation in slab geometry drives strong trench-parallel flow beneath the slab and a complex flow pattern above the slab. This new mechanism for driving trench-parallel flow provides a good explanation for seismic anisotropy observations from the Middle and South America subduction zones, where both slab dip and overriding plate thermal state are strongly variable and correlated, and thus may be an important mechanism in other subduction zones. The location and strength of trench-parallel flow vary with the time-dependent evolution of the slab, suggesting that the global variability in seismic anisotropy observations in subduction zones is in part due to the non-steady-state behavior of these systems. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Olivine grains subject to shear deformation through the dislocation creep mechanism have a non-Newtonian rheology (Hirth and Kohlstedt, 2003) (i.e., a viscosity that depends on the strain rate through a power-law) and develop a lattice-preferred orientation (LPO) (e.g., Karato et al., 2008 and references therein). Shear waves traveling through the mantle can be split into two orthogonally polarized waves with different travel times due to the olivine LPO. Under most conditions in the mantle, the fast seismic direction is aligned with the mantle flow (or shear) direction (Karato et al., 2008). Observations of seismic anisotropy suggest a good global correlation between plate motion and mantle flow far from subduction zones (Becker et al., 2003; Conrad et al., 2007), but near subduction zones the inferred mantle flow does not always align with the plate motion (Long and Becker, 2010; Long and Silver, 2008) (Fig. 1). Observations of seismic anisotropy in subduction zones are derived from measurements of shear wave splitting (time delay and

*

Corresponding author. E-mail address: juan.rodriguez@fis.ucm.es (J. Rodríguez-González).

http://dx.doi.org/10.1016/j.epsl.2014.06.013 0012-821X/© 2014 Elsevier B.V. All rights reserved.

fast direction) from different types of shear waves (e.g., SKS, ScS, local S), each of which have different paths through the subduction zone depending on the source-receiver locations. SKS, ScS and local S waves are only sensitive to anisotropy on the path from the core to the receiver or local source to receiver (Long and Becker, 2010). However, one may also use source-side splitting measurements to determine sub-slab anisotropy for stations known to have no splitting for local phases (Eakin and Long, 2013). In general, each set of measurements can be modeled to better determine the orientation, thickness and strength of multiple anisotropic layers, however large data sets at each station are needed to determine multi-layered structure, while smaller data sets are often averaged and may miss complex behavior or give incorrect results (Eakin and Long, 2013). Seismic anisotropy below the slab is usually trench-parallel in the central region of the slab (Anderson et al., 2004; Polet et al., 2000; Rabble et al., 2011; Russo and Silver, 1994) and trench perpendicular near the slab edges (Abt et al., 2010; Hoernle et al., 2008; Soto et al., 2009). In contrast, shear-wave splitting patterns above the slab, in the mantle wedge, are substantially more complicated, often showing abrupt transitions between trench-parallel

228

J. Rodríguez-González et al. / Earth and Planetary Science Letters 401 (2014) 227–235

Fig. 1. Observations of the Middle America and South America subduction systems. Colors indicate regions with increased depth of the lithosphere–asthenosphere boundary (LAB; solid patterns) and surface heat flow (hashed patterns); double arrows indicate the fast polarization directions from shear wave splitting below (red arrows) and above the slab (black arrows) with the associated delay time (in seconds) between orthogonally polarized waves (Long and Becker, 2010; Long and Silver, 2008, 2009) red lines are slab contours (Hayes et al., 2012) every 40 km. Areas of thick lithosphere in South America correspond to the Amazonian (north) and Rio de la Plata (south) Cratons, and in Middle America correspond to the Paleozoicage continental margin. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and trench-perpendicular fast directions and sharp changes in LPO intensity (Long and Silver, 2008; Rabble et al., 2011). Observations of seismic anisotropy in the Middle and South America subduction zones illustrate these general observations of trench-parallel anisotropy below the slab and multi-directional anisotropy in the mantle wedge (Fig. 1). More recently, studies in these regions have also shown that seismic anisotropy beneath the flat portion of the Peruvian slab is aligned trench perpendicular near the slab (Eakin and Long, 2013). Similarly, the subslab anisotropy beneath the Nazca slab in southern South America includes some regions of trench-perpendicular and oblique anisotropy below the Chilean flat-slab, in addition to the trenchparallel alignment beneath the steeper slab portion further south (MacDougall et al., 2012). There has been long-lived subduction in both Middle America and South America since at least the Early Cretaceous (140 Ma; Lawyer et al., 2003). However plate tectonic reconstructions along both subduction zones indicate a time-varying subduction history, in particular for the past 60 Myr (Sdrolias and Müller, 2006). In Middle America, the Cocos plate is currently subducting at a rate of 5 cm/yr, having decreased from convergence rates of up to

15 cm/yr at 10 to 20 Myr ago (Sdrolias and Müller, 2006). Similarly, in South America the convergence rate has been steadily decreasing from a high of 12 to 14 cm/yr about 15 Myr ago to present-day convergence rates of only 5 cm/yr (Sdrolias and Müller, 2006). In addition, in Middle America there is a possible episode of slab break-off at 4 to 10 Ma evidenced by gaps in the tomographic image of the slab and surface uplift (Rogers et al., 2002). Notably, in both of these regions there is a strong correlation between lithospheric thickness (thermal or elastic) in the overriding plate and slab dip variations (Manea et al., 2012; O’Driscoll et al., 2009; Pérez-Gussinyé et al., 2008). The alongstrike change in overriding plate thermal thickness is summarized in Fig. 1 and is caused by the juxtaposition of old cratons adjacent to younger mobile belts (Artemieva, 2006). Artemieva (2006) calculated lithospheric thickness from heat flow measurements and xenolith geotherms and, for Middle America, obtained a transition from values between 100 and 125 km in the Yucatan Peninsula to values between 50 and 100 km in the Chortis Block. Lithospheric thermal thickness estimates obtained modeling the Bouguer anomaly in South America (Tassara and Echaurren, 2012; Tassara et al., 2006) show a transition from thickness greater than 130 km in the region of the Peruvian flat slab to 60 km further south. Surface heat flow data (Blackwell and Richards, 2004; Hamzaa and Muñoz, 1996) further illustrate significant thermal thickness variations in Middle and South America. Previous models have shown that overriding plate thermal state influences slab dip (Rodríguez-González et al., 2012; Schellart and Moresi, 2013) and that variations in slab dip can cause trenchparallel flow above the slab (Kneller and van Keken, 2007). Therefore, these results suggest a causal link between overriding plate structure, slab geometry and mantle flow in subduction zones. Trench-parallel flow below the slab and toroidal flow around its edges are commonly related to roll-back of the slab in laboratory and numerical models of subduction (Faccenda and Capitanio, 2012; Funiciello et al., 2006; Kincaid and Griffiths, 2003; Schellart and Moresi, 2013), but roll-back does not explain trenchparallel flow below the central region of the slab or variable flow patterns in the mantle wedge. Some 3D models of subduction show that trench-parallel flow in the mantle wedge can be induced by along-strike changes in the geometry of the slab (dip, curvature and length) (Jadamec and Billen, 2010; Kneller and van Keken, 2007), or by a subducting plate with laterally varying density (Capitanio and Faccenda, 2012). These models provide valuable knowledge of factors governing induced mantle flow, but they are limited in some aspects (e.g., use of Newtonian rheology, kinematic boundary conditions, no overriding plate, or are only solved instantaneously) and/or do not address the origins of the variations in the slab geometry that are responsible for the induced mantle flow. In this study we implement generic 3D time-dependent thermomechanical numerical models of buoyancy-driven subduction to study the influence of a non-uniform overriding plate thermal thickness on slab geometry and how this in turn modifies the induced mantle flow. We compare the resulting mantle flow with the flow pattern inferred from observations of shear wave splitting in the abovementioned subduction zones. 2. Methodology We use the finite-element code CitcomS (Tan et al., 2006; Zhong et al., 2000) to solve for the time-dependent solid-state flow of the lithosphere and mantle. The equations of conservation of mass, momentum and energy are solved for a 3D incompressible fluid. Time-dependent subduction is driven by thermal density anomalies prescribed by the initial thermal structure.

J. Rodríguez-González et al. / Earth and Planetary Science Letters 401 (2014) 227–235

Fig. 2. Schematic view of full 3D model domain. Geometry and thermal state of the subducting, lateral and overriding plates (red color for warm portion and blue for cold portion); plate boundary shear zone (green line); location of the trench (solid triangles) and location of the cross sections in Figs. 3, 5 and 6 (dashed grey lines). The depth of the 1000 ◦ C isotherm is z1000 = 99 km under the cold portion and z1000 = 58 km under the warm portion, with each portion extending ∼1300 km along strike. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2.1. Model set up The model domain spans 45◦ by 35◦ by 1500 km in latitude, longitude and depth, respectively (Fig. 2), and is discretized by a finite element grid with 480 by 488 by 160 nodes (37,478,400 elements total). Mesh resolution ranges from 0.04◦ to 0.1◦ in the longitudinal direction, 5 to 15 km in depth, and is uniform in the latitudinal direction with an element size of 0.1◦ . The highest resolution coincides with the plate boundaries and the slab. The modeled domain includes lateral plates to allow the mantle to flow around the edges of the slab. The subducting plate is surrounded by a plate boundary shear zone (PBSZ) modeled as a narrow low viscosity zone (Jadamec et al., 2013) that extends to a depth of 110 km. To the eastern boundary of the subducting plate the PBSZ (subduction hinge) dips 30◦ and decouples the overriding and subducting plates to allow for subduction. To the north and south the PBSZ is vertical and decouples the subducting plate from the lateral plates to allow for transform motion along these boundaries. To the west, the PBSZ is vertical and decouples the subducting plate from the western boundary of the mesh to simulate an oceanic spreading center. The reduced viscosity in the PBSZ accounts for the different weakening mechanisms taking place at the plate boundary (e.g., shear heating, faulting, fracture or fluid release). There are no direct rock mechanics constraints on the effective viscosity and width of this shear zone; therefore, we have tested different values to obtain plate-like motion and realistic plate velocities. Specifically, the amplitude of plate velocities in the models ranges between 3 and 14 cm/yr. These rates are controlled by the plate boundary shear zone and yield strength of the slab, more so than the viscosity of the mantle. Two- and three-dimensional tests were carried out to determine the necessary box dimensions in order to minimize boundary condition effects on the flow within the subduction zone. All the boundaries have a free-slip condition, except the overriding and lateral plates, which are fixed, the top boundary temperature is set at 0 ◦ C, the bottom boundary is set at 1400 ◦ C and the lateral boundaries are thermally insulated. 2.2. Initial temperature distribution In the models presented here subduction is dynamically-driven, that is, the models start with a partially subducted proto-slab that is negatively buoyant and provides a pull force that drives subduction. To obtain the temperature distribution of this proto-

229

slab we perform a simulation where subduction is kinematicallydriven by applying a convergence velocity of 5 cm/yr on the top of the subducting plate (Model 0). The initial thermal structure for the simulation of kinematically-driven subduction is built according to a half-space cooling model (HSCM) (Parker and Oldenburg, 1973). The age of the subducting plate varies with the distance from the oceanic spreading center at the western boundary. Given the extent in longitude of the subducting plate (25◦ ) and the convergence velocity, the age of the subducting plate varies from 0 Myr at the western boundary to 50 Myr at the subduction hinge. A uniform age of 50 Myr is used for overriding and lateral plates, for which the depth of the 1000 ◦ C isotherm, the temperature at which olivine behavior changes from quasi-rigid to viscous (Parsons and McKenzie, 1978) is z1000 = 58 km. Therefore this depth represents the thickness of the mechanical lithosphere, whereas the thickness of 50 Myr old thermal lithosphere is about 92 km (given by 2.52 times the diffusion depth, following Turcotte and Schubert, 1982) and the surface heat flow is about 60 mW/m2 , in agreement with observed values for unperturbed lithosphere in South America (Hamzaa and Muñoz, 1996). The kinematic simulation is run for 9 Myr, at which time the slab has reached a depth of ∼350 km and therefore has enough slab pull to drive self-sustained subduction. The proto-slab has a uniform dip of ∼75◦ measured along the slab surface from 120 to 220 km depth. The final temperature structure from this simulation is used as the initial condition for the self-driven subduction models, where a free-slip condition is applied on top of the subducting plate. In models with non-uniform overriding plates we replace the northern portion of the overriding plate with colder lithosphere (z1000 = 99 km). 2.3. Rheology We use a temperature, pressure and strain-dependent composite (dislocation and diffusion creep mechanisms) viscosity for the upper mantle, whereas we assume a linear viscosity (diffusion) for the lower mantle (Jadamec and Billen, 2010). The parameters used for the viscosity law (Hirth and Kohlstedt, 2003) are listed in Suppl. Table 1. With these parameters the viscosity predicted by each mechanisms (diffusion and dislocation creep), is equal to 1020 Pa s at 250 km depth for a strain rate of 10−15 s−1 . The choice of grain size in the lower mantle gives a tenfold increase in viscosity between the upper and lower mantle in agreement with constraints from the geoid (Hager, 1984). We set maximum and minimum cutoff viscosities of 1024 Pa s and 1019 Pa s respectively. In order to account for plastic deformation occurring at low temperatures and pressures, the effective viscosity is limited by a depth-dependent yield stress given by σ y = 15z + 0.1 (σ y in MPa and depth z in km) for σ y < σy max and σ y = σy max for σ y > σy max (see Suppl. Table 1 and Burkett and Billen, 2009 for details). The viscosity in the PBSZ is smoothly blended into the background viscosity, using a similar method to that in previous studies (Jadamec and Billen, 2010). The minimum allowed viscosity in the weak zone is 1021 Pa s. The viscosity in the weak zone is overwritten if the computed composite viscosity is lower. Within the multi-resolution finite-element mesh, the PBSZ spans a specified minimum number of elements (usually 15 elements, or ∼65 km) rather than a specific width, with the minimum viscosity occurring within the middle 2–3 elements. This allows us to constrain the viscosity change across the elements, which has been shown to be important for convergence in models with large viscosity variations (Jadamec et al., 2012; Moresi and Gurnis, 1996). We found good convergence by limiting the viscosity change between two adjacent nodes to a factor of three.

230

J. Rodríguez-González et al. / Earth and Planetary Science Letters 401 (2014) 227–235

Fig. 3. Results of viscosity structure and slab dip variations for Model B (with a non-uniform overriding plate) after 4.4 Myr. a) and b), cross-sections taken at 7.5◦ N (AA ) and 7.5◦ S (BB ) showing the viscosity structure (colors), temperature (contours every 200 ◦ C) and vertical projection of the velocity field (arrows). c) Time evolution of the slab dip. d) Slab dip variations along the trench at five different time steps. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3. Results To illustrate the effects of thermal heterogeneities in the overriding plate on slab geometry and induced mantle flow, two models are presented: Model A with a uniform overriding plate and Model B with a non-uniform overriding plate that includes a 12.5◦ wide portion of colder lithosphere (Figs. 2 and 3), with a HSCM age of 150 Myr (z1000 = 99 km). The proto-slab is identical in both models and has initially uniform dip of ∼75◦ . Characteristics and resulting along-trench slab dip variations of eight representative 3D simulations are listed in Suppl. Table 2. 3.1. Spatial variation of slab-induced mantle flow Model A (Figs. 4–6; Suppl. Figs. S1–S4; Videos 1, 2) has a nearly uniform slab dip along strike throughout the simulation. During the initial stage of the simulation the induced mantle flow in Model A shows the typical pattern of retreating slabs (Funiciello et al., 2006; Schellart and Moresi, 2013). Two toroidal cells appear on the edges of the slab (Fig. 4a). This return flow generates a locally intense trench-parallel flow extending about 1000 km away from the edges of the slab. Beneath the slab, the trench-parallel component of the velocity, v TP is very small, especially below the central portion of the slab. The trench-parallel flow is even weaker in the mantle wedge, especially in the central region of the slab. Therefore, this component only dominates over the radial and longitudinal components (velocity ratio v TP / v > 0.5) in small regions below and above the slab edges (Figs. 4a, 5a). In contrast, Model B (Figs. 3–6; Suppl. Figs. S2–S4; Videos 1, 2) evolves with a non-uniform dip. The higher viscosity in the mantle wedge beneath the cold portion of the overriding plate in Model B causes greater hydrodynamic suction than beneath the warm portion. This differential suction causes the initially uniform slab to

Fig. 4. Slab-induced mantle flow for a) Model A with a uniform overriding plate and b) Model B with a non-uniform overriding plate. Map view at 300 km depth after 4.4 Myr of evolution. Trench-parallel component of the velocity (colors: red and blue indicate flow directed northward and southward, respectively) and horizontal projection of the velocity field (arrows) are shown. Color-scale is truncated at ±2 cm/yr to mask fast but localized velocities (up to 3 cm/yr) in the mantle wedge. The location of the slab is indicated by the 1000 ◦ C isotherm (gray line). Regions in which the trench-parallel component dominates are outlined by green contours (v TP / v > 0.5), while orange contours outline regions where mantle flow is almost entirely trench-parallel (v TP / v > 0.9). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

J. Rodríguez-González et al. / Earth and Planetary Science Letters 401 (2014) 227–235

231

Fig. 5. Vertical trench-parallel sections for a) Model A and b) Model B taken at 29◦ E (CC ; top) and 22◦ E (DD ; bottom). Color arrows and contours are the same as in Fig. 4.

Fig. 6. Vertical trench-perpendicular sections for a) Model A and b) Model B taken at 2.5◦ N (EE ; top) and 2.5◦ S (FF ; bottom). Color arrows and contours are the same as in Fig. 4.

deform as it sinks. After 4.4 Myr, when the slab deformation is at a maximum, slab dip beneath the cold portion of the overriding plate diminishes from the initial value of 75◦ to 60◦ and slightly increases beneath the warm portion to 80◦ (Fig. 3). When the slab reaches the lower mantle, at about 5 Myr, the overall slab dip decreases (reaching 55◦ in the cold region and 60◦ in the warm region), due to the partial support of the slab by the more viscous lower mantle. Between 7 Myr and the end of the simulation at 20 Myr the slab dip remains approximately constant in both regions, but the slab dip variations along the trench are still significant.

The induced mantle flow is much more complex in Model B, due to the variable geometry of the slab (Jadamec and Billen, 2010; Kneller and van Keken, 2007, 2008) induced by the along-strike variations in the thermal state of the overriding plate (Figs. 4b, 5b). The flow beneath the slab resembles a cork-screw shape aligned along the strike of the slab: trench-parallel components couple to trench-perpendicular components of flow (Figs. 4b, 6b). After 4.4 Myr there is a uniform and intense trench-parallel component of flow below the slab. The region of high trench-parallel flow (v TP > 1 cm/yr) is more extensive than in Model A, extending 2000 km in the latitudinal direction (reaching the central region

232

J. Rodríguez-González et al. / Earth and Planetary Science Letters 401 (2014) 227–235

Fig. 7. Time evolution of the trench-parallel flow for Model A (top) and Model B (bottom). Trench-parallel component of the velocity computed at latitudes 2.5◦ N (a and c) and 2.5◦ S (b and d) at points indicated by green and orange crosses in Fig. 6, corresponding to positions below (green lines) and above (orange lines) the slab. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of the slab), 800 km in the longitudinal direction and 500 km in the radial direction. In addition, in Model B there are wide regions in which the trench-parallel component of the velocity dominates over the other components (v TP / v > 0.5). The induced mantle flow in the mantle wedge in Model B also shows a complex pattern, with sharp changes in intensity and direction in agreement with observations of seismic anisotropy (Long and Becker, 2010; Long and Silver, 2008). First, variations in the thermal state of the overriding plate induce small-scale vertical convective cells, which lead to two narrow regions of high trenchparallel flow (v TP > 3 cm/yr) at the edges of the colder region (Figs. 4b, 5b). Beneath the colder portion of the overriding plate, there is a wide region of trench-parallel flow extending into the central portion of the slab (v TP > 0.5 cm/yr): this flow is generated by variations in the slab geometry. This region of trench-parallel flow extends horizontally up to 500 km from the slab surface. Close to the slab the vertical component of the velocity dominates over the trench-parallel flow due to the influence of the sinking slab, which drags material down with it as it subducts. On the contrary, further away, where the influence of slab is lower, the trench-parallel component of the flow dominates (v TP / v > 0.5; Fig. 6). 3.2. Temporal variation in slab-induced mantle flow In both models, the slab geometry varies significantly during the first 7 Myr, which causes the mantle flow to vary with time. After this stage the slab dip in the upper mantle, in both models, remains approximately constant and the induced mantle flow stabilizes (Figs. 3 and S1). By computing the trench-parallel velocity at two points, one above the slab and another one below it (see location in Fig. 6a), we consistently find that the trench-parallel flow is more significant during the first 7 Myr of the evolution in Model B (Fig. 7). In Model A the trench-parallel velocity below and above the slab is consistently less than 0.5 cm/yr in the central region (Fig. 7a, b). In this model there are three different stages of evo-

lution of the flow. During the first 5 Myr the slab is steepening into the upper mantle, which induces return flow around the slab edges. Once the slab has reached 670 km, at about 5 Myr, the subduction velocity decreases due to the increased viscous resistance of the lower mantle. Between 5 and 7 Myr the slab dip decreases because it is easier for the slab to migrate laterally to accommodate the decrease in sinking rates between the lower viscosity upper mantle and the higher (Newtonian) viscosity lower mantle (Billen and Hirth, 2007). The return flow is inverted during this second stage in response to a decreasing slab dip. From 7 Myr until the end of the simulation the slab dip remains stationary and there is no toroidal flow around the slab edges. In Model B the overall mantle flow is more vigorous throughout the simulation and also shows three stages of evolution (Fig. 7c, d). For the first 5 Myr the trench-parallel flow increases rapidly to values between 1 and 2 cm/yr, and is directed towards the cold overriding plate region (COPR) below the slab, and directed towards the warm overriding plate region (WOPR) above it. Between 5 and 7 Myr, the southern portion of the slab rolls forward rapidly due to penetration into the lower mantle, causing the inversion of the flow pattern, which is now directed towards the COPR above the slab and reaches values higher than 3 cm/yr beneath the COPR. After 7 Myr the slab dip varies less with time and the induced mantle flow decreases, although it is still more significant than in Model A due to the contorted geometry of the slab. The spatial pattern of flow induced in the mantle does change significantly with time as the slab shape changes. First, as the part of the slab that is already deformed (by interaction with the overriding plate at shallow depths) sinks deeper into the mantle, trench-parallel flow is induced at greater depths. Second, because the variations in slab geometry in Model B decrease with time the trench-parallel velocity also eventually decreases. However, the induced mantle flow is always more complex in Model B than in Model A, showing wide regions of trench-parallel flow (reaching v TP ∼ 3 cm/yr) with a depth extent of a few hundred kilometers and larger regions of dominant trench-parallel flow (v TP / v > 0.5) (Figs. 4–6, S2–S4). In Model A, the trench-parallel flow is restricted

J. Rodríguez-González et al. / Earth and Planetary Science Letters 401 (2014) 227–235

to regions near the edges of the slab and has lower values (v TP < 2 cm/yr) throughout the simulation. 4. Discussion The results of the numerical simulations are in agreement with mantle flow inferred from seismic anisotropy (Long and Silver, 2008, 2009) reproducing the inferred trench-parallel flow below the slab and the complex flow pattern in the mantle wedge, as illustrated by the observations in the Middle and South America subduction zones. In these two subduction zones, regions of shallow slab dip correspond to regions with a cooler (thicker) overriding plate (Fig. 1). Globally variations in overriding plate thickness result from the juxtaposition of lithosphere of different ages or different thermo-tectonic histories. Therefore, contrary to the premise that slab dip primarily responds to variations of subducting plate buoyancy or radial variations in mantle viscosity, overriding plate structure is an additional factor that may have a dominant effect on slab evolution in some subduction zones. More generally we demonstrate that the spatial and temporal evolution of the slab geometry in three dimensions, regardless of the underlying mechanism, is responsible for driving strong trench parallel flow. This result is consistent with several recent studies showing that mantle flow around the slab can be decoupled from the flow driven by the dominant sinking forces of the slab in the presence of strong along-strike variations in either slab geometry (Jadamec and Billen, 2012; Kneller and van Keken, 2008) or slab buoyancy (Capitanio and Faccenda, 2012), or strong roll-back of the slab (Faccenda and Capitanio, 2012). The models presented here take these previous results a step further by linking the spatial and temporal changes in slab shape, and providing a mechanism that can drive the non-steady-state subduction required to produce strong trench-parallel flow. Other mechanisms may also lead to strong spatial and temporal variation in slab shape, including along-strike variation in subducting plate age, subduction of buoyant crust, interaction with other nearby slabs, having multiple overriding plates with different relative plate motions, interaction with upwellings from the lower mantle (Gurnis et al., 2000) or relict slabs in the transition zone and upper part of the lower mantle (Hall and Spakman, 2002). In order to isolate the effect of overriding plate heterogeneities without further complications we imposed a fixed plate boundary. This limitation is less restrictive in wide subduction zones (e.g., Middle and South America subduction zones), where observed trench motion is small (Schellart et al., 2007). However, further work is needed to evaluate the effects of overriding plate thermal variations on trench migration (e.g., Capitanio and Faccenda, 2012; Capitanio et al., 2011; Yamato et al., 2009). In addition, by fixing the trench plate boundary we have prevented slab roll-back, which could help to sustain and accentuate shallowing of the slab under the colder overriding plate (Cızkova et al., 2002). In future studies, removal of this and other simplifications, such as omitting the positive compositional buoyancy of the craton and the ridge-push induced motion of the overriding plate, would accentuate shallowing of the slab; the combination of these effects may be responsible for the flat (<20◦ ) dip of the Peruvian slab. The composite rheology used in these models is key factor for developing the trench-parallel flow regions both above and below the slab. In models with Newtonian viscosity the mantle is strongly coupled to the adjacent sinking slab leading to the classic corner flow geometry near the center of the slab. In models with non-Newtonian viscosity, localized lower viscosities occur in high strain-rate regions allowing the flow to decouple from the slab and respond to the along-strike pressure gradients created by the variations in slab dip (Jadamec and Billen, 2012). In the models presented here, the regions with significant trench-parallel

233

flow have gradients in flow velocity from a fast core to a slower edge, creating significant trench-parallel shear strain, and are deforming through the dislocation creep mechanism (non-Newtonian rheology), which likely leads to the development of LPO (Karato et al., 2008). While the strong trench-parallel flow is a transient feature lasting a few million years in these models, for several reasons strong LPO is still likely to develop and persist as trenchparallel flow is present in Model B for the entire model run time. First, an estimate of how well LPO will align with a flow field can be determined by comparing the rate at which the LPO fabric rotates toward the infinite strain axis (Ω ISA ), to the rate of rotation of the infinite strain axis within a flow field, Ω flow (Conrad et al., 2007; Kaminski and Ribe, 2002). Conrad et al. (2007) found that if Π = Ωflow /ΩISA is less than 0.5, then the ISA is a good indication of the orientation of LPO. Ω ISA is proportional to the strain-rate, while Ω flow depends on how the angle θ between the flow and ISA direction change both in time (∂θ/∂ t) at a particular place in the flow and space (u × dθ/dx) as the material moves through the flow field. For the models presented here, the regions of strong trenchparallel flow, with velocities of 5 to 10 cm/yr, are deforming by dislocation creep with strain-rates of 10−13 to greater than 10−11 s−1 (Figs. S5–S6). At a strain-rate of 10−13 s−1 , as long as the flow at a particular location changes by less than 5◦ /Myr and flow orientation changes by less than 1◦ /km along the flow path, then Π < 0.5, and the ISA, or the shear direction, is a good indicator of the LPO direction. In regions where the strain-rate reaches 10−12 s−1 , the LPO re-aligns so quickly that even changes in the flow field of over 30◦ /km and 30◦ /Myr result in Π < 0.5. Therefore, the trench-parallel induced mantle flow below the slab in Model B is expected to result in more easily detected and larger regions of trench-parallel seismic anisotropy. Second, two simplifications in our model design likely limit the development of trench-parallel flow: a fixed trench and a minimum viscosity cut-off. Fixing the trench was necessary for computational reasons but, as discussed above, a migrating trench may roll back adding to the along-strike variations in slab dip and driving stronger trench-parallel flow (Faccenda and Capitanio, 2012). Using a minimum viscosity cut-off accelerates the solver convergence, however allowing for even lower mantle viscosities would lead to larger strain-rates and more decoupling of mantle and slab flow (Jadamec and Billen, 2010). In addition to the above implications, the non-steady-state evolution of the slab calls for the interpretation of present-day observations to take into account how the subduction system, and in particular how the slab shape, is changing in time. First, interpretation of seismic anisotropy patterns will depend on how the slab shape (and induced mantle flow) is changing along strike, not just the present day slab dip (e.g., a presently flat slab could be steepening). Second, the stress-state in the overriding and subducting plates will depend on whether the slab is locally flattening or steepening, and could lead to along strike variation in trench motion and back arc spreading. This creates a potential for positive feedback, for example, between regions with cold cratons to drive slab shallowing, which in turn drives compression and thickening of the overriding plate, which then drives further shallowing of the slab. The central portion of the South America subduction zone may be an example of a location where this type of positive feedback has led to the shallow Nazca slab beneath the Andes (Allmendinger and Gubbles, 1996; Allmendinger et al., 1997; Arriagada et al., 2008; Isacks, 1988). Finally, the location and petrology of volcanism within the overriding plate (Behn et al., 2007; Grove et al., 2012; Honda and Mitsunobu, 2003) will also depend on the location of both over-

234

J. Rodríguez-González et al. / Earth and Planetary Science Letters 401 (2014) 227–235

riding plate thickness variations and along-strike variations in slab shape. 5. Conclusions Strong along-trench variations of slab dip and overriding plate thermal and mechanical structure are commonly found in subduction zones. Type localities for these variations are the South America and Middle America subduction zones. Seismic anisotropy observations in these zones suggest there is trench parallel flow beneath the slab and a complex pattern of flow in the mantle wedge. While these and similar observations elsewhere have been challenging to explain using steady-state models of subduction, the additional observations of a correlation between slab dip and overriding plate thermal structure suggest a causal link between the flow in the mantle, the overriding plate and slab dynamics. Using 3D, time-dependent numerical simulations we have shown that significant trench-parallel flow is generated by time-dependent along-trench changes in slab dip driven by along-strike variations in overriding plate thermal structure. We demonstrate that these variations produce highly contorted slab geometries, and therefore represent a plausible alternative to buoyancy heterogeneities in slabs to explain along-trench slab dip variations. In contrast to models with a homogeneous overriding plate, which predict trench parallel flow restricted to slab edges, models accounting for thermal variations in the overriding plate predict extensive regions of dominant trench parallel flow below the slab, extending more than 2000 km in the slab-parallel direction and reaching the central region of the slab. In the mantle wedge, a complex flow with wide regions of trench-parallel flow emerges from both the spatial variation in slab dip and overriding plate thermal structure. This mechanism for generating spatially and temporally variable flow in subduction zones provides new insight for interpreting seismic anisotropy observations at other subduction zones. In order to gain better insights into these and other subduction-related observations, we further emphasize that the lithospheric structure of the overriding plate and the time evolution of the slab geometry should be taken into account in analysis and modeling studies of subduction zones. Acknowledgements We would like to thank two anonymous reviewers for thorough and constructive reviews that helped to improve the manuscript. This work was funded by the Spanish Plan Nacional del MCINN project CGL2012-37222. This is a contribution of the ConsoliderIngenio 2010 team CSD2006-00041 (TOPO-IBERIA). J. RodríguezGonzález acknowledges the support from Spanish Ministry of Education for a Ph.D. grant and additional funding for a 4-month stage at UC-Davis. M.I. Billen acknowledges support from NSF grants 6877321 and 0748818. Appendix A. Supplementary material Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2014.06.013. References Abt, D.L., Fischer, K.M., Abers, G.A., Protti, M., González, V., Strauch, W., 2010. Constraints on upper mantle anisotropy surrounding the Cocos slab from SK(K)S splitting. J. Geophys. Res. 115, B06316. Allmendinger, R.W., Gubbles, T., 1996. Pure and simple shear plateau uplift, Altiplano–Puna, Argentina and Bolivia. Tectonophysics 259, 1–13. Allmendinger, R.W., Jordan, T.E., Kay, S.M., Isacks, B.L., 1997. The evolution of the Altiplano–Puna plateau of the Central Andes. Annu. Rev. Earth Planet. Sci. 25, 139–174.

Anderson, M.L., Zandt, G., Triep, E., Fouch, M., Beck, S., 2004. Anisotropy and mantle flow in the Chile–Argentina subduction zone from shear wave splitting analysis. J. Geophys. Res. 31, L23608. Arriagada, V., Roperch, P., Mpodozis, C., Cobbod, P.R., 2008. Paleogene building of the Bolivian orocline: tectonic restoration of the central Andes in 2-D map view. Tectonics 27, TC6014. Artemieva, I.M., 2006. Global 1◦ × 1◦ thermal model TC1 for the continental lithosphere: implications for lithosphere secular evolution. Tectonophysics 416, 245–277. Becker, T., Kellogg, J.B., Ekstrom, G., O’Connell, R.J., 2003. Comparison of azimuthal seismic anisotropy from surface waves and finite strain from global mantlecirculation models. Geophys. J. Int. 155, 696–714. Behn, M.D., Hirth, G., Kelemen, B., 2007. Trench-parallel anisotropy produced by foundering of arc lower crust. Science 317, 108–111. Billen, M.I., Hirth, G., 2007. Rheologic controls on slab dynamics. Geochem. Geophys. Geosyst. 8, Q08012. Blackwell, D., Richards, M., 2004. Geothermal Map of North America. American Assoc. Petroleum Geologist. Burkett, E., Billen, M.I., 2009. Dynamics and implications of slab detachment due to ridge-trench collision. J. Geophys. Res. 114. Capitanio, F.A., Faccenda, M., 2012. Complex mantle flow around heterogeneous subducting oceanic plates. Earth Planet. Sci. Lett. 353–354, 29–37. Capitanio, F.A., Faccenna, C., Zlotnik, S., Stegman, D.R., 2011. Subduction dynamics and the origin of Andean orogeny and the Bolivian orocline. Nature 480, 83–86. Cızkova, H., van Hunen, J., van den Berg, A.P., Vlaar, N.J., 2002. The influence of rheological weakening and yield stress on the interaction of slabs with the 670 km discontinuity. Earth Planet. Sci. Lett. 199, 447–457. Conrad, C.P., Behn, M.D., Silver, P.G., 2007. Global mantle flow and the development of seismic anisotropy: differences between the oceanic and continental upper mantle. J. Geophys. Res. 112. Eakin, C.M., Long, M.D., 2013. Complex anisotropy beneath the Peruvian flat slab from frequency-dependent, multiple-phase shear wave splitting analysis. J. Geophys. Res. 118, 4794–4813. Faccenda, M., Capitanio, F.A., 2012. Development of mantle seismic anisotropy during subduction-induced 3-D flow. Geophys. Res. Lett. 39, L11305. Funiciello, F., Moroni, M., Piromallo, C., Faccenna, C., Cenedese, A., Bui, H., 2006. Mapping mantle flow during retreating subduction: laboratory models analyzed by feature tracking. J. Geophys. Res. 111, B03402. Grove, T.L., Till, C.B., Krawczynski, M.J., 2012. The role of H2 O in subduction zone magmatism. Annu. Rev. Earth Planet. Sci. 40, 413–439. Gurnis, M., Ritsema, J., van Heijst, H.-J., Zhong, S., 2000. Tonga slab deformation: the influence of a lower mantle upwelling on a slab in a young subduction zone. Geophys. Res. Lett. 27, 2373–2376. Hager, B.H., 1984. Subducted slabs and the geoid: constraints on mantle rheology and flow. J. Geophys. Res. 86, 6003–6015. Hall, R., Spakman, W., 2002. Subducted slabs beneath the eastern Indonesia–Tonga region: insights from tomography. Earth Planet. Sci. Lett. 201, 321–336. Hamzaa, V.M., Muñoz, M., 1996. Heat flow map of South America. Geothermics 25, 599–646. Hayes, G.P., Wald, D.J., Johnson, R.L., 2012. Slab 1.0: a three-dimensional model of global subduction zone geometries. J. Geophys. Res. 117, B01302. Hirth, G., Kohlstedt, D.L., 2003. Rheology of the upper mantle and the mantle wedge: a view from the experimentalists. In: Eiler, J. (Ed.), Inside the Subduction Factory. AGU, Washington, DC, pp. 83–105. Hoernle, K., Abt, D.L., Fischer, K.M., Nichols, H., Hauff, F., Abers, G.A., van den Bogaard, P., Heydolph, K., Alvarado, G., Protti, M., Strauch, W., 2008. Arc-parallel flow in the mantle wedge beneath Costa Rica and Nicaragua. Nature 405, 1095–1097. Honda, S., Mitsunobu, S., 2003. Small-scale convection under the back-arc occurring in the low viscosity wedge. Earth Planet. Sci. Lett. 216, 703–715. Isacks, B.L., 1988. Uplift of the Central Andean Plateau and bending of the Bolivian Orocline. J. Geophys. Res. 93, 3211–3231. Jadamec, M.A., Billen, M.I., 2010. Reconciling surface plate motions with rapid threedimensional mantle flow around a slab edge. Nature 465, 338–342. Jadamec, M.A., Billen, M.I., 2012. The role of rheology and slab shape on rapid mantle flow: three-dimensional numerical models of the Alaska slab edge. J. Geophys. Res. 117, B02304. Jadamec, M.A., Billen, M.I., Kreylos, O., 2012. Three-dimensional simulations of geometrically complex subduction with large viscosity variations. In: Proceedings of the 1st Conference of the Extreme Science and Engineering Discovery Environment: Bridging from the eXtreme to the Campus and Beyond. ACM, Chicago, Illinois. Jadamec, M.A., Billen, M.I., Roeske, S.M., 2013. Three-dimensional numerical models of flat slab subduction and the Denali fault driving deformation in south–central Alaska. Earth Planet. Sci. Lett. 376, 29–42. Kaminski, É., Ribe, N.M., 2002. Timescales for the evolution of seismic anisotropy in mantle flow. Geochem. Geophys. Geosyst. 3, 1051. Karato, S., Jung, H., Katayama, I., Skemer, P., 2008. Geodynamic significance of seismic anisotropy of the upper mantle: new insights from laboratory studies. Annu. Rev. Earth Planet. Sci. 36.

J. Rodríguez-González et al. / Earth and Planetary Science Letters 401 (2014) 227–235

Kincaid, C., Griffiths, R.W., 2003. Laboratory models of the thermal evolution of the mantle during rollback subduction. Nature 425, 58–62. Kneller, E.A., van Keken, P.E., 2007. Trench-parallel flow and seismic anisotropy in the Mariana and Andean subduction systems. Nature 450, 1222–1224. Kneller, E.A., van Keken, P.E., 2008. Effect of three-dimensional slab geometry on deformation in the mantle wedge: implications for shear wave anisotropy. Geochem. Geophys. Geosyst. 9, Q01003. Lawyer, L.A., Dalziel, I.W.D., Gahagan, L.M., Martin, K.M., Campbell, D.A., 2003. The PLATES 2003 atlas of plate reconstructions (750 Ma to present day). PLATES progress report, p. 97. Long, M.D., Becker, T.W., 2010. Mantle dynamics and seismic anisotropy. Earth Planet. Sci. Lett. 297, 341–354. Long, M.D., Silver, P.G., 2008. The subduction zone flow field from seismic anisotropy: a global view. Science 319, 315–318. Long, M.D., Silver, P.G., 2009. Mantle flow in subduction systems: the subslab flow field and implications for mantle dynamics. J. Geophys. Res. 114, B10312. MacDougall, G.J., Fischer, K.M., Anderson, M.L., 2012. Seismic anisotropy above and below the subducting Nazca lithosphere in southern South America. J. Geophys. Res. 117, B12306. Manea, V.C., Pérez-Gussinyé, M., Manea, M., 2012. Chilean flat slab subduction controlled by overriding plate thickness and trench rollback. Geology 40, 35–38. Moresi, L.N., Gurnis, M., 1996. Constraints on the lateral strength of slabs from three-dimensional dynamic flow models. Earth Planet. Sci. Lett. 138, 15–28. O’Driscoll, L.J., Humphreys, E.D., Saucier, F., 2009. Subduction adjacent to deep continental roots: enhanced negative pressure in the mantle wedge, mountain building and continental motion. Earth Planet. Sci. Lett. 280, 61–70. Parker, R.L., Oldenburg, D.W., 1973. Thermal models of ocean ridges. Nat. Geosci. 242, 137–139. Parsons, B., McKenzie, D., 1978. Mantle convection and the thermal structure of the plates. J. Geophys. Res. 83, 4485–4496. Pérez-Gussinyé, M., Lowry, A.R., Phipps Morgan, J., Tassara, A., 2008. Effective elastic thickness variations along the Andean margin and their relationship to subduction geometry. Geochem. Geophys. Geosyst. 9, Q02003. Polet, J., Silver, P.G., Beck, S., Wallace, T., Zandt, G., Ruppert, S., Kind, R., Rudloff, A., 2000. Shear wave anisotropy beneath the Andes from the BANJO, SEDA and PISCO experiments. J. Geophys. Res. 105, 6287–6304.

235

Rabble, W., Koulakov, I., Dinc, A.N., Jakovlev, A., 2011. Arc-parallel shear deformation and escape flow in the mantle wedge of the Central America subduction zone: evidence from P wave anisotropy. Geochem. Geophys. Geosyst. 12, Q05S31. Rodríguez-González, J., Negredo, A.M., Billen, M.I., 2012. The role of the overriding plate thermal state on slab dip variability and on the occurrence of flat subduction. Geochem. Geophys. Geosyst. 13, Q01002. Rogers, R.D., Ka’rason, H., van der Hilst, R.D., 2002. Epeirogenic uplift above a detached slab in northern Central America. Geology 30, 1031–1034. Russo, R., Silver, P.G., 1994. Trench-parallel flow beneath the Nazca plate from seismic anisotropy. Science 263, 1105–1111. Schellart, W.P., Freeman, J., Stegman, D.R., Moresi, L., May, D., 2007. Evolution and diversity of subduction zones controlled by slab width. Nature 446, 308–311. Schellart, W.P., Moresi, L., 2013. A new driving mechanism for backarc extension and backarc shortening through slab sinking induced toroidal and poloidal mantle flow: results from dynamic subduction models with an overriding plate. J. Geophys. Res. 118, 1–28. Sdrolias, M., Müller, R.D., 2006. Controls on back-arc basin formation. Geochem. Geophys. Geosyst. 7, Q04016. Soto, G.L., Ni, J.F., Grand, S.P., Sandvol, E., Valenzuela, W., Speziale, M.G., Gómez González, J.M., Domínguez Reyes, T., 2009. Mantle flow in the Rivera–Cocos subduction zone. Geophys. J. Int. 179, 1004–1012. Tan, E., Choi, E., Thoutireddy, P., Gurnis, M., Aivazis, M., 2006. GeoFramework: coupling multiple models of mantle convection within a computational framework. Geochem. Geophys. Geosyst. 7, Q06001. Tassara, A., Echaurren, A., 2012. Anatomy of the Andean subduction zone: threedimensional density model upgraded and compared against global-scale models. Geophys. J. Int. 189, 161–168. Tassara, A., Götze, H.-J., Schmidt, S., Hackney, R., 2006. Three-dimensional density model of the Nazca plate and the Andean continental margin. J. Geophys. Res. 11, B09404. Turcotte, D.L., Schubert, G., 1982. Geodynamics: Applications of Continuum Physics to Geological Problems. Cambridge University Press, New York. Yamato, P., Husson, L., Braun, J., Loiselet, C., Thieulot, C., 2009. Influence of surrounding plates on 3D subduction dynamics. Geophys. Res. Lett. 36, L07303. Zhong, S., Zuber, M.T., Moresi, L.N., Gurnis, M., 2000. The role of temperaturedependent viscosity and surface plates in spherical shell models of mantle convection. J. Geophys. Res. 105, 11063–11082.