Abhishek Prakash,Caleb W. Holyoke,Peter B. Kelemen,Stephen H. Kirby,Andreas K. Kronenberg,William M. Lamb
Carbonates and intermediate-depth seismicity: Stable and unstable shear in altered subducting plates and overlying mantle
intermediate-depth earthquakes,double Wadati–Benioff zones,thermal shear instability,magnesite in mantle,carbonates
Significance Earthquakes at subduction zones extend to great depths where slip on faults cannot occur by conventional frictional processes. Earthquakes at depths greater than 15 km are commonly attributed to high pore pressures, low effective pressures, and dehydration embrittlement of serpentine and other hydrous silicates. However, recent evaluations of temperatures and pressures of intermediate-depth earthquakes do not always correspond to breakdown reactions of antigorite serpentine, amphiboles, or talc. Altered mantle rocks of the upper lithosphere and mantle wedges also contain magnesian carbonates, which are stable to great depths, pressures, and temperatures. Carbonates of downgoing slabs are weaker than anhydrous mantle silicates. Our model suggests that intermediate-depth earthquakes may occur by strain localization and thermal shear instabilities within carbonate horizons.
Carbonates and intermediate-depth seismicity: Stable and unstable shear in altered subducting plates and overlying mantle

Earthquakes at subduction zones extend to great depths where slip on faults cannot occur by conventional frictional processes. Earthquakes at depths greater than 15 km are commonly attributed to high pore pressures, low effective pressures, and dehydration embrittlement of serpentine and other hydrous silicates. However, recent evaluations of temperatures and pressures of intermediate-depth earthquakes do not always correspond to breakdown reactions of antigorite serpentine, amphiboles, or talc. Altered mantle rocks of the upper lithosphere and mantle wedges also contain magnesian carbonates, which are stable to great depths, pressures, and temperatures. Carbonates of downgoing slabs are weaker than anhydrous mantle silicates. Our model suggests that intermediate-depth earthquakes may occur by strain localization and thermal shear instabilities within carbonate horizons.

This paper introduces a mechanism for earthquakes, from ~50 to 400 km depth (Fig. 1) in subduction zones, roughly corresponding to intermediate-depth events (between 50 and 300 km), which are not readily explained by conventional brittle mechanisms or phase changes under deviatoric stress (1–3). At shallower depths, rapid rupture rates and stress drops of megathrusts and earthquakes that nucleate in the upper lithosphere are commonly explained by brittle failure and frictional sliding on faults subject to tectonic stresses at low mean stresses (e.g., refs. 4 and 5). This process, facilitated by high pore pressures that reduce the effective normal stress on faults, is widely accepted as the mechanism for subduction zone earthquakes to depths of ~50 km. On the contrary, where subducting plates extend into the mantle transition zone (and below) from ~410 to 700 km, deep-focus earthquakes may be explained by localized transformation plasticity of olivine and pyroxene to wadsleyite, ringwoodite, and other high pressure phases under deviatoric stress (6–8). Dehydration reactions of hydrous minerals such as serpentine and lawsonite can extend the range of conditions where effective pressures are low by generating high pore pressures locally over narrow temperature intervals (e.g., refs. 9–11). However, the dehydration embrittlement hypothesis is not universally accepted for all intermediate-depth earthquakes (12–16). Experiments on localized, dehydrating shear zones within serpentine exhibit velocity hardening and stable shear rather than unstable softening (17). Serpentine dehydration has also been invoked for double Wadati–Benioff zones (15, 18–21). In the lower band of seismicity of double Benioff zones, sometimes 40 km below the top of the paleo-seafloor of the subducting plate (22), there is little or no evidence for the presence of hydrous phases other than the earthquakes themselves (12). The brittle failure-permeability model of Faccenda et al. (23) predicts penetration of the uppermost lithosphere by seawater and formation of serpentine and other hydrous silicates to a maximum depth of 30 km, corresponding to seismic constraints on upper plate alteration (24). Neither this modeling nor imaging predicts fluid penetration where bending stresses become compressive and inelastic deformation occurs by crystal plasticity. While some studies have suggested that the lower band of earthquakes of double Wadati–Benioff zones closely follows isotherms that coincide with antigorite serpentine dehydration (1, 2, 11), more detailed comparisons of predicted temperature and depth conditions of intermediate-focus earthquake locations in the Hokkaido, Tohoku, and Tonga subduction zones (15, 16, 25–27) demonstrate significant departures from the conditions for antigorite breakdown (SI Appendix, section 1 and Fig. S1). Many events fall at temperatures well below antigorite breakdown, and some exceed the maximum temperatures of antigorite stability. These events can be explained by dehydration and phase transitions that lead to stress transfer to neighboring unaltered peridotite (14, 15, 28, 29). However, intermediate-depth earthquakes at depths greater than 200 km require some other mechanism of nucleation. Seismicity continues to great depths, after dehydration and hydrous fluid expulsion may be complete along normal steady-state subduction zone geotherms (30). Earthquakes at depths of 50 to 400 km are observed in many subduction zones (Fig. 1) (e.g., refs. 1, 15, and 31) and the focal mechanisms of these events indicate failure within slabs. Such intraslab earthquakes require rapid-shear-failure mechanisms at high pressures and temperatures; in the absence of reactions that generate high pore fluid pressures, brittle failure and frictional sliding cannot explain the occurrence of these intraslab seismic events (3, 4). Another proposed explanation for intermediate-focus earthquakes involves viscous shear instabilities arising from positive feedback of viscous shear heating and temperature-dependent strain rates (12, 32, 33). In the model of Kelemen and Hirth (12), strain localization is initiated within fine-grained zones of peridotite by grain size-sensitive flow. This mechanism is an alternative to dehydration embrittlement, particularly in the lower band of seismicity in double Benioff zones, and it could also explain deeper earthquakes in fully dehydrated subducting plates, at depths between 200 and 400 km. Generalizing, the strain localization and thermal runaway shear predicted by the Kelemen and Hirth (12) model could result from any weak horizon within a stronger mantle host.

Ocean-floor carbonates and fluids may enter accretionary wedges at subduction zones, leading to reduced frictional strengths at high pore pressures and low effective pressures (34). Just as mantle peridotites react with aqueous fluids to form hydrous silicates, they may react with CO2 to form magnesian carbonates during hydrothermal alteration and weathering (35–37) or at high temperature and pressure conditions of the lower lithosphere (38–41). In oceanic plates, extensive low temperature carbonate alteration occurs, as documented for exhumed peridotites that have reacted at the seafloor (42, 43) or at the base of ophiolites (44–46). In many cases, field observations indicate that carbonates are distributed inhomogeneously as lenses and veins (46, 47) (SI Appendix, Fig. S2). There is also evidence for mantle-sourced CO2 in magmas (e.g., refs. 48–50). The presence of this CO2 is consistent with experimentally determined phase relations in peridotite C–O–H systems; CO2-bearing magmas can be generated, for example, via decompression melting (e.g., refs. 51 and 52). Some peridotite and pyroxenite xenoliths in basalts host fluid inclusions that are CO2 rich (53, 54). In some cases, high-temperature fluid inclusion formation is indicated by the presence of coexisting melt inclusions (55, 56). The presence of melt does not require great depths, and some CO2-rich fluids trapped in inclusions may have formed via magma outgassing at depths less than 30 km (57, 58). On the contrary, microstructures associated with the CO2-rich inclusions indicate formation during depressurization while high temperature dislocation creep processes were active at depths near or below the base of the lithosphere (59, 60). Deep CO2 inclusion generation is also supported by the presence of nanocrystalline diamonds in certain CO2 (+H2O + H2S) fluid inclusions (61). CO2-rich fluids that migrate through the mantle lithosphere are likely to react with peridotite to form carbonate veins and lenses at depths ≤100 km in oceanic plates (48, 49). This may be most common where plume-derived, undegassed primary magmas, such as those beneath Hawaii, have molar CO2:H2O ratios close to 1, and CO2-rich fluids begin to exsolve preferentially from these magmas near the base of the lithosphere (48, 50) (SI Appendix, section 2). Magnesian carbonates in the lithosphere will decompose to form enstatite +CO2-rich fluid or melt at mantle depths where temperatures are high (>1,000° to 1,200 °C; ref. 38), but not at the lower temperatures of subducting slabs (SI Appendix, section 2 and Fig. S3). At depths less than approximately 60 to 70 km, the upper temperature stability of dolomite will be limited by a decarbonation reaction in which dolomite and orthopyroxene (Opx) react to form clinopyroxene (Cpx), olivine, and CO2. At greater depths, the stability of magnesite is limited only by the carbonated peridotite solidus. Magnesite in cold subducting plates and in the relatively low temperature base of overlying mantle wedges may be stable well beyond the depth of intermediate-depth earthquakes, depending on oxygen fugacity (62–65). In and around subduction zones, strain may localize within weak bands of hydrous silicates where these bands are stable and within carbonates where hydrous silicates have reacted and fluid products have escaped (66). Experimentally constrained magnesite and dolomite flow strengths (67) indicate that magnesian carbonates are stronger than antigorite serpentinite (68–70) where subducting plates are hydrated, but weaker than wet or dry olivine (e.g., refs. 71 and 72) at all subduction zone conditions. Thus, for altered downgoing plates that contain both serpentine and carbonate veins, strains will be localized within serpentinites. Where hydrous silicates are not stable, or are absent due to low H2O contents, carbonates will focus viscous deformation within altered downgoing plates. Earthquakes arising from viscous shear heating instabilities during grain size-sensitive creep localized within thin, fine-grained, olivine-rich shear zones hosted in coarser grained peridotite were modeled by Kelemen and Hirth (12). This paper uses a similar approach to evaluate the potential for viscous earthquakes due to runaway shear heating in carbonate horizons within mantle peridotite. We have adapted the model of Kelemen and Hirth (12) using experimentally determined flow laws for magnesite, dolomite, and olivine (Methods), to evaluate the potential for unstable shear at the pressure and temperature conditions of intermediate-depth earthquakes reported for the Tonga and NE Japan subduction zones (2, 25, 27). Earthquakes of the Tonga subduction zone occur at all depths from very shallow trench events to the base of the mantle transition zone (25, 27), without a gap in seismicity at intermediate depths. Focal depth determinations of the NE Japan subduction system reveal double Wadati–Benioff zones with well-developed seismicity near the top and base of the downgoing slab (e.g., ref. 2).

Model parameters and results for the mechanical behavior of magnesian carbonate horizons within an olivine host are given in SI Appendix, with the range of parameters listed in SI Appendix, Table S1, flow laws given in SI Appendix, Table S2, stable and unstable shear results given in SI Appendix, Table S3 for magnesite, and stable and unstable shear results given in SI Appendix, Table S4 for dolomite. The one-dimensional thermal–mechanical model applied to olivine-hosted carbonate layers predicts strain localization in the carbonate layers at all conditions of subducting plates beneath Tonga and NE Japan at depth-temperature conditions inferred by Wei et al. (27) and Yamasaki and Seno (2), respectively (SI Appendix, section 3). Both stable, steady-state creep and unstable, runaway mechanical behaviors are predicted, depending on the relative rates of heat generation and diffusion. Stable, localized shear of a magnesite layer within altered peridotite is illustrated in Fig. 2 for an initial temperature Ti of 700 °C, pressure of 4 GPa, initial differential stress (σ1−σ3)i of 10 MPa, and far-field velocity (ffv) of 0.1 m/y, with a carbonate layer width wsz of 0.200 m and carbonate grain size d of 1,000 µm. Model stresses, strain rates, and temperatures are calculated, rising from the initial values, followed by complex behavior associated with the nonlinear flow laws of magnesite and olivine. Predicted magnesite layer temperatures rise to ~780 °C after 1,500 y of deformation, and magnesite strain rates reach nearly constant values of ~10−8 s−1, while the viscous strain rate of the peridotite host remains much slower (~10−16 s−1). Differential stress rises, from 10 MPa to a peak value of ~275 MPa, followed by a gradual decay to ~125 MPa (after 1,500 y) associated with gradual heating. Accelerating shear heating and strain rate, with a sudden stress drop, are not predicted even after 25,000 y. Unstable shear of a magnesite layer within peridotite is illustrated in Fig. 3 for an initial temperature Ti of 600 °C, pressure of 4 GPa, and initial stress of 10 MPa, for a carbonate shear zone width wsz of 0.200 m and carbonate grain size d of 1,000 µm. The model predicts a nearly linear rise in stress prior to an accelerating shear strain event at about 275 y, with complex nonlinear changes in viscous strain rate and temperature without imposing an ad-hoc mechanism of triggering. The shear zone temperature for each localized shear event rises slowly at first, and then increases rapidly (SI Appendix, Fig. S4), reaching ~1,200 °C, the solidus of carbonated peridotite at a nominal pressure of 4 GPa (as shown in SI Appendix, Fig. S3). Strain rates within the magnesite layer accelerate to ~10/s within seconds as temperature rises (SI Appendix, section 4 and Figs. S4 and S5), yielding rates of displacement across the carbonate layer that are equivalent to earthquake slip rates. Differential stress reaches a peak value of ~1.3 GPa during early stages of the model, when viscous displacement lags behind the far-field, largely elastic displacement. Runaway temperatures and strain rates are accompanied by a rapid stress drop, all of which define a single unstable shear event, governed by the positive feedback of viscous shear heating to the temperature-dependent strain rate. Once conditions of melting are reached and solid-state flow laws of the model no longer apply, the stress is reset, and strain rate and shear heating rate decrease rapidly. This model remains quasi-static and we do not treat the dynamic slip that is expected once melting occurs. Once stress and strain rate decrease, temperature of the shear zone declines via conductive heat loss. Once temperature drops and viscous displacement in the shear zone drops below the far-field displacement, stresses once again increase. Following each slip event, the shear zone temperature remains higher than the initial model temperature, establishing the new beginning point for each subsequent cycle. Running the model forward, we find that viscous shear heating produces repeated earthquake events (Fig. 3). Repeating intermediate-depth earthquakes may occur along the same carbonate layer or they may represent events on multiple carbonate layers. Shear heating varies as stresses and strain rates of the shear zone evolve and heat is conducted from the shear zone to surroundings by diffusion. For the example shown here, a magnesite shear zone 0.200 m wide with an initial temperature of 600 °C exhibits a recurrence interval of unstable shear events of ~250 y. Peak stresses in subsequent events decline as the minimum temperature of the carbonate layer between events rises. Ultimately, when residual temperatures of the magnesite layer and surrounding peridotite at 4 GPa approach 700 °C, strain rates, stress, and temperature within the magnesite layer reach steady values. In addition to the initial temperature, the transition between stable and unstable behavior depends on far-field strain rate, shear zone width, and grain size, as analyzed for fine-grained olivine shear zones by Kelemen and Hirth (12). Strain rates in the magnesite layer can reach stable, steady-state values over a wide range of conditions, for example at high temperatures over a wide range of far-field strain rates (Fig. 4A). Strain rates become unstable at other conditions, for example at low temperatures, with runaway strain rate events favored by high far-field velocities and narrow magnesite layers. The transition between stable and unstable shear strain rates is illustrated in Fig. 4B by plotting peak strain rates achieved in the shear zone as a function of ffv normalized by the magnesite shear zone width, or time-averaged fictive shear zone strain rate (as defined in ref. 12). Model predictions for stable and unstable shear have also been made for dolomite horizons for depths up to 90 km (SI Appendix, section 5 and Table S4). Although unstable shear is also predicted for dolomite lenses for altered downgoing lithosphere, its field of stability is unlikely to extend beyond that of antigorite, and dehydration embrittlement may continue to be the primary earthquake mechanism.

Magnesite and dolomite deformed by dislocation creep and diffusion creep have experimentally constrained effective viscosities that are many orders of magnitude lower than effective viscosities of peridotite (67). Strain rates predicted in carbonate layers in our model are at least five to seven orders of magnitude larger than the maximum viscous strain rates of the host olivine under all conditions tested (SI Appendix, Table S1), even during stable, localized deformation in carbonate layers at high temperatures (greater than 700 °C to 900 °C at P = 3 GPa and 15 GPa, respectively). At low temperatures (500 to 650 °C at 4 GPa), peak strain rates in carbonate layers (with grain sizes of 50 to 1,000 µm) accelerate and lead to repeated shear instabilities, which can explain intermediate-focus earthquakes of cold subducting plates at high plate convergence rates. As examples, we apply our one-dimensional model to earthquakes at depths of 200 to 400 km and temperatures inferred for the top, bottom, and core of downgoing lithosphere below Tonga and NE Japan (27 and 2, respectively). Parametric testing of localized viscous deformation in carbonate layers, carried out over 350 ≤ T ≤ 950 °C, 3 ≤ P ≤ 15 GPa, 0.0001 ≤ ffv ≤ 0.1 m/y, 0.02 ≤ wsz ≤ 0.2 m, wez = 1,000 m, 10 ≤ initial stress ≤ 100 MPa, 10 ≤ grain size ≤ 1,000 µm (SI Appendix, Table S1), reveals conditions over which stable and unstable viscous deformation may occur as subducting plates reach increasingly higher temperatures and pressures and temperatures of the overlying mantle change. Model parameters that govern stable versus unstable viscous deformation include initial temperature, pressure, far-field strain rate, shear zone width, and grain size. Model predictions are not strongly dependent on initial stress or the melting temperature at which shear heating is assumed to stop, since runaway strain rates are achieved only at the high stress and temperature conditions of each rapid, unstable event. Model results indicate that shear instabilities of magnesite horizons in altered subducting plates and overlying mantle wedge can occur at conditions that correspond to intermediate depths and double Wadati–Benioff zones where seismic activity is observed (Fig. 5). Earthquakes may be due to dehydration embrittlement where antigorite of the altered slab reaches its upper temperature limit at depths up to ~200 km (reaction curve in green after ref. 73). Fig. 5A shows the temperature–depth paths inferred for the top, near bottom (~65 km below paleo seafloor), and core (~35 km below paleo seafloor) of the subducting oceanic plate beneath Tonga (ffv = 0.24 m/y), with instabilities predicted for carbonate horizons at low temperatures and our highest ffv (shown in pink) and stable shear at high temperatures (above pink curve). For the thermal models of Wei et al. (27), unstable shear events may occur within the top, bottom, and core of the subducting plate and within the base of the overlying mantle beneath Tonga to depths of 450 km. Fig. 5B shows the temperature–depth paths inferred for the top and bottom (~55 km below paleo seafloor) of the subducting plate and the adjacent mantle wedge beneath NE Japan (ffv = 0.08 to 0.09 m/y), with conditions of dehydration embrittlement shown at depths less than 200 km (antigorite breakdown in green) and those of viscous shear heating instabilities in carbonate lenses (in pink). Based on the earthquake locations and thermal model of Yamasaki and Seno (2), unstable shear events in carbonate horizons are predicted at deeper levels of both the base and the top of the subducting plate and for the base of the mantle wedge. While seismic events in both upper and lower bands of double Benioff zones have been ascribed to dehydration embrittlement (2, 10), there is no evidence that seawater alteration is present at the depth of the lower band beneath NE Japan. Most seismic events in the lower band occur at depths less than 250 km (2), 40 km below the paleo-seafloor in the subducting plate. Our model of viscous shear heating events of carbonate layers within a peridotite host predicts thermal instabilities in the subducting plate over this depth range. Differential stresses predicted by this one-dimensional model just prior to unstable shear are significantly larger than seismological constraints for stress drops of intermediate-focus earthquakes (74). Yet, our model input parameters come from laboratory measurements of flow laws, elasticity, and conductive loss of heat, at conditions inferred for well-studied subduction systems. More generally, high stresses are a feature of many experimentally derived yield envelopes for lithosphere (75). Lower stress drops may be calculated by fully coupled, 2- or 3-dimensional thermomechanical models. Alternatively, magnitudes of stress drops predicted by our model and by seismic moments may differ because our model applies to the initiation of unstable shear over a small nucleation area, while seismological stress drops correspond to dynamic rupture following nucleation, averaged over much larger shear zone areas (76). Our model of viscous strain localization within carbonates of downgoing altered lithosphere may explain seismicity by thermal shear instabilities at intermediate focal depths of other cold subduction zones (e.g., Sumatra subduction zone; ref. 30). In contrast, stable localized shear is predicted for carbonate lenses in hot subduction zones (at temperatures above the transition between stable localized creep and unstable strain softening, Fig. 5) where seismicity does not extend to depths greater than 150 km (e.g., Cascadia and Nankai subduction zones; ref. 26). For subduction zones in which unstable shear is predicted for low initial temperatures, shear heating during each runaway viscous shear event raises the local temperature. As temperatures in the subducting plate and adjacent mantle wedge exceed conditions for viscous earthquakes, stable deformation localized within carbonate layers will continue, with shear heating raising local temperatures above average values for the overlying mantle wedge and the top, core, and bottom of the subducting plate at a given depth. If such a process is ubiquitous and prolonged, it may change subduction zone geotherms based purely on diffusive heating of the slab as it penetrates increasing mantle depths in the absence of other heat source terms. The global occurrence of earthquakes at subduction zones correlates with temperature distributions (e.g., ref. 26); intermediate-depth earthquakes are more frequent for colder downgoing plates. At the same time, the extent of peridotite alteration and related occurrences of heterogeneously distributed carbonates and their grain sizes may vary. Real variations in temperature, alteration products, and distributions of weaker phases may explain why some subduction zones show wide ranges in focal depths between 50 and 400 km, why some show well-defined double Wadati–Benioff zones, and why intermediate-depth earthquakes are absent in others (77). The primary uncertainties in our proposed mechanism for intermediate-depth earthquakes are the extent to which subducting plates and the overlying mantle wedge are altered and the spatial distributions of carbonate alteration products. The variable frequency of earthquakes at depths of 200 to 400 km in different subduction zones (77) may be related to different proportions and spatial distributions of carbonate lenses in altered plates. However, this model predicts unambiguously that shear heating and thermal softening lead to seismic displacement rates if and when carbonates (or other low viscosity materials) are heterogeneously distributed in subducting plates. Model inputs for carbonate layers, their geometry, and dimensions come from direct observations of dolomite and magnesite vein networks in altered peridotites in the mantle section of ophiolites, formed both by low temperature near-surface weathering (46, 47, 78) and by higher temperature reactions of the mantle wedge with CO2-rich aqueous fluids derived from subducting sediments and altered basalts (79–81). However, no field observations are available for the deeper parts of the mantle lithosphere in subducting oceanic plates. We can only infer that carbonates formed by reaction with mantle plume-derived CO2 within the lithosphere could be distributed heterogeneously (38, 41). In the absence of direct field observations, the presence of carbonates in altered mantle rocks may be inferred from seismological observations, such as mid-lithospheric low-velocity horizons (82) and seismic anisotropy and fabric in subducting plates (83). Earthquakes at depths >50 km are not readily explained by brittle and frictional processes without abnormal pore pressures, and a number of exotic earthquake mechanisms have been proposed. Localized viscous shear heating of peridotite-hosted carbonate layers, as modeled here, adds a mechanism that can produce rapid stress drops over a larger depth range than many proposed processes. This mechanism can occur in regions where significant proportions of fluid at lithospheric pressure are absent, and it may also explain some earthquakes in the depth range where other processes are responsible for some of the seismic events. Localized viscous shear heating may also act in concert with other proposed mechanisms, for example facilitating very rapid heating and devolatilization, or leading to strain rates and localized differential stresses sufficient for brittle failure at high confining pressure. The most frequent earthquakes of subduction zones occur within the uppermost 50 km (Fig. 1B) where brittle and frictional processes are most important, facilitated by elevated pore pressures and low effective pressures. Deep earthquakes from ~410 to 700 km are probably due to localized transformation plasticity of olivine to wadsleyite and ringwoodite (6, 7). Earthquakes from ~50 to 400 km are more difficult to explain. Earthquakes between ~50 and 200 km are commonly explained by serpentine-dehydration embrittlement (2, 3, 10). However, there are several observations and lines of reasoning that suggest this single mechanism cannot explain all earthquakes of this depth range: i) experiments on localized, dehydrating shear zones within serpentine exhibit velocity hardening and stable shear rather than unstable softening (17), ii) in the lower band of seismicity of double Benioff zones, sometimes 40 km below the top of the subducting plate (22), there is little or no evidence for the presence of hydrous phases (12), and, iii) the temperatures and depths predicted for the locations of intermediate-depth earthquakes generally do not correspond to conditions of antigorite serpentine dehydration (SI Appendix, Fig. S1) (25–27). Our model, can explain intermediate-depth earthquakes in regions where dehydration reactions are slow or absent, and high pore pressures are not present. Below ~200 km, thermal models suggest that temperatures in the subducting plate and overlying mantle will lead to complete serpentine dehydration, even for cold subduction zones (12) (Fig. 5). Thus, dehydration embrittlement is unlikely at depths of ~200 to 400 km in most subduction zones. Our model offers a plausible explanation for earthquakes in this depth range. The model presented here can be generalized to encompass viscous shear heating instabilities in any thin, weak layer within a stronger host rock, provided that the viscosity contrast remains present over the interval of heating. Candidates for such weak layers include fine-grained mylonites (or fault zones that originally deformed as cataclasites) that deform by grain size-sensitive creep within mineralogically similar, coarser host rocks, carbonate veins in altered basalt, foliated serpentinite shear zones in anhydrous peridotite, and wet opal/chalcedony/quartz bands within carbonate or greywacke sedimentary sections. All might produce viscous earthquakes over appropriate ranges of temperature, layer width, and other environmental parameters.

1. Unstable viscous shear heating instabilities can result from localized viscous deformation of thin carbonate layers within peridotite, in subducting plates and overlying mantle. 2. In cold, thick subducting plates and overlying mantle, such instabilities may explain intermediate-focus earthquakes at conditions that do not favor brittle failure or dehydration embrittlement. 3. Earthquakes of double Wadati–Benioff zones may result from a number of mechanisms, including dehydration embrittlement at low effective pressures, particularly in the upper band of earthquakes near the altered top of the subducting plate, while shear heating instabilities may give rise to sudden stress drops and high rupture rates in both the upper and lower bands of earthquakes. 4. Viscous shear-heating instabilities, such as those modeled here for weak carbonate layers within a stronger peridotite host, may provide a general explanation for earthquakes beyond the limits of dehydration reactions, at subduction zone depths from ~200 to 400 km.

The MATLAB code developed by Kelemen and Hirth (12) was modified to model stable and unstable viscous deformation via dislocation and diffusion creep in a carbonate layer within a uniform olivine host, subject to a constant layer-parallel ffv (far-field velocity) (Fig. 6). Carbonate layer widths, wsz, were varied from 0.020 to 0.200 m, based on field observations of magnesite veins of the Koniambo massif (47) formed by near-surface alteration of mantle peridotite and far-field velocities were varied from 0.0001 m/y to 0.1 m/y. The stronger peridotite host rock, with a width wez (of 1 or 10 km), mainly deforms elastically and represents a mechanical half-space adjacent to the carbonate layer. Strain rates within carbonate layers were determined using dislocation and diffusion creep laws of magnesite and dolomite from Holyoke et al. (67, 84). The sum of these strain rates determines the viscous strain rate in the shear zone, : where σ is differential stress in MPa, d is grain size in µm, exponents n and m are dimensionless, E* is the activation energy (kJmol−1), P is pressure (Pa), V* is the activation volume (m3 mol−1), R is the gas constant (8.314 Jmol−1K−1), and T is the temperature (K). Flow law parameters for dislocation creep (Eq. 1) of magnesite and dolomite are given in SI Appendix (SI Appendix, Table S2), as well as parameters for diffusion creep (Eq. 2) of both carbonates. The pressure effect on dislocation and diffusion creep of magnesite is given by the PV* terms, using V* values reported by Millard et al. (85). We did not include the PV* terms for dolomite creep since its activation volumes for dislocation and diffusion creep are not known, and the combined PV* terms for dolomite creep are not likely to be much larger, in its field of stability, than PV* of olivine dislocation creep. The sum of strain rates calculated for carbonate dislocation and diffusion creep was used to determine the total viscous response of the carbonate shear zone. Unlike the model of Kelemen and Hirth (12), we held the grain size of the carbonate shear zone constant, making no assumptions about grain growth or reduction by recrystallization during deformation. Reductions in grain size of coarse-grained carbonates could lead to further evolution of strain rates as diffusion creep becomes more important with strain. However, the recrystallized grain size–stress relation of neither magnesite nor dolomite is known, and physically based models of grain size reduction have not been developed. Instability is nevertheless predicted for shear zones consisting of both carbonates within an olivine host when conditions favor greater mechanical heat generation than diffusive heat loss. We model deformations of carbonate shear zones of widely varying grain sizes to evaluate conditions of stable and unstable shear due to dislocation and diffusion creep laws of magnesite and dolomite, without following changes in grain size and strain rate evolution over changing proportions of dislocation and diffusion creep. The mainly elastic peridotite zone undergoes the same stress history as the carbonate shear zone. The small strains of the surrounding peridotite were determined using the elastic properties and inelastic creep laws for wet olivine (71, 86), with V* determinations for olivine (SI Appendix, Table S2). Flow law parameters of dislocation creep (Eq. 1), diffusion creep, and dislocation-accommodated grain boundary sliding (both of the form of Eq. 2) for wet olivine were used to calculate strain rates and summed to determine the total viscous response of the peridotite host (SI Appendix, Table S2). We ran models with far-field strain rates of 10−12 to 10−15 s−1, based on the ffv and elastic zone width, wez. Total strain was determined by the product of total time elapsed and ffv divided by (wez + wsz). Stress was taken to be constant everywhere in the model at each time step, given by σ = [Golivine (wez + wsz) εtotal – Gcarbonate wsz εviscous]/(wez + wsz), where shear moduli, G of the carbonate shear zone and of the olivine host, are given in SI Appendix, Table S1. Heat Q generated by dissipation of viscous work (W = σ εviscous) in the model is lost from the shear zone by one-dimensional conductive heat flow H, and shear zone temperatures are determined from dT/dt = −κ[d2T/dx2] + σ /(Cp ρ) (Fig. 6), for all T < ~1,000 to 1,600 °C (below the solidus of carbonated peridotite), where κ is thermal diffusivity of the olivine, Cp is heat capacity, and ρ is density of magnesian carbonates. The Crank–Nicholson method is used to calculate temperature of the host olivine over a 10-km-wide region around the shear zone, holding temperature constant at the far boundaries. The time step is the smaller of 1°/dT/dt, wsz2/κ, or 0.1/ . Once conditions of a carbonated peridotite solidus are reached, partial melting will rapidly decrease viscosity, stress, and shear heating within the shear zone. Thus, we impose an upper bound to temperature in the model, at the solidus of carbonated peridotite. We used initial model temperatures ranging from 350 °C to 950 °C, as predicted within subducting plates and at the base of the overlying mantle wedge in the Tonga and NE Japan subduction zones at depths from the seafloor to more than 400 km (SI Appendix, Fig. S3) (25–27). Modeling has been performed for pressures of 3 to 15 GPa; representative results for stress, strain rates and temperature as a function of time (Figs. 2 and 3 and SI Appendix, Figs. S4 and S5) are shown for a nominal pressure of 4 GPa. Initial differential stresses of 10 MPa to 100 MPa are assumed and grain sizes of the carbonate shear zones are assumed to be 10, 50, and 1,000 µm for different model runs. The full range of parameters tested is given in SI Appendix, Table S1. All models are run for time intervals of 5,000 y to predict stable or unstable events in the shear zone. During each unstable shear heating event, strain rate and temperature increase rapidly, reaching melting conditions of carbonate-bearing peridotite (51, 64). Melt viscosities many orders of magnitude lower than those of solid-state viscous laws are important to the dynamic shear event, but are not treated here. Nevertheless, the stress with increasing melt fraction must drop precipitously and is assumed to reach ~1 MPa, as assumed in the Kelemen and Hirth model (12). Accordingly, viscous strain rate and the rate of shear heating are reduced, leading to a period of cooling following the runaway shearing event.

Click here for additional data file.