Spatially Consistent Small‐Scale Stress Heterogeneity Revealed by the 2008 Mogul, Nevada, Earthquakes

The relatively large numbers of small earthquakes have the potential to provide detailed information about active fault geometries and their stress state. Precise relocations of earthquake sequences, including swarms, foreshocks, and aftershocks, are revealing complex subsurface geometries and seismicity migration in both tectonic (e.g., Ruhl et al., 2016; De Barros et al., 2020; Shelly et al., 2023) and induced settings (e.g., Shapiro et al., 2002; Bourouis and Bernard, 2007; El Hariri et al., 2010). Systematic analyses of earthquake stress drops provide additional information regarding earthquake behavior, identifying spatially coherent stress‐drop variations at both regional (e.g., Shearer et al., 2006, 2022; Oth, 2011; Zhang et al., 2022) and local scales (e.g., Abercrombie et al., 2017; Ruhl et al., 2017; Pennington et al., 2021). However, in most cases it is unclear what is causing these stress‐drop variations and what they may reveal about seismicity evolution, fault interactions, and hazard estimation. The large uncertainties, both random and systematic, in spectral estimates of stress drop remain a challenge (e.g., Abercrombie, 2021; Bindi et al., 2023a,b), making multimethod confirmation of any reported variations valuable (e.g., Pennington et al., 2021).

Earthquake swarms are of particular interest, because they are indicators of physical driving mechanisms that promote earthquake triggering (e.g., Shapiro et al., 2002; Bourouis and Bernard, 2007; El Hariri et al., 2010; Ruhl et al., 2016; De Barros et al., 2020; Shelly et al., 2023). Here, we focus on an exceptionally well‐recorded (station spacing 1–2 km) shallow tectonic swarm beneath Mogul, Nevada (Ruhl et al., 2016). The previous work (Ruhl et al., 2017) revealed intriguing spatial variations in stress drop on the main fault plane using an empirical Green’s function (EGF) method applied to about 150 of the largest earthquakes (⁠$ML>2.1$⁠). Applying a new spectral decomposition approach (Shearer et al., 2022) that reduces uncertainties associated with attenuation and other path corrections, we compute stress drops for over 4000 earthquakes in the Mogul swarm. We confirm the main results of Ruhl et al. (2017), and obtain additional details regarding the temporal and spatial variations of small‐magnitude earthquake stress drops during the Mogul sequence, including foreshocks and aftershocks of the swarm’s largest event (⁠$Mw$ 4.9, $ML$ 5.1).

The 2008 Mogul, Nevada, earthquake sequence was an unusually shallow, energetic foreshock–mainshock–aftershock sequence that occurred in the transtensional Walker Lane region just west of Reno, Nevada (Fig. 1a). Two months of shallow (<5 km; Ruhl et al., 2016) foreshocks (up to $ML$ 4.2) were widely felt and led to the deployment of a dense network of temporary broadband stations, producing exceptional recording of the sequence. The largest earthquake, an $Mw$ 4.9 strike‐slip event, occurred on 26 April 2008 at ∼3.5 km depth (Fig. 1a,b). It ruptured the previously unknown (despite its shallow depth) northwest–southeast‐trending strike‐slip fault in a region of mapped, mostly southwest–northeast‐trending, normal faults. The total geodetic deformation observed was approximately twice that of the cumulative seismic moment, implying a large amount of aseismic slip (Bell et al., 2012). The timing of the aseismic slip is largely postseismic but could also have occurred during the most intense foreshock phase in the week before the largest earthquake.

Ruhl et al. (2016) relocated 7500 earthquakes using waveform cross correlation and computed over 1000 focal mechanisms (subset shown in Fig. 1a), which revealed significant spatiotemporal complexity and clustering within the sequence. Foreshock hypocenters spread out spatiotemporally (Fig. 1c), first along a fault parallel to and southwest to that of the main event, then moving to the eventual mainshock fault plane, ultimately highlighting the bulk of the mainshock fault zone prior to mainshock rupture. The well‐resolved hypocenters and mechanisms reveal a complex fracture mesh of interacting strike‐slip and normal faults; the migration rates are broadly consistent with pore‐pressure diffusion, but individual clusters exhibit much higher migration rates consistent with aseismic slip (Ruhl et al., 2016).

To explore further the complexity of the fault system and stress field, Ruhl et al. (2017) exploited the high‐resolution data recorded by the dense station coverage, and estimated source dimensions and static stress drops (⁠$Δσ$⁠) for ∼150 of the largest earthquakes (Fig. 1d–g) using a spectral ratio (EGF) approach. The results revealed distinct spatial and spatiotemporal patterns that were consistent for independent measurements obtained separately for P and S waves. The clearest relation between seismicity location, timing, and stress drop was on and around the mainshock fault plane, for which the foreshocks revealed a region of relatively high‐stress drop to the north of the mainshock hypocenter and an area of relatively low stress drop to the southeast (Fig. 1d,e). There were few aftershocks in the region of high‐stress‐drop foreshocks, and the aftershocks in the region of lower stress‐drop foreshocks had similarly low values (Fig. 1f,g). These observations suggested that the spatial patterns of stress drop were unaffected by the mainshock, and that the seismicity void surrounding the high‐stress‐drop foreshocks may represent a distinct fault patch that ruptured either aseismically or during the mainshock. Ruhl et al. (2017) also found no relation between stress drop and depth or focal mechanism but did find the events off the mainshock fault plane tended to have lower stress drop. However, the relatively small number of events well‐resolved by the EGF method limited resolution of such patterns.

Here, we further explore the complex spatiotemporal stress‐drop variations within this sequence by calculating source parameters for over 4000 of the Mogul earthquakes using a spectral decomposition method. In contrast to the EGF approach, this method analyzes larger numbers of events using the modified spectral decomposition approach developed by Shearer et al. (2022) to minimize bias and resolve reliable relative spatial variations. We compare the resulting corner frequencies of common events between this study and Ruhl et al. (2017), and find the new results confirm and provide further resolution of the previously reported patterns.

We start with the catalog of earthquakes relocated by Ruhl et al. (2016) and compute P‐wave displacement spectra at all available stations within 80 km using the method described in Shearer et al. (2022). We also compute a noise spectrum for each record from a time window of the same length, immediately before the P‐wave arrival. We only use P‐wave spectra with an average signal‐to‐noise ratio (SNR) greater than five over the frequency band 2.5–6 Hz. To avoid selection bias effects, we do not apply an SNR cutoff at higher frequencies (Shearer and Abercrombie, 2021). Next, we apply spectral decomposition to solve for best‐fitting event, station, and travel‐time terms (e.g., Shearer et al., 2006; Oth et al., 2011). We assume that the catalog magnitude is equal to moment magnitude (⁠$Mw$⁠) at M 3.5, and use the best‐fitting linear relationship between catalog magnitude and the spectral amplitudes of the event terms at 2.7 Hz to estimate seismic moment (⁠$M0$⁠) for each event term.

We stack the event terms at intervals of 0.25 in relative log moment, as plotted in Figure S1, available in the supplemental material to this article. For comparison, we also plot stacks of the pre‐event noise obtained using the same processing scheme. The event terms are well above the average noise levels out to frequencies of at least 35 Hz, and that the noise increases sharply beyond 40 Hz. The event stacks plotted in Figure S1 describe only the relative shapes of the averaged spectra as a function of event size, reflecting the nonuniqueness in spectral decomposition resulting from trade‐offs among the event, station, and path terms. To obtain true source spectral estimates, the event terms must be corrected by an empirical correction spectrum (ECS). However, as documented in Shearer et al. (2019) and Abercrombie et al. (2021), many previous methods for estimating the ECS function suffer from large uncertainties because of the limited bandwidth and dynamic range of the data, the number of free parameters involved, and the consequent simplifying assumptions required.

To avoid these problems, we apply the approach described by Shearer et al. (2022), in which the average corner frequency of small earthquakes is fixed to a constant value. This method stabilizes the problem, and avoids the parameter tradeoffs and large uncertainties in empirical path corrections inherent to many previous approaches. Fixing the small earthquake properties removes the possibility of resolving spatial variations in their stress drops, but ensures that any variations seen in larger events are real and not the result of inaccurate path corrections. More details about the method are described in Shearer et al. (2022), who applied it to estimate stress drops across southern California. However, because the Mogul sequence is compact (less than 8 km in extent) and contains high‐quality records to magnitudes as low as zero, we adjusted some details of the method for our dataset.

Following Shearer et al. (2022), we estimate relative moment from the low‐frequency (2.67 Hz) part of the event terms, and find a roughly linear relationship with unit slope between catalog magnitude (⁠$Mcat$⁠) and relative moment (Fig. S2). This implies that $Mcat=1.5Mw$ (moment magnitude) for earthquakes below M 3, as prior studies have shown (e.g., Hanks and Boore, 1984). Assuming $Mcat=Mw$ at $Mcat$ 3.5, we then convert our relative moment estimates to true moment. Shearer et al. (2022) set $fc=30Hz$ for a reference catalog magnitude (⁠$Mref$⁠) near 1.5, based on observations from the deep Cajon Pass borehole seismometer (Abercrombie, 1995; Shearer and Abercrombie, 2021), located below the highly attenuating near‐surface layers. Because the Mogul sequence contains high‐quality records to magnitudes as low as zero, we used $Mref∼0.4$⁠. Assuming self‐similarity and that $Mw$ is proportional to 2/3 $Mcat$⁠, we fix $fc$ to 70 Hz at $Mcat$ ~0.4 to yield the same stress drop as $fc=30Hz$ at $Mcat$ ∼1.5. Shearer et al. (2022) searched for calibration events (⁠$M∼Mref$⁠) within 2 km in depth and 5 km in horizontal distance for every earthquake in their database. Because the Mogul sequence is compact (extending only about 8 km horizontally and 3 km in depth), and we were particularly interested in possible depth variations in stress drop, we defined our calibration events as within 1 km in depth and did not implement any horizontal distance cutoff.

For each calibration event (0.3 < M < 0.5), we compute a Brune spectrum (Brune, 1970) based on our assumption that $fc=70Hz$ at M ~0.4. We estimate the ECS for each target event as the average difference between the event‐term log spectra and these Brune model predictions, following smoothing with a fourth‐order polynomial. Next, we subtract the ECS from the log event term and fit this corrected spectrum with a Brune model (with an assumed high‐frequency fall‐off rate of 2) to obtain an estimate of $fc$⁠. From our $M0$ and $fc$ estimates for each event, we compute earthquake stress drops assuming the circular rupture model of Kaneko and Shearer (2014) and a fixed S‐wave velocity of 3.46 km/s. Requiring at least three contributing spectra per event, we obtain moment and stress‐drop estimates for 4175 earthquakes in the Mogul sequence.

Figure 2 shows a comparison of seismic moments estimated and/or used by the two studies as well as the corner frequencies for 163 common events (P‐wave estimates of Ruhl et al., 2017). We have generally good agreement between corner frequency estimates, especially at the lower magnitudes (with generally higher corner frequencies). The seismic moment in the spectral decomposition is underestimated for the largest events (M > 3.5), as expected given the frequency range of analysis (Fig. S1 and Shearer et al., 2022), resulting in likely overestimation of their corner frequencies (Fig. 2). Because there are relatively few M > 3.5 events, retaining them in the figures does not affect our interpretations, which are based on spatial and temporal averages over many events. The average or median stress drop in the complete data set depends strongly on the ECS (e.g., Shearer et al., 2019) and the assumption that average $fc=70Hz$ for M ~0.4 events. Thus, we focus our interpretations on the relative values of stress drop, which are more reliably determined. Variations in stress drop among individual M ~0.4 events can be resolved even though their average $fc$ is fixed.

The foreshock hypocenters began in a well‐defined fault fracture mesh to the west of the mainshock fault zone (Figs. 1a,3a) and spread diffusively from the initiation point. Activity on the mainshock fault plane initiated at a depth of about 4 km, and spread out first in the northwest part of the fault and then extended to the southeast (Ruhl et al., 2016). The earliest foreshocks on the mainshock fault plane form a dense ring of relatively high‐stress drops outlining an area of lower seismicity (red double ellipse in Fig. 3a). The foreshocks to the southeast of the mainshock hypocenter tend to have lower stress drops (purple ellipse in Fig. 3a) in an overlapping depth range. Few aftershocks occur in the region of the earliest foreshocks, but the aftershocks on and around the mainshock fault plane show a similar pattern with higher stress drops to the northwest and lower to the southeast.

This distinct and temporally stable spatial pattern of stress drop is similar to that previously observed by Ruhl et al. (2017), see Figure 1 and Figure S5. The spatiotemporal patterns obtained here confirm the observations presented by Ruhl et al. (2017), and support the interpretation that the source complexity and spatiotemporal variations are real. The high spatial variation of stress‐drop estimates is independent of timing and fault orientation. We also find that the spatial variations are robust with respect to the number of stations required to contribute spectra for each event. Imposing a more stringent requirement of five stations compared to three reduces the number of events, but the spatial pattern of stress drop remains unchanged (see Fig. S3). The overall median value of stress drop depends on the assumed value of $fc=70Hz$ at M ~0.4. A higher or lower value than 70 Hz will shift the stress drops to correspondingly higher or lower values (see Fig. S4) for lower magnitude events. However, the results converge at higher magnitudes, and spatial variations in relative stress drop are robust with respect to changes in the fixed $fc$ value.

Figure 4 shows the complex spatial and temporal evolution of the Mogul sequence in three depth and time intervals, together with the available focal mechanisms. The region of high‐stress drops on the mainshock fault plane is the highest in the entire sequence. Clear smaller scale patterns are also revealed. For example, the small northeast–southwest‐trending fault to the southwest of the mainshock fault plane shows significant variability in stress drop, despite having consistent strike‐slip focal mechanisms (longer term aftershocks, >3 km depth; Fig. S8). Parallel to this fault, some of the other small structures tend to be mainly high‐ or low‐stress drop. Two structures activated within the foreshock sequence (>4 km depth) show contrasting stress drops, with higher values on a normal fault and lower on a strike‐slip fault. To determine whether there is any consistent relationship between fault orientation, stress field, and resulting stress drop, we calculate the slip tendency of the optimal nodal plane in each focal mechanism following Jansen et al. (2019) and Qin et al. (2022). We observe a very weak decrease in spectral stress drop with higher slip tendency, but it is orders of magnitude smaller than the total stress‐drop variability (Fig. S6).

We observe strong and consistent spatial variations in spectral stress drop for the Mogul sequence that confirm the previous findings of Ruhl et al. (2017). The spatial variability is stable with time, and no strong temporal variation is observed. Neither these spatial variations show any simple correlation with fault orientation or sense of slip or slip tendency (under average stress field), nor do we observe any systematic variation with depth (Fig. S7), consistent with the previous findings of Abercrombie et al. (2021). The scale of the variation we resolve is of the order of a few kilometers, similar to that seen by Shearer et al. (2019) for aftershocks of the 1992 Landers earthquake, by Folesky et al. (2021) for aftershocks to the Iquique earthquake, northern Chile, and by Allmann and Shearer (2007) and Zhang et al. (2022) for the San Andreas fault at Parkfield. In general, these and multiple other studies are findings that the size of small‐scale variations in stress drop is of equal or greater magnitude than any larger scale differences in the same region (e.g., Abercrombie, 2021; Shearer et al., 2022).

A key question is whether these variations are stable in time or evolve in response to changes in background stress resulting from moment release on faults. Because seismicity patterns evolve with time, it can be difficult to separate temporal and spatial stress‐drop variations (e.g., Chen and Shearer, 2013; Folesky et al., 2021). At Parkfield, the high rate of seismicity and excellent recording networks enabled both Allmann and Shearer (2007) and Zhang et al. (2022) to resolve significant spatial variations in stress drop at Parkfield that were largely unchanged after the 2004 M 6 Parkfield earthquake. They concluded that spatial variations at Parkfield are significantly larger, and thus better resolved than any temporal variations in agreement with our observations for the Mogul sequence. Temporal variation has been seen in detailed analysis of some repeating earthquake sequences (e.g., Cauchie et al., 2020; Chaves et al., 2020). Yoshida et al. (2017) were more confident that the temporal changes were real in the Yamagata–Fukushima swarm following the 2011 M 9 Tohoku–Oki earthquake as the greater density of earthquakes eased comparison of collocated earthquakes from different time periods.

In general, our Mogul results are consistent with the previous studies that have found persistent small‐scale stress‐drop heterogeneity in and among fault zones that does not noticeably evolve with time. There are various possible causes for such small‐scale, persistent spatial variability. Mogul is an unusually shallow sequence with complex, discrete fault structures between 2 and 5 km depth. Despite the clear fault‐plane variation, we do not observe any strong dependence on fault orientation or slip sense. We find that significantly different stress drops occur among and along parallel fault structures (e.g., Fig. S8), indicating that inherent fault‐zone properties (rather than fault orientation) may be responsible for variations in average rupture characteristics. Other structural explanations for spatially varying earthquake stress drop include fault complexity (e.g., Kemna et al., 2021), fault roughness, and width of the fault damage zone (e.g., Moyer et al., 2018). In laboratory experiments, higher stress drops were reported on relatively smoother fault surfaces (Blanke et al., 2021). The part of the Mogul mainshock fault plane with higher stress drops may be a narrower, simpler structure (Fig. 4), consistent with this interpretation.

Variation in pore pressure is another possibility to explain source parameter variation. The seismicity migration patterns and nearby hot springs (Ruhl et al., 2016; Fig. 1a) suggest the involvement of pore fluids in the Mogul earthquake sequence (Jansen et al., 2019). However, strong dependence of stress drop on pore pressure has only been found very close to injection sites at higher pore pressures than we have any evidence for here (e.g., Goertz‐Allmann et al., 2011).

Detailed analyses of tectonic and induced earthquake sequences are demonstrating the importance of small‐scale heterogeneity in their development. Comparative studies of spatiotemporal seismicity patterns are providing insight into the relative roles of pore pressure and aseismic slip as driving mechanisms (Ruhl et al., 2016; Danré et al., 2022; Fischer et al., 2023). Adding robust source parameter measurements to these characterizations will enable a better understanding of the principal factors governing the progress of dynamic failure and potentially leading to improved temporarily varying seismic hazard estimation.

The fine spatial resolution of earthquake properties that we obtain for the Mogul sequence was made possible by the very dense station coverage immediately above the earthquakes, highlighting the advantages of rapid field campaigns to monitor seismicity episodes.

All waveform data were obtained through the Nevada Seismological Laboratory (NSL) and are also available for download from the Earthscope Consortium Incorporated Research Institutions for Seismology Data Management Center (IRIS‐DMC): University of Nevada, Reno (1971; doi: 10.7914/SN/NN); University of Nevada, Reno (1992; doi: 10.7914/SN/SN). Earthquake relocations, clusters, and focal mechanisms are published in Ruhl et al. (2016), and the previous stress drops in Ruhl et al. (2017). We provide our estimates of corner frequency and stress drop in Table S1 (available in the supplemental material to this article).

The authors acknowledge that there are no conflicts of interest recorded.

Funding for this research was provided by the National Earthquake Hazards Reduction Program/U.S. Geological Survey (NEHRP/USGS) Award Numbers G21AP10125, G21AP10124, and G20AP00027 and by the Southern California Earthquake Center (SCEC) Award Numbers 21082 and 20097. This research was supported by the SCEC (Contribution Number 12735). The SCEC is funded by the National Science Foundation (NSF) Cooperative Agreement EAR‐1600087 and USGS Cooperative Agreement G17AC00047. The authors thank two anonymous reviewers, and TSR Editor‐in‐Chief Keith D. Koper and Associate Editor Ruth Harris for their help in improving the article for publication.