Research ArticleGEOPHYSICS

A laboratory nanoseismological study on deep-focus earthquake micromechanics

See allHide authors and affiliations

Science Advances  21 Jul 2017:
Vol. 3, no. 7, e1601896
DOI: 10.1126/sciadv.1601896


Global earthquake occurring rate displays an exponential decay down to ~300 km and then peaks around 550 to 600 km before terminating abruptly near 700 km. How fractures initiate, nucleate, and propagate at these depths remains one of the greatest puzzles in earth science, as increasing pressure inhibits fracture propagation. We report nanoseismological analysis on high-resolution acoustic emission (AE) records obtained during ruptures triggered by partial transformation from olivine to spinel in Mg2GeO4, an analog to the dominant mineral (Mg,Fe)2SiO4 olivine in the upper mantle, using state-of-the-art seismological techniques, in the laboratory. AEs’ focal mechanisms, as well as their distribution in both space and time during deformation, are carefully analyzed. Microstructure analysis shows that AEs are produced by the dynamic propagation of shear bands consisting of nanograined spinel. These nanoshear bands have a near constant thickness (~100 nm) but varying lengths and self-organize during deformation. This precursory seismic process leads to ultimate macroscopic failure of the samples. Several source parameters of AE events were extracted from the recorded waveforms, allowing close tracking of event initiation, clustering, and propagation throughout the deformation/transformation process. AEs follow the Gutenberg-Richter statistics with a well-defined b value of 1.5 over three orders of moment magnitudes, suggesting that laboratory failure processes are self-affine. The seismic relation between magnitude and rupture area correctly predicts AE magnitude at millimeter scales. A rupture propagation model based on strain localization theory is proposed. Future numerical analyses may help resolve scaling issues between laboratory AE events and deep-focus earthquakes.


Deep earthquakes, that is, those with hypocenter depths greater than ~70 km, constitute about a quarter of all recorded events, with moment magnitudes greater than 5 in the International Seismological Centre catalog (1). They occur in association with convergent margins, defining planar regions known as the Wadati-Benioff zones (2), which delineate the cold, down-going cores of subducting slabs. Deep earthquakes pose significant seismic hazards and have played a key role in illuminating the structure of Earth’s mantle and core (3). The physical processes that permit the occurrence of deep earthquakes remain poorly understood. Brittle-frictional processes, which produce shallow seismogenic faults through open cracks and frictional sliding via mode-I and mode-II/III fractures, are suppressed by conditions at depths greater than ~70 km, because high pressure inhibits opening cracks and frictional sliding on existing fractures (4, 5). Furthermore, higher temperatures at depths render rocks more ductile (68).

Several mechanisms have been proposed to explain deep faulting. High fluid pressure resulting from dehydration reactions may reduce the effective normal stress and allow movement on preexisting faults (911). A number of dehydration reactions have been shown to trigger rock failure in the laboratory (1214). Transformational faulting, in which a metastable mineral, such as olivine, transforms into a denser phase, such as wadsleyite and/or ringwoodite, triggering instability by strain localization and stress concentration in the mantle transition zone (MTZ), has been suggested to explain deep-focus earthquakes (DFEQs; depths, >350 km) (1522). Once instability is triggered, shear failure may propagate as ductile faulting or plastic instabilities (2325) rather than as brittle failure, and rupture propagation may be facilitated by melt production at the tip of the fracture zone (26) or failure of the surrounding mantle materials (27). In all of these proposed mechanisms, there is a marked lack of knowledge as to how fracture self-organizes and propagates, leading to ultimate failure. Without these details, it is difficult to develop realistic physical models for the failure process and scale laboratory observations to geological activities.

Recently, we conducted transformational faulting experiments on Mg2GeO4 olivine with controlled deformation and phase transformation by combining in situ synchrotron x-ray diffraction (XRD) and imaging along with acoustic emission (AE) monitoring (see Materials and Methods and fig. S1 for details) (28). The onset of the olivine-spinel transformation was observed to correlate closely with an abrupt jump in AE occurring rate and a marked stress drop. All the AE events were located within the sample. X-ray microtomography (XMT) on the recovered samples revealed conjugated faults. Scanning and transmission electron microscopy (SEM and TEM) confirmed the XMT results and allowed more details to micrometer scales. These results showed that the AE events were linked to the final failure. However, locations of the AE events were too poor to yield clear spatial correlation with observed fault zones by XMT and SEM (28). Here, we adapt a suite of techniques developed for seismology for quantitative data analysis. We apply waveform cross-correlation (CC) and a double-difference relocation algorithm [hypoDD (29)] for accurate hypocenter relocation of the recorded AE events. Together with first P-wave amplitude analysis, we determine moment tensors of the events and improve both event location and detection sensitivity by a factor of ~10. The relocated events show excellent spatial correlation with faults observed in XMT and SEM, permitting an examination of temporal distribution of events, with unprecedented details for the faulting process. During deformation, the b values in the Gutenberg-Richter relation decrease from >2 to ~1 shortly before the ultimate failure of the specimen. The empirical relation between seismic moment magnitude and rupture area (30) is shown to be applicable down to laboratory millimeter scales. These observations strongly suggest that faulting is a self-affine phenomenon and that rupture initiation, nucleation, and propagation are scalable from laboratory samples to large tectonic units, over ~19 orders of magnitude in elastic energy release. On the basis of well-known shear localization theories (31, 32), we propose a new model for the rupture nucleation and propagation process.


XMT on fault distribution in the sample

Figure 1 shows two XMT images of sample D1247, which failed at 4-GPa effective mean stress. The original 2-mm-diameter and 3-mm-long sample underwent significant axial deformation to become ~2.5 mm in diameter and ~2 mm in height. Figure 1A is a view with the sample tuned to 85% transparent, to show faults (the linear dark traces) with various orientations and locations. Faults form conjugated configurations. The three-dimensional (3D) fault distribution is shown in movie S1. Because of the spatial resolution, tightly closed faults (such as those shown in fig. S2) cannot be resolved. Figure 1B is a vertical cut-through view, where parts of two conjugated faults are visible. A 3D display of the “ortho-cut” view (that is, three sections orthogonally intersecting at the midpoint of the sample) is given in movie S2. The two fault planes (hereafter referred to as faults 1 and 2) are digitized to determine their precise locations and orientations, to compare with relocated AE events described in the following section.

Fig. 1 XMT images of faults and correlation with relocated AE events.

(A) View of 3D tomographic reconstruction of the recovered sample, with the bulk tuned to 85% transparent, revealing sets of conjugated fault traces (dark lines). (B) Cut-away view of the same image showing conjugated faults. The faults are digitized and plotted in (C) and (D). For 3D representations of (A) and (B), see movies S1 and S2, respectively. (C) Group 1 events superimposed with digitized faults to show their spatial correlation with fault 1. (D) Group 2 events superimposed with digitized faults showing correlation with fault 2. 3D representations of spatial correlation between group 1 and 2 events and the faults are shown in movies S3 and S4, respectively.

AE event relocation

A total of 593 events were recorded in sample D1247 in triggered mode. We measure differential arrival times of all the events at each of the six transducers (or channels) by waveform CC (transducer locations are given in fig. S1), using a time window of 6 μs starting at 1 μs before the first P-wave arrival. With the requirements described in Materials and Methods, relative locations of 523 events are determined by hypoDD (29) in 55 separate groups. Selected waveforms of the two largest event groups, which contain 121 and 99 events, respectively, are shown in fig. S3 (hereafter referred to as groups 1 and 2). A comparison between the two groups shows that the general features of the waveforms are all similar for channels 1 and 2 (on the bottom and top anvils) but differ noticeably for the horizontal channels (fig. S3).

Results of relocation for events in groups 1 and 2 are shown in Fig. 1 (C and D, respectively), with symbol size proportional to their scalar moment M0 (in Newton meter). All the events are located within the sample, with estimated location uncertainties of about 10 μm, more than one order of magnitude improvement over the technique used in our previous study (28). Faults 1 and 2 are overplotted with the events, showing distinct spatial correlation of group 1 and 2 events with faults 1 and 2, respectively (Fig. 1, C and D). The 3D relations between faults 1 and 2 and groups 1 and 2 are displayed in movies S3 and S4, respectively.

Event source parameters from nanoseismological analyses

We use a grid search method to determine moment tensors of AE events based on amplitudes of the first P-wave pulses at transducers 1 to 4 (transducers 5 and 6 are excluded because of their abnormal instrument responses). Specifically, we determine source parameters of scalar moment M0, isotropic (ISO) component strength ζ, compensated-linear-vector-dipole (CLVD) component strength χ, and strike φ of the fault plane (see Materials and Methods). We assume normal faulting; thus, fault slip rake λ = −90° (see Materials and Methods), as deformation was axial shortening in the vertical direction (z axis) in the experiment. We find that the data have no sensitivity to the dipping angle of fault planes, so we fix it to 45°. Figure S4 shows grid search results for AE event 277. Both ζ (=0.20 ± 0.02) and χ (=−0.01 ± 0.02) are well resolved. Strike φ is also well resolved (24 ± 2°) except for the ambiguities of ±φ and π ± φ due to the z-axial symmetry.

Figure S5 shows distributions of ζ, χ, and φ of all 593 events. Because of the quadrant ambiguity in φ, only values in the first quadrant are used. Whereas the CLVD parameter has a broad distribution from −0.3 to 0.5, most of the ISO parameter values concentrate in a narrow range from 0 to 0.2. Distribution of strike exhibits two peaks, one from 30° to 40° and the other from 50° to 70°. The first peak coincides with the strike distribution peak of events in group 1, and the second one matches the strike distribution of events in group 2. On the other hand, neither the distribution of ISO parameter nor the CLVD parameter shows differences between groups 1 and 2.

Synthetic seismograms of the events are then made using the whole-space displacement solution (33) and the moment tensors from the grid search. They are compared with observations at transducers 1 to 4 after low-pass–filtered at 3 MHz in Fig. 2, where observed displacements are calculated on the basis of transducer piezoelectric constant and amplification (see Materials and Methods). Synthetic waveforms fit the first P-wave pulses very well. The long-period near-field and intermediate terms are strong in synthetics but are difficult to identify in the data because of the narrow bandwidth of the sensors and wave reverberations in the instrument. The main difference between the two events is the relative amplitudes of channels 3 and 4, resulting in different strikes of fault planes (also demonstrated in fig. S5).

Fig. 2 Examples of synthetic waveforms (red curves) based on moment tensors from the grid search.

Events 277 and 288 are from groups 1 and 2, respectively. Only the first P-wave arrivals are modeled, which show good agreement with observed waveforms after low-pass–filtered at 3 MHz (black curves). Complex waveforms beyond the first P-wave arrival are most likely due to the narrow bandwidth of the sensors and wave reverberations in the instrument. The relative amplitude difference between the horizontal channels 3 and 4 (CH3 and CH4, respectively) is controlled by the strike of fault planes (also see fig. S5).

Moment magnitudes (MW) of the 593 events range from −9.5 to −8 (Fig. 3); events with MW > −8 are clipped due to the predefined amplitude limit in the detection system. The number of events with MW < −9 is incomplete due to the prescribed threshold in event detection to avoid false triggering by noise.

Fig. 3 Relation of AE event counts (in log scale) versus moment magnitude.

Dark gray histograms are from the 593 events above the detection threshold (MW between −9.5 and −8.0). Light gray histograms are based on the refined event detection algorithm based on waveform cross-correlation for events below the threshold and scaling of saturated events based on duration (see text for details). Detectable magnitude range has doubled. This yields a more tightly constrained b of 1.5. The relation for b = 1.0 is also plotted for reference (dashed line).

Event detection by waveform CC

To capture events missed in the triggered event catalog, we apply the sliding CC (SCC) detection technique (34) to the continuous waveform data recorded in channel 2. We use the waveform of event 437 (MW = −8.3) as a template. A CC coefficient threshold of 0.75 is found best for eliminating false detections caused by noise and glitches in the continuous data. We detect 5818 events by SCC, about 10 times more than the number of triggered events. Waveforms of representative detected events of different magnitudes are shown in fig. S6. For events with unsaturated waveforms, their moment magnitudes are calculated using the amplitude ratios with respect to the template event. The smallest event thus detected has been extended down to MW = −9.8.

There are nine large events whose waveforms were clipped by the recording system (fig. S7), so their magnitudes cannot be estimated from amplitudes. We measure the widths of direct P-wave pulses, ranging from 0.5 to 1.5 μs, as the source rupture durations. Their moment magnitudes are then estimated using ratios of their source duration to that of the template event according to the scaling relationship log(duration) ∝ 0.5MW (30, 35). The largest events are found to have MW = −7.1. Between MW of −7.0 and −9.5, event distribution defines a b value of 1.5 (Fig. 3) in the Gutenberg-Richter statistics log10(N) = ab MW, where N is the number of earthquakes with magnitudes ≥MW, and a and b are constants (36).

Spatial-temporal variations of AE events

With event occurrence times also refined by hypoDD, it is straightforward to examine them in both space and time. Figure 4 shows spatial-temporal evolution for events in groups 1 and 2. Two plots are provided for each group. One is a time series showing distances of events from the fault plane (Fig. 4, A and C), and the other is a view normal to the XMT fault planes (Fig. 4, B and D). Movies S5 and S6 are time sequences of events in groups 1 and 2, respectively. Events occurring in the early stage of the deformation/transformation process show no clear spatial correlation (Fig. 4, A and B): They appear to occur at all distances from the final fault plane. With increasing deformation and transformation, more events occur near the final fault planes with increasing magnitude. Shortly after the 2000-s mark, several very large events took place, with their waveforms saturated. As a result, these large events cannot be located precisely and therefore are not shown in Figs. 3 and 4 (A to D). These event sequences, revealed for the first time in laboratory experiments, provide a unique opportunity to follow rupture nucleation and propagation process.

Fig. 4 Event evolution throughout the deformation/transformation process.

(A) History of group 1 events projected along the strike of fault 1, whose average trace is shown as the red dashed line. Most of the large events after 1700 s are located in fault 1. (B) Group 1 events projected onto fault 1. Size of the events is proportional to their MW to minimize overlap. (C and D) Similar plots for group 2 events and fault 2. (E) Moment magnitudes of 5818 events detected from the continuous record in the time sequence as they occurred. Superimposed are pressure (blue), strain (red), and differential stress (green). Stress peaked near 1700 s and then dipped when large AE events stared bursting. In addition, plotted are the b values over 600-s time windows during the deformation/transformation process. A continuous decrease trend in b value is evident. After the 1700-s mark, b decreases rapidly from 1.5 to 1.0, suggesting that later events are preferentially distributed in a planar fashion.

Figure 4E is a time series of all detected events on channel 2, showing that larger events occur shortly after the 2000-s mark during the deformation/transformation period. The occurrence of the largest events is accompanied by a sudden increase in strain rate and a ~100-MPa drop in differential stress. We use a moving time window of 600 s with 50% overlap with the neighboring windows to estimate the b values defined by the events in each window. The results show that b stays in a range of 2 to 3 for the initial two-thirds of the duration and then decreases rapidly to 1 as large AE events begin to occur.

A closer examination of the triggered events reveals that variation of ζ exhibits a tightening trend with event magnitude (Fig. 5A). On the other hand, event rate appears correlated with proximity of the criticality point, which occurs at approximately 1650 s (vertical arrow in Fig. 5B). Between 400 and 1500 s, the average event rate is 10.3 per min; between 1500 and 2250 s, event rate more than doubles at 24.5 per min, followed by a sudden drop. The ζ parameter again shows a tightening trend with time after the criticality point (Fig. 5B).

Fig. 5 Variation of ISO strength (ζ) and event rate as a function of MW and time.

(A) Events below (red circles) and above MW ≈ −8.5 (blue circles) display different spread in ζ. Dashed lines are a guide to the eye. (B) The deformation/transformation process is divided into four regions based on event rates (given on top of each colored region, in number of events per minute). Event rate reaches maximum (24.5/min) between 1500 and 2100 s, presumably due to the nucleation and propagation of NSBs. The critical point is indicated by the vertical arrow. Note the rapid tightening of ζ values after ~2000 s (dashed lines are a guide to the eye), followed by a “quiescent period” between ~2150 and 2550 s (event rate drops to 4.0/min), before the sample failed shortly after ~2700 s.

Faulting microstructure

As shown in fig. S2, both recovered samples contain numerous visible faults. The samples are predominantly olivine, which appears dark gray in the backscattered electron (BSE) images. In all of the BSE images, the cylindrical axis of the samples is vertical. Grains of pyroxene (MgGeO3) are dispersed throughout and appear bright against the olivine matrix. In sample D1247, large conjugated faults are apparent (fig. S2A), with smaller visible fault segments. Sample D1253 also contains multiple visible fault segments along the northwest-southeast (NW-SE) diagonal (fig. S2B). Some fault segments appear discontinuous when viewed at large scales (several gaps between the fault segments are indicated by red arrow pairs in fig. S2B). At slightly high magnification, numerous smaller faults are observed, most of which are under 0.2 to 0.3 mm long [thin dark traces along SE-NW and NE-SW (northeast-southwest) in Fig. 6A, corresponding to the white box in fig. S2A], subparallel to the conjugated main faults. These “faultlets” segments and the larger fault segments are often associated with long and narrow bands, which appear slightly brighter than olivine, but darker than pyroxene (Figs. 6B and 7). XRD, SEM, and TEM indicate that these bands consist of nanograined Mg2GeO4 spinel, that is, the transformation product (28). In what follows, we refer to these bands as nanoshear bands (NSBs).

Fig. 6 Low-magnification SEM BSE images showing the distribution of faultlets and fault segments.

Maximum compressive stress is approximately vertical. Pyroxene grains are the bright and round crystals against the dark gray olivine matrix. Spinel has ~8% higher density than olivine and appears slightly brighter. (A) Area corresponding to the white box in fig. S2A. The two macroscopic faults (faults 1 and 2) are those indicated in Figs. 1 and 4. Numerous faultlets (thin and dark traces; a few are indicated by white arrows) are present, with lengths of <0.2 to 0.3 mm primarily subparallel to either fault 1 or fault 2, in conjugated relations. Their orientations are consistent with AE event distribution with time, as shown in Fig. 4. (B) Area corresponding to the white box in fig. S2B. Tips (white arrows) of two fault segments (faults a and b) are connected by a long NSB. Near horizontal cracks are due to decompression.

Fig. 7 BSE micrographs showing the presence of long NSBs connecting visible fault segments.

(A) Image corresponding to the white box in Fig. 6B, showing the long and continuous NSB in the NW-SE direction (red arrows). Note a bifurcated of NSB (white arrows), about 20° away from the main branch. Near-horizontal cracks were due to decompression. (B) Image corresponding to the red rectangle in fig. S2B. The “gap” between the upper left and lower right fault segments (white arrows) is shown in (C) and (D), each consisting of two SEM images. (C) An NSB runs from upper left to lower right (indicated by the red arrows). One pyroxene grain has been sheared into two halves, with an offset of about 30 μm (white double-headed arrow). (D) The same NSB in (C) continues into this image. Overall, this NSB is more than 150 μm long, connecting the two fault segments in (B). (E) Composite TEM image from an FIB foil cut at the white line in (D), showing that the NSB is curvy at depths, and passing through olivine grain G2 and merging into the boundary between G3 and G4.

Numerous long NSBs are observed to connect adjacent fault segments. Figure 7A, the area corresponding to the white box in Fig. 6B, shows a long (>50 μm) NSB (red arrows) bifurcating into a new branch (white arrows), at an angle of about 20° from the main branch. Figure 7B shows the area corresponding to the red box in fig. S2B, with two higher-magnification SEM images shown in Fig. 7 (C and D). The very long and continuous NSB is present (red arrows). A pyroxene grain has been sheared so severely that it broke into two halves, with the top half offsetting by ~30 μm (white double-headed arrow in Fig. 7C). This NSB continues from Fig. 7C to Fig. 7D, eventually connecting the visible fault at the lower right corner of Fig. 7D. In total, this NSB is more than 150 μm long. Because the average grain size of the samples is ~50 μm, it must have cut through multiple grains.

Figure 7E is a composite TEM image taken from a focused ion beam foil extracted from the location indicated by the white line (“FIB”) in Fig. 7D. The imaged surface of Fig. 7D corresponds to the left side of Fig. 7E, with increasing depth to the right. Multiple olivine grains are present, some of which are labeled (G1 to G5). For clarity, several grain boundaries are marked by the dashed white arrows. The NSB (indicated by red arrows) is at least 7 to 8 μm deep from the surface of Fig. 7D. It passes through grain 2 and curves at the boundary between G2 and G3 and then continues into the boundary between G3 and G4. Electron backscatter diffraction (EBSD) (37) confirms that NSBs propagate across boundaries of grains with significantly different orientations (fig. S8).

Thus, the NSBs are ribbon-like in shape. The band angle θ (the angle between the plane normal of an NSB and the maximum compressive stress) can be estimated on the basis of the traces of the NSBs in the SEM images. In general, the NSB traces are inclined from the imaged surface; hence, only an apparent band angle can be observed, which is lower than the true band angle. The apparent angles we observed are all <50° (Fig. 7 and fig. S9). When the “strike” of an NSB is horizontal and parallel to the imaged surface, the apparent band angle is close to 0° (for example, the horizontal NSBs indicated by green arrows in fig. S9B). Only when an NSB is viewed edge-on will the true band angle be equal to the apparent angle. Therefore, we estimate the maximum true band angle to be approximately 50°.

Lenticular-shaped nanospinel pockets are also present (fig. S9). These are referred to as “anticracks” in previous studies (15, 22), in an analogy to the closure of voids in rocks due to pressure solution (38) or compaction (39). However, these anticracks only have ~8% volume reduction and do not fully collapse, as in the case of classical definition of an anticrack. Therefore, we will refer to them as localized volume reduction (LVR). LVRs are predominantly isolated small pockets in the metastable olivine matrix, with their shortest dimensions approximately aligned with maximum compressive stress (vertical). They are less frequently observed and, when present, are not closely associated with the visible faults.


Mechanism of fault self-organization and dynamic propagation

In previous studies of transformational faulting in Mg2GeO4 at pressures of 1 to 2 GPa (15, 22, 40), two types of precursors are reported preceding macroscopic faults. The first are anticracks, which nucleate on grain boundaries, intracrystalline flaws, and phase boundaries (15, 22). They are believed to dominate at high temperatures and/or low strain rates (40). The second precursors, believed to be dominant at low temperatures and/or high strain rates, are slip bands that form as a result of dislocation pileups in the metastable olivine in the (110) and (010) planes (40). Both precursors are filled with nanograin spinel. These results provide important information on the onset of rupture nucleation process before failure but do not offer details as to how initial ruptures coalesce and self-organize into larger faults, leading to macroscopic failure.

The microstructure of our recovered samples is dominated by NSBs, whose lengths are several micrometers and above, and faultlets at somewhat greater length scales. The range of AE moment magnitude (M0), from ~10−6 to ~0.05 N·m, provides useful estimates on rupture areas (S) and slip distance (D), according to M0 = G × S × D, where G is the shear modulus. With G ≈ 100 GPa, M0 ~ 10−6 N·m corresponds to S several micrometers in linear dimension and D on the order of 1 μm, whereas M0 ~ 0.05 N·m corresponds to S of hundreds of micrometers in size and D of tens of micrometers. These estimates agree with microstructural observations, suggesting that NSBs and faultlets are the cause of recorded AE activities (also recall that essentially no triggered AEs were recorded before the sample temperature was jumped from 973 to 1173 K, for example, Fig. 4).

The NSBs and faultlets are likely the growth products of the precursory slip bands and somehow remained planar in shape throughout deformation. Most of the faultlets are subparallel to NSBs (fig. S9B), and all of the NSBs observed have more or less the same thickness of ~100 nm, regardless of their lengths. In a sample deformed in the Griggs apparatus (41), even when a fault reached several millimeters in length, the nanospinal-filled fault zone still remained to below ~100-nm thickness. What is the mechanism controlling the propagation and organization of the NSBs into macroscopic faults with this limited thickness?

In situ AE monitoring again provides critical firsthand information. With increasing deformation, AEs spread in subparallel plane-like distributions, forming conjugated groups (Fig. 4). Magnitude of the AEs grows larger as the sample approaches and shortly passes the criticality point. This corresponds to increasing slip distance D and explains the large number of faultlets and fault segments: The faultlets are a result of propagating NSBs, which connect heads and tails of adjacent faultlets and fault segments, reorganizing them into larger faults (Fig. 7). Most of the NSBs (Fig. 7 and Fig. 8) and faultlets (Fig. 6) are subparallel to the macroscopic conjugated faults (Fig. 6A) and to the near-planar AE distributions (Fig. 4). Frictional resistance of the NSBs has been estimated (28, 41). The effective friction coefficient of the nanograined spinel fault gouge is <0.01, and the strain rates across the narrow fault gouge zone can be as high as 103 to 104 s−1. Therefore, the propagation of NSBs is seismogenic and is responsible for the observed AEs.

Fig. 8 BSE micrographs showing NSBs initiating from stress concentrations.

(A) Two LVRs with NSBs radiating from the corners (white arrows). The NSB originated from the LVR on the right has sheared the pyroxene (pyx) grain with an offset of ~3 μm. (B) Two NSBs radiating from two pyroxene grains. The originating points are indicated by the white arrows.

Why do the NSBs remain planar with a constant thickness and certain orientation relations to the stress field as they pass across individual olivine grains? We note that orientations of the NSBs are consistent with theoretical predictions of strain localization for common pressure-sensitive materials (31, 42). Under axisymmetric stress conditions, angle θ of a shear band normal relative to the maximum compressive stress is (31, 42)Embedded Image(1)whereEmbedded ImageHere, ν is the Poisson’s ratio, β is the ratio of increments of inelastic volume strain to inelastic increments of shear strain during deformation, and μ is the friction coefficient of the shear band. According to the theory, compaction bands with θ = 0° form when Embedded Image, and the material hardens on further incremental deformation. When Embedded Image, shear bands form with 0° < θ < 90°, and the material weakens with continued incremental deformation (42).

The treatment of strain localization only requires a pressure-sensitive material (volume compaction or dilation upon compression) and is quite general. Volume reduction due to a phase transformation is one of the special cases. Because μ < 0.01 (28, 41) in the NSBs, θ may be estimated at −0.07 ≤ β + μ ≤ 0.01 (as −0.08 ≤ β ≤ 0), yielding θ ≈ 40° to 41° (for ν = 0.3), if the material is ISO. Olivine is anisotropic, and microfaults in our samples generally have apparent θ angles less than 50° (Fig. 7). The very high strain rates estimated in the NSBs suggest extremely steep strain gradients across the thickness of the NSBs. According to strain localization theories (43, 44), under these conditions, the thickness of shear bands approaches a constant value, which is approximately 15 times the grain diameter in the shear band. For Mg2GeO4, the spinel “fault gouge” has average grain size on the order of 5 to 10 nm (28, 41). Thus, shear band theories give thickness of NSBs, consistent with our observations (~100 nm).

Even without invoking any phase transformation, shear localization theory predicts that, within the regime of shear band formation [Embedded Image], rocks weaken with increasing deformation, leading to instability. In the case of exothermic olivine-spinel transformation, once an NSB is formed, the latent heat released promotes further transformation. This is a positive feedback process similar to adiabatic shear runaway process, driving the propagation of NSBs. Thus, NSBs are likely the growth products of slip bands, as observed in samples preceding macroscopic faulting (40). NSBs may also develop from the tips of LVRs (Fig. 8A) and from the interfaces between olivine and pyroxene grains (Fig. 8B) due to stress concentration. Recorded AE activities (Fig. 4) suggest that the “avalanche” of NSB propagation most likely took place around the critical point (1700 s). Beginning around 1900 s, increasingly larger events took place preferentially in conjugated planes, presumably formed by dense arrays of NSBs. These NSBs are superplastic weak zones that behave like highly lubricated surfaces, leading to macroscopic shear offsets shortly after 2300 s, as seen in the micrographs (Fig. 7).

Microstructural interpretation of the nanoseismic behavior

Analyses based on P-wave amplitude provide important constraints on source parameters. It is now feasible to obtain true moment magnitudes of AE events. With six transducers located in three orthogonal directions, the diagonal elements of moment tensor of any event are measured directly. Therefore, the percentage of isotropic component of the event is well resolved. Most events are found to have a dilatational isotropic component of about 1% to 4% (=ζ2) (fig. S5) (45). There are two possible sources of error, which introduce systematic bias in measuring ζ. The first is of instrument origin. Unfortunately, before the experiments, the transducers were not quantitatively calibrated; hence, their individual responses might differ slightly due to imperfect coupling to the anvils. The second is the uniform whole-space assumption used in moment tensor inversion. Any departure from this simplified model will surely affect the moment tensor results. Both may contribute to the misfit of synthetic waveforms to the observations (Fig. 2). The resultant source parameters shown in fig. S5 may also contain systematic errors and should not be taken as quantitative at this point. In future AE experiments, careful calibration of transducers will be performed to obtain more reliable information on source parameters.

Although most of the ζ values fall within ~0 and 0.2, smaller events (MW < −8.5) exhibit large variations, from ~−0.3 to +0.4. This range is significantly above the estimated system bias. A tightening trend is observed for events with MW > −8.5 (Fig. 5A). This is further confirmed by examining the average event rate (in events per minute) with time (Fig. 5B). This jump occurs at ~200 s before the criticality point, after which event rate drops sharply to 4.0 between 2250 and 2550 s, before increasing again to 6.4 prior to the ultimate failure.

On the basis of the observed microstructure (Figs. 6 and 7 and fig. S8), we postulate that MW < −8.5 events with large |ζ| are the result of precursory slip bands or smaller NSBs growing and cutting through individual olivine grains with associated grain boundary sliding, causing instantaneous creation and annihilation of small voids due to geometric mismatch. At the early stage (below 1500 s), these events may not be self-similar, thus the large apparent b values before the critical point (Fig. 4E). Above 1500 s, NSBs become increasingly more active and dominant, hence the increased event rate. This may be an important step in “preparing” the propagation NSBs and faultlets. Shortly after the criticality point, ζ values tighten rapidly before the ultimate failure. Accordingly, b values decrease rapidly toward 1.

The CLVD components are reasonably resolved and, overall, significantly large. The complex arrays of subparallel and conjugated faults, as shown in fig. S9, are expected to propagate contemporaneously. For P- and S-wave velocities of 7.4 and 4.2 km/s (46), respectively, dynamic triggering of NSBs by the adjacent microfault (ca. 10 μm away) requires a time lag on the order of 1 and 2 ns, respectively, too short to be detected with our current detectability. The shortest source duration of 0.25 μs that can be recorded by our AE system may easily contain multiple dynamically triggered sub-events. Therefore, the observed fault zones are not singular fault planes but several segments propagating together, or one triggering another, resulting in the large CLVD components observed.

It has been shown that, on average, deep earthquakes tend to have larger CLVD components than shallow earthquakes (47). This has been interpreted as due to relative motion of complex shear plane networks (48). For DFEQs, the large CLVD components are attributed to the dynamic triggering by the passage of the seismic waves from the initial rupture (49). Recently, it has been shown that many DFEQs contain multiple sub-events and display great variability of focal mechanisms, resulting in significant CLVD components (50). The closeness in time of the ruptures on separate planes and the distance between the planes suggest dynamic triggering, where seismic waves from the first rupture initiate rupture on the second plane, resulting in significant non–double-couple components in moment tensor solutions (51). These interpretations are consistent with the microstructures observed in our laboratory specimens. The subparallel faults in Fig. 7C, for example, exhibit various offsets. Superposition of multiple events with different rupture displacements may be responsible for the CLVD components. This may also explain observations of near- or supersonic source propagation speeds in some DFEQs (52). Finally, the AE history we observe (Fig. 5) suggests that DFEQs may have microseismic precursors, just like their shallow interplate counterparts (53).

Can the laboratory rupture behavior be linked to natural earthquakes?

The much improved detection sensitivity allows events with much lower magnitudes to be identified. This, in turn, better constrains b values in the Gutenberg-Richter statistics. The well-defined b value suggests that rupture events in our experiments are scale-invariant up to a few millimeters. Our observed b-value evolution before macroscopic failure also generally agrees with AE detections on large testing pieces of rocks and concrete components up to several meters. In particular, it has been noted in testing large rock and concrete blocks that b value decreases from ~1.5 to 1.0 during deformation, from the critical point (stress peak) to the “final collapse” (54). This is consistent with our observations (Fig. 4E). Assuming the failure process to be self-affine, the b values have been argued to relate to fractal dimensions, Df, of the fracture network: Df = C × b, where C is a constant typically ranging from 4/3 to 2 (55, 56). Thus, b = 1.5 implies a 3D fracture network, and b = 1.0 corresponds to a 2D fault system. On the other hand, laboratory samples are size-limited, making it difficult to extrapolate results to larger scales (57).

A recent analysis (58) based on more than 40 years of observations reveals a bimodal pattern for DFEQs at a depth of 500 to 700 km: In the cold Tonga slab, b is close to, or slightly greater than, 1, whereas in warmer slabs (for example, South America, Japan-Kuril, and Izu-Bonin-Mariana), b is close to 0.5 for intermediate magnitudes (MW = 5.3 to 6.5) and increases to ~1 for large magnitudes (MW > 6.5).

Linking b values to different faulting mechanisms is not straightforward (57, 59). The simplest explanation for the above observations appears to be the availability of nucleation sites for rupture to propagate. In western Pacific and South America, seismogenic layers in the MTZ are only a few tens of kilometers thick (60). These layers, presumably the slab “cores” enriched in metastable olivine, may be approximated by 2D sheets (61, 62). On the other hand, the largest DFEQ swarm of the world in Tonga occupies a substantial 3D volume in the MTZ, with lateral extensions of >200 km (63). This swarm of DFEQs is interpreted as failure of stagnant cold slab materials accumulated in the MTZ, due to the unusually fast slab descending rate in this region (63, 64). Thus, within the framework of transformation-induced instability, the different b values for DFEQs may be simply related to the extent of available metastable material: In western Pacific and South America, failure occurs within thin and planar metastable slab cores, which are essentially 1D for earthquake ruptures (hence, b ~ 0.5), whereas in Tonga, more potential nucleation sites are available due to the large volume of metastable slab materials (therefore, b ~ 1 to 1.1) (58), similar to what is observed in laboratory samples after reaching the criticality point. On the other hand, the greater b values of ~1 for the large DFEQs in warm slabs may be due to propagation of ruptures outside the metastable olivine wedge, thereby suggesting a different mechanism (58).

Finally, the following empirical relation between rupture area (A) and MW is shown to hold for earthquakes regardless of faulting styles (30, 35)Embedded Image(2)where c1 and c2 are constants. Two statistical analyses give identical c1 = 1.02, with c2 ranging from −4.15 (35) to −4.01 (30). Recent analyses (65, 66) show that these relations are valid for A up to ~500 km2. Theoretically, if rupture processes are scale-invariant, c1 = 1 (67). Equation 2, determined from earthquakes MW ≥ 5, predicts that the rupture length (≈Embedded Image) corresponding to MW = −7 is 2.2 to 2.7 mm. The maximum rupture length corresponding to the largest AE events (MW = −7) is approximately the sample diameter (~2.5 mm). Considering the more than 12 orders of magnitude difference in energy release from MW = −7 to MW > +5, the agreement is remarkable, suggesting that rupture is generally a scale-invariant process from millimeters to kilometers.


High-pressure, high-temperature deformation with AE monitoring

High-pressure, high-temperature deformation experiments were conducted in the deformation DIA apparatus (D-DIA) at 13-BM-D Beamline of GESCARS Sector, Advanced Photon Source, Argonne National Laboratory. Details of the apparatus and operation were described elsewhere (68). The D-DIA consists of six anvils, each with a square truncated tip of 6 mm in linear dimension, compressing a cubic pressure medium, in which a sample with initial dimensions of 2 mm in diameter and 3 mm in length was loaded. Figure S1 shows the conceptual setup with the four horizontal anvils. Top and bottom anvils are not shown. Two anvils were made of sintered diamond (SD; Versimax from Sandvik Hyperion). The x-ray transparent SD anvils allowed complete diffraction Debye rings to be recorded with a Rayonix area detector (MAR-165) for the determination of differential stress (69). The other two anvils were made of tungsten carbide (WC). Top and bottom anvils were also made of WC. A graphite tube heater was placed around the sample, separated by a soft hexagonal boron nitride liner. Fully densified alumina pistons were placed at the upper and lower ends of the sample, with 5-μm gold foils marking the interface. At high pressure and temperature, the differential rams built in the D-DIA module were advanced to deform the sample. Sample axial strain was monitored and measured by the shadows of the gold foils with x-ray imaging.

A PZT (lead-zirconate-titanate) piezometric acoustic transducer was mounted on the back side of each anvil. Acoustic signals from the transducers were recorded through six channels. Channels 1 and 2 corresponded to the bottom and top anvils, respectively. Channels corresponding to the horizontal anvils are given in fig. S1. All the transducers were from Boston Piezo-Optics Inc., with a center frequency of 3.6 MHz, except for channel 5, which was from OMEGA Engineering Inc., with a center frequency of 2.0 MHz. Both transducer types have narrow bandwidths of <1 MHz.

The acoustic signal that reached the transducers was converted into an electric signal, which was then amplified and recorded with an 8-channel AE and 18-bit ultrasonic acquisition system manufactured by Itasca Microseismic and Geomechanical Evaluation (IMaGE) ( with a maximum sampling frequency of 50 MHz. In this system, amplifier circuitry has selectable gains between 30 and 70 dB and a broad frequency response with a plug-in filter circuit that analytically tailors the unit to the desired frequency range. Signals were registered by an OMNIBUS trigger and hit count (THC) unit. THC was equipped with transient recorders that continuously recorded information in temporary memory, which was discarded unless signals were energetic enough to exceed a user-defined threshold. Such a unit is crucial because it markedly reduces the amount of memory required for an experiment, as only the most important information is recorded, while ensuring the capture of entire waveforms whose onset takes place before threshold crossing. Therefore, invaluable information associated with the precise onset of the signal can be used for location analysis and event characterization. The system was also equipped with a pulser-amplifier that allowed switching between pulsing and receiving modes for all the transducers, allowing each to act as both a transmitter and a receiver. This was used in our experiment to establish a velocity model for the cell assembly that accounts for faster P-wave velocities in the SD anvils (relative to WC anvils) as well as across the vertical alumina deformation pistons in the cell that act as waveguides. In the experiment, triggered AE signals were captured on six oscilloscopes and six Cecchi data acquisition units. One channel (channel 2) was split into two so that continuous data acquisition was recorded with a Richter system, with a sampling rate of 10 MHz. This allowed for the detection of nonimpulsive and/or long-period signals otherwise missed.

We used the software package InSite-Lab (also by IMaGE) for data collection and processing. With this software, signals can be treated using standard high- and/or low-pass filters. The triggered waveforms can be automatically processed to pick acoustic wave arrival times, maximum amplitude, and first amplitude. P-wave arrival time differences can then be interpreted in terms of travel time differences within the cell. Location of each event can be obtained by inversion calculations using various algorithms based on the velocity model. Focal mechanisms can also be retrieved from recorded AE data.

The AE data analyzed in this study were first reported by Schubnel et al. (28). The sample of run D1247 was chosen in the present study. In total, 593 AE events were recorded in the triggering mode in this experiment; all occurred during the olivine-spinel transformation at 4 GPa and 1173 K during deformation, when sample axial strain exceeded 8%. All of the events have identical record length of 163.8 μs. Sample D1253, which was deformed and transformed at 5 GPa and 1173 K with 73 AE events recorded, was also reanalyzed for microstructure and compared with D1247. Both samples contained ~15 volume % pyroxene (MgGeO3). This ratio of olivine to pyroxene is similar to harzburgite, which is considered the major component of subducting slabs. Pyroxene grains also served as convenient strain markers.

Microtomography on recovered samples

The recovered samples with their surrounding pressure medium were cast in epoxy resin. The entire assemblies were then examined by XMT. XMT data were collected in the high-pressure XMT apparatus (70) but under ambient condition. A monochromatic beam of 31 keV was used, with linear pixel dimensions of 0.0029 mm. Estimated spatial resolution was about 5 to 6 μm. 3D XMT images were processed using the GSECARS software package tomo_display ( Thresholding, smoothing, and display of the image were conducted using the open-source image processing program ImageJ ( Special care was taken during the data collection so that the coordination system of the XMT image and that of the AE events were aligned to with approximately 5° rotation error about the vertical axis. This made it possible to compare locations of the faults with AE event hypocenters.

AE event relocation using hypoDD

We used the hypoDD algorithm (29) to relocate the AE events. We first cross-correlated waveforms between events for each channel, using a 6-μs-wide time window, starting at 1 μs before the first P-wave arrival. The CC results were visually inspected. For channels 1 to 4, CC coefficients greater than 0.8 could identify similar waveforms and align them correctly. For channels 5 and 6, a higher CC coefficient of 0.9 was needed to avoid misalignment of waveforms due to cycle skipping and noise due to higher noise levels. The CC algorithm determines relative location and origin time between any two events. 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 the trade-off 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. Differential arrival times of event pairs measured by the waveform CC were then used in the hypoDD algorithm to relocate the AE events. In the analysis, P-wave velocity of Mg2GeO4 olivine [7.4 km/s (46)] was used because the amount of spinel was less than 5 volume % according to XRD. The 593 events could not be located in one group to satisfy the above requirements. The algorithm divided them in groups of different sizes and determined the relative locations of the events within each group.

Source parameter determination

We assume that each AE event can be treated as a point source. This assumption was justified for most of the detected AE events, whose rupture lengths were estimated to be less than tens of micrometers at most, much smaller than the wavelength of signals (~2 to 3 mm) and the distances between the events and the sensor (~25 mm away from the sample). Large events with saturated amplitude were not used in the moment tensor analysis.

Using far-field displacement solution of a point source in a homogeneous whole space (33)Embedded Image(3)where M is the source moment tensor, and α and ρ are whole-space P-wave velocity and density, respectively, the area under the first P-wave pulse at each of the six transducers can be expressed asEmbedded Image(4)where i = x, y, and z, corresponding to transducer pairs 3-5, 4-6, and 1-2, respectively. A factor of 2 was introduced to take into account the free surface reflection effect at the transducer locations. It is clear that with the six transducers in the three orthogonal directions, only the three diagonal elements of M can be determined. Because we were more interested in the partition of isotropic and deviatoric components of the source and fault plane orientation, we rewrote the moment tensor M in terms of the scalar moment M0 and a “normalized” tensorEmbedded Image(5)and decompose D into the ISO, double-couple (DC), and CLVD componentsEmbedded Image(6)whereEmbedded Imageζ is the nondimensional ISO contribution between −1 (implosion) and 1 (explosion), χ is a nondimensional CLVD contribution between −0.5 and 0.5, and Embedded Image and Embedded Image are fault plane’s normal and shear slip directions, respectively, expressed in terms of strike φ and dip δ of the plane and rake λ of the slip (33), Embedded Image.

For every set of possible values (ζ, χ, and φ), we first calculated the three diagonal elements Dii using Eq. 6 and then calculated M0 by using the L2-norm ratioEmbedded ImageThe best solutions of ζ, χ, and φ were found that minimize the L2-norm misfitEmbedded Image

For “absolute” displacement amplitude estimation, we converted the voltage amplitude of the AE events to displacement using the manufacturer’s piezoelectric constant (586 × 10−12 m/V), after taking into account the 60-dB signal amplification in the recording system.

It is possible that each recorded AE event may consist of several smaller sub-events; therefore, the estimated moment tensors represent some kind of average or, more precisely, a sum of source process in the event record. However, most events’ source durations were less than 0.3 μs. In such a case, the sub-events within an event must be highly correlated both temporally and spatially. They most likely occurred on planes of similar orientations with similar slip directions. Therefore, the double-couple component of the obtained moment tensor will represent the overall orientation of the fault plane and slip direction.

In the reported experiments, we only measured amplitudes of direct P waves in three orthogonal directions. Because there were six independent parameters for a general moment tensor, our source inversion was an underdetermined problem. Given the axial-symmetric compression setup, it was reasonable to assume that the rake angle was −90° (that is, normal faulting). With this assumption, the amplitude data had no sensitivity to the dip angle of the fault planes, so we fixed it at 45°. These two assumptions will certainly affect the moment tensor solutions, with the most affected parameter being the strike of the fault plane. The moment magnitude and isotropic component estimates were relatively robust.


Supplementary material for this article is available at

fig. S1. Top view of the experiment setup.

fig. S2. Low-magnification SEM images showing cross sections of the entire samples of the experiments.

fig. S3. Example waveforms of AE events in groups 1 and 2.

fig. S4. An example of grid search result.

fig. S5. Histograms of ISO parameter ζ, CLVD parameter χ, and strike φ of 593 moment tensor solutions (light gray).

fig. S6. An example of cross-correlation template searching results for smaller AE events.

fig. S7. Waveforms of events with different moment magnitudes.

fig. S8. Orientation relations between olivine grains and NSBs.

fig. S9. BSE micrographs of recovered samples.

movie S1. XMT image of the recovered sample D1247.

movie S2. XMT image of the recovered sample D1247.

movie S3. Relocated AE events in group 1 superimposed with the digitized faults, as shown in movie S2.

movie S4. Relocated AE events in group 2 superimposed with the digitized faults, as shown in movie S2.

movie S5. Group 1 events viewed along the strike of fault plane 1, throughout the deformation/transformation history.

movie S6. Group 2 events viewed along the strike of fault plane 2, throughout the deformation/transformation history.

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


Acknowledgments: We thank H. Green for providing sintered Mg2GeO4 polycrystals for this study. The EBSD analysis was conducted at the Centre Commun de Microscopie at Université de Lille. We are grateful to Z. Zhan and an anonymous reviewer for the constructive comments that significantly improved this manuscript. Funding: Portions of this work were performed at GeoSoilEnviroCARS (University of Chicago, Sector 13), Advanced Photon Source, Argonne National Laboratory. GeoSoilEnviroCARS is supported by the NSF Division of Earth Sciences (EAR-1128799) and U.S. Department of Energy (DOE), GeoSciences Division (DE-FG02-94ER14466). This research used resources of the Advanced Photon Source, a DOE Office of Science User Facility operated for the DOE Office of Science by Argonne National Laboratory under contract no. DE-AC02-06CH11357. Y.W. acknowledges NSF support EAR-1361276 and 1661489. L.Z. and Z.L. were supported by NSF grants EAR-1249701 and 1661519. Author contributions: Y.W. and L.Z. conceived the idea of this project. A.S., N.H., J.G., F.B., and Y.W. conducted the AE experiments. T.Y., M.R., and Y.W. collected the XMT data. A.S., N.H., F.B., and Y.W. conducted the microstructural analyses. L.Z., F.S., and Z.L. conducted the seismological analyses. A.A. performed the TEM EBSD analyses. D.D. conducted the SEM analyses. Y.W., L.Z., and F.S. wrote the manuscript, with all authors providing inputs. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article