Subduction-transition zone interaction: A review

The mantle transition zone controls the material and heat exchange between upper and lower mantle. The upper-lower mantle boundary has long been realized to coincide with an endothermic phase transition (Ringwood, 1975; Ito and Yamada, 1982) and to be the likely site of a viscosity increase (Hager, 1984), both of which will hamper flow across the boundary. Mechanisms of Benioff earthquakes and the morphologies of upper-mantle slabs inferred from seismicity and seismic tomography provide direct evidence that the transition zone indeed hinders flow, because many slabs are deformed, and several flattened at the base of the upper mantle (e.g., Lay, 1994; Li et al., 2008). It has been debated (1) what causes some slabs to flatten while others do not; (2) on what time scales slabs are stagnant in the transition zone if they have flattened; (3) how slab-transition zone interaction is reflected in plate motions; and (4) how lower-mantle fast seismic anomalies can be correlated with past subduction.

Over the past 20 years since the last Subduction Top to Bottom volume and several reviews of subduction through the transition zone (Lay, 1994; Christensen, 2001; King, 2001; Billen, 2008; Fukao et al., 2009), information about the nature of the transition zone and our understanding of how slabs dynamically interact with it has increased substantially. In this review, we will show how, when evaluated together, this past work allows answering the questions above as well as narrowing the plausible properties of the transition zone and slabs. These insights can help the next generation of global-evolution models of mantle convection and can aid further unravelling of what factors may have been responsible for changes in plate motions, plate boundaries, and surface tectonics in the past.

In the first part of the paper (Section 2), we review observational constraints on slab morphology from seismic tomography and slab seismicity, and constraints on subduction dynamics from the plate motion record. Subsequently, in Section 3, we discuss what has been learned from modeling slab dynamics in the transition zone. Figure 1 summarizes the main forces that affect this interaction. Section 3.1 discusses the role of mantle resistance, i.e., the viscosity jump, phase changes, and a potential density jump between the upper and lower mantle; while Section 3.2 covers the role of slab properties, i.e., density and strength, and the mobility of the trench. Finally, in Section 4, we combine the observations with the modeling insights to answer the questions posed above.

The most detailed information on the structure of slabs as they interact with the transition zone has been obtained from seismic tomography. Before tomographic inversions, Benioff seismicity (Isacks and Molnar, 1971) and analyses of slab structure, for example using the residual sphere technique (Lay, 1994), already indicated that there are different styles of subduction-transition zone interaction. The first tomographic studies mapped the structure of slabs subducting below the western Pacific, along the Japan-Kurile and Izu-Bonin-Mariana trench (Zhou and Clayton, 1990; Van der Hilst et al., 1991; Fukao et al., 1992) and revealed that some slabs flatten at the base of the transition zone (Izu-Bonin and Japan–southern Kurile), while others appear to pierce straight through to larger depths (Marianas and northern Kurile).

Several generations of tomographic models later, the structure of most subducting slabs in the upper mantle and transition zone has been mapped. In Figure 2, we compiled the range of slab morphologies in the transition zone from a wide range of tomographic studies (references used are numbered and listed in the caption). Most of the studies we used consist of tomographic inversions of P-wave travel times (using direct as well as surface- or core-bounced phases); a few used S-wave travel times. The studies differ in their data selection and parametrization and regularization of the inversion. Nonetheless, they agree on many aspects of transition-zone slab morphology, which is what we try to capture by our compilation. In the few cases where studies significantly disagree on the interpretation of slab extension or continuity (e.g., below Tonga-Kermadec and North America) or resolution is very limited (e.g., below the Scotia plate), parts of the slab shapes are shaded light blue. To aid discussion of slab volumes and time scales below, the slabs are labeled with the approximate length of the flat segments perpendicular to the present-day trench (in white) or with the approximate depth to which lower-mantle slabs are continuous with those in the upper mantle (in black).

Stagnant slabs are prevalent under the western Pacific, including, besides those listed above, the Ryukyu, Aleutian, and Tonga slabs. Those that do penetrate into the lower mantle below the northwest Pacific are generally not continuous below ∼1000 km depth, although there are also significant volumes of fast anomalies deeper in the lower mantle below the region (Van der Hilst et al., 1997; Grand, 2002; Simmons et al., 2012). Lower-mantle anomalies below Tonga and Kermadec extend from ∼700 km down to at least 1700 km, but it is not clear whether they are connected to upper-mantle slabs (Van der Hilst, 1995; Schellart and Spakman, 2012).

The classic example of a deeply penetrating plate is the Cocos plate, subducting along the eastern side of the Pacific, below Central America (Grand et al., 1997). However, this is clear only along a few cross sections through a rather complex slab, which changes its strike from more or less west-east in the upper mantle below Mexico, parallel to the Trans-Mexican Volcanic Belt, to more or less north-south in the lower mantle below the Gulf of Mexico and the United States (Gorbatov and Fukao, 2005). Along other profiles below the Americas, significant thickening is observed in the transition-zone slab below Peru, where it penetrates the lower mantle down to only ∼1500 km depth, while slab flattening appears to occur in the transition zone below Cascadia, Mexico and southern South America (Sigloch et al., 2008; Simmons et al., 2012; Fukao and Obayashi, 2013). Slab fragments that lie in the transition zone below the United States are generally associated with previous (Laramide) subduction (Van der Lee and Nolet, 1997; Sigloch, 2011), and their connection to the prominent lower-mantle anomaly below eastern North America is debated (Bunge and Grand, 2000; Schmid et al., 2002; Ren et al., 2007; Sigloch and Mihalynuk, 2013).

Fukao and coworkers (Fukao et al., 2001; Fukao et al., 2009; Fukao and Obayashi, 2013) proposed that there is a significant amount of material trapped below 660 km and above 1000 km depth, in particular below Indonesia, Tonga-Kermadec, and South America. In tomographic models from other groups, there is little evidence of ponding above 1000 km, apart from the material below Borneo, Indonesia, which is found in many studies (Hafkenscheid et al., 2001; Replumaz et al., 2004; Simmons et al., 2012). The recent P-velocity model of Simmons et al. (2012) instead finds even more extensive amounts of material sequestered above 660 km, including below Kamchatka and Hellenic arc, where most previous studies have suggested the slabs penetrate. Some of the slabs (e.g., below South America) may appear flat when cross sections are taken obliquely through relatively complex slab geometries. For Kermadec, the interpretations of the lower-mantle anomalies range from flattened (Fukao and Obayashi, 2013) to just thickened (Schellart and Spakman, 2012). French and Romanowicz (2014) note that their recent global shear-speed model confirms the flattened slab material that Fukao et al. (2009) identified in the uppermost lower mantle. However, they only show a cross section for Kermadec, which could be interpreted as either flattening over ∼500 km below 660 km or thickening. Rather than being flat, the Indonesia deep flat slab has been suggested to represent a cluster of slab material subducted at different zones and not connected directly to the steep upper slab below Java (Hall and Spakman, 2015). Thus it appears that very few, if any, slabs actually flatten below 700 km depth.

A few previous studies (Loiselet et al., 2010; King et al., 2015) have noted a broad correlation between stagnant slabs and their estimated thermal structure (based on the age of the subducting plate and its subduction velocity). Indeed, stagnant slabs are common in the western Pacific, where very old lithosphere is currently subducting (90–145 m.y., Müller et al., 2008), while some deeply penetrating slab sections are found along the eastern rim of the Pacific, where the lithosphere currently subducting is relatively young (0–50 m.y., Müller et al., 2008). In more detail, there may be further complexity. For example, the penetrating segment of the Izu-Bonin-Marianas slab corresponds to the oldest lithosphere along this trench, while at the Tonga-Kermadec and the Japan-Kuriles-Kamchatka trenches, the youngest part of the lithosphere penetrates. However, all lithosphere currently subducting along these trenches has an age (>90 m.y.) that exceeds that above which plates appear to no longer thicken and increase their negative buoyancy (Stein and Stein, 1992; Davaille and Jaupart, 1994). Additionally, it has long been observed that there is no simple correlation between present-day slab dip and age at the trench (Lallemand et al., 2005; Sdrolias and Müller, 2006). Along each trench there has been a large variation in age of the subducting plate during the Cenozoic, which may be one reason that there are no simple correlations with slab geometry (e.g., Garel et al., 2014).

Also noted (in white text) on Figure 2 are the approximate lengths of the flat slabs in the transition zone. These range from several hundred km (below Calabria) to up to 2000 km (below North America and, if this is indeed a flat slab, below Borneo). At subduction rates of 5 cm/yr (i.e., around the global average; Lallemand et al., 2005; Sdrolias and Müller, 2006), and without significant stretching or thickening, 1000 km corresponds to ∼20 m.y. of subduction. Comparison of tomography with plate reconstructions has led to estimates of accumulation times, i.e., stagnation times, around such a number, e.g., for the ∼1000-km-long, flat Izu-Bonin slab, 20–30 m.y. (Van der Hilst and Seno, 1993; Miller et al., 2006; Sdrolias and Müller, 2006); for the ∼1600 km Honshu flat slab, more than 15–20 m.y. (Miller et al., 2006; Goes et al., 2008); for the Tonga flat segment of ∼1300 km, ∼30–40 m.y. (Van der Hilst, 1995; Sdrolias and Müller, 2006; Schellart and Spakman, 2012); and for the ∼500-km-long, Calabrian flat slab, 8–10 m.y. (Wortel and Spakman, 2000; Faccenna et al., 2001). Schmid et al. (2002) estimated that the Farallon material in the transition zone corresponds (due to buckling) to subduction since 50–60 m.y. ago. The deeper flat slab below Indonesia has been associated with opening of the South China Sea 30–40 m.y. ago (Replumaz et al., 2004). These numbers imply that material can indeed stagnate in the transition zone for a few tens of millions of years. However, there is no seismic evidence of older material in the upper mantle and transition zone. Because the time scales to diffuse the thermal anomaly of a cold slab and hence make it seismically invisible are on the order of hundreds of millions of years, material subducted in the early Cenozoic and Mesozoic has apparently sunk into the deeper mantle.

Focal mechanisms provided the earliest evidence that slabs encounter resistance as they reach the base of the upper mantle. Most slabs exhibit either downdip extension or compression, where globally, extension is most common above 300–350 km depth, while below this depth all slabs that exhibit deep seismicity are dominantly under compression (Isacks and Molnar, 1971; Apperson and Frohlich, 1987; Alpert et al., 2010). Richter (1979) showed that, in addition to increasing downdip compression in earthquake mechanisms, there is a bulge in seismic energy release above the base of the upper mantle, between 500 and 700 km depth, which can be reconciled with increased resistance at this depth. Initially, it was thought that these observations and the absence of seismicity below 700 km depth meant slabs stopped at this depth, until several studies established that the depth extent of seismicity is linked to how deep the slab core remains cooler than a potential temperature of 600–800 °C (Molnar et al., 1979; Wortel, 1982), and studies of seismic structure showed that many slabs extend beyond their Benioff zones (e.g., Van der Hilst et al., 1991; Fukao et al., 1992; Lay, 1994; Grand et al., 1997).

Vassiliou and Hager (1988) showed that earthquake stress directions could be matched with models of viscous sinking slabs encountering a viscosity increase of, at least, an order of magnitude at 660 km depth, the mechanism previously proposed by Isacks and Molnar (1971). More recent studies (Alpert et al., 2010; Bailey et al., 2012), using global-flow models driven by plate motions and slabs with viscosities 1–2 orders of magnitude higher than the background upper mantle that track the Benioff zones, can match the global pattern of downdip slab extension and compression with an upper-lower mantle viscosity jump of 30–100.

There is regional variation on the global pattern: while most slabs exhibit extension at intermediate depths (100–350 km), some, in particular in the western Pacific, are compressional throughout their depth range (including Japan, Kurile, Izu-Bonin, and Tonga; see Fig. 2) (Isacks and Molnar, 1971). A few studies have modeled the difference in stress patterns as the result of a net westward rotation of the lithosphere relative to the lower mantle, which enhances compression in westward-dipping slabs and extension in eastward-dipping ones (Alpert et al., 2010; Carminati and Petricca, 2010). However, another study was able to match the difference in Benioff stress patterns in the Tonga and South American slabs in a global mantle flow model as a result of differences in slab geometry (Alisic et al., 2010). Earlier studies (Isacks and Molnar, 1971; Vassiliou and Hager, 1988) had proposed instead that the slabs with compression throughout were more supported by the lower mantle than those with intermediate depth extension. The compilation in Figure 2 provides support for the earlier interpretation because (1) not all westward-dipping slabs are in compression throughout (e.g., Marianas and Kermadec), and not all eastward-dipping slabs exhibit intermediate depth extension (e.g., Calabria); and (2) the slabs with compression throughout tend to be those that flatten in the transition zone (with the exception of the Chilean slab). This includes slabs that are too warm to sustain deep seismicity, such as the Aleutian, Ryukyu, and Scotia flattened slabs, where the former two have previously been identified as outliers of the usual trend of downdip extension in slabs without Benioff seismicity below 350 km (Isacks and Molnar, 1971; Bailey et al., 2012).

Others have used seismicity to estimate the amount of deformation that occurs in upper-mantle slabs, a measure that can be used as a constraint on slab strength. Already from Benioff-zone geometries, it is clear that there is slab bending and contorting (e.g., Alpert et al., 2010; Hayes et al., 2012). Moment tensors provide further evidence of bending, folding, and shearing on top of the overall downdip stress patterns (Giardini and Woodhouse, 1984; Apperson and Frohlich, 1987; Bailey et al., 2012; Myhill, 2012). Moment tensor summations for the seismicity in the Tonga slab (by far the most seismically active Benioff zone) show that ∼50% of the sinking is accommodated by internal deformation (Fischer and Jordan, 1991; Holt, 1995), implying effective upper-mantle slab viscosities of no more than 2 orders of magnitude higher than the surrounding mantle (e.g., review by King, 2001; Alisic et al., 2010).

Lower-mantle seismic anomalies that are attributed to past subduction come in a wide range of shapes, from relatively linear and dipping (Farallon under Central America, northern Kurile, Aegean), to vertical and short (Marianas), flattened in top of lower mantle (Java-Borneo) to broad anomalies that do not have a slab shape, e.g., below the Indian plate and in the deepest mantle around the Pacific (Karason and Van der Hilst, 2000; Grand, 2002; Hafkenscheid et al., 2006). Many studies have concluded that the lower-mantle anomalies correspond to thicker volumes of material than the upper-mantle slabs, even considering limitations of seismic resolution. Several studies have estimated that some slabs, e.g., the Cocos, Java, and Hellenic slabs, thicken by a factor of 4–5 upon entering into the lower mantle, to widths of ∼400 km (Ribe et al., 2007; Loiselet et al., 2010), or even 400–700 km for the slabs below North America (Sigloch and Mihalynuk, 2013); while thickening in the lower-mantle slab anomalies below the Indian plate has been estimated to be about a factor of 3 (Hafkenscheid et al., 2006). Such strong thickening has commonly been attributed to slab buckling that starts in the transition zone (Ribe et al., 2007; Běhounková and Čížková, 2008; Lee and King, 2011). Thickening has been inferred to be less for some of the older penetrating slabs below the western Pacific (Marianas, Kermadec, Loiselet et al., 2010) but may still be a factor of two, an amount which can also be achieved by simple compressional thickening (Gurnis and Hager, 1988; Běhounková and Čížková, 2008).

At mid-mantle depths, there are two prominent belts of high-velocity anomalies (Van der Hilst et al., 1997; Bijwaard et al., 1998; Grand, 2002). One stretches below the Americas from the northern tip of Canada to Bolivia and is generally attributed to Farallon subduction (Grand, 1994; Van der Hilst et al., 1997; Bunge and Grand, 2000; Ren et al., 2007), although recently an interpretation involving additional oceanic plates has been advanced (Sigloch and Mihalynuk, 2013). The other anomaly more or less underlies the Tethyan suture zone stretching from the Aegean to Indonesia (Van der Hilst et al., 1997; Van der Voo et al., 1999b). These Tethys- and Farallon-attributed anomalies extend to depths from 700 to 1000 km down to 1700–2000 km depth (Van der Hilst et al., 1997; Bijwaard et al., 1998; Grand, 2002). Each of these belts is made up of a number of segments, and within the Tethyan anomalies several parallel belts, attributed to separate subduction zones, have been identified (Van der Voo et al., 1999b; Hafkenscheid et al., 2006; Nerlich et al., 2016). By comparing these anomalies with plate reconstructions, they have been estimated to represent subduction of the past 100–140 m.y. (Grand, 1994; Hafkenscheid et al., 2006; Sigloch and Mihalynuk, 2013). These ages imply lower-mantle slab sinking rates of 1–2 cm/yr (Grand, 1994; Hafkenscheid et al., 2006; Van der Meer et al., 2009; Sigloch and Mihalynuk, 2013; Butterworth et al., 2014; Domeier et al., 2016). Deeper anomalies may correspond to subduction going back as far as 200 m.y. ago and in a few locations even 300 m.y. (Richards and Engebretson, 1992; Ricard et al., 1993; Van der Voo et al., 1999a; Van der Meer et al., 2009), although it has been questioned how robust correlations between plate-motion history and anomalies below 2300 km depth are (Domeier et al., 2016).

As discussed above, slabs appear to thicken as they transit from the upper to the lower mantle. A further change in morphology occurs toward the deep mantle, where bulky fast regions rather than linear belts of anomalies are found. High-velocity anomalies in the deepest 500 km of the mantle more or less encircle the Pacific with a particularly large pond of fast material below eastern Asia, as well as others below the Americas and the southern Indian Ocean (Van der Hilst et al., 1997; Grand, 2002; Simmons et al., 2012). These changes in morphology with depth have been suggested to be a simple consequence of slab deformation during sinking (Grand, 2002; Jarvis and Lowman, 2007).

In sum, the tomographic evidence shows that slabs experience significant resistance at the base of the transition zone leading to flattening, buckling, and thickening, and further deformation to less linear forms on their journey through the lower mantle. Comparisons with plate reconstructions constrain the time scales of slab sinking into the lower mantle to 10–50 m.y. and through the rest of the mantle to 150–200 m.y. Even if many subduction zones have been long lived (e.g., along the circum-Pacific), the slabs are generally not continuous over the full mantle depth range, indicating that the sinking process is discontinuous, as previously noted by Grand (2002), which may reflect the influence of stagnation in the transition zone.

Most Cenozoic downgoing plate motions range between ∼5 and 10 cm/yr (Sdrolias and Müller, 2006; Müller et al., 2008; Goes et al., 2011; Zahirovic et al., 2015), and there is a correlation of increasing velocity with the age of the plate at the trench (Carlson et al., 1983; Lallemand et al., 2005; Goes et al., 2011). This range and the age trend have been interpreted as an expression of upper-mantle slab pull as driving force of the plates (Conrad and Lithgow-Bertelloni, 2002; Faccenna et al., 2007; Goes et al., 2011). Some segments of the Farallon plate moved notably faster (≥15 cm/yr) (Müller et al., 2008), as did the Indian plate before 50 m.y. ago (Lee and Lawver, 1995; Müller et al., 2008; Capitanio et al., 2010a). It has been suggested that these fast motions are due to: plate avalanching into the lower mantle (Goes et al., 2008), to the setting up of a whole mantle convection circulation, as opposed to upper-mantle confined cells (Becker and Faccenna, 2011; Faccenna et al., 2013), and for India, the push of the Reunion plume and consequent lowering of asthenospheric viscosities below the plate (Cande and Stegman, 2011; van Hinsbergen et al., 2011), a pull by double subduction (Jagoutz et al., 2015) or the removal of the Indian cratonic keel (Kumar et al., 2007), or the interplay of several of these factors (Zahirovic et al., 2015).

The large majority of trenches (∼70% of the segments throughout the Cenozoic) retreat (Sdrolias and Müller, 2006; Müller et al., 2008; Schellart et al., 2008; Goes et al., 2011; Williams et al., 2015). This is true irrespective of the hotspot reference frame used, but an Indo-Atlantic or global hotspot reference frame is most consistent with slab shapes and dynamic considerations (Schellart et al., 2008; Goes et al., 2011; Williams et al., 2015). Most trench motions are between –1 cm/yr (advancing) and +2 cm/yr (retreating) (Williams et al., 2015), and, on average, trench retreat accounts for ∼10% of the convergence velocities (Goes et al., 2011). Some trenches retreated significantly faster (motions up to 10 cm/yr) during part of the Cenozoic (Tonga, Izu-Bonin, Japan, and Calabria). A few trenches advance(d) (Marianas and Kermadec since ∼10 m.y. ago, and the Sumatran end of the Sunda trench and Alaskan end of the Aleutian trench at slow rates, <1 cm/yr). Throughout the Cenozoic, trench motions that exceed 1–2 cm/yr (both in advance and retreat) are only observed for older plates (>50 m.y.) (Fig. 3; Molnar and Atwater, 1978; Sdrolias and Müller, 2006; Goes et al., 2011).

It has been noted that although there is some correlation between absolute downgoing-plate motions and age, there is no correlation between present-day trench motions and slab morphology or dip (King, 2001; Lallemand et al., 2005; Billen, 2008; Fig. 3). On the other hand, it is also well established that there is a correlation between trench-motion history and plate morphology (e.g., Van der Hilst and Seno, 1993; Faccenna et al., 2001) (Fig. 4). In Figure 4, we compiled estimates of the total trench retreat (references in the caption), measured approximately perpendicular to the present-day trench during the past 50 m.y. Uncertainties in the measurements are probably ∼±200 km, due to uncertainties in the reconstructions and reference frame. It should be noted that the total trench motions were generally achieved in only part of these 50 m.y. (Clark et al., 2008), e.g., for Tonga during the past 35 m.y. in two phases (Sdrolias and Müller, 2006), for Scotia during the past 40 m.y. during several phases (Barker, 2001), for Calabria during the past 15 m.y. (Faccenna et al., 2001), and corresponded to trench velocities from 2 to ∼10 cm/yr. Most flattened slabs seem to have undergone trench retreat at a rate of at least 2 cm/yr at some point in their Cenozoic history (Sdrolias and Müller, 2006). When integrated over the past 50 m.y., all trenches, except India, which advanced by ∼1500 km (Replumaz et al., 2004; Hafkenscheid et al., 2006) (not included in Fig. 4), show a net retreat, over distances ranging from a few hundred km to 2000 km.

The history of trench retreat generally correlates with the morphology of the slab in the transition zone (Fig. 4). The Indian slab (not shown) is the only one that is rolled over and the only one that underwent a net trench advance. All other slabs are reclining, as expected from a net trench retreat. Most slabs that penetrate into the lower mantle subducted at relatively stable trenches, while flat slabs all have a history of significant trench retreat (>800 km). For example, along the Tonga-Kermadec, Izu-Bonin-Marianas, Japan-Kurile-Kamchatka trenches, the penetrating parts of the slab subducted where the trench moved the least over the past 50 m.y. The trenches bordering the west coast of the Americas overall moved westward by ∼1000 km (Grand, 1994; Ren et al., 2007; Müller et al., 2008); however, due to the shape of the trench and relative North and South American motions, trench motion reached a minimum near Peru and was almost zero off Central America, the only place where a slab appears to penetrate directly and deeply into the lower mantle.

For a few of the penetrating slabs, under the Antilles, Java, and Kermadec, the trenches retreated as much as for some of the flat slab cases. The Antilles slab seems to penetrate even though most plate tectonic studies (starting with Ross and Scotese, 1988) have inferred that its trench has undergone >1000 km of retreat. However, some studies have challenged this view and concluded the trench was formed only ∼300 km westward of its current position (James, 2009). It is also possible that interaction with the nearby Cocos plate has played a role in forcing the penetration of the Antillean slab with the transition zone. In several reconstructions (Replumaz et al., 2004; Müller et al., 2008), the Java trench retreated over ∼2000 km, largely achieved while the South China Sea opened, and the flat slab segment at ∼1100 km depth has been linked to this retreat (Replumaz et al., 2004). On the other hand, a recent study attributed both the extension phase and the (inferred to not be flattened) slab segments in the lower mantle to different subduction systems (Hall and Spakman, 2015). Kermadec is another location where a lower-mantle flat slab segment, of up to ∼800 km length, similar to the reconstructed amount of trench retreat, has been suggested (Fukao and Obayashi, 2013), although others interpret the deeper slab structure as thickened rather than flattened (Schellart and Spakman, 2012).

The lengths of the flat slabs are similar (most in the range 800–1600 km) to the estimates of trench retreat (Fig. 4), where it should be borne in mind that uncertainties in flat-slab length are likely larger than those in the trench motion estimates, due to variable seismic resolution and subjectivity in defining the extent of the slab.

The many modeling studies done, on how slabs interact with the mantle transition zone, show that there are two main sets of factors (Fig. 1) that can affect whether slabs stagnate or penetrate: (1) the mantle resistance to flow through the upper-lower mantle boundary; and (2) the shape and strength of the slab as it starts interacting with this boundary, where the ability of the trench to retreat plays a crucial role in the former. We will discuss the relative importance of each of these below.

Several factors have been proposed to contribute to the resistance to sinking at the base of the upper mantle, most prominently a jump in viscosity between upper and lower mantle (Section 3.1.1) and the deflection of the endothermic phase transition from ringwoodite to the postspinel phases bridgmanite and magnesiowüstite in the cold slab, which locally reduces the slab’s negative buoyancy (Section 3.1.2). Some studies have suggested that additionally a (small) change in intrinsic density due to variations in chemistry between upper and lower mantle may contribute (Section 3.1.3).

Modeling of the geoid, postglacial rebound, free-air gravity, and dynamic topography observations constrains the average viscosity of the lower mantle to be a factor 10–100 higher than the average upper-mantle viscosity (e.g., Chen and King, 1998; Panasyuk and Hager, 2000; Mitrovica and Forte, 2004). Similar increases in viscosity of the lower relative to the upper mantle have been inferred from models that reproduce lower-mantle slab sinking velocities estimated from the match of plate reconstructions with tomography (Čížková et al., 2012).

Because experiments show that viscosity varies exponentially with pressure (and exponentially with the inverse of temperature), part of this increase may be a gradual pressure effect (in particular in the shallow lower mantle, e.g., Marquardt and Miyagi, 2015; King, 2016), a trend that will, to some extent, be offset by the decrease in viscosity due to adiabatically increasing temperature with depth (Ranalli, 1995; King, 2016). On the other hand, the phase transition in most of the mantle minerals at the base of the upper mantle is a logical place for a step in viscosity to occur (Karato, 1989; Faccenda and Dal Zilio, 2017). Some joint inversions also require a jump close to this depth (Mitrovica and Forte, 2004). Furthermore, experiments show that bridgmanite, which comprises ∼80% of the lower-mantle composition, is a relatively strong phase (Yamazaki and Karato, 2001). Additionally, geoid modeling of the signature above subduction zones indicates that the resistive effect of an endothermic phase transition has a low signal and an increase in viscosity at the boundary, by a factor 30–100, is required to explain smaller-scale geoid structure (King, 2002; Tosi et al., 2009), where it should be borne in mind that dynamically, a strong viscosity-depth gradient will have a similar effect as a sharp jump.

Global-scale models of mantle convection find that the range of upper- to lower-mantle viscosity jumps inferred from observations does not lead to layered convection or lead to a significant accumulation of slab material in the transition zone (Van Keken and Zhong, 1999; Yanagisawa et al., 2010). Regional models, of a single subduction system, do generate a range of slab morphologies in response to such a viscosity increase but only if the trench can move freely (Kincaid and Olson, 1987; Gurnis and Hager, 1988; Guillou-Frottier et al., 1995; Christensen, 1996; Olbertz et al., 1997; Garel et al., 2014). The regional-model slab morphologies resemble the diversity of geometries imaged tomographically (Fig. 5).

However, in models with only a viscosity jump, even if the slabs flatten, lower-mantle sinking velocities are not negligible, and as a result, they are not able to generate flat-slab segments of the >1000 km length inferred from tomographic images (Garel et al., 2014). On the other hand, if there is no viscosity jump, it is difficult to develop any flat slabs (Tagawa et al., 2007; Torii and Yoshioka, 2007; Yanagisawa et al., 2010; and our own models for Fig. 6). At the high end of the geoid and rebound-based viscosity jumps—a factor of 100 or higher—slab sinking into the lower mantle is strongly impeded (Běhounková and Čížková, 2008; Loiselet et al., 2010), trench retreat decreased (Čížková and Bina, 2013; Garel et al., 2014), and the range of slab morphologies is limited to flattening and folding, even if slab cores are 2–3 orders of magnitude more viscous than the surrounding mantle (Garel et al., 2014).

The downward deflection of the endothermic ringwoodite to bridgmanite + magnesiowüstite transition within the cold slabs leads to a local increase in buoyancy of the slabs, which hampers their further descent (e.g., Bina et al., 2001; Faccenda and Dal Zilio, 2017). In a pyrolitic composition, along a background (1300 °C potential) mantle temperature profile, the density jump at the ringwoodite to postspinel transition is 220 kg/m³ (calculated with the database from Stixrude and Lithgow-Bertelloni, 2011) (a 5.5% density increase relative to the density above the transition). Compared with the excess density of slabs, which ranges from ∼30 kg/m³ for young plates to up to 100 kg/m³ for very old plates (Cloos, 1993), it is clear that a downward deflection of the transition makes the slabs locally buoyant, thus providing a very strong stagnating force.

The slab’s buoyancy due to the deflected phase transition in the downgoing slab is only a localized effect. Once part of the slab manages to sink below the transition, its negative buoyancy is regained. The joint influence of the density jump and phase boundary deflection is described by the phase buoyancy parameter, a measure of the local buoyancy due to the phase transition compared with overall thermal convective driving force (Christensen and Yuen, 1985). The most recent experimental estimates for the Clapeyron slope, Γ, of the ringwoodite to postspinel transition are between –0.5 and –2.0 MPa/K (Hirose, 2002; Fei et al., 2004; Katsura et al., 2004; Litasov et al., 2005; Stixrude and Lithgow-Bertelloni, 2011). At the lowest temperatures (below 900–1400 °C [Stixrude and Lithgow-Bertelloni, 2011; Holland et al., 2013]), the also strongly endothermic transition (Γ between –2 and –5 MPa/K) from ringwoodite through akimotoite to bridgmanite and magnesiowüstite may further contribute to the slab’s phase buoyancy (Faccenda and Dal Zilio, 2017).

Many modeling studies have investigated how slab material may be sequestered in the transition zone by the effect of a negative Clapeyron slope or a combined viscosity jump and Clapeyron slope effect, starting with Richter (1973), who showed that it could somewhat hinder convection. Christensen and Yuen (1984) did a detailed investigation of the interaction of a vertically sinking slab across a boundary with an endothermic phase change and/or chemical density changes of variable strength. They showed that for very high slopes, more than –4 to –6 MPa/K, the phase buoyancy effect is so strong that convection becomes layered. Many subsequent models have shown that for intermediate Clapeyron slopes, above –2 to –3 MPa/K (i.e., at the high end of current experimental constraints), partial layering of mantle convection is possible, with catastrophic penetration events (“avalanches”), which affect the global flow field, once sufficient material has accumulated above the phase boundary (Machetel and Weber, 1991; Tackley et al., 1993; Solheim and Peltier, 1994). At lower Clapeyron slopes, the phase boundary has little effect in these global-scale models.

However, such global models generally do not fully capture the effects of plate strength and trench mobility. Many regional subduction models have documented that when trench retreat is possible or imposed, partial stagnation can already occur at relatively low Clapeyron slopes (from –0.5 to –1 MPa/K, i.e., well within the experimentally expected range) (Christensen, 1996; Torii and Yoshioka, 2007; Agrusta et al., 2017). This will be discussed further in Section 3.2.

The experiments from Agrusta et al. (2017) (Fig. 6) illustrate that in models where trenches are allowed to migrate freely, changing the Clapeyron slope of the postspinel transition between –0.5 and –3 MPa/K has a stronger effect on whether stagnation or partial stagnation occurs than changing the mantle viscosity jump between 5 and 30, although without a viscosity jump, stagnation is not easy to achieve either. These results agree with those from Yanagisawa et al. (2010), who found in a 3-D spherical model with mobile trenches that a mix of stagnant and penetrating slabs is only created through the joint effect of a viscosity jump and endothermic phase boundary, where the latter is key in allowing stagnation and hence intermittent sinking of slabs into the lower mantle.

The exothermic phase transition expected at ∼410 km depth when olivine transforms to its wadsleyite phase can encourage penetration, because it enhances slab sinking and piling in the transition zone (Christensen, 1996; Čížková et al., 2002; Běhounková and Čížková, 2008). However, even with the effect of such an exothermic 410, slabs that can induce trench retreat can still get stalled by Γ as low as –1 MPa/K (Agrusta et al., 2017).

On the other hand, it has been proposed that the transition from olivine to wadsleyite may be kinetically delayed in the core of cold, i.e., old, slabs, to depths 200–300 km below the equilibrium depth. This would reduce the negative buoyancy of the oldest plates and aid their stagnation in the transition zone (Schmeling et al., 1999; Bina et al., 2001; Tetzlaff and Schmeling, 2009; Agrusta et al., 2014). The olivine component (in olivine, wadsleyite, or ringwoodite phase) makes up only ∼60% of a peridotitic composition above 660 km. Experimental studies indicate that the pyroxene to garnet phase transitions in the remaining 40% of the composition, as well as in the crustal part of the slab, may be also be kinetically inhibited at typical slab temperatures (Hogrefe et al., 1994; Nishi et al., 2008; Van Mierlo et al., 2013). The delay of this transition, from its equilibrium depth between 200 and 300 km to possibly as deep as the base of the transition zone, adds so much buoyancy that it could aid the stagnation of both younger and older plates, although the effect should be stronger in the latter (Agrusta et al., 2014; King et al., 2015).

Near 660 km depth, it has been suggested that the effect of the ringwoodite to postspinel transition with a negative Clapeyron slope and the transition from garnet to bridgmanite, which also occurs near this depth and has a strongly positive Clapeyron slope, may largely cancel out any dynamic effects (Tackley et al., 2005). The depth range of the garnet transition and its Clapeyron slope are uncertain (Hirose, 2002; Xu et al., 2008; Stixrude and Lithgow-Bertelloni, 2011; Holland et al., 2013). In any case, the transition of the garnet components is rather diffuse and may spread over a depth interval as large as 100–250 km. The Clapeyron slope of the garnet transitions may also be low below ambient mantle temperatures. It is furthermore predicted to have a milder buoyancy effect than the transition of ringwoodite (Bina et al., 2001; Faccenda and Dal Zilio, 2017). It remains to be tested in regional subduction models with mobile trenches what the dynamic effect is of the interplay of these transitions. However, if the garnet transition spreads over a larger depth than the thickness of the slab, its effect would likely be significantly less than that of the olivine component transitions (Christensen and Yuen, 1985; Bina et al., 2001).

The modeling done so far, with a viscosity jump and a single endothermic phase transition at 660 km depth, strongly indicates that phase buoyancy plays a key role in causing slab stagnation. In our study (Agrusta et al., 2017) (where we used a Δρ_ph of 350 kg/m³, as appropriate for a 100% olivine composition), we found that at Δη = η_LM/η_UM = 30, old plates can stagnate for Γ stronger than –0.7 to –0.8 MPa/K. For a pyrolitic density jump, this translates into a critical Clapeyron slope of –1.1 to –1.3 MPa/K. Thus partial stagnation is quite plausible for properties of the ringwoodite to postspinel transition well within the recently determined experimental range.

Although the most common assumption is that the whole mantle has a pyrolite-like composition (e.g., McDonough and Sun, 1995; Lyubetskaya and Korenaga, 2007)—with some amount of heterogeneity from recycled subducted slabs (Xu et al., 2008)—some studies have questioned this assumption. Discrepancies between the velocity and density jumps in spherical seismic models and estimates from experiments on phase transitions in olivine or pyrolite, as well as Earth evolution and geochemical and cosmochemical arguments have been used to advocate that there is an intrinsic chemical difference between upper and lower mantle (e.g., Stixrude et al., 1992; Khan et al., 2008; Javoy et al., 2010). If there is any chemical difference, it would need to result in a density jump at 660 km depth that is less than 6%, or else it would lead to fully layered mantle convection (Christensen and Yuen, 1984), which seismic tomography, plate evolution modeling, and considerations about the convective and thermal evolution of the Earth have shown is implausible (Silver et al., 1988; Ricard et al., 1993; Van der Hilst et al., 1997; McNamara and Van Keken, 2000). Furthermore, if it does not lead to layered convection, a chemical interface cannot survive over many convective cycles, although chemical gradients may persist (Van Keken and Zhong, 1999; Tackley et al., 2005; Brandenburg and Van Keken, 2007). For example, Ballmer et al. (2015) invoked compositional density gradients, due to a progressive enrichment in basaltic components with depth in the lower mantle, as a mechanism to stagnate slabs within the top of the lower mantle. Our compilation of transition-zone slab morphologies in Figure 2 shows, however, that slabs very rarely stagnate below 660, if at all.

From the discussion above, it is already clear that trench retreat plays an important role in allowing some slabs to stagnate in the mantle transition zone, as was first proposed by Van der Hilst and Seno (1993). Trench motion turns out to be intimately linked to the strength of the subducting plate as it bends into the mantle. Furthermore, slab strength within the transition zone has been suggested to play a role in governing how easily slabs penetrate (Karato et al., 2001). Below we will discuss the role played by these slab-intrinsic factors and by trench mobility.

Most slabs actually do not just move in an along-dip direction but also exhibit retrograde motion (seaward away from the upper plate) that results in trench retreat (Elsasser, 1969; Chase, 1978; Garfunkel et al., 1986). Such retrograde motion is the consequence of vertical sinking of the plate, as expected for the most basic form of free subduction, driven by the slab’s negative buoyancy and not significantly hampered by the upper plate or mantle convective forcing (Kincaid and Olson, 1987; Zhong and Gurnis, 1995; Becker et al., 1999; Funiciello et al., 2003a; Funiciello et al., 2003b; Schellart, 2008). Motivated by the free-subduction models of Kincaid and Olson (1987), Van der Hilst and Seno (1993) proposed that trench retreat could be responsible for the slab flattening as tomographically observed in the transition zone below the western Pacific.

Zhong and Gurnis’s (1995) model that incorporates mobile plate boundaries clearly established that trench migration can produce slabs that stagnate over significant distances in the transition zone. Many other, regional-scale models prescribed trench motions and determined how much motion is required for flattening (Griffiths et al., 1995; Guillou-Frottier et al., 1995; Christensen, 1996; Olbertz et al., 1997; Čížková et al., 2002; Torii and Yoshioka, 2007). In this way, they inferred critical trench motions or dips. For example, Christensen (1996) found a critical retreat rate of 2–4 cm/yr. Further work indicates that the critical rate may scale with sinking velocity, i.e., a critical dip may be more appropriate, and furthermore, the cutoff depends on transition-zone resistance and likely also on slab strength (Torii and Yoshioka, 2007; Agrusta et al., 2017), so that there is not a single critical retreat rate.

Models in which trench motion develops dynamically show that variable trench mobility governed by slab internal and external forcing is the most likely cause for variable modes of subduction-transition zone interaction (e.g., Ribe, 2010; Stegman et al., 2010b; Čížková and Bina, 2013; Garel et al., 2014).

Models in which subduction and retreat velocities develop self-consistently have been used to study what factors control trench retreat, at first without any interface at the base of the upper mantle or with an impenetrable lower mantle, all without any forcing from an overriding plate. In such “free” subduction models, subduction usually involves trench retreat even if there is no boundary that hampers sinking at the base of the mantle (Ribe, 2010; Fourel et al., 2014).

The speed of trench retreat is controlled by subducting-plate strength and buoyancy (Bellahsen et al., 2005; Capitanio et al., 2007; Di Giuseppe et al., 2008; Schellart, 2008; Ribe, 2010; Stegman et al., 2010b; Capitanio and Morra, 2012; Fourel et al., 2014) (Fig. 7). Both high subducting-plate density and high resistance to bending at the trench encourage trench retreat. This can be understood as follows. The time that the plate has available to bend is controlled by its sinking rate, which again is governed by its density (slab pull). If the sinking time is too short to allow the plate to bend from a horizontal to a vertical geometry, the trench will retreat while the slab sinks. The stronger the resistance to bending, i.e., the higher the time that would be required to viscously bend the plate at the trench, the more retreat will occur.

In some models, slabs with high viscosity and low densities lead to trench advance (Bellahsen et al., 2005; Di Giuseppe et al., 2008; Funiciello et al., 2008; Stegman et al., 2010b). The advance occurs when plates are able to (i.e., have sufficient sinking time to) fully bend over but do not have a chance to unbend before reaching the base of the upper mantle (Ribe, 2010). Interaction of an overturned plate (dip >90°) with a resisting upper- to lower-mantle interface triggers trench advance (Bellahsen et al., 2005; Funiciello et al., 2008; Schellart, 2008; Ribe, 2010). Weaker plates will unbend in response to the viscous mantle flow driven by the sinking and moving plate and thus arrive at the base of the upper mantle at an angle less than 90°. Strong and fast plates never bend beyond 90. For a thermal aging trend of oceanic plates, density and viscosity will both increase, and retreating modes are expected to dominate. Furthermore, if the elastic effects on rheology are considered, plates unbend more easily, in particular when they have a higher viscosity (i.e., higher Maxwell time, which equals viscosity over the elastic Young’s modulus) and hence behave more elastically (Farrington et al., 2014; Fourel et al., 2014). Therefore, the case of trench advance driven purely by the downgoing plate is probably not common on Earth. More likely, advance is the result of additional forcing either at the trailing end of the downgoing plate (like ridge push [Capitanio, 2013]) or from the upper plate or the arrival of buoyant features at the trench (Capitanio et al., 2010a; Magni et al., 2012; Fourel et al., 2014; Čížková and Bina, 2015; Goes et al., 2014).

Thermomechanical models (Garel et al., 2014; Agrusta et al., 2017) confirm that old plates, with high density and high viscosity have a higher propensity to retreat than younger plates. In principle, older plates should be less hampered by phase buoyancy, and indeed this is the behavior found in models of sinking plate segments of different thermal buoyancy (Ballmer et al., 2015). However, self-consistent subduction models show that their stronger tendency for retreat leads to older plates flattening and stagnating in the transition zone, while, for intermediate Clapeyron slopes and viscosity jumps less than 100, young plates tend to penetrate at relatively steep dips (Goes et al., 2008; Garel et al., 2014; Agrusta et al., 2017) (Fig. 6). Variations in plate age are thus a likely contributor to the observed variation in styles of slab transition zone interaction.

The insights into the role of subducting plate strength and buoyancy as the primary controls of trench motion have mostly come from 2-D models. In 3-D, additional variations in trench motion occur because of along-strike plate bending in response to flow around the plate edges (Funiciello et al., 2006; Morra and Regenauer-Lieb, 2006; Schellart et al., 2007; Loiselet et al., 2009; Stegman et al., 2010a; Li and Ribe, 2012), leading to lower retreat of the plate edges and, for very wide plates, a trench stagnation point can develop in the center. Lateral changes in slab buoyancy can also lead to variations in trench motions, which depend on the scale of the features of different buoyancy and their relative contribution to overall slab pull (Martinod et al., 2005; Morra et al., 2006; Goes et al., 2008; Mason et al., 2010; Magni et al., 2014; Goes et al., 2014).

Subduction zones are in reality never completely “free” to move, and it is a challenge to disentangle the respective roles and feedbacks of subducting versus neighboring plates (upper and side plates) dynamics. While the subducting plate generally drives trench motion, interplate coupling and upper-plate forcing provide resistance. In some cases, the upper plate may provide an additional driving force (Van Hunen et al., 2002; Arcay et al., 2008; Van Dinther et al., 2010; Čížková and Bina, 2015). Stronger plate coupling limits trench mobility and can, when high enough, even preclude subduction (De Franco et al., 2006; Běhounková and Čížková, 2008; Androvičová et al., 2013) (Fig. 8). Upper-plate mobility depends not only on the forces that the plate experiences at its other boundaries, but in addition, thicker, more buoyant and longer overriding plates provide more resistance to trench motion (Zhong and Gurnis, 1997; Capitanio et al., 2010b; Van Dinther et al., 2010; Capitanio et al., 2011; Garel et al., 2014; Holt et al., 2015).

The effect of the upper plate is often studied by using kinematic conditions, most commonly by prescribing a fixed versus mobile trench (e.g., Čížková and Bina, 2013; Agrusta et al., 2017). The choice of kinematic conditions can lead to quite different slab stress patterns (Čížková et al., 2007), and it needs to be borne in mind that kinematic forcing may not allow the energetically most favorable modes of subduction (Han and Gurnis, 1999). In models where the trench is held fixed, higher slab strength is found to encourage penetration (Zhong and Gurnis, 1994; Arredondo and Billen, 2016), contrasting with dynamic trench-motion models where higher slab strength aids trench retreat and slab flattening (e.g., Zhong and Gurnis, 1995; Capitanio et al., 2007).

Trench retreat is not a consequence of the slab’s interaction with an interface that hampers sinking. However, trench motion is enhanced and modulated by interaction with a viscosity and/or phase change (Becker et al., 1999; Bellahsen et al., 2005; Capitanio et al., 2010b; Čížková and Bina, 2013; Garel et al., 2014; Agrusta et al., 2017). As the slab accumulates in the transition zone, alternating phases of somewhat higher and lower retreat velocities (with accompanying changes in dip) tend to develop (Capitanio et al., 2010b; Lee and King, 2011; Čížková and Bina, 2013; Garel et al., 2014). These fluctuations are an expression of slab buckling, because even where the slab flattens at the base of the transition zone, the deformation is essentially a buckling response (Houseman and Gubbins, 1997; Ribe et al., 2007), as is seen most clearly in models with weaker slabs or ones where sinking is enhanced (Lee and King, 2011; Čížková and Bina, 2013; Cerpa et al., 2014). Trench-motion fluctuations due to buckling can be enhanced if the upper plate can break and heal (Clark et al., 2008). Models where the slab buckles in its interaction with the transition zone produce variations in upper-plate stress on time scales similar to those observed in, for example, the Andes or backarc spreading phases for Tonga (Clark et al., 2008; Capitanio et al., 2010b; Lee and King, 2011).

Karato et al. (2001) proposed that slab weakening in the transition zone associated with small grain size around a wedge of metastable olivine in the slab’s core may be key for trapping slab material in the transition zone and that this would be most effective for intermediate age slabs, which are hot and thin enough to bend but with a core that is cold enough to induce a sluggish olivine-to-wadleyite phase transition. Čížková et al. (2002) tested the effect of weakening the slab in the transition zone and found it was secondary to that of trench motion. King and Ita (1995) also found that slab strength does not exert a major influence on whether slabs penetrate or not, when the trench is free to move. Others further tested the effect of local slab weakening in the transition zone and confirmed that the effect is not dominant, but it may be sufficient to stagnate a slab that is otherwise marginally penetrating (Tagawa et al., 2007; Agrusta et al., 2017). The buckling of weaker slabs can actually increase the slab’s Stokes’ sinking velocity (which scales with diameter squared) and hence help it penetrate a high-viscosity lower mantle rather than stagnating it. These results imply that strength of the subducting plate at the trench is much more important in controlling how slabs interact with the transition zone than any changes in slab strength within the transition zone.

Several studies did find that over the longer term, weakening helps stagnant slabs destabilize and flush into the lower mantle (Nakakuki et al., 2010; Agrusta et al., 2017). This flushing is akin to a Rayleigh-Taylor instability, for which it is well established that it is facilitated by lower viscosities. Several studies have found that such stagnant slab destabilization can affect mantle flow and overlying plate velocities (Pysklywec et al., 2003; Pysklywec and Ishii, 2005; Motoki and Ballmer, 2015) (see further in Section 4.2).

Aging of the subducting plate and the consequence of this on trench mobility can result in a change in transition-zone dynamics in time, from penetration to stagnation (Agrusta et al., 2017). This would lead to a flat-slab morphology with a leading end that has entered the deep mantle. The opposite switch is not as easily induced, but decreasing plate age at the trench does lead to decreased trench motion and can under certain conditions trigger lower-mantle sinking of a previously stagnated slab. Detachment of such increasingly buoyant slabs can facilitate the start of lower-mantle sinking (Agrusta et al., 2017). This switch would lead to slabs that start entering the lower mantle at the hinge between dipping and flattened slab segments as shown in Figure 9C. Changing the upper-plate resistance to trench motion can be a quite efficient mechanism to induce switches from penetration to stagnation as well as vice versa (Agrusta et al., 2017). However, regional models usually test kinematic extremes of upper-plate forcing (either fixed, free, or forced upper plates), and the effects may be more subtle or different in a global dynamic system.

Dynamic modeling identifies variable trench motions as the most probable cause for the different ways in which slabs interact with the transition zone (Section 3). This proposed link is consistent with the observation that slab morphology correlates well with the history of trench motions (Fig. 4). Furthermore, the models predict the largest trench motions for the oldest densest plates, and indeed, although trench motions do not simply follow an age trend (Sdrolias and Müller, 2006; Goes et al., 2011) (Fig. 3), the largest trench velocities have been found to occur for plates older than ∼50 m.y. (Molnar and Atwater, 1978; Sdrolias and Müller, 2006; Goes et al., 2011), and slabs of old plates seem to stagnate more often than young plate slabs (King et al., 2015; Section 2.1). This indicates that the forcing by downgoing plates older than ∼50 m.y. is generally sufficiently large to overcome other forces (e.g., by the upper plate or the convective mantle) that may impede trench motions, consistent with the long-standing conclusion that the downgoing plate pull is the primary force driving plate motions (Forsyth and Uyeda, 1975; Zahirovic et al., 2015). A recent global dynamic modeling study also concluded that subducting plates provide the main forces that drive formation and migration of plate boundaries (Mallard et al., 2016).

Variations in slab transition zone interaction in space and time are thus likely the consequence of variations in trench mobility, which can be caused by either changes in the downgoing-plate driving force (e.g., due to variations in plate age or arrival of buoyant features on the subducting plate at the trench) or by changes in the forces that constrain or force trench motion, e.g., by the upper plate (Agrusta et al., 2017). For example, the arrival of continental material at the trench has been shown to be able to induce trench advance (Magni et al., 2012), like that observed for the trench north of India. Thus the start of continental collision may be responsible for a switch from a retreating subduction mode, which has been proposed to explain lower-mantle slab anomalies below India (Hafkenscheid et al., 2006), to the Cenozoic advancing mode (Capitanio et al., 2010a). Similarly, arrival of the Australian continent at the Banda trench may have halted its trench retreat and induced the sinking of a slab segment that previously flattened in the transition zone. On the other hand, the aging of the Pacific plate as it subducted along the Tonga and Japan trenches during the Cenozoic may have led to the change from penetrating to stagnating subduction that has been inferred from seismic tomography (Van der Hilst, 1995; Goes et al., 2008) (Fig. 2).

Several previous studies have noted that transition-zone slabs will only remain stagnant for a limited amount of time (Zhong and Gurnis, 1995; Christensen, 1996; Pysklywec and Ishii, 2005; Nakakuki et al., 2010; King et al., 2015). As slabs warm up, phase boundary topography reduces and hence phase buoyancy decreases, while the slab’s density anomaly diffuses, although the overall integrated thermal density anomaly decreases only very slowly (Turcotte and Schubert, 1982). The diffusion time scale over which the slabs warm sufficiently to make the phase buoyancy effect too small to keep a slab stagnant has been estimated to ∼250–300 m.y. (Christensen, 1996). As a slab warms up, it will also weaken, which further facilitates its deformation and destabilization (Nakakuki et al., 2010).

As soon as part of the slab sinks below the phase boundary, that part will be gravitationally unstable and can eventually pull the rest of the slab along into the lower mantle. Various regional-scale models, with or without free trench motion, show that if slabs have a similar strength to that of the surrounding (lower) mantle, they can destabilize in a few tens of millions of years (Griffiths et al., 1995; Zhong and Gurnis, 1995; Nakakuki et al., 2010; Motoki and Ballmer, 2015), while slabs that have viscosities that exceed those of their surroundings by about two orders of magnitude take between 100 and over 300 million years to destabilize (Zhong and Gurnis, 1995; Nakakuki et al., 2010; Agrusta et al., 2017). It would be expected that, due to the Arrhenius temperature sensitivities of mantle flow laws (Ranalli, 1995), the colder slab has a higher viscosity than its surroundings. The occurrence of some earthquakes in the flat parts of slabs (Lundgren and Giardini, 1994) is an indication that they are indeed stronger than the surrounding mantle, because the seismicity requires that at least some flat-slab segments retain sufficient strength to allow elastic energy accumulation. Furthermore, although grain-size reduction at a metastable phase transition may lead to weakening of the slab core (Karato et al., 2001), this is unlikely to persist into the stagnant slab. This implies that the time scales of destabilization are probably closer to several hundred million years.

Such time scales are much longer than the time scales of stagnation estimated from the seismically imaged volume of material in the transition zone, which range from ∼20–60 m.y., even if some amount of buckling is accounted for (Section 2.1) (Tajima and Grand, 1998; Schmid et al., 2002). This would indicate that penetration is commonly triggered before the stagnant slab becomes negatively buoyant (Agrusta et al., 2017).

It has been suggested from models (Zhong and Gurnis, 1995; Tackley et al., 2005) and observations (Goes et al., 2008) that the sinking of previously stagnant slabs leads to large increases in subducting-plate velocities. However, as discussed in Section 3.2, neither global nor regional models (e.g., Christensen, 1996; Yanagisawa et al., 2010) confirm that for the postspinel Clapeyron slopes most likely in the mantle, lower-mantle sinking of previously stagnant slabs is associated with large peaks in plate velocities. Some regional plate acceleration by up to a factor of 2 can occur (Christensen, 1996; Pysklywec and Ishii, 2005; Nakakuki et al., 2010; Agrusta et al., 2017), and this may explain some of the plate motions that are relatively high for their slab’s buoyancy (as identified by Goes et al. [2008]), such as that of the subduction below Japan in the early Cenozoic. However, the high Mesozoic velocities for the motion of India and early Cenozoic motions for the Farallon plate exceed current average velocities by a factor of 2–4 and, hence, probably require an additional driving mechanism, e.g., an increase convective vigor due to the lowering of viscosity by the influx of hot material from depth and/or the development of a large whole-mantle–scale convective cell (Becker and Faccenna, 2011; Cande and Stegman, 2011; van Hinsbergen et al., 2011).

As discussed in Section 1.3, once slabs enter the lower mantle, they appear to thicken, by as much as a factor 2–5. Models predict significant buckling and thereby thickening when the viscosity contrast encountered by the slabs at the base of the transition zone is at least a factor of 20 or 30 (Gurnis and Hager, 1988; Gaherty and Hager, 1994; King, 2002; Běhounková and Čížková, 2008; Lee and King, 2011; Čížková and Bina, 2013) The increase in slab pull due to a shallowing olivine to wadsleyite phase transition in the slab ∼400 km depth may enhance buckling (Běhounková and Čížková, 2008). Finally, some amount of transition-zone slab buckling is consistent with focal mechanisms (e.g., Myhill, 2012), as well as with high seismic strain rates inferred for the Tonga slab (Giardini and Woodhouse, 1984; Bevis, 1988; Holt, 1995).

The ability to buckle and thicken puts a limit on transition-zone slab strength, which should not be stronger than ∼2–3 orders of magnitude higher viscosity, or a yield stress of 100–300 MPa (Houseman and Gubbins, 1997; Běhounková and Čížková, 2008; Lee and King, 2011). Effective slab viscosities of ∼2 orders of magnitude higher than the surrounding mantle are also most consistent with the geoid signal of subduction zones (King and Hager, 1994; Moresi and Gurnis, 1996; Zhong and Davies, 1999) On the other hand, the existence of deep Benioff zones requires that those slabs retain sufficient strength to store elastic energy. This may require that slabs have an effective viscosity at least 2–3 orders of magnitude higher than that of the surrounding mantle (Fourel et al., 2014), i.e., a strength near the upper bound defined by the other constraints.

Upon sinking, all slabs will develop a more rounded form because they try to achieve a shape that sinks more efficiently than the planar geometry they start subducting with (Jarvis and Lowman, 2007; Loiselet et al., 2010). For simple vertical sinkers, this is achieved by rounding and widening of the leading edge (Jarvis and Lowman, 2007; Loiselet et al., 2010). For flattened slabs, this may be achieved by folding. Some examples of how slabs evolve their shapes while sinking in the transition zone are shown in Figure 9. A combination of buckling, folding, rounding, and diffusive widening leads to a significant change in slab shape as they sink through the lower mantle. Note that stagnant slab shapes may be preserved for a while as such slabs start sinking into the lower mantle (Fig. 9B). This could be an explanation for the apparent flat slab at 1100 km depth below Borneo. Folding slabs have not been proposed previously; however, it may be worth re-examining some complex or broad high-velocity anomalies below the transition zone in this light (e.g., those from Sigloch [2011] below North America).

In comparing plate reconstructions with tomography, it is often assumed that the mantle anomalies reside more or less (within ∼500 km) below the location of the trench where the slab originally subducted (Van der Meer et al., 2009; Sigloch and Mihalynuk, 2013). However, models show that upon sinking into the lower mantle, previously stagnant slabs may end up as far as 1000–2000 km displaced laterally from where they originally subducted (Billen, 2008; Garel et al., 2014; models in Fig. 9). Both the shape evolution and lateral displacement that accompany it will make it more challenging to relate lower-mantle anomalies to plate-motion histories and increasingly so with depth, which may explain the difficulty of correlating plate-motion histories with the deepest mantle structures (Domeier et al., 2016).

Seismic tomography and Benioff seismicity highlight that subducting slabs encounter resistance to sinking in the transition zone and that they can respond to this resistance by either buckling and thickening before they sink into the lower mantle or by flattening and apparently stagnating at the base of the transition zone over distances of 500–2000 km for 20–60 m.y. From modeling and observational studies over the past 50 years, we now have a quite good understanding about what governs this behavior and what its implications are for plate motions and time scales of mantle mixing. A summary of the main insights:

• 1. When slabs stagnate, they mainly do so above 700 km depth (Fig. 2). Very few, if any, are flat in the upper part of the lower mantle. Slabs that appear flattened in the shallow lower mantle are sometimes only apparently flat (when the cross sections studied are oblique to the actual structures) or may represent slabs that previously flattened in the transition zone and subsequently started sinking through.

• 2. An increase in viscosity at the base of the transition zone by a factor 20–50 provides sufficient resistance to sinking to induce slab buckling and/or trench retreat. Such an increase in viscosity also explains the dominance of downdip compression below 300 km depth seen in Benioff seismicity globally (Vassiliou and Hager, 1988; Alpert et al., 2010), as well as the local geoid signature around subduction zones (Chen and King, 1998; King, 2002; Tosi et al., 2009). It does not, however, sufficiently stagnate slabs to accumulate flat slabs of up to 2000 km long (Garel et al., 2014).

• 3. The endothermic ringwoodite to postspinel phase change takes place at a larger depth in the cold slab than in the surrounding mantle, thus providing the additional buoyancy that can make slabs flatten over distances >1000 km and stagnate in the upper mantle over time scales of several tens of millions of years or more (Yanagisawa et al., 2010; Agrusta et al., 2017). However, for the likely values of the Clapeyron slope of this transition (between –0.5 and –2 MPa/K), stagnation only occurs when trench retreat is possible (Yanagisawa et al., 2010; Agrusta et al., 2017).

• 4. The variable style of slab transition zone interaction on Earth is governed by variations in trench mobility. Dynamic subduction models show that subduction accompanied by trench retreat is the most basic form of subduction (Kincaid and Olson, 1987). Trench motion is controlled by downgoing plate buoyancy and strength (Bellahsen et al., 2005; Capitanio et al., 2007; Tagawa et al., 2007; Ribe, 2010) and by upper-plate deformability and mobility (Garel et al., 2014; Holt et al., 2015). As a result, old strong slabs have a stronger tendency to retreat than young weak slabs (Agrusta et al., 2017). Indeed, most trenches on Earth retreat, and when integrated over the Cenozoic, all of the major trenches retreated, with the exception of the Himalayan front, which advanced in response to the continental collision. That downgoing plate buoyancy plays a key role is further borne out by the observation that significant trench motions are only observed for relatively old (>50 m.y.) subducting plates (Molnar and Atwater, 1978; Sdrolias and Müller, 2006).

• 5. The upper-mantle slabs imaged by Benioff seismicity and seismic tomography correspond to subduction of the past 15–60 m.y. (Fig. 2). Because thermal assimilation of slabs takes on the order of several hundred million years, this is a strong indication that earlier subducted slabs have sunk into the lower mantle. Indeed lower-mantle, high seismic-velocity anomalies can mostly be correlated with Mesozoic to early Cenozoic plate-motion histories (Ricard et al., 1993; Grand, 1994; Van der Meer et al., 2009).

• 6. Modeling shows that slabs once stagnated above a viscosity and phase transition take 100–300 m.y. to destabilize (Christensen, 1996). This implies that in natural subduction, lower-mantle penetration is generally triggered. Such triggering can be induced by changes in slab buoyancy (e.g., due to changing age or the arrival of continental lithosphere at the trench) or in resistance to upper-plate motion (Agrusta et al., 2017). Irrespective of whether lower-mantle sinking is triggered or evolves independently, models do not predict that the deeper sinking of previously flattened slabs is accompanied by large fluctuations in plate velocities (Christensen, 1996; Pysklywec and Ishii, 2005; Yanagisawa et al., 2010; Agrusta et al., 2017). Hence periods of much higher subduction velocities (e.g., India and Farallon before 50 m.y.) are likely the result of a different mechanism, e.g., deep subduction in conjunction with forcing by plumes (Becker and Faccenna, 2011; Cande and Stegman, 2011; van Hinsbergen et al., 2011).

• 7. Once slabs sink into the lower mantle, they will slowly lose their linear shape, either by folding and/or bulk thickening and rounding. This and the fact that slabs may end up as far as 1000–2000 km from the position where they originally subducted, makes associating lower-mantle slabs with past trenches more difficult with increasing depth in the mantle. This has indeed been noted in studies that have tried to correlate plate-motion history with tomographic anomalies in the lower mantle (Domeier et al., 2016).

• 8. These insights in slab transition zone interaction indicate that, at present, mixing between upper and lower mantle is at most delayed by a few tens of millions of years. On the time scale of 100–200 m.y., most subducted material is expected to sink into the lower mantle. It has been suggested that mixing across an endothermic phase boundary may have been less efficient in an earlier hotter Earth, due to the higher Rayleigh number (Christensen and Yuen, 1985; Davies, 1995). However, hotter plates would likely have been younger and thinner and thus less able to drive trench retreat, and hence the role of the transition zone in hotter Earth mixing may need re-evaluation (Agrusta et al., 2016).

Thanks to all those we have had discussions with over the years, which helped shape our ideas in this paper, in particular, Hana Čížková, Fabio Capitanio, Claudio Faccenna, Francesca Funiciello, Loic Fourel, Rhodri Davies, and Huw Davies. We also thank Scott King and an anonymous reviewer for their thoughtful comments on the manuscript. This work was funded by Natural Environment Research Council (NERC; grants NE/J008028/1 and NE/I024429/1) and made use of the facilities of N8 High Performance Computing provided and was funded under the embedded CSE program of the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). RA and JvH also acknowledge funding from the European Research Council (grant ERC StG 279828).