Earthquake rupture below the brittle-ductile transition in continental lithospheric mantle.

Earthquakes deep in the continental lithosphere are rare and hard to interpret in our current understanding of temperature control on brittle failure. The recent lithospheric mantle earthquake with a moment magnitude of 4.8 at a depth of ~75 km in the Wyoming Craton was exceptionally well recorded and thus enabled us to probe the cause of these unusual earthquakes. On the basis of complete earthquake energy balance estimates using broadband waveforms and temperature estimates using surface heat flow and shear wave velocities, we argue that this earthquake occurred in response to ductile deformation at temperatures above 750°C. The high stress drop, low rupture velocity, and low radiation efficiency are all consistent with a dissipative mechanism. Our results imply that earthquake nucleation in the lithospheric mantle is not exclusively limited to the brittle regime; weakening mechanisms in the ductile regime can allow earthquakes to initiate and propagate. This finding has significant implications for understanding deep earthquake rupture mechanics and rheology of the continental lithosphere.

INTRODUCTION

The role of temperature in the rheology of the lower crust and lithospheric mantle is now well understood, especially in oceanic lithosphere. Earthquakes occur throughout the oceanic crust and upper mantle, the latter with olivine-dominated lithology; the brittle-ductile transition is associated with a threshold temperature of about 600°C (–). Upper-mantle seismicity is less common in the continental lithosphere. Earthquakes tend to occur mostly in the upper crust and, in some cases, in the uppermost mantle (). For the latter case, olivine-dominated lithology predicts brittle failure at temperatures below ~600°C (). It has been proposed (, ) that the rare continental lithospheric earthquakes occur in cold and anhydrous mantle rocks (T < ~700° ± 100°C). Earthquakes located below the Moho may thus indicate a strong, seismogenic upper mantle that can sustain large stresses, which are later released during brittle rupture (–).

An isolated earthquake with a moment magnitude (MW) of 4.8 occurred deep in the continental lithosphere on 21 September 2013 near the Wind River Range in central Wyoming at a depth of 75 ± 8 km (). The focal mechanism indicates a dominant strike-slip faulting mechanism with a small reverse-slip component (). It was followed within 2 hours by a single M 3.0 aftershock at a similar hypocentral depth. The Wyoming earthquake and its aftershock are unusual because they occur deep in the lithospheric mantle, far from any convergent plate boundary (Fig. 1). A few previous earthquakes have occurred in the region at a similar depth (), including a M 3.8 earthquake in 1979 beneath Randolph, UT (), but the 2013 Wyoming earthquake is by far the best recorded. The complexity of the velocity structure in this transitional region and the large uncertainties in the temperature in the hypocentral regions meant that previous studies (, , ) could exclude the possibility that these earthquakes represent brittle rupture in cold mantle.

Fig. 1

Earthquake location and P- and S-source time functions.

(A) Map of the seismic stations (triangles) used here and location of the 21 September 2013 Wyoming earthquake (star) and its focal mechanism (beach ball). Colored triangles stand for seismic stations used in the EGF procedure. Black triangles show the extended set of stations used in the radiated energy estimate. Inset shows the study area (red box) as well as the epicenter of the earthquake (red star). Single-station STFs obtained from the EGF procedure on P wave (B) and S wave (C). STFs are sorted as a function of station azimuth from the fault strike. Waveform colors correspond to those of the stations in (A). The gray dashed lines are the predicted widths of the P- and S-STFs (see main text).

Here, we analyze broadband seismograms at regional distances to determine the static and dynamic source parameters of the 2013 Wyoming earthquake and for clues concerning the rupture conditions. Estimates of the seismic radiation, the source dimension, and the amount of slip of an earthquake can be used to investigate the stress drop and fracture energy (). We also perform detailed thermal modeling to ascertain whether the Wyoming earthquake occurred below or above the brittle-ductile transition temperature.

RESULTS

Earthquake source analysis

We use the M 3.0 aftershock as our empirical Green’s function (EGF) to remove propagation effects from the recorded seismograms and to obtain the relative source time functions (STFs) and spectral ratios (see Materials and Methods). The STFs at most of the stations show a clear source pulse, with consistent shapes between nearby stations (Fig. 1). Figure 1 (B and C) shows a clear azimuthal dependence of the width of single-station P- and S-STFs, which reflects predominant rupture propagation along the subvertical northwest-oriented plane; stations aligned along the rupture direction have narrower pulses, whereas stations away from the rupture direction have broader pulses. This angular dependence of apparent duration can be quantified using a directivity function D, which depends on the rupture velocity VR, the wave velocity around the source c, and the angle between the ray takeoff angle and the rupture direction φ at station iwhere Tapp and TR are the apparent and the actual rupture durations, respectively. Because we have both P and S measurements, we can combine them to obtain the length of the rupture without assuming the rupture duration. We stretch the STFs in the time domain according to Eq. 1, varying the rupture velocity, fault length, and rupture direction, and compute the cross-correlation between all STFs to find the model that best explains the azimuthal variation and has the highest cross-correlation (see Materials and Methods). This approach eliminates uncertainties inherent in picking start and end times of individual STFs or in modeling spectral parameters. We include the possibility of bilateral rupture () but find that the data are best fit by a slow unilateral and subhorizontal rupture (Fig. 2A), with VR = 1.3 ± 0.3 km/s (1-σ uncertainty) (~0.28VS) (Fig. 2, B and C). Even the maximally allowed value of VR at the 95% confidence level (~0.41VS) is significantly lower than those commonly assumed (~0.8VS to 0.9VS) in source studies (Fig. 2D). We sum the single-station STFs corrected for directivity effects using our estimated VR (Fig. 2C) to improve the signal-to-noise ratio (SNR) and to obtain clear and unbiased, average P- and S-stacked STFs (Fig. 2B). From these STFs, we calculate a model-independent rupture duration TR = 0.47 ± 0.04 s and hence the length of the rupture L = VR × TR ≈ 610 m.

Fig. 2

Rupture directivity and velocity estimate.

(A) Contours of normalized variance reduction as a function of rupture velocity VR and percent unilateral rupture values e. The variance of the best-fit model in this plot is set to 1. Dark blue contours indicate variability within 10%, and the dashed line shows the best-fit rupture velocity. (B) Stacked P- and S-STFs resulting from the correction using our best-fit rupture velocity VR = 1.3 km/s. (C) Single-station S-STFs corrected for directivity effects using the stretching method with VR = 1.3 km/s. (D) Same as (C), with VR = 3.8 km/s. Red waveform at the bottom of each figure is the stack of all the corrected STFs. This correction involves stretching the STFs at station i with a factor of 1/D.

Calculation of the earthquake stress drop (proportional to slip/length or seismic moment/area2/3) is source model–dependent. We use the global centroid moment tensor seismic moment M0 = 2.17 × 1016 N·m (www.globalcmt.org). If we assume a square fault with a length of 610 m, we obtain Δσ ≈ 60 MPa, while a circular fault with a radius of 305 m yields Δσ ≈ 330 MPa (see Materials and Methods). The average corner frequency (~2 Hz) would give a stress drop of ~50 to 150 MPa, depending on source model assumptions (see fig. S2). All these values are relatively high compared to those of shallow earthquakes () but similar to those of other intermediate-depth earthquakes (, ) of similar size.

We calculate the earthquake-radiated energy ES by averaging frequency-domain estimates from P and S waves at all available stations (see Materials and Methods; Fig. 1A). We obtain a total radiated seismic energy ES = 3.9 × 1012 J. The scaled energy (ES/M0) is relatively high (1.7 × 10−4), and the apparent stress, defined as σa = μES/M0, is equal to ~12 MPa, where μ is 7 × 1010 N m−2. The stress drop obtained from the directivity modeling (60 to 330 MPa) is at least five times larger than the apparent stress and the radiation efficiency ηR ≤ 0.4. We compare this value to the radiation efficiency expected from our estimated VR using crack theory, resulting in our preferred radiation efficiency ηR ≈ 0.2. This value is smaller than those of most shallow earthquakes, similar to those of intermediate-depth earthquakes in Japan and Mexico (, ), and larger than those of other deep earthquakes (, ).

In summary, the deep Wyoming earthquake has an anomalously low rupture velocity, a high stress drop and scaled radiated energy, and a low to average radiation efficiency. We interpret this event to represent a dissipative source mechanism that was slow to propagate with a significant portion of the available potential energy used in or around the source region.

Temperature modeling

Knowledge of the temperature at the hypocentral depth of the Wyoming earthquake is necessary to determine whether the earthquake nucleated above or below the expected brittle-ductile transition. Surface heat flow measurements in the Wyoming Craton range from 40 to 60 mW m−2. Temperature estimates at the hypocentral depth of the earthquake vary from ~600° () to more than 800°C (), but large uncertainties remain. The earthquake is located along a strong lateral gradient of shear wave velocities, similar to other lithospheric mantle earthquakes (, ) as well as in a transition of lithospheric thickness (, ). In addition, estimated temperatures at depth based on xenoliths in the Wyoming area suggest a fairly warm lower crust (, ), and recent volcanism in Leucite Hills (about 100 km to the south and <3 million years ago) may require high temperatures in the upper mantle ().

We estimate the depth-dependent temperature profile using the grain-size viscoelastic relaxation model of olivine in the mantle () to predict the shear wave velocities at depth and compare to tomographic models [(, ); see Materials and Methods]. The best-fitting geotherms (Fig. 3) predict temperatures at the hypocentral depth of the Wyoming earthquake that are substantially higher than the 600°C transition from velocity weakening to velocity strengthening of olivine extrapolated from laboratory experiments (). In common with some other studies (, , ), we favor the warmer geotherms in Fig. 3 because those that predict T < 750°C at the earthquake hypocenter also predict low temperatures down to 400 km, in disagreement with lithosphere-asthenosphere boundary depth estimates (, ). A similar analysis of the hypocentral region of the 1979 Randolph, UT earthquake predicts even higher temperature ranges (see fig. S1).

Fig. 3

Temperature modeling of the Wyoming earthquake.

(A) The VS model () in the western United States at a depth of 75 km and east-west cross sections showing the location of the Wyoming earthquake (star) in the strong velocity gradient. (B). Geothermal models (left; inset shows a zoom around the hypocenter) and predicted shear wave velocities using the approach of Faul and Jackson () that best fit the VS tomography models (right) of Schaeffer and Lebedev () (blue area) and Shen et al. () (red area). The geothermal models assume a crustal thickness of 50 km and produce surface heat flows between 40 and 60 mW m−2. The range of temperatures at the hypocentral depth predicted is from 750° to 850°C in the upper range of the brittle-ductile transition.

The temperature modeling implies that the Wyoming earthquake occurred in the upper range of the brittle regime or in the ductile regime. The source modeling found a slow rupture with a dissipative rupture process, which is more consistent with a ductile regime. Our estimates of rupture velocity are extremely low, similar to that of the Bolivia earthquake () or the initial rupture of intermediate-depth earthquakes (, ), and are lower than oceanic lithospheric mantle earthquakes (–).

DISCUSSION

Previous studies of the 2013 Wyoming earthquake focused on its location and orientation (, ). The focal mechanism shows similar stress orientation to that of the shallow seismicity, suggesting that they are responding to the regional stress regimes (, ). The location of the earthquake within a high-velocity region (), combined with the large uncertainties in the temperature profile, meant that previous studies were unable to exclude the possibility that it was caused by brittle failure in cold, anhydrous mantle. More detailed modeling of the thermal structure of the lithosphere suggests that temperatures in the hypocentral region are quite high (), regardless of the presence of the high-velocity region (), and our analysis suggests that T > 750°C (Fig. 3). Our source analysis found that the earthquake had low radiation efficiency, implying that a large amount of energy was dissipated during rupture, perhaps as frictional heating (, ), in contrast to brittle-like, fast rupture of other deep earthquakes (). The additional lines of evidence presented here raise questions about the brittle nature of the 2013 Wyoming earthquake and potentially of other continental lithospheric mantle earthquakes elsewhere.

We hypothesize that an ongoing ductile deformation in the mantle, for example, along a midlithospheric discontinuity (), may lead to shear band formation and slip concentration along subvertical zones. The earthquake may have nucleated along one of these shear zones through a thermal shear runaway mechanism or aided by dehydration reactions, triggering the MW 4.8 earthquake and its only aftershock (, , ). The 1979 deep Randolph earthquake probably nucleated in a similar manner (see fig. S1 for temperature modeling).

Another possible explanation is that compositional or mineralogical heterogeneity in the mantle [for example, grain-size heterogeneity (–)] can promote brittle behavior or instability in a small volume, wherein the earthquake nucleates and then propagates into the mostly ductile mantle rocks. This mechanism has been proposed for the Newport-Inglewood Fault (), where deep earthquakes are attributed to a system of seismic asperities in a ductile fault zone. Nucleation may additionally be encouraged by the presence of fluids or fluid migration in the mantle, either by reducing the effective normal stress or promoting strain localization along the shear zones (, ).

Source analysis of oceanic lithospheric earthquakes off the shore of Sumatra show deep slip in the hotter, velocity-strengthening region (, ), but this deep slip is only observed in the greatest earthquakes (). Aseismic creep is also observed at shallow depths, and the small repeating earthquakes at Parkfield, CA are thought to involve slip in both velocity-weakening and velocity-strengthening regions (). In contrast to both of these cases, the nature of rupture in the Wyoming earthquake implies that the rupture was dominated by slip in the more ductile, aseismic region.

The presence and mechanisms of earthquakes in the continental lithospheric mantle have been a highly controversial topic in geophysics (, –, , ). Our results on the 2013 Wyoming earthquake shed new light on the rupture mechanism of these rare earthquakes (, , ) and provide important constraints on rheological properties of the lithospheric mantle (, , , ).

MATERIALS AND METHODS

EGF method

We studied the source parameters using the EGF method (). This approach used a smaller event that is collocated with the earthquake of interest to remove the path, site, and instrumental effects from the earthquake source signal in the seismic recording (). Fortuitously, the single aftershock was an ideal EGF candidate. Although its location is less well constrained than that of the main shock, their waveform similarity suggests a similar location, orientation, and depth (), making these events ideal candidates for a good EGF earthquake pair ().

We used raw seismograms recorded in a 10° × 10° area around the epicenter zone of the two earthquakes, shown as colored triangles in Fig. 1A. We performed the analysis on vertical (P waves) and radial and transverse (S waves) seismograms. After the EGF procedure (), we computed the earthquake STF using the multitaper deconvolution algorithm () at stations for which the cross-correlation coefficient between the two event waveforms is ≥0.6 for P waves and ≥0.5 for S waves.

Directivity from STFs: Cross-correlating stretched waveforms

Because some earthquakes exhibit bilateral rupture, we used the directivity function ()where 0 ≤ e ≤ 1 is the “percent unilateral rupture” that reflects the degree of unilateral rupture, and the notation ζ = VR/c cos(φ) is introduced to simplify the expression.

Static source parameters are usually estimated by modeling the source signal with simple source models (). However, earthquake ruptures may be much more complex than models suggest, and significant azimuthal variation of the spectral characteristics of the recorded signals is expected (, ). We proposed here to use the directivity observations to estimate the rupture velocity and duration without assuming any particular source model. Our technique was based on the idea that directivity effects appear as stretching or compression of the waveforms. In view of the different velocities of P and S waves, directivity effects were expected to be stronger on S wave recordings. As a result, the S-STF at a given station i turns out to be a stretched (or compressed) version of the P-STF. Following Eq. 1, we could write the ratio between single-station S- and P-STFs apparent durations T as

This ratio removed any dependence on TR and could be used to estimate the rupture velocity VR.

In practice, we measured the ratio between S- and P-STFs durations T at each station as the factor ε, where the time axis of S-STF has to be stretched or compressed to obtain the best correlation with P-STF. We then performed a grid search on VR to look for the rupture velocity that best explains our measured ε over the set of stations, with the model in Eqs. 1 and 2. Our results showed that a slow unilateral rupture is required to explain the data (Fig. 2A). The best fit between the predicted and the measured T was obtained for e ≈ 1 and VR = 1.3 km/s (~0.29VS).

This measurement of VR can then be used to correct the directivity effect on the individual STFs, as expressed by the actual source pulse observed at each station. Following our directivity model in Eq. 1, this correction came down to stretching the STFs at station i with a factor of 1/D using the values that best fit the data (that is, e = 1 and VR =1.3 km/s). Figure 2C shows the resulting single-station S-STFs corrected for directivity effects. In comparison, Fig. 2D shows the corrected single-station S-STFs, assuming the widely used value of VR = 0.9VS (, , , ). In the latter case, we clearly saw that STFs at stations aligned in the direction of the rupture are overcorrected, implying that the actual rupture velocity was much lower than 0.9VS. This test definitively confirmed that a slow rupture velocity is required to explain our data.

We determined that rupture duration is 0.44 and 0.51 s on the stacked P- and S-STFs, respectively. Estimation of the width of the resulting stacked STFs was straightforward and did not rely on a particular source model, as is usually the case when interpreting source spectra. The rupture duration was 0.44 and 0.51 s on the stacked P- and S-STFs, respectively, which led to an estimate of TR = 0.47 ± 0.04 s. Note that we also investigated the time resolution limit by computing the delta function (that is, the deconvolution of the earthquake by itself) at each station. All the resulting delta functions were narrower than 0.1 s. The pulse duration of single-station STFs are always >0.2 s, confirming that our STFs are all well above the time resolution limit.

Radiated seismic energy

To improve the stability of our energy estimate, we extended the number of stations used in the estimate, including colored and black triangles in Fig. 1. We corrected the instrumental response at each station and the propagation effect using a frequency-independent quality factor (Q = 1000). We estimated the quality factor Q that best explains the EGF spectral ratio in the previous step, with values ranging from 500 to 1500. Finally, we performed the radiated energy analysis on almost 70 stations compared to 30 single-station EGF estimates (see Fig. 1A).

We calculated the radiated seismic energy at each station in the frequency domain by integrating the velocity spectrum from both P and S wave ground motions for the three components, which are summed to obtain the total seismic energy (). We took radiation pattern corrections into account but only used stations with radiation pattern coefficients of 0.2 or larger in our energy estimates () to limit the influence of stations near the nodal plane.

Estimates of radiated energy required the integration over a wide frequency band, but the SNR limited our observable band to about 0.5 to 20 Hz; thus, we extrapolated the velocity spectrum, assuming a Brune spectral model (ω−2 falloff) at both low and high frequencies. When not all of the components were available, we set them equal to the average energy value of the available ones. Not all stations had P and S energy estimates because of the SNR threshold limitation (we chose stations for which SNR > 2 over at least 50 samples, that is, a bandwidth of ≥7 Hz), resulting in 44 P and 17 S wave station estimates. The total energy radiated by the earthquake was calculated by summing the contributions of the P and S wave energies E and E after being logarithmically averaged over all the stations. We assumed VS = 4500 m/s and obtained a total radiated seismic energy E = 3.9 × 1012 J and an S-to-P energy ratio of 22.

Calculation of stress drop

We used multiple approaches to estimate the rupture area, and hence the stress drop, because they were model-dependent. If we assumed a unilateral strike-slip rupture with a length of 610 m (calculated from the directivity analysis) and that fault width (W) was equal to length (L), we would obtain (, )

This should be considered a minimum estimate because W could be much smaller than L. If we assumed instead a circular fault with radius R = 610/2 m, we would obtain (, )

Alternatively, we could use the average S wave corner frequency of 2 Hz (fig. S2) to estimate a stress drop (), obtainingwhere VS = 4500 m/s from the PREM (Preliminary Reference Earth Model) model, and k is determined by the source model. In the simple circular model of Madariaga () where k = 0.21, Kaneko and Shearer () calculated k = 0.26 for a circular source and k = 0.19 for an asymmetrical ellipse, both with a lower rupture velocity (0.7VS). These values gave stress drop estimates of ~90, 50, and 120 MPa, respectively. However, note that our rupture velocity estimates were lower, and thus the stress drop listed here should be considered as minimum.

Although it may be preferable to work with more directly observable parameters, such as potency and strain drop (), we followed common seismological practice and used seismic moment and stress drops. Stress drop variations as a function of depth may be explained by the larger shear modulus at greater depth, but the equivalent strain drop (~0.02) of the 2013 Wyoming earthquake was also larger than the strain drop of shallow crustal earthquakes (). The stress drop difference reported here cannot be explained solely by the rigidity differences at depth.

Calculating geotherms

We estimated the depth-dependent temperature profile using the grain-size viscoelastic relaxation model of olivine in the mantle () constrained by surface heat flow and crustal thickness. Our approach was that of a forward calculation in which we generated a large number of geotherms and then, using the viscoelastic relaxation model, calculated the corresponding shear wave velocity as a function of depth and compared it with published velocity models in the mantle (, ). We generated random geothermal models with physical constraints, including a surface heat flow between 40 and 60 mW m−2, a crustal thickness of 50 km (), and thermal conductivity in the crust of 2.5 W m−2, and a temperature-dependent conductivity in the mantle (). For each geotherm, a velocity profile was predicted. The range of these predicted shear wave models that best fit the data and their corresponding geotherms are shown in Fig. 3 (and fig. S1) for both velocity models used [blue () and red () areas]. Our modeling was not capable of reproducing the strong positive velocity gradient between depths of 50 and 80 km, suggesting that shear wave velocities are affected by factors other than temperature ().