Research ArticleGEOPHYSICS

Deep slab seismicity limited by rate of deformation in the transition zone

See allHide authors and affiliations

Science Advances  27 May 2020:
Vol. 6, no. 22, eaaz7692
DOI: 10.1126/sciadv.aaz7692


Deep earthquakes within subducting tectonic plates (slabs) are enigmatic because they appear similar to shallow earthquakes but must occur by a different mechanism. Previous attempts to explain the depth distribution of deep earthquakes in terms of the temperature at which possible triggering mechanisms are viable, fail to explain the spatial variability in seismicity. In addition to thermal constraints, proposed failure mechanisms for deep earthquakes all require that sufficient strain accumulates in the slab at a relatively high stress. Here, I show that simulations of subduction with nonlinear rheology and compositionally dependent phase transitions exhibit strongly variable strain rates in space and time, which is similar to observed seismicity. Therefore, in addition to temperature, variations in strain rate may explain why there are large gaps in deep seismicity (low strain rate), and variable peaks in seismicity (bending regions), and, possibly, why there is an abrupt cessation of seismicity below 660 km.


Deep earthquakes occur within cold (<900° to 1000°C) slabs at pressures of 10 to 25 GPa (350 to 680 km). At these high pressures, brittle failure by frictional sliding or fracture, as occurs near the surface of Earth, is inhibited. Therefore, brittle failure at these high pressures requires a different mechanism to overcome either the large normal stresses (e.g., embrittlement through high pore fluid pressure) or a weakening mechanism. Therefore, much of the research on deep earthquakes has focused on the conditions at which sufficient fluids are present for embrittlement to occur (1, 2) or on the conditions at which weakening mechanisms can trigger failure. Other proposed failure mechanisms include thermal shear instability (35) and transformational faulting of metastable olivine (MO) (6, 7) or pyroxene (8), both of which trigger faulting through localized weakening. However, at both low and high pressure, the rock must be at conditions where strain energy accumulates and can be released through the failure processes. This study is focused on better understanding the physical conditions that determine where sufficient strain energy accumulates and how those conditions are related to the observed spatial distribution of seismicity.

Despite the differences in pressure-temperature conditions and triggering mechanisms, seismological observations of deep earthquakes suggests that these events are similar to shallow events in several ways. Deep earthquakes have (mostly) double-couple mechanisms indicating that they occur as a shear failure, aftershock sequences follow a standard Omori law for aftershock decay, energy/moment ratios are similar to shallow events, and long-range triggering of deep earthquakes has been observed [for comprehensive reviews, see (911)]. Analysis of source-time functions also shows remarkable similarity between shallow and deep earthquakes when depth-dependent differences in rigidity are taken into account (12).

However, compared with shallow events, deep earthquakes tend to have shorter rupture duration for a given earthquake size, exhibit more rupture complexity (i.e., multiple subevents), have a larger range in stress drop (1 to 100 s MPa) and rupture velocities [0.2 to 0.9 times shear velocity, with some earthquakes exhibiting supershear rupture velocity; see references in (11)], and exhibit depth-dependent aftershock productivity (e.g., almost absent at intermediate depths but present deeper than 550 km). In addition, some very deep events have very low radiation efficiency (<0.1), which has been interpreted to indicate melting during the rupture process (13). These seismic observations indicate that the rupture process for shallow and deep earthquakes is similar (i.e., a shear failure) despite different failure mechanisms and physical conditions (pressure, temperature, deviatoric stress), which affect the details of the rupture process and strain energy release.

Seismicity in the subducting lithosphere is often presented as the number of earthquakes per year versus depth for the world’s subduction zones (i.e., the global slab seismicity depth profile; Fig. 1A). This profile has been interpreted as indicating that there are likely two mechanisms for slab earthquakes, with dehydration embrittlement occurring at intermediate depths (50 to 300 km) where fluids are being released from the slab, and transformational faulting occurring deeper, where the slab is likely drier and metastable olivine and pyroxene may be present. There is growing seismological evidence for a MO wedge in the Japan slab, although it has not been definitively detected elsewhere due to a lack of dense seismic networks above other slabs (11). In addition, Gutenberg-Richter statistics (i.e., b values) also differ for seismicity above and below 350 km, supporting the hypothesis that there are two different mechanisms (9).

Fig. 1 Subduction zones exhibit strong variability of seismicity and strain rate in depth.

(A) Global seismicity in subduction zones for earthquakes greater than 100-km depth for Mw > 4 (blue) and Mw > 5 (gray) for the period 1964–2014 from the ISC-EHB catalog (52). (B to I) Regional seismicity rate and strain rate versus depth. Tonga (B), Kermadec (C), and Java-Sumatra (D) all have regional seismicity depth profiles that mimic the global profile in (A). The Kuriles (F), Japan (G), and Marianas (H) do not exhibit the characteristic increase in seismicity rate within the transition zone (400 to 660 km). Chile (E) and Peru (I) have a distinct lack of seismicity from 300 to 500 km but do have earthquakes in the transition zone from 500 to 660 km.

This explanation for the depth distribution of seismicity, based on a depth-controlled failure mechanism, needs to be reevaluated in light of new observations. First, there is growing evidence that the deep slab may not be dry and fluids can be transported and released from the slab in the deep transition zone (14). Second, regional differences in b values suggest that rupture may occur by a combination of mechanisms including transformational faulting and thermal shear instability for deep earthquakes (15). Third, recent studies suggest that shear instability is also a viable mechanism for intermediate-depth earthquakes (4, 5, 16).

All three proposed mechanisms for triggering deep earthquakes are, to first order, controlled by the requirement that the slab remains sufficiently cold to the base of the transition zone. The thermal structure of slabs is primarily controlled by the age of the slab at the time of subduction (tsub; older slabs are colder and thicker) and the rate at which slabs sink into the mantle (Vs; slower subduction allows more time for the slab to heat up). These two variables have been combined to define the thermal parameter, ϕ = Vstsub, and compared with the maximum depth of seismicity in the slabs, zmax. This comparison has been shown to be consistent with the predicted depth extent of an MO wedge in kinematic thermal models (17) and in dynamic models (18), although thermal models accounting for the variability of thermal conductivity in the slab (19) or water effects on reaction kinetics (20) suggest the MO wedge could be substantially shorter.

However, perhaps more problematic for any possible failure mechanism that is thermally limited is the observation that many slabs exhibit large gaps in seismicity below 410 km (see figs. S2 and S3). Most notably, the Peru and Chile slabs have earthquakes at 550 to 650 km, but no deep earthquakes from 410 to 550 km, and the same is true for sections of the Java-Sumatra slab. Similarly, the shallow-dipping section of the Japan slab (between profiles 3 and 4 in fig. S2Da) has a large aseismic region (~150 km across from 200- to 660-km depth).

These observations appear to be inconsistent with any mechanism for deep earthquakes that is primarily controlled by temperature. This is because temperature contours in the slab are elongated concentric surfaces, and therefore, the critical temperature for the failure mechanism will be present at all depths up to the maximum depth. If earthquakes are occurring at 660 km because the slab is cold enough for the failure mechanism to operate at that depth, then there is a portion of the slab at all depths above 660 km that is also cold enough for earthquakes to occur. Large gaps in deep earthquake seismicity indicate that being cold enough for the deep earthquake failure mechanisms to be viable, while necessary, is not a sufficient condition to explain the distribution of deep earthquakes. Therefore, some other physical factor, in addition to temperature, must control the depth distribution of deep earthquakes.

An alternative explanation for the peak in seismicity in the transition zone is that it is due to higher stresses in the slab at this depth caused by the viscous resistance in the lower mantle (21), buoyancy forces, or stresses related to volume contraction associated with both equilibrium or metastable phase transitions [see (22) and references therein]. These studies demonstrate that the combination of available forces and slab rheology predicts high stress magnitudes (>500 MPa). However, they are instantaneous calculations, which cannot show if the resulting deformation of the slab is consistent with observations or occurs at high enough strain rates. They also do not address the spatial variability in seismic strain rate.

At shallow depths (<100 km) in the crust and lithosphere, the interior of tectonic plates is aseismic, while seismicity occurs at plate boundaries where deformation is localized and strain rates are high. Localized deformation at plate boundaries occurs through a feedback between tectonic forces and the rock rheology [e.g., (23)]. High strain rate is also known to be a factor affecting failure by thermal shear instability (4) and transformational faulting (6). This suggests that the discontinuous distribution of seismicity in slabs may also be determined by feedbacks between the rheology and forces acting on the slab, which leads to a discontinuous distribution of high strain rate regions in the cold slab.

Early consideration of the causes of deep earthquakes recognized that slab rheology is also an important factor in determining where deep earthquakes occur. Wortel (24) considered the requirement that the rheology of slabs must allow for the accumulation of stresses that are released in the process of an earthquake and used this to determine a critical temperature for deep seismicity of <900° to 1000°C. Above this temperature, the slab strength (viscosity) is too low and the stresses will be accommodated by viscous flow. This temperature is consistent with more recent estimates for the maximum temperature at which slabs deform through low-temperature plasticity or yielding (25). Brodholt and Stein (26) later argued that slabs are rheologically strong enough to support deep earthquakes beyond 660 km but assumed uniformly low strain rate of 10−18 s−1 in the slab. At this low strain rate, the slab would be essentially rigid and would sink through the mantle without deforming internally, which is inconsistent with the occurrence of earthquakes and the deformed shapes of slabs inferred from seismicity (27, 28). Therefore, both the occurrence of earthquakes in slabs and the geometry of slabs require that slabs support relatively high stresses but are able to deform internally. This is an important constraint not only for the generation of deep earthquakes but also for the rheology used in long-term subduction models.

The requirement that slabs support high stress and deform internally is met by dynamic models of subduction that use a rheology with a strong temperature dependence from a composite diffusion-dislocation creep viscosity for olivine and either yielding [e.g., (29, 30)] or some other approximation (31) of low-temperature plasticity (25). Here, I show that the strain rate pattern within deforming slabs in these models varies spatially and temporally and mimics the strong variability in seismicity within the world’s subduction zones. I use this correspondence to argue that observed seismic strain rate, while only a fraction of the total strain rate, reflects the actual variations in strain rate in slabs. That is, in addition to the thermal constraints on the possible failure mechanisms, earthquakes occur where the strain rate is high, and the gaps in seismicity within cold slabs reflect regions of the slab that are deforming more slowly.

This is not an entirely new concept, Tao and O’Connell (32) showed that the peak in seismicity in the transition zone could be explained by the high strain rate in a weak slab (same viscosity as the mantle) due to a jump in viscosity into the lower mantle. However, such a weak slab could not support the stresses required for earthquake generation nor explain the magnitude or orientation of stress in the slab (21). More recently, others have shown that earthquakes preferentially align in regions of high curvature [e.g., (27)] and that weakening the slab leads to strain rates in the deep slab that are sufficiently high to drive thermal shear instability (33). Here, I revisit this largely overlooked factor in the process of deep earthquake generation to explicitly argue that the seismicity distribution directly reflects the spatial variation in strain rate within strong but deforming slabs. Note that this is different from arguing that seismicity occurs where the stresses are higher [e.g., (21, 22)]: In the models presented, the stress within the cold slab is high everywhere because it is deforming at the yield stress (1 GPa), but a reduction in viscosity through plastic yielding allows the slab to deform at a higher strain rate locally.


The hypothesis presented here was first motivated by examining the strain rate evolution in two-dimensional (2D) dynamic models of subduction and recognizing the similarity to the seismicity pattern for intermediate to deep earthquakes. Therefore, I first present observed seismic strain rate estimated from the moment release rate in the slab (see Methods) and then compare this to the strain rate distribution in slabs from simulations (see Methods and in Supplementary Materials).

Seismic strain rate depth profiles

Although much consideration of the mechanisms for deep earthquakes have been motivated by the distinctive global seismicity depth profile (Fig. 1A), there is considerable variability both between and within individual subduction zones. Figure 1 (B to I) shows the regional seismicity and strain rate profiles for eight subduction zones. The strain rate curves largely follow the seismicity pattern, except in regions with large events relative to the regional average, which exhibit strain rate peaks (e.g., Kuriles at 600 km). There are three types of regional depth profiles. First, Tonga, Kermadec, and Java-Sumatra are all similar to the global profiles with seismicity present at all depths and a seismicity/strain rate peak in the transition zone. Second, both Chile and Peru are distinctive due to the large gap in seismicity from 300 to 500 km (except for a few events in Chile). Last, the Kuriles, Japan, and Marianas lack the peak in the transition zone. In Japan, there are only a few earthquakes deeper than 600 km, and the Marianas appears to have two peaks centered at 450 and 600 km.

The variability in seismicity and strain rate observed between subduction zones is also present within individual subduction zones. Figure S2 (A to H) shows profiles spaced 100 or 200 km apart for all eight subduction zones. Figure 2 shows the seismicity and strain rate profiles for the central section of the Tonga-Kermadec slab from 32°S to 21°S. Even in this central region of an old (85 to 100 Ma) and cold slab, far from the effects of slab edges, the seismicity pattern and strain rate vary substantially from one profile to the next. Some of the profiles have the characteristic peak in seismicity within the transition zone, but the depth and width of the peak vary; still, other profiles do not have this peak at all. In adjacent profiles, the depth and width of the peak change continuously, suggesting that the processes or conditions that determine the location of seismicity are also changing over length scales of less than 100 km along-strike.

Fig. 2 Seismicity and strain rate exhibit continuous changes in peak depth and width along adjacent profiles.

Examples from the Tonga (1 to 5) and Kermadec (7 to 11) slabs. From profiles 1 through 5, there is a continuous change in the depth and width of the transition zone peak and the emergence of a narrow strain rate peak at 400-km depth. From profiles 11 to 7, the transition zone peak disappears, being replaced by two smaller peaks, and then a shift to a broad peak centered at 400-km depth. Locations of profiles are shown in fig. S2 (A and B).

The global average, while illustrative, does not capture the observed variability in seismicity and strain rate. More importantly, ignoring this variability is ignoring important information about the conditions at which deep earthquakes occur and the physical requirements for any failure mechanism. If strain rate is a key factor controlling the variable distribution of seismicity in slabs, then the seismic strain rate should reflect the actual strain rate in the slab (just as seismicity reflects strain rate within the plates at the surface). In this case, multiple failure mechanisms for deep earthquakes are possible, but all require sufficient strain accumulation, and therefore, earthquakes are limited to cold and high strain rate regions. Alternatively, if the strain rate is uniform within the slab, then the strain rate is not an important factor, and in this case, some other factor controlling appropriate conditions for the various failure mechanisms will determine the seismicity pattern. Because it is not possible to directly measure the strain rate in slabs, I use numerical simulations to estimate the strain rate distribution.

2D dynamic models of subduction

The strain rate (dε/dt) within the slab depends both on the rheology used in the models and the effect of phase transitions on the time-dependent evolution of the slab. For the strongly temperature-dependent viscosity of olivine, the slab rheology is primarily determined by the yield strength and maximum viscosity allowed in the models. This means that the slab interior is stiff and resists internal deformation. Therefore, the stresses caused by sinking of the slab are primarily accommodated by flow in the surrounding mantle. However, during episodes of slab folding, yielding within the slab leads to regions of localized internal deformation with higher strain rate.

Model 1 shows the evolution of an 80-Ma-old slab in which no phase changes are included (Fig. 3A and movie S1). At the start of subduction, the stress due to negative slab buoyancy is small, and the main area of deformation is the bending region at the outer rise, which exhibits a classic hourglass-shaped yielding region (Fig. 3A, a). As the slab lengthens, it sinks rapidly through the upper mantle and experiences high strain rates throughout; the slab sinks faster than the trailing plate and stretches due to low viscous resistance from the surrounding mantle (Fig. 3A, b). Once the slab reaches the viscosity jump at the top of the lower mantle, it slows down as the buoyancy of the slab becomes partially supported by higher viscosity in the lower mantle. At the same time, the internal strain rate decreases markedly, and the effective viscosity increases (Fig. 3A, c). Subsequently, the slab dip shallows slightly, but there is otherwise little change to the slab shape or sinking rate for the remainder of the simulation (Fig. 3A, d).

Fig. 3 In simulations of subduction, phase transitions cause folding and buckling of the slab, leading to localized peaks of high strain rate that change in time.

Cross sections of subducted slab with strain rate (color) within the 1000°C contour. Red bars indicate the pattern of shortening (compression) directions. (A) Model 1 with no phase transitions shows little strain rate variation in the slab. (B) Model 2 with phase transitions has strong variations in strain rate both spatially and temporally as the slab changes shape. (B, e to f) Depth profiles of the maximum strain rate occurring below a specified temperature in the slab interior for model 2 at the times shown in (B, a to d). Colder temperature limits show the strain rate at the interior of the slab. These strain rate depth profiles are similar to the observed strain rate profiles calculated from slab seismicity. (C) Model 3, a younger slab (40 Ma), is a warmer and therefore weaker slab and deforms at higher strain rates but shows similar folding and buckling behavior to model 2.

In model 2, the evolution of the slab is quite different, owing to the effect of the phase transitions (Fig. 3B and movie S2). The evolution of the slab is similar to model 1 until the slab starts to interact with the more viscous lower mantle. At this point, there is a strongly time-dependent strain rate pattern associated with the bending and buckling of the slab. As shown previously, the density anomalies associated with the phase transitions cause folding and buckling of the slab as well as forward and retrograde motion of the trench (30, 34). High strain rate regions form in regions of bending with an hourglass pattern characteristic of bending with a neutral plane. In addition, there is a region of high strain rate at 550- to 650-km depth, which occurs with or without slab bending. This high–strain rate region occurs between the garnet-to-ilmenite (g-i; elevated) and garnet-to-bridgmanite (depressed) transitions in the harzburgite layer.

For model 2, the maximum strain rate occurring below a temperature limit (e.g., 1000°C) is plotted as a function of depth (Fig. 3B, e to h). The depth profiles show high strain rates at shallow depths (<100 km) corresponding to the high rates of deformation along the slab surface and in the bending region at the trench. Below this depth, the strain rate magnitude generally decreases with a minimum near 300 km. In the transition zone, the shape of the profile is strongly time variable as the slab folds and buckles. During sizeable bending events, the maximum strain rate does not depend strongly on temperature, while at other times, the maximum strain rate is higher at temperatures of 900° to 1000°C. The depth and width of strain rate peaks, as well as the number of strain rate peaks, change in time. Note also that the strain rate drops sharply, crossing into the higher-viscosity lower mantle. These strain rate profiles have similar characteristics to the strain rate profiles calculated from observed seismicity. It is these similarities in strain rate between the long-term subduction models and the observations that argue in favor of strain rate as an important environmental variable determining the distribution of deep earthquakes.

Model 3 is the same as Model 2 except that the initial subducting plate age is younger (40 Ma). Here, the evolution of the slab and strain rate pattern is similar to model 2, exhibiting peaks in strain rate during bending and folding (Fig. 3C and movie S3). However, because the slab is younger, it is also warmer with the 700° to 800°C isotherms restricted to depths less than 400 km and periods of time when the 1000C contour does not extend past 660 km or becomes broken. Also, because the slab is warmer, it has a smaller integrated strength, and it deforms at higher strain rates overall. Note that in the deeper, warmer regions of the slab, the stress levels are lower because the slab is not deforming at the yield stress, and therefore, despite the higher strain rates, seismicity might not be expected to occur in these regions.

In all three models, the stress orientations indicate that the slab exhibits down-dip compression (DDC) along the top/central part of the slab (~200 to 600 km), except in regions of folding, consistent with earlier studies (21, 35). In the folding regions, the location of the DDC shifts to the high–strain rate region on the underside of the slab, and the stretching orientation is parallel to the folded slab surface. This suggests that the envelope of seismicity defining the location of deep slabs may shift from the top of the slab to the bottom of the slab, depending on the orientation of folding at that depth.

Comparison of observations and models

The models show a dynamic view of continuously changing slab shape and peaks in strain rate. For Earth, however, we have only a single snapshot in time. One way around this is to recognize that for 3D slabs, their shape evolves in time and in space so that adjacent profiles capture the time evolution of the changing shape. This is a commonly used space for time substitution from structural geology. Using this approach, comparison shows that both the observations and the models exhibit peaks in strain rate of variable magnitude, depth, and width. In the models, peaks follow regions of bending in time (Fig. 4), while in observations, there are many examples of the depth and shape of peaks systematically changing along-strike (Fig. 2 and fig. S3) and of seismicity following lines of high curvature along-strike (27). In the models, there are times when strain rate is low throughout much of the slab, just as there are profiles in the observations that have very low seismicity adjacent to regions with higher seismicity (e.g., Japan).

Fig. 4 Time evolution of the maximum strain rate in the slab for model 2 demonstrates spatial and temporal variability.

The color is the maximum strain rate occurring in the slab at temperatures less than 900°C as a function of depth and simulation time. The depth and width of strain rate peaks migrate in time following bending regions and folds. There is also a peak centered at 600-km depth that occurs at three times, independent of a major fold in the slab (see Fig. 5).

In addition to the time-variable evolution of the strain rate profiles, there are two more general characteristics of the strain rate profiles to explain. First, note that in the models, the strain rate decreases to a minimum value beneath 660 km, for all times, regardless of the shape of the slab at this depth (Figs. 3 and 4). This drop in strain rate within the slab is caused by the increased viscous support provided by the higher viscosity of the lower mantle. This higher viscosity slows down the overall rate of deformation of the slab (by a factor of 100) and the surrounding mantle. The increase in viscous support by the surrounding mantle also means that much less stress must be supported by the slab itself. And, the slab is no longer bending and buckling but rather sinks passively (<1 cm/year) into the lower mantle. This suggests that the cessation of seismicity at 660-km depth could be controlled by the change in rheology, which causes an overall lower strain rate in the slab and surrounding mantle. However, it is also true that transformational faulting from MO also shuts off beyond 660 km because the transformation to bridgmanite plus ferropericlase is endothermic (36).

Second, there is a peak in strain rate just below 600 km that persists even when there is no appreciable bending of the slab at this depth (Fig. 4). Figure 5 shows that this peak occurs at the garnet-to-bridgmanite transition in the harzburgite layer, which lies between the elevated garnet-ilmenite (g-i) phase transition and the depressed ilmenite-to-bridgmanite (i-B) and ringwoodite-to-bridgmanite plus ferropericlase (R-B + f) phase transitions. The region between these two phase transitions is subject to net compression caused by the negative buoyancy above g-i and positive buoyancy below i-B and R-B + f. This compression causes a higher localized stress, which due to yielding leads to a reduction in viscosity and higher strain rate. This result suggests that accurately accounting for the phase transitions in both the pyroxene and olivine components of the slab, as is done in these simulations, is essential for fully accounting for the buoyancy forces affecting the overall deformation of the slab and causing deep slab seismicity.

Fig. 5 Strain rate peak below 600-km depth is associated with garnet-to-bridgmanite transition in the harzburgite layer.

(A) Zoom-in on the strain-rate in the slab for model 2 at 13.9 Ma (Fig. 3E) showing the location of the phase transitions (gray), and temperature contours (black). (B) shows the profiles of the strain rate peak at 1000°C (black) and 900°C (blue). This strain rate peak occurs for periods of 5 to 10 Ma (see Fig. 4) in the absence of appreciable folding or bending of the slab. gt, garnet; il, ilmenite; rw, Ringwoodite; fp, ferropericlase; brg, bridgmanite.

Because the numerical simulations are fully dynamic, the evolution of subduction rates and trench motion, and thus slab shape, is determined by the time-evolving balance of forces. Therefore, the models do not correspond to any particular subduction zone on Earth. However, because the time steps in the models correspond to thousands of years, and the strain rates are determined by the instantaneous balance of forces and the rheology, snapshots from the models with similar geometry can be compared with observed profiles.

First is a comparison between a model snapshot and a profile from Chile (Fig. 6A). The Chile slab has a fairly planar shape below the shallow flat-slab segment (28) and is thought to be sinking directly into the lower mantle (37), similar to the model snapshot. Both the model and observations exhibit low strain rates above the transition zone with a peak around 600 km associated with the phase transitions at this depth. Similar profiles are also seen in Peru and portions of Java-Sumatra (see fig. S3, F and G).

Fig. 6 Comparisons of model snapshots and observations show similarities between slab shape and strain rate distribution.

Comparison for (A) the Chilean slab, (B) the Marianas slab, and (C) the Bonin slab. (a) Snapshots from model 2 at times indicated. Colors and contours are the same as in Fig. 3. (b) Maximum strain rate profiles. Line colors are the same as in Fig. 3. (c) Earthquake histogram and strain rate for profile 7 in Chile, profile 7 in the Marianas, and profile 9 in Bonin. (d) Cross section (depth versus distance) of earthquakes for the same profile as the histograms. Colors indicate depth and are the same as used in fig. S2.

Second is a comparison for the Mariana slab (Fig. 6B). The model snapshot shows an overturned slab with high strain rate regions between 350 and 500 km and a second peak associated with the g-i transitions near 600 km. While typical profiles of slabs do not exhibit this overturned shape, sections of the Mariana and Java-Sumatra slabs are clearly overturned (see fig. S2). The distance from the westernmost limit of the slab and the slab tip at 700 km is about 250 km in the model and about 200 km in the observations. The high strain rate region in the model at 350- to 500-km depth corresponds to a cluster of events in the Marianas slab with no seismicity above or below this bending region (except for one very deep event). However, note that the high stain rate region at 350- to 500-km depth in the model occurs on the bottom of the plate with DDC, while the deeper strain rate peak occurs across the slab width. If the location of deep earthquakes is controlled by strain rate, then these results require that some events initiate in the lithospheric portion of the plate, rather than the crust or harzburgite layers: an important observation for determining which failure mechanisms may be active in different locations in the slab.

Third is a comparison between a model snapshot and a profile from the Bonin slab (Fig. 6C). The profile is from the region of the slab between the shallow-dipping, planar slab to the north in Japan and the steeply dipping, curved slab to the south in the Marianas (see map in fig. S2E). Note that the shape of the profile seen here exists over 200 to 250 km along-strike. The seismicity indicates that the slab flattens to the west just below 500-km depth [the 1982 moment magnitude (Mw) 6.7 outboard event is noted]. However, there is also an event at 680-km depth behind (east) of the shallower seismicity and separated by a gap. This deep earthquake is the 2015 Mw 7.8 Ogasawara (Bonin) Islands event (38). One possibility to explain the change in geometry across such a narrow region is a tear in the slab (39). A second possibility is that the slab is folded (38), as shown in the model snapshot. Note that in this snapshot, there is a region of high strain rate in the outboard, top portion of the fold, but the strain rate in the bottom of the fold is low. The fold in the model snapshot is also broader than the fold needed to explain the earthquakes, indicating that more intense weakening of the slab is required during folding.


The numerical models show that high strain rate regions occur where the slab is bending or folding, and locally between phase transitions with opposite Clapeyron slopes. Both the modeled strain rate and observed strain rate profiles exhibit the following: (i) peaks of variable magnitude, depth, and width; (ii) regions of very low strain rate (gaps in seismicity); and (iii) a sharp drop-off in strain rate below 660 km. The similarity in the strain rate profiles from the models and from the seismicity shows that seismic strain rate directly reflects the actual strain rate of the slab. Therefore, the spatial variation in deep earthquake seismicity is determined by the spatial pattern in strain rate within strong, deforming slabs.

This conclusion requires that the slab rheology is such that the slab is viscously strong (strong temperature-dependent viscosity) and that it can yield in response to localized higher stresses. If a lower yield stress (or lower maximum viscosity) were used, the whole slab would deform at a higher rate. However, in these models, a lower yield stress causes the slab to break off because of the added negative buoyancy from the phase transitions (34) compared with other models that use a lower yield stress [e.g., (40)]. Also, a lower yield stress may not be consistent with the large stress drops (up to 100s of MPa) estimated for some deep earthquakes (9, 41) or the differential stress at which low-temperature plasticity would occur in cold slabs (25).

The yield stress, or use of a strong power law exponent (31), is an approximation to low-temperature plasticity (i.e., Peierls creep). At shallow depth, the yield stress is also used to approximate brittle failure through frictional processes (i.e., Byerlee’s law). A better approximation of the low-temperature plasticity may be achieved using a temperature-dependent power law exponent (42) and would likely lead to overall higher strain rates in the cold interior of the slab, but would still exhibit peaks in strain rate in regions of bending or folding. Similarly, including the effects of elasticity with a viscoelastic rheology would also decrease the magnitude of stress in the slab and can result in a lower apparent slab viscosity (43). It is also important to note that slab temperature remains an important factor for deep seismicity. Young and slowly sinking slabs will not have deep seismicity because they are too warm and therefore deform through diffusion and dislocation creep. In this case, the stress in the slab is relaxed viscously.

The idea that strain rate is an important factor in understanding deep earthquakes is not new, but it has been largely ignored or forgotten in the literature relating potential failure mechanisms to the physical state of the slab. The mechanism of shear instability explicitly relies on having a high enough strain rate (and stress) to cause shear heating (3); however, localized grain size reduction is likely necessary to reach sufficient strain rates for this mechanism to be viable (5). Transformational faulting of olivine also requires a sufficient strain rate, and the strain rate affects the window of temperatures at which this mechanism occurs in the laboratory [see Figure 4 in (6)]. However, while many previous models explore the state of stress in the slab [e.g., (21, 22, 44)] and evolution of an MO wedge, they do not also assess the strain rate requirements [e.g., (18)]. Last, the correspondence of high strain rate with regions of bending and buckling in the models also agrees with the observation that earthquakes appear to align with regions of high slab curvature (27) and the higher rates of seismicity in strongly deformed slabs [e.g., Tonga; (45)].

While several studies have used seismic strain rate as a minimum constraint on the slab strain rate (and maximum viscosity) required in dynamic models of subduction [e.g., (40, 46)], these models have not considered the spatial variability in strain rate and how this may be related to the rheology of the slab. In contrast, many past studies of subduction have used simplified rheologies and inherently weak slabs [see (32) and review in (47)]. While these weak slab models may meet the minimum average strain rate requirement, they are inconsistent with laboratory constraints on rheology, the requirement that slabs be strong enough to store stresses that are released seismically (not through viscous flow), and the large observed stress drops for some deep earthquakes.

While the models presented here provide compelling evidence that strain rate is an important factor, together with the thermal structure, in determining the distribution of deep earthquakes, there are limitations that must be explored in future models, including low-temperature plasticity, elasticity, tectonic overpressure from phase transitions, inclusion of an MO wedge, compressibility, and 2D geometry (this list is addressed in the Supplementary Materials). Systematically addressing these limitations will surely affect the details of the slab deformation, and the magnitudes of the strain rates and stresses, but will unlikely affect the primary conclusions presented here because the slabs will still bend and buckle with spatially and temporally variable strain rate. To further link the dynamical models to possible failure mechanism, running identical subduction models, both with and without MO, for a range of thermal parameters would facilitate analyzing how stress magnitudes and orientations are related to portions of the slab that have appropriate temperature and strain rate conditions for the possible deep earthquake failure mechanisms.

Incorporation of better approximation of low-temperature plasticity, followed by elasticity and compressibility, should be primary targets for further study because this will allow more direct comparison between the strain rate and stress states in the models and the conditions required for each deep earthquake failure mechanism. It is important to note that Farrington et al. (43) show that the mode of subduction and slab morphology are not affected by inclusion of elasticity in 3D models of free subduction; therefore, the distribution of high–strain rate regions associated with bending in the slab would also not be affected. However, the location of the maximum stress (and its magnitude) and the stress orientations within the bending region is shifted because the stress and strain rate are not generally coaxial in viscoelastic materials. It is for this reason that the stress orientations from the slab models have not been analyzed in detail or compared directly with moment tensor solutions.

Progress in understanding the triggering mechanisms for deep earthquakes has been slowed by an insufficient physical framework to adequately demonstrate or refute the viability of proposed failure mechanisms. Taking into account the added constraint of strain rate should help to resolve which of these mechanisms are active in the subducting lithosphere, with the possibility that multiple mechanisms may be required to explain the variability within and between slabs of varying age and state of deformation. Considering the rate of deformation and the process through which strain accumulates in the slab more explicitly may also foster reevaluation of seismic observations related to rupture behavior (moment release, rupture time, b values, aftershock occurrence, and stress drop), and how these are related to different triggering mechanisms. For example, considering the detailed orientations of stress and strain in bending regions in conjunction with a physical model of the failure mechanism may help to explain the preference for near-horizontal failure planes for deep earthquakes (48) and the large depths (near 660 km) for the largest deep earthquakes (11). In addition, it is possible that the observed correlation of b values with thermal parameter (15) may also be related to the strain rate of the slab. This is because (i) faster sinking slabs are expected to undergo more internal deformation to accommodate the viscous resistance to sinking into the lower mantle, and (ii) the correlation between b value and thermal parameter is primarily controlled by the sinking rate. This is why, for example, Tonga has a higher thermal parameter than Japan, even though the subducting plate in Tonga is younger than that in Japan. Last, the spatially variable strain rate can be used to further constrain the appropriate rheology for the lithosphere by providing a direct link between the short time scale phenomena of strain release through earthquakes and the long-term deformation of slabs.


The spatial distribution of deep earthquakes exhibits marked variability between different subduction zones and along-strike within individual subduction zones. Existing explanations for deep earthquakes assume that the distribution of seismicity is primarily determined by whether the appropriate thermal conditions exist in the slab for a variety of special triggering mechanisms. Here, I have shown that the strain rate distribution from simulations of strong but deforming slabs has the same variability as the observed strain rate: peaks in strain rate at variable depths, regions of low strain rate, and a sharp drop-off in strain rate at 660 km. The results presented here cannot distinguish between possible failure mechanisms for deep earthquakes. However, they do suggest a new approach for testing these mechanisms that would combine the thermal and strain rate constraints with appropriate rheological models. With these models, it may be possible to more directly link the required conditions (e.g., pressure, temperature, and strain rate) for viable failure mechanisms with seismic observations of deep earthquakes (e.g., stress drop, radiation efficiency, and b values) and to better constrain the rheology of the lithosphere and mantle.


I compared the time-dependent evolution of deformation within subducting lithosphere to the seismically accommodated strain rate observed in present-day slabs. For the numerical simulations, I used the second invariant of the strain rate tensor to quantify the magnitude of the strain rate. The seismic strain rate was calculated following the analysis of Bevis (49). While these are both measuring deformation within the slab, the seismic strain rate is, by definition, the deformation that is not accommodated by viscous deformation. Therefore, the spatial pattern (depth dependence) of the strain rates was compared, but not the magnitudes.

Dynamic subduction models

The subduction models are fully described in Billen and Arredondo (30). The time-dependent evolution of the sinking lithosphere (slab) was modeled in a 2D slice of a spherical shell extending from the surface to the core-mantle boundary and 61° in longitude (fig. S1). Simulations were run using the CitcomS finite element code (50). CitcomS solves the conservation equations for mass, momentum, and energy using the extended Boussinesq approximation, which assumes incompressibility, but includes an initial adiabatic gradient, shear heating, and latent heat from phase transitions (51). The model setup allows for fully dynamic simulations in which only buoyancy forces drive subduction, plate, and trench motions (free subduction): All boundaries have a zero normal velocity and no tangential stress (free slip). To allow the plates to move freely toward or away from the sidewalls, we imposed a boxed region at the trailing end of both plates that has a fixed thermal profile and low viscosity. Subduction was initiated with a proto slab extending to a depth of 200 km.

Key features of the model setup include (fig. S1) the following: (i) a layered compositional density structure for the subducting and overriding plates, (ii) a composite viscoplastic rheology based on laboratory experiments for olivine, and (iii) compositionally dependent phase transitions. In addition, the basaltic crustal layer is modeled as a weak layer (maximum viscosity of 1020 Pa s), which allows the subducting plate to slide past the overriding plate. The maximum viscosity reverts to the global maximum of 1024 Pa s as the basalt transitions to eclogite (at 80 to 100 km). All of the parameters for the composite viscosity and phase transitions are documented in Billen and Arredondo (30).

Seismic strain rate calculation

The strain rate associated with seismicity in the subducted lithosphere was calculated following Bevis (49). This calculation relates the moment released within in a volume of the slab to the down-dip strain rateϵ̇s=TMo2μVT(1)where ∑TMo is the total moment released within a volume, V, during a time period, T, and μ is the rigidity. For this calculation, a volume of slab material, V = WHL is considered, where W is trench parallel width, H is slab thickness, and L is the down-dip slab length. Equation 1 results from considering that the slip in each event, D, is related to the seismic moment by Mo = μAD, where the fault area is A. The average slip accumulated within a time period, T, isDˆ=TMo/μWH(2)where the average fault area is taken as A = WH. Assuming this slip accommodates down-dip deformation, the change in down-dip length isδ=TMo/2μWH(3)

The down-dip strain is then given by ϵ = δ/L and the strain rate by ϵ̇=δ/LT.

The strain rate is calculated along evenly spaced trench-perpendicular profiles located every 200 km (W) along the trench, in 10-km-depth intervals (dz), assuming a constant seismogenic thickness of 80 km (H). For Tonga, I used W = 100 km because the seismicity rate is much higher than in other subduction zones. A value of 60 GPa was used for the rigidity (49). In the Preliminary Reference Earth Model (PREM), the rigidity increases by a factor of 2 with depth in the upper mantle. However, at the same time, the seismogenic width is also expected to decrease with depth as the slab warms. Therefore, these changes will largely cancel out. In addition, the changes in strain rate of interest vary by 10× to 100×. Therefore, for simplicity, I used a constant value for rigidity and seismogenic thickness. The down-dip length of the slab within each depth bin depends on the average dip of the slab within each depth bin, α, as L = dz/ sin (α). Using the Slab 2.0 geometry model (28), the average slab dip was calculated from all points in the depth bin and located within 100 km (0.5 W) of the profile.

Earthquake data for a 50-year time period (1964 to 2014) are downloaded from the ISC-EHB (International Seismological Centre-Engdahl-vanderHilst-Buland) catalog for depths of 100 to 700 km and magnitudes of 4.0 and greater (52). A map and the earthquake profiles for each region are shown in figs. S2 and S3. For this time period, all the earthquakes have been relocated using the EHB algorithm, which provides better earthquake locations and, in particular, better depth estimates than previous ISC determinations (53). Moment magnitudes are used when available; otherwise, body-wave magnitude is converted to moment magnitude using the relationship Mw = 0.85mb + 1.03 (54). The moment for each event is then given by Mo = 101.5(10.7 + Mw) (note that this gives the moment in dyne-cm; 1 N/m is 107 dyne-cm). Last, for each subduction zone, the regional strain rate as a function depth is determined by summing the profiles. Plots of the seismicity rate and strain rate as a function of depth along each of the profiles are included in the Supplementary Materials (fig. S3, A to H).

For any study using seismicity to constrain rates of deformation, one must be aware that the relatively short (50 years) duration over which data are available may not capture seismicity that occurs on a longer time scale. In particular, since the rate of seismicity in deep slabs is quite low compared with rates along plate boundaries at the surface, apparent gaps in seismicity may be filled in by longer observation times. Similarly, an isolated, but rare, large events can result in an apparent spike in strain rate. For example, the Mw 8.3 Okhotsk event in the Kuriles is the largest deep earthquake recorded and appears as a strain rate spike at 600-km depth in Fig. 1F and fig. S3C (profile 8).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: I would like to thank two anonymous reviewers and Z. Zhan for the constructive comments. All earthquake data are publicly available at the websites cited. Figures were made with the GMT (Generic Mapping Tool) plotting software. The Computational Initiative for Geodynamics (CIG) supports the CitcomS software used for the numerical simulations. Funding: This project was supported by a fellowship of the Alexander von Humbolt Foundation. The numerical simulations research was funded by the NSF, award #1246864. Author contribution: M.I.B. conceived of this study using the existing simulations from Billen and Arredondo, 2018 (30). M.I.B. analyzed the earthquake data and model simulations, and wrote the manuscript. Competing interests: The author declares that she has no competing interests. Data and materials availability: An archive for the numerical models is included with previous publications cited. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the author.

Stay Connected to Science Advances

Navigate This Article