Metamorphism-facilitated faulting in deforming orthopyroxene: Implications for global intermediate-depth seismicity.

SignificanceThe exothermic metamorphic reaction in orthopyroxene (Opx), a major component of oceanic lithospheric mantle, is shown to trigger brittle failure in laboratory deformation experiments under conditions where garnet exsolution takes place. The reaction product is an extremely fine-grained material, forming narrow reaction zones that are mechanically weak, thereby facilitating macroscopic faulting. Oceanic subduction zones are characterized by two separate bands of seismicity, known as the double seismic zone. The upper band of seismicity, located in the oceanic crust, is well explained by dehydration-induced mechanical instability. Our newly discovered metamorphism-induced mechanical instability provides an alternative physical mechanism for earthquakes in the lower band of seismicity (located in the oceanic lithospheric mantle), with no requirement of hydration/dehydration processes.

Intermediate-depth (∼60 to 300 km) earthquakes usually delineate two subparallel seismicity bands, collectively known as the double seismic zones (DSZs). The upper band of seismicity (UBS) is concentrated mostly in the basaltic crust, and the lower band of seismicity (LBS) is ∼10 to 40 km below and well within the lithospheric mantle (1, 2). Fig. 1 shows nine subduction corridors with well-characterized DSZs. An example of seismicity distribution is taken from NE Japan (Fig. 1). The two bands merge at various depths, but seismicity often continues to ∼300 km in depth, mostly as a deeper extension of the LBS (1, 2). When we refer to LBS in this paper, we implicitly include this deeper extension.

Fig. 1.

Subduction corridors, DSZ characteristics, and expected phase relations in the LBS. (A) Nine subduction zone corridors with DSZs analyzed in this study. These subduction corridors have both thermal structure models (orange patches; reported by ref. 18, except for New Britain [corridor] modeled in the present study) and LBS earthquake distribution (reported in ref. 2). (B) Earthquake distribution under northeastern (NE) Japan (corridor 7) (74), from 39.5°N to 40.5°N, as a representative example of DSZ, together with the thermal model of ref. 18. The earthquake data used by ref. 74 are from the Japan Meteorological Agency catalog and National Research Institute for Earthquake Science and Disaster Resilience between 2004 and 2008. Red dots denote LBS earthquakes and are the primary focus of this study. (C) T-d trajectories of LBSs in the nine subduction zones in comparison with phase boundaries of antigorite and the CMAS system. Red (slabs younger than 100 Ma at trench) and blue (older than 100 Ma) solid lines with numerals are the T-d paths of the LBSs in A (). The dashed portion of these lines corresponds to the deeper extension of LBS where UBS no longer exists. Gray dashed lines delineate simplified phase relations in antigorite (Atg) (20). The four phase boundary segments are (a) Atg ↔ phase A (phA) + enstatite (En), (b) Atg ↔ Fo + En, (c) Atg ↔ Fo + talc, and (d) En + phA ↔ Fo. Black solid lines mark major phase boundaries of the CMAS system (65). Solid gray lines labeled “Al-n” are equilibrium pyrope content (in mol%) in Opx. The solid gray line “Py-99” is the pyrope content (99 mol%) in Grt (with 1 mol% majorite). Pyrope content decreases with increasing pressure. HCpx, a high-pressure polymorph of Opx. If 10% Fe is added to the system, the boundary between Fo + Cpx + Opx + Sp and Fo + Cpx + Opx + Grt will be lowered in pressure by roughly 0.2 GPa. The pyrope isopleths will also be lowered by ∼0.2 GPa (32).

Given adequate elastic strain energy, the occurrence of an earthquake requires three conditions: 1) large enough differential stress to cause failure, 2) localization of the failure into a thin zone of slip, and 3) rapid weakening of the slip zone to cause slip instability and propagation. In the cold and brittle part of the lithosphere, either at shallow depths or in a subducting slab, there is no shortage of differential stress to meet the first condition. At shallow depths, frictional faulting can bring about the other two conditions. In the slab, however, as confining pressure increases with depth, frictional faulting is increasingly prohibited (3), and standard rock mechanics can only predict brittle failure to occur in the form of cataclastic flow. Therefore, other mechanisms that promote localized friction-like failure must be sought to explain the observed occurrence of earthquakes at greater depths.

Dehydration-triggered embrittlement is by far the most prevailing hypothesis to explain both UBS and LBS earthquakes (4–11). It is assumed that extremely high pore fluid pressure due to the dehydration of hydrous minerals acts to offset the effect of high confining pressure so that the frictional faulting mechanism commonly seen at shallow depths can also operate at large depths, although there are still unresolved issues regarding building sufficiently high pore fluid pressure in open systems typical of subduction zone environments (12). Dehydration-driven mechanisms remain reasonable hypotheses for the UBS in the basaltic crust (13) and may apply to the LBS to some extent (10), but it is more difficult to envision their extensive operation in the LBS that is confined within the subducted mantle, which, similar to harzburgite or lherzolite, consists primarily of olivine (Mg,Fe)2SiO4 and orthopyroxene (Opx), here simplified as (Mg,Fe)2(Si,Al)2O6 (14). The physical mechanism for H2O to penetrate several tens of kilometers below the slab surface before subduction (to enable dehydration embrittlement after subduction) is still a matter of debate. Exhumed oceanic mantle rocks (>60 km in depth) from ancient subduction zone segments along the Western Alps contain dry pseudotachylyte and minerals devoid of hydration (15, 16). A recent seismological study reported no anomalous Vp/Vs ratios (where Vp and Vs are seismic compressional and shear velocities, respectively) in the LBSs under northern Chile and northeastern Japan, suggesting that the oceanic lithospheric mantle is dry (17).

Fig. 1 shows temperature-depth (T-d) trajectories along the lengths of the LBSs (), based primarily on a slightly refined version of slab thermal models by Wada and Wang (18, 19) (). The model for New Britain was constructed in this study using the same approach ( and ). In young slabs (corridors 1 to 4; ), the separation between UBS and LBS is narrow (≤15 km). LBSs are cold at shallow depths (red curves) with temperature rising with depth. By ∼150 km, temperatures in these LBSs reach the antigorite dehydration phase boundary (20), where, however, seismicity has largely ceased in LBSs in most young slabs (). In older slabs (corridors 5 to 9), LBSs are further away from the UBS (separation of ≥∼30 km), with little thermal interaction with the overlying mantle. Therefore, LBSs are largely isothermal instead of following dehydration boundaries (Fig. 1), with seismicity continuing to more than 300-km depths ().

The difference in T-d trajectories shown in Fig. 1 is primarily controlled by heat conduction from the overlying mantle into the slab. It is convenient to characterize the thermal states in LBS and UBS using the so-called thermomechanical parameter , where and are heating rate and strain rate, respectively (5). While typical values of in UBSs fall in the range of 102 to 104 K (9), for LBSs, the values are on the order of 10 K for young slabs and essentially 0 for old slabs. The different T-d trajectories and thermomechanical parameters of LBSs in young versus old slabs demonstrate the limitation of the dehydration hypothesis for LBS earthquakes.

Another hypothesis to explain the LBS is thermal runaway, based primarily on theoretical and numerical analyses (21, 22), although some experimental results began to emerge (23). In this hypothesis, positive feedback between shear localization and shear heating turns viscous shear into runaway slip to produce an earthquake (24, 25). Grain size reduction is proposed to play an important role in this process (24). Thermal runaway does not require specific mineralogy and/or water content and has also been suggested to operate in deep-focus (>300-km) earthquakes (26–28). However, the mechanism tends to predict stress drops much larger than observed in real earthquakes (29). The absence of a preexisting brittle fault plane also presents challenges to building a connection with contemporary theories for the dynamics of rupture propagation.

In addition to the difference in thermomechanical state, different lithology should also be considered. While UBSs are within the basaltic crust, LBSs are in the oceanic mantle, which consists predominantly of olivine and Opx (30). To seek an additional or alternative mechanism for LBS earthquakes, we examine metamorphic transformations in aluminum-bearing Opx, one of the major minerals in the oceanic lithospheric mantle. Some background information about the phase relations in the mantle lithosphere is given in the next section. We examine mechanical behavior of Opx aggregates during metamorphic reactions in controlled deformation under pressure and temperature conditions corresponding to deep subducted slabs. Using a combination of in situ acoustic emission (AE) monitoring and postmortem microstructural analyses, we confirm the presence of metamorphism-facilitated shear localization, weakening, and brittle faulting in deforming natural Opx. Based on the correlation among metamorphic phase relations, kinetics arguments, reported slab thermal modeling, and observed seismicity distribution in the LBSs, we suggest that metamorphism-facilitated mechanical instability in Opx is a viable mechanism for earthquakes in subducted mantle.

Results

Petrological Background: Phase Relations in Dry Mantle Lithosphere.

Also shown in Fig. 1 is the phase diagram of the CaO-MgO-Al2O3-SiO2 (CMAS) system (31), which forms the basis for understanding the metamorphic reactions in a harzburgitic or lherzolitic mantle lithosphere. From the surface to ∼20 to 30 km in depth, thermodynamically stable phases are forsterite (Fo), Opx, clinopyroxene (Cpx), and plagioclase (Plg), the latter being the main host for Al. Below ∼20 to 30 km, Plg decomposes, forming spinel (Sp), another major Al-hosting mineral, with Ca contents primarily entering Cpx. Below ∼50 to 60 km, garnet (Grt), the high-pressure Al host, exsolves from Opx through substitution , where Me stands for Mg and Fe, both occupying the same crystallographic site (e.g., MeSiO3 = [Mg,Fe]SiO3) via the following simplified reactions:and

Reaction breaks down the initially Al-rich Opx into a fine-grained mixture of Al-rich Grt and Al-poor Opx. Within the Fo + Cpx + Opx + Grt field, as pressure increases, Opx sheds Al continuously by exsolving Grt, whose Al content also decreases with increasing pressure. Consequently, Grt volume fraction also increases. Even after the remaining Opx transforms to its high-pressure polymorph known as high-pressure clinopyroxene (HCpx) in the Fo + Cpx + HCpx + Grt field, Grt continues to exsolve with increasing Si content (reaction ). The pressure (depth) required for the onset of reaction is sensitive to both temperature and Al content in Opx. For example, in a lithospheric mantle containing Opx initially with 4% Al per 12 oxygen, the equilibrium reaction pressure is ∼2 GPa (60 km) and 2.5 GPa at temperatures of 700 K (approximate conditions for Alaska and northern Sumatra LBSs) and 900 K (condition for sourthern Mariana LBSs), respectively, as indicated by the Al isopleths in Fig. 2. For an Opx with 2% Al, the reaction pressure is ∼2 GPa (60 km) at 900 K.

Fig. 2.

Mechanical behavior of deforming Opx during the metamorphic reaction. (A–C) An example of experimental observations in D2071. (A) Differential stress and strain versus time. The two vertical dashed lines highlight the time window of accelerating stress decrease. Note an acceleration in strain, as indicated by the red dashed lines and arrows. (B) AE rate (blue bars) and cumulative event counts (black curve) during deformation. The AE burst shortly before the 1,800-s mark corresponds to the stress peak in A. (C) Relative AE event amplitude (blue dots with sticks) and energy release rate (black curve). Most events have similar magnitudes, corresponding to moment magnitude (Mw) around −8.5 to −9.8. Energy release rate peaking around 3,300 s marks the beginning of the final stress drop stage in A.

Under thermodynamic equilibrium conditions, once the reaction begins, Grt continuously exsolves from Opx with increasing depth over the entire Fo + Cpx + Opx + Grt field, according to reaction . While Opx becomes increasingly Al depleted, Al content in Grt also decreases slightly from 100% pyrope to 99% (the Py-99 isopleth) following reaction . After Al content is exhausted in Opx, reaction becomes dominant. With further increases in pressure, both Cpx and HCpx begin to transform to Grt. Volume fraction of Grt continues to increase with depth, with increasing majorite content. The above description follows strictly the CMAS phase diagram (Fig. 1) for simplicity. Available experimental data show that replacing 10 mol% Mg with Fe in Opx changes the major phase boundaries in Fig. 1 only by ∼0.2 GPa without affecting the general phase topology (32).

Within subducting slabs, the onsets of reaction and reaction are expected to be delayed to greater depths due to the significantly colder environments (Fig. 1). The effects of these kinetically delayed metamorphic reactions on the mechanical behavior of subducted oceanic lithosphere are the main subject of this study.

Embrittlement of Metastable Opx during Syn-Deformation Metamorphism.

We conducted triaxial deformation experiments on presintered, nominally dry natural Opx, with FeO and Al2O3 contents of ∼14.5 mol% and ∼2 mol%, respectively (; ). This composition does not represent realistic components in the oceanic lithosphere; rather, it is considered an analog for us to understand the physics of mechanical response to metamorphic reactions. Ferrous iron (Fe2+) is predominant, which behaves similarly to Mg2+ in reaction and reaction . Experimental conditions (Table 1) correspond to depths from ∼50 to 250 km. The sample was under maximum compressive stress σ1 vertically and confining stress σ2 = σ3 horizontally. By experimental necessity, strain rates used in the laboratory (ca. 10−5 s−1) have to be many orders of magnitude higher than those in subducting slabs (ca. 10−15 s−1). There is a well-known tradeoff between temperature and strain rate, to be discussed to some detail in . The temperature values in the description of the results below are those in the laboratory experiments, but they should be understood to represent much lower temperatures in real slabs. The deformation experiments were conducted under isothermal conditions, corresponding to a thermomechanical parameter of zero.

Table 1.

Mechanical behavior and AE activities in polycrystalline Opx samples under various conditions

Experiment ID	P (GPa)	T (K)	Opx stability	Strain rate (s−1)	Peak stress (MPa)	Maximum strain (%)	No. of triggered AE events	Mechanical behavior	Phases identified
D2070	1.8	1,273	Stable	6.2 × 10−5	512	24.1	0*	Ductile	Opx
D2071	5.9	1,123	Metastable	6.8 × 10−5	1,750	29.8	261	Micro- and macrofaults	Opx + Cpx + Grt
D2202	5.3	1,173	Metastable	7.0 × 10−5	1,400	36.2	337
D2204-2†	5.1	1,273	Metastable	9.8 × 10−5	415	24.1	2,223
D2081	4.8	1,373	Metastable	7.6 × 10−5	285	34.5	2,698
D2069	5.1	1,473	Metastable	8.3 × 10−5	385	34.8	2,547	Microfaults only
D2077	5.2	1,523	Metastable	8.3 × 10−5	330	25.6	0*	Ductile	Opx + Cpx + Grt

*AE activity was monitored, but no events were detected.

†Experiment with two deformation segments. At 673 K (D2204-1), the sample was deformed up to ∼12% strain. AE was monitored, but no events were detected. Then, temperature was ramped while deformation continued. By ∼25% strain, sample temperature reached 1,273 K, where the second isothermal deformation segment (D2204-2) began (). During temperature ramping, AE was monitored, but no events were detected until ∼1,173 K, where the first AE appeared ().

We conducted one experiment to simulate the relatively low-pressure condition in which Opx is thermodynamically stable, corresponding to the Grt-absent field in previous studies on Fe-containing systems (33). In this experiment (D2070), at 1.8 GPa and 1,273 K, not a single AE event was recorded with axial strain up to 24.1%. It serves the purpose of verifying the ductile behavior of stable Opx reported previously (34, 35). All other experiments were conducted at higher pressures (∼4 to 6 GPa) and were designed to investigate how exsolution of Grt from metastable Opx affects deformation behavior (Table 1). In these experiments, samples failed consistently by macroscopic faulting while producing AEs, except for runs at very low (D2204-1) or very high (D2069 and D2077) temperatures.

Fig. 2 uses experiment D2071, conducted at 5.9 GPa and 1,123 K, to illustrate typical sample behavior during deformation. Differential stress peaked at 1.75 GPa around 3% strain (Fig. 2), where the first burst of AEs was recorded just before the 1,800-s mark (Fig. 2 ). Differential stress decreased afterward, and, after reaching ∼10% strain, AEs began to burst abundantly with several intense episodes. The largest single AE event occurred shortly before the 3,600-s mark (Fig. 2), after which the decrease in differential stress accelerated, with an abrupt strain jump around 4,500 s (Fig. 2), followed by maximum AE activity (Fig. 2). In total, 261 AE events were recorded. Relocation of AEs showed that all of them occurred inside the sample (). shows similar high-strength brittle behavior in other samples, although details differ from Fig. 2.

Between 4 and 6 GPa, no brittle behavior as exemplified by D2071 (Fig. 2) was observed when the temperature was lower than 1,123 K. In experiment D2204-1 (5.1 GPa, 673 K), the sample produced no AEs despite a very high differential stress of 2.6 GPa, indicating that Opx is strong and deforms in a ductile fashion at the low-temperature condition. However, when the same sample was heated further to 1,273 K (experiment D2204-2), more than 2,200 AEs were recorded with a constant differential stress of ∼0.4 GPa (). At 1,523 K (experiment D2077), not a single AE was detected, indicating a ductile behavior, although the peak differential stress was identical to those at 1,373 K (D2081) and 1,473 K (D2069).

The mechanical behavior of Opx with increasing temperature as seen in our experiments, excluding the low-pressure D2070, which represents stable Opx, is summarized in Fig. 3. Because all these experiments were at similar pressure conditions, we can focus on how temperature controls the modes of deformation in metastable Opx. The mechanical behavior can be roughly divided into three thermal regimes. Within regime 1, starting at least from 673 to ∼1,100 K, Opx is strong and ductile without AEs, with measured strength consistent with previous creep experiments (36). Thermal regime 2 is between ∼1,100 and 1,400 K, where strength decreases by a factor of >4 (solid circles), and the deformation is brittle with numerous AEs during deformation and macroscopic faults in recovered samples. In thermal regime 3 (above ∼1,500 K), Opx is weak and ductile, with no AEs detected, and the recovered sample contains neither faults nor intragranular fractures. A brittle-to-ductile transitional region may exist between regimes 2 and 3, as indicated by our observations at 1,473 K, where, although many AEs were recorded, the recovered sample contains no macroscopic faults but with numerous intragranular fractures (see details in the next section).

Fig. 3.

A summary of peak stresses observed in all the experiments between 4 and 6 GPa. Solid circles represent samples emitting numerous AEs during the experiment and contain macroscopic faults. Open circles represent samples with no AE activity nor macroscopic faults. The half-filled circle represents an experiment with some AE activity but no macroscopic faults. B-D transition: brittle-ductile transition.

Microstructure of Deformed Samples.

We characterized microstructures of recovered samples using scanning electron microscopy (SEM) and transmission electron microscopy (TEM). shows several recovered samples at low magnification. All SEM images are orientated such that the compression axis is vertical. The stable Opx sample that was deformed at 1.8 GPa and 1,273 K () contains no faults, except a horizontal decompression crack, consistent with the ductile behavior. The microstructure of samples at the Grt exsolution condition further verifies the mechanical behavior summarized in Fig. 3. For example, at ∼5 GPa, samples that were deformed in the “brittle” thermal regime 2 at 1,123 K (D2071; Fig. 2) and 1,373 K (D2081; ) contain numerous macroscopic faults, whereas those that were deformed at the “weak and ductile” thermal regime 3, 1,473 and 1,523 K (), contain no macroscopic faults, except for some decompression cracks.

We focused on the sample of experiment D2071 (Table 1, Figs. 4 and 5) for detailed microstructural characterization. The sample of D2081 () shows similar microstructure features. Both samples show well-developed conjugated macroscopic faults, inclined from the σ1 direction. In D2071, large conjugate faults cut through the entire sample (Fig. 4), with fault zones filled with fine-grained reaction products, predominantly Grt and Al-poor Opx (through reaction ) and occasionally Cpx (Fig. 4). Fine-grained Grt decorates the boundaries of Opx grains and forms thin lamellae within, forming continuous shear zones (Fig. 4). TEM shows that the shear zones are characterized by extremely fine-grained (<50-nm) reaction products (Fig. 4). Together, intergranular and intragranular Grt form elongated arrays subparallel to the macroscopic faults, whose displacements may be estimated using the needle-like ilmenite (Ilm; FeTiO3) impurities as markers. Fault offsets are common. The fault shown in Fig. 4 continues into Fig. 5, where a fault offset up to 50 μm can be identified by the sheared Ilm needle. Electron back-scatter diffraction analysis on large Opx grains shows a low degree of grain orientation spread in Opx (≤ 2°), suggesting grain boundary sliding to be a dominant deformation mechanism (). This apparently is facilitated by the presence of fine-grained Grt in and around Opx grains (Fig. 3; ). Intragranular Grt lamellae are primarily parallel to (100)Opx, consistent with those observed in natural samples (37). Cpx lamellae () are also present in the Opx grains and parallel to (100)Opx. Additional details and identification of reaction products are given in the .

Fig. 4.

Multiscale images of microstructure in D2071. (A–C) SEM images oriented with the macroscopic maximum compressive stress vertical. (A) Backscattered electron (BE) image of the entire recovered sample showing conjugate faults (NE–SW and NW–SE; near horizontal cracks are due to decompression). The two white boxes correspond to images of B and Fig. 5. (B) Enlarged view of the box in A. Grains appear in three gray levels. The darkest gray grains are Opx, and the brightest are Ilm (FeTiO3), which is from the original hand specimen. Grains with intermediate gray are fine-grained reaction product, which is predominantly Grt. The white dashed line with label “Fault” highlights the small fault cutting through the box in A. (C) Expanded view of the box in B. Grt-rich reaction product forms elongated lamella along Opx grain boundaries and within Opx grains. Overall, the lamella can be seen to form thin reaction bands in two general directions, approximately NE–SW (green arrows) and NW–SE (yellow arrows). (D) Bright-field TEM image of an ultrafine-grained reaction shear band (between yellow dashed lines) ∼0.5 μm thick. The inset is an enlarged area around the intersection of the shear bands, showing equant grains as small as <50 nm in linear dimension.

Fig. 5.

Examples of microfaults indicated by Ilm inclusions. BE images. The location of this image is shown in Fig. 4. One fault (with the upper side of the fault zone marked by the red dashed line) running NE–SW is joined by another fault (yellow dashed curve), offsetting a long Ilm needle by 50 μm (the fault continues into Fig. 4). Both fault zones are filled with extremely fine-grained reaction products. The inset is an enlarged view of the white dashed box, a better-preserved fault segment, showing details of the fine-grained material. Elsewhere, within large Opx grains, Grt-rich reaction products are also present in forms of either extremely thin lamellae or small patches. The image is oriented with the macroscopic maximum compressive stress vertical.

The sample that was deformed at ∼5 GPa and 1,473 K (D2069, transition between thermal regimes 2 and 3) contains no macroscopic faults, despite the numerous AEs recorded. SEM shows extensive Grt precipitation (> 10 vol%), essentially engulfing all the Opx grain boundaries (). We attribute the AEs detected to intragranular cracking, which is widespread in large Opx grains in this sample. We interpret this sample to be in the brittle-ductile transitional regime (Fig. 3).

Discussion

Tradeoff between Temperature and Strain Rate When Considering the Mechanical Behavior of Opx.

In our experiments, strain rate is ∼10−4 s−1, and the differential stresses range from 290 to 1,700 MPa in thermal regime 2 (Fig. 3). To compare the temperature range of LBS activities (∼700 to 900 K; Fig. 1) with that in laboratory thermal regime 2 (∼1,100 to 1,400 K), where brittle behavior of Opx is observed, we need to examine experimentally determined Opx flow laws in the form of , where n is the stress exponent, H is the activation enthalpy, and R is the gas constant. A recent creep study on Opx single crystals with similar Fe contents to our samples reported n∼4 at pressures of 3 to 6 GPa and strain rates of 10−4 to 10−5 s−1 (35). The temperature range (1,300 to 1,600 K) partially overlaps with our thermal regime 2 (1,100 to 1,400 K). The reported values of H are 290(60) and 490(60) kJ/mol, respectively, for the [001](100) and [001](010) slip systems. Based on these results and comparing with the microstructure and fabrics of naturally deformed Opx samples, these authors concluded that [001](100) slip has high activity under lithospheric conditions (35). For example, the dislocation microstructure they reported is similar to natural Opx samples from Balmuccia Massif (NW Italy), which underwent deformation at temperatures as low as ∼1,000 K (38). These last authors also concluded that microstructure of naturally deformed Opx from Balmuccia Massif is similar to that of laboratory deformed samples. Assuming that Opx deforms by the same mechanism in the lithospheric mantle, a laboratory shear stress of ∼0.2 GPa at ∼10−4 s−1 and T∼1,200 K (Table 1) would imply subduction zone stress of 10 MPa at ∼10−15 s−1 and T∼800 K, for an H of ∼270 kJ/mol, which is in general agreement with laboratory-determined activation enthalpy.

Mechanism of Metamorphism-Facilitated Embrittlement in Polycrystalline Opx.

The dramatic decrease in strength with increasing temperature in the brittle thermal regime 2 (Fig. 3) in our experiments is clearly due to metamorphic reaction, which produces ultrafine-grained Grt and Al-poor Opx along intragranular flaws (e.g., twin boundaries and stacking faults) and grain boundaries (Figs. 4 and 5, ). With increasing deformation, the reaction products form thin reaction bands, which then coalesce to form macroscopic shear zones in the sample. The ultrafine-grained product in the reaction bands is extremely weak and potentially an excellent fault gouge material for the following reasons. Most materials follow a similar d grain size dependence in their grain size sensitive flow laws, with p between 2 and 3 (39). For Grt, creep experiments show that p = 2.5 (40). Assuming that the nanomaterial in the reaction bands follows the same grain size dependence, for a factor of ∼1,000 grain size reduction in our experiments (from 50 μm to 50 nm), strain rates in the reaction bans are expected to be approximately ∼108 times as fast, capable of generating bursts of AEs to release elastic strain energy. In materials science, it has been known for decades that multiphase nanomaterials are remarkably weak and can flow to >1,000% strain even in tensile tests, while the same materials with larger grain size would fracture at strains less than ∼10%. A range of multiphase oxide ceramics with average grain size of ∼100 nm has been shown to deform plastically at very high strain rates up to 1 s−1 in tensile tests (41). Ductility should be even more pronounced under compression. Indeed, in friction experiments under a normal stress of 1 MPa, a fault gouge made of MgO nanoparticles (average grain size of 50 nm) has been shown to be able to reach seismic slip speed of 1.3 m s−1, most likely due to rolling of the nanoparticles during shear deformation, thereby dramatically reducing the dynamic friction coefficient (42).

Based on the mechanical data and microstructural observations, we propose a grain-scale model of metamorphism-facilitated embrittlement for polycrystalline Opx. During deformation of a metastable Opx rock, Grt exsolution takes place. This exothermic reaction (43) releases latent heat and, along with the resultant ultrafine-grained reaction products, induces a decrease in strength in Grt-rich reaction zones, triggering shear localization through enhanced grain boundary sliding and intragranular deformation. As the reaction progresses, strain localization becomes more pervasive in the reaction zones (Fig. 4). Grain boundaries that are “lubricated” by the ultrafine-grained reaction product connect with intragranular lamellae to form large faults in our samples (Fig. 5). As argued in ref. 44, shear heating is much faster under high confining pressure. For the same total shear, the thinner the shear zone, the more intense the shear heating will be. These thin slip zones preconditioned through the Grt exsolution process are expected to be greatly weakened by shear heating once slip is initiated.

It is difficult to estimate the required Grt volume fraction quantitatively for the embrittlement. Qualitatively, this embrittlement mechanism requires low volume fraction of reaction product, because the Grt exsolved from Opx is extremely fine-grained and expected to be superplastic. Our SEM and TEM observations show that such ultrafine-grained Grt tends to form very thin layers along grain boundaries, analogous to liquid. The highly selective spatial distribution of the weak Grt is the key in controlling the mechanical behavior: little volume faction of Grt is required to lubricate grain boundaries and cause stress concentration to produce large-scale shear localization.

This mechanism works only in a limited temperature range. In our experiments, when the temperature is below 1,100 K (thermal regime 1), kinetics are so sluggish that no metamorphic reactions occur within our experimental time scale. When the temperature is too high (>1,373 K; thermal regime 3), Grt exsolution becomes so widespread that shear localization is no longer possible. Only within thermal regime 2 can the exsolution kinetics trigger instability. Under such temperatures, the Al content in metastable Opx varies depending on the over-pressure of the Opx relative to its equilibrium Al content and the dwelling time. The “most metastable” Opx (i.e., the Opx that is farthest away from its equilibrium Al content) first overcomes kinetic barriers and exsolves Grt, whereas somewhat less metastable Opx requires more time and/or higher pressure to overcome kinetic barriers. This process takes place continuously, providing a continuous source of weak Grt to promote instability throughout the Fo + Cpx + Opx + Grt field.

Implications for LBS Seismicity.

The micromechanical model presented above is for pyroxenite rocks, while oceanic lithospheric mantle is polymineralic, consisting mostly of olivine and Opx. Both laboratory experiments and field observations show that, while making up no more than 30% of the total volume, Opx plays a major role in controlling mechanical properties of the lithospheric mantle (45–47). Before discussing failure models in subduction zone settings, it is necessary to consider ductile deformation and ductile shear zone development in oceanic lithosphere prior to subduction, because they prepare the rock texture eventually needed by metamorphism-facilitated seismic faulting after subduction.

The lithospheric mantle undergoes continuous basal shear as it travels toward the trench (48, 49). Ductile shear zones develop throughout the journey. Microstructure observations on field samples show that shear zones are characterized by significant grain size reduction, from coarse-grained porphyroclastic near the edge to fine-grained mylonite (grain size of ∼0.1 to 0.5 mm) toward the center (47). Elongated ultramylonite patches may be dispersed within the very center of the ductile shear zones, with grain size down to a few tens of microns as a homogeneous multiphase mixture (47).

Laboratory experiments show that the development of ductile shear zones is closely related to chemical interaction between olivine and Opx (50). While remaining to be separate phases in the upper mantle, the two minerals contain one common species (Mg,Fe)O (or MeO). Under conditions that render Opx thermodynamically stable, deformation in olivine and Opx aggregates enhances MeO interdiffusion between the two phases via the following reaction (50):where, for simplicity, the Al-containing component in Opx has been ignored. Under appropriate temperature conditions, MeO, the fastest diffusing species, diffuses out of olivine–olivine grain boundaries under stress, leaving behind new Opx grains (reaction proceeds from right to left). The newly formed equant Opx grains are substantially smaller due to limited diffusion distance, which is controlled mainly by temperature. The excess MeO resulted from forming new Opx grains diffuses through Opx–Opx grain boundaries to form new olivine grains (reaction from left to right). The newly formed grains of one phase pin grain boundaries of the other phase, significantly retarding grain growth. With increasing deformation, the nucleation and growth of the small grains produce more olivine–Opx boundaries, resulting in intense phase mixing. With sufficient deformation, most grains in the center of the shear zone are replaced with newly nucleated grains, resulting in a layer of fine-grained aggregate. In the Opx-rich shear bands, Opx grain size is smaller than olivine, typically by a factor of ∼10 (47). Such texture is consistent with experimentally deformed synthetic samples (50). Exsolution of Grt in Opx, once initiated after subduction, would further drastically weaken the ductile shear zone.

The deformation mechanism controlled by Me interdiffusion is expected to be dominant only within a narrow temperature range. Indeed, microstructures of exhumed mantle xenolith samples suggest that ductile shear localization occurs at temperatures below 1,000 to 1,100 K (47). This may be understood by considering the competing deformation mechanisms at different temperatures. At low temperatures, where the mantle is strong and ductile, interdiffusion is too slow to play a major role in deformation. At higher temperatures, dislocation creep dominates, making interdiffusion irrelevant in controlling the mechanical behavior.

Based on the above results, we propose the following model for intermediate-depth earthquakes in LBS (Fig. 6).

Fig. 6.

A proposed model for intermediate-depth earthquakes in LBS. (A) Schematics of a typical subduction zone. As the lithospheric mantle travels from the midocean ridge to the subduction zone, shear localization zones develop due to basal shear. The inset schematically shows ductile shear zones (black lines) developed prior to reaching the trench. Earthquakes in the shallow part of the system (<60 km) are not shown in this figure. (B–D) Schematics of microstructure development in preexisting ductile shear zones. Panels B, C, and D correspond to the small tilted black rectangles b, c, and d, respectively, in A. Green grains are olivine, and yellow grains are Opx. (B) Microstructure of mylonite within the ductile shear zone, showing grain-size gradient at ∼centimeter scale. The concentration of Opx in the center is due to MeO diffusion-induced grain size reduction of a previous large Opx grain surrounded by olivine, resulting in a banded structure (modified from ref. 48, with permission from John Wiley and Sons). (C) As the material descends to a certain depth, reaction is activated. Ultrafine-grained Grt forms (red dots), mostly concentrating in the center of the shear band. (D) With continued Grt exsolution, the center of the shear band loses strength to produce brittle slip (gray shaded band). Continuous shear deformation forces the shear bands to coalesce, forming larger-scale brittle shear zones, which are further enhanced by shear-induced heating, to facilitate seismic faulting. Opposing white arrows in C and D schematically represent the maximum principal stress.

Condition 1) for earthquake production, as summarized earlier, is readily satisfied because the relatively cold and strong slab mantle is stressed due to various loading mechanisms, such as gravitational pull and plate flexure, and stores large amounts of elastic strain energy.

Condition 2) for earthquake production is also met and is the main point of this paper. Preexisting ductile shear zones contain finer-grained olivine and Opx near their center, where grain size of the secondary phase Opx has been reduced to ∼0.5 mm (47). At the centimeter scale, a banded texture has been repeatedly reported in field observations (46, 47, 51). Fig. 6, modified from ref. 47, schematically shows such a band structure in mylonite. Laboratory experiments show that the directions of the banding are generally within 10 and 15° of the shearing direction (50). As the lithospheric mantle subducts to the depth where metamorphic reaction is kinetically activated, deformation mechanism undergoes a profound change, triggering embrittlement as described in the following. The metamorphic reaction occurs preferentially within the center of the preexisting ductile shear zone, where the fine-grained materials possess additional energy due to their much higher surface-to-volume ratio as well as greater strain concentration. Opx breaks down into an ultrafine-grained mixture of Grt and (Al-poor) Opx (Fig. 6). Throughout the Opx breakdown reaction, nucleation of new phases is promoted, but the nanograined nature of the multiphase mixture prevents grain growth, especially during deformation. Nanometric Grt and Opx grains engulf Opx grains as well as along the Opx–olivine phase boundaries, forming extremely thin nanometric reaction bands (e.g., Fig. 4) within the center of the preexisting ductile shear zone. Following the arguments presented for the Opx failure model above, the thin reaction bands are extremely weak and hence facilitate slip.

LBS seismicity is restricted to temperatures below ∼1,000 K (Fig. 2). Such temperature conditions are consistent with the transition from ductile behavior to metamorphism-induced embrittlement in Opx (Fig. 3) observed in our experiments, considering the tradeoff between strain rate and temperature. The LBS temperatures are also similar to the conditions for ductile shear zone formation in the oceanic lithospheric mantle before subduction, supporting the notion that seismogenic faulting tends to nucleate in preexisting shear zones in the lithospheric mantle that are ductile at geological timescales. Finally, the Grt exsolution microstructure observed in our partially metamorphosed Opx samples is compatible with microstructures in many upper-mantle harzburgitic and lherzolitic xenoliths (37, 52–54). Fine-grained Grt-rich reaction products form thin layers along Opx grain boundaries (38) and in cleavage and kink bands (52). Grt lamellae within Opx grains are preferentially parallel to {001}, {011}, and {0-11} planes of the host Opx (53). Based on the observed microstructures, it has been argued that shear stress enhances Grt exsolution in lithospheric mantle (52, 55).

All laboratory experiments are limited in sample size, which is a few millimeters in our case. We have demonstrated that syn-deformational metamorphism in Opx samples causes slip instability at such a length scale. Whether this process can operate in LBS settings and propagate to form large-scaled instability must be discussed based on fundamental physics. We reason that condition 3) will also be met in the LBS. Assuming that grain size of the reaction products is on the order of 50 nm as observed in the laboratory, a ∼104 grain size reduction corresponds to a factor of 1010 viscosity decrease for p = 2.5. As the metamorphic reaction continues, the weak nanograined materials coalesce to form a large-scale continuous gouge with inherent heterogeneity, allowing stress concentration, providing a structural condition for the initiation of large-scale slip instability. Although there are currently no experimental data on the mechanical behavior of this gouge material, we note that a three-phase mixture of olivine, enstatite (an Al-free Opx), and diopside with grain size of <500 nm displays a remarkable superplastic behavior in tensile tests under ambient conditions (56), indirectly attesting the potential of ultrafine-grained silicate minerals to behave as a lubricating gouge. Upon the initiation of slip, frictional heat along the very thin slip zone will enhance further exothermic metamorphic reaction, and the positive feedback with the consequent release of latent heat is expected to cause unstable seismic slip (Fig. 6). The key difference of this mechanism from the hypothesis of thermal runaway is the presence of localize shear zone (fault) before shear is initiated, which not only allows frictional heating to be more efficient in causing dynamic weakening but also enables propagation of the rupture front.

As shown in Fig. 1, the initiation depth of reaction depends on the Al content in Opx before subduction begins. Kinetic onset of Grt exsolution in LBSs is a complex process with tradeoffs among temperature, Al content, and strain rate. This may explain the wide range of LBS age, depth, temperature, and separation from UBS in subduction zones around the globe. As mentioned earlier, the UBS can be explained by dehydration embrittlement in the subducting crust and is not addressed in our model. As slabs descend beyond the depth where UBS and LBS merge, our proposed metamorphic-facilitated mechanism continues to operate insofar as the slab is cold enough to maintain slow kinetics for the activation of reaction , with increasing Grt volume fraction. Since the equilibrium Grt composition varies with depth, Grt will remain metastable, and new Grt with different composition will continue to exsolve. This may explain the presence of seismicity to 300-km depths in colder slabs.

Materials and Methods

LBS Seismicity in Various Subduction Zones.

Numerous studies have reported earthquake distributions in DSZs. However, most hypocenter data in public data repositories have large uncertainties, making it difficult to distinguish LBS from UBS events (57). Here, we use primarily the earthquake location data of LBSs from ref. 2, which has addressed issues that may affect determination of hypocenters and developed criteria and procedures to locate earthquakes in DSZs (57). The drawback is that the events plotted are only a small subset of LBS earthquakes, because the criteria, technique (based on the time difference between the depth phase pP and the direct P), and station locations used by these authors have significantly limited the number of earthquakes that could be confidently identified. Seismic activity patterns (rise and fall in event frequency), however, should be more or less preserved. displays all nine subduction zone corridors individually in terms of LBS earthquake distribution and temperature profiles (red and blue curves). The temperature profiles are along the centerlines of the LBS, mostly based on the slightly updated slab thermal models of Wada and Wang (18). For New Britain, thermal models are from this study () using the approach identical to that of ref. 18. Considering the finite thickness of the LBSs (10 to 20 km), temperature variations in each LBS are estimated to be on the order of ±100 K (the shaded bands along the red curves). summarizes some of the key parameters for the subduction zones.

Experimental Starting Materials.

The starting material was from a fresh, coarse-grained, granulite xenolith specimen from Chifeng, Inner Mongolia, China (58). The granoblastic specimen, whose protolith was likely an orthopyroxene gabbro (58), contained Opx, diopside, and Plg as major phases, with biotite, amphibole, and Ilm as minor phases. The rock was disaggregated, cleaned, and separated into monomineralic grains. Opx grains were collected and ground to an average grain size of 23 to 38 μm and resintered into cylinders at 3 GPa and 1,273 K for 4 h using the 10 MN multianvil apparatus at GeoSoilEnviroCARS (GSECARS), Advanced Photon Source, Argonne National Laboratory (59). The resulting nominally dry Opx rocks were ∼2.1 mm in diameter and ∼3 mm in length. Chemical analyses on these sintered samples yielded an FeO content of ∼27 mol% and an Al2O3 content of ∼2.07 mol% (). This iron content is higher than what is expected in Opx (∼10 mol%) in subduction slabs and will shift major phase boundaries and isopleths shown in Fig. 1 somewhat (60), but neither the amount of Fe nor the Fe3+/ΣFe ratios affect the general trend of metamorphic reactions (61, 62). The average grain size was 23 to 38 mm, and water content, determined by infrared spectroscopy was about 63 parts per million (). The samples were therefore considered dry.

Experimental Procedures.

The experimental setup in the deformation device known as the D-DIA (63), with in situ AE and X-ray monitoring at beamline 13-BM-D of the GSECARS facilities, is shown in . The D-DIA uses six anvils, four of which are made of tungsten carbide (WC) and two are made of sintered diamonds (SD), which are transparent to X-rays, allowing complete diffraction Debye–Scherrer rings to be recorded on an area detector. The cell assembly used is illustrated in . Thermocouple was not used during deformation. Temperature estimate was based on temperature–power relations obtained from several calibration runs before the experiments. Variation in these calibration runs suggested that uncertainties based on power calibrations were between 50 and 100 K. Gold foils were placed on the top and bottom of the sample as strain markers so that strain can be monitored in situ using X-ray radiography. The cell assembly was similar to that of ref. 35, with no special efforts to buffer the oxygen fugacity in the cell, as it has been demonstrated that for Fe-bearing Opx, the creep behavior exhibits no dependence on oxygen fugacity over 8 log units, with log fO2 from −2 to −10 (64). Pressure and temperature conditions were within the general stability field of Opx (Table 1), but Al content in Opx would change according to the phase relations (60, 65). In each experiment, hydrostatic pressure was first applied to the sample and maintained at a desired value, and then the sample was heated to the targeted temperature. Deformation was carried out after annealing at the target temperature for 30 min to minimize effects of cold compression and possible water absorption during cell preparation. Heating power was monitored continuously throughout the experiment.

The incident X-ray beam was collimated to 200 × 200 μm by WC slits. A Rayonix SX-165 area detector was used for X-ray diffraction (XRD), with monochromatic radiation at 45 keV (wavelength λ = 0.2755 Å). Typical XRD data collection time was 300 s. Detector orientation relative to the incident beam was calibrated using a CeO2 powder standard using the data collection and analysis software Dioptas (66). Stress analysis was carried out on the two-dimensional XRD patterns using the Multifit–Polydefix software (67). Possible detector center drift during the experiments was corrected in the stress analysis by adjustment of the diffraction ring center in Multifit–Polydefix.

Sample axial strains were measured using radiography with a large beam (2.0 × 3.5 mm) by opening up the incident beam WC slits. Each radiograph was acquired for 8 s using a Ce-doped yttrium aluminum garnet (YAG) scintillator and a Point Gray charge-coupled device (1.3 μm/pixel). Sample bulk axial strain was defined as ε = ΔL/L, where ΔL = L – L, L is the sample length during deformation, and L is the length just before starting the deformation in each run under pressure and temperature. Positive strain corresponds to axial shortening. Typical uncertainty in strain measurement was ∼0.2%. During a deformation cycle, diffraction and imaging data were collected alternatively by switching the diffraction/imaging optics. The process was automated with minimal human intervention. The total time required for a pair of strain and stress measurements was about 6 min, defining the time resolution in real-time faulting strain detection. During deformation, the differential rams were kept at a constant speed. It usually took 10 to 15 min from the starting point for stress to build up on the sample. As a result, sample strain rate increased slowly until about 2 to 3% strain, beyond which the strain rate was more or less constant, as exemplified by Fig. 2.

Stress Analysis.

Under hydrostatic pressure, the d-spacing for a given set of lattice planes (hkl) in a random powdered sample gives perfectly circular Debye–Scherrer rings on an area detector. Under a macroscopic differential stress, the lattice strain is a function of orientation of the diffracting planes relative to the stress geometry (i.e., the azimuth angle), resulting in pseudoelliptical Debye–Scherrer rings. The pseudoellipticity of powder diffraction rings recorded by XRD is a linear elastic function, a proxy for the differential stress the polycrystalline sample supports. The diffracting plane azimuth angle ψ was defined as the angle between the minimum compression axis and the diffracting plane normal (ψ = 90° corresponds to the maximum compression axis).

Each XRD pattern was integrated over 5° slices of azimuthal angle, resulting in 72 2θ-versus-intensity spectra, using macros in Fit2D (68). At any azimuth angle, simple elastic theory, assuming no lattice preferred orientation and under an axial-symmetric stress field, predicts the relation between measured d-spacing d and lattice strain Q to bewhere d() is the d-spacing under hydrostatic pressure. For each plane hkl, d and ψ were measured from the XRD pattern (sin ψ = cos θ × cos η), where η is apparent azimuth angle on the detector and θ the diffraction angle. Q and d() were extracted by fitting d versus ψ according to Eq. .

Differential stresses t were then calculated for plane hkl according towhere the “effective modulus” G was calculated for the elastic stiffness tensor (C) for orthoenstatite (69) with experimentally determined pressure and temperature derivatives (69, 70). We used the Reuss (iso-stress) assumption to calculate differential stress for each hkl, believed to be closer to reality than the Voigt (iso-strain) assumption. The flow stress quoted in this work is the unweighted average of differential stresses, t, calculated from the lattice strains Q of Opx. This approach may be biased by the availability of diffraction lines and has no physical ground; however, it is used in many studies as there is still no complete theory for evaluation of macroscopic stress from powder diffraction data.

Unit-cell volumes of Opx were calculated using four crystallographic planes of 321, 610, 421, and 531, near the “magic angle” sin2ψ = 1/3. Confining pressures were calculated following the third-order high-temperature Birch–Murnaghan equation of state, with 1 atm thermal expansion data for Opx (69) and pressure and temperature derivatives of the bulk modulus (69, 70). Estimated uncertainties in pressure were about 0.2 GPa. Stress data are summarized in ref. 71.

AE Monitoring.

The AE setup () was similar to that reported in an earlier study (72). Six piezoelectric lead–zirconium–titanate transducers (10 mm in diameter, 0.5 mm in thickness) were mounted on the back side of the D-DIA anvils to monitor AEs arising from the sample. Resonant frequencies of the transducers were ∼2 MHz. The back end of the anvils was mirror polished to create a perfectly planar interface for transducer mounting. For electrical insulation, each transducer was glued onto a 0.25-mm thick Al2O3 disk. Six amplifiers filtered analogically (0.3 to 5 MHz) were used to amplify the AE signals (at 60 dB). AE waveforms were acquired in triggered mode simultaneously from all six transducers using a digital multichannel oscilloscope. When the signal was higher than a prescribed trigger level (typically 150 mV) on a given channel, the waveform was recorded. A typical triggered waveform contained 2,000 to 8,000 points at a sampling rate of 50 MHz. shows an example of the AE waveforms recorded. Due to faster wave speeds of diamond, arrivals on the two SD anvils were significantly earlier than in the other four WC anvils. To calibrate travel times through the anvils, a 500-V pulse was sent from one sensor, with P-wave signals recorded by the opposing transducer. Travel times along the three pairs of anvils were then used to build a velocity model, which was used to locate AE events based on differential arrival times.

AE Location using HypoDD.

AE events were relocated using the seismological algorithm hypoDD (73). We first cross-correlated waveforms between events for each channel by using a time window of 6 μs, starting 1 μs before the first P-wave arrival. The cross-correlation (CC) results were visually inspected. For all six channels, CC coefficients greater than 0.8 could identify similar waveforms and align them correctly. This procedure yielded multiple groups of events. Within each group, events had similar waveforms. Differential arrival times of event pairs within each group were then used in the hypoDD algorithm to relocate the AE events. The hypoDD algorithm determines relative location and origin time between any two events for each group. In the analysis, a P-wave velocity of 8.4 km/s was used for Opx, and the amounts of Grt and Cpx were neglected, as they were less than 5 vol% according to SEM.

For each pair of events, differential P-wave arrival times of at least four channels were needed. To avoid the situation where all four channels were located in a plane, which would lead to a tradeoff between the origin time and relative location in the direction perpendicular to the plane, we required differential P-wave arrival times of at least five channels for the pair. shows the relocation results of faulted specimen D2071. AE location data are summarized in ref. 71.

Microstructural Characterization.

Electron microscopy was conducted at Argonne Laboratory Microscopy Center. SEM was carried out using a Hitachi S-4700-II SEM. TEM specimens were prepared using a focused-ion beam from selected areas in the samples. High-resolution TEM (HRTEM) images were obtained using a chromatic aberration-corrected TEM equipped with an image corrector to correct both spherical and chromatic aberrations. Correction of spherical and chromatic aberrations enabled the microscope to reach a spatial resolution less than 0.1 nm (measured by Young’s fringes) at 80 kV. A low-dose (∼5 eÅ−2s−1) HRTEM technique with an accelerating voltage of 80 kV was used to minimize electron beam damage to the sample. High-angle annular dark-field imaging and X-ray energy-dispersive spectroscopy mapping near cracks were performed by an FEI Talos F200X TEM. Microstructure data are summarized in ref. 71.