Seasonal water storage, stress modulation, and California seismicity

See allHide authors and affiliations

Science  16 Jun 2017:
Vol. 356, Issue 6343, pp. 1161-1164
DOI: 10.1126/science.aak9547

Saving earthquakes for the wet season

Earthquakes can be triggered by changes in crustal stress, such as variations in fluid pore pressure. As a result, the alternating wet and dry cycles in earthquake-prone California should affect the earthquake rate. Johnson et al. asked whether this is indeed the case by combining detailed earthquake records with high-resolution GPS data from the past 9 years. Slight changes in stress did indeed influence the timing of earthquakes, which confirms that the annual hydrological loading cycle modulates microseismicity in California.

Science, this issue p. 1161


Establishing what controls the timing of earthquakes is fundamental to understanding the nature of the earthquake cycle and critical to determining time-dependent earthquake hazard. Seasonal loading provides a natural laboratory to explore the crustal response to a quantifiable transient force. In California, water storage deforms the crust as snow and water accumulates during the wet winter months. We used 9 years of global positioning system (GPS) vertical deformation time series to constrain models of monthly hydrospheric loading and the resulting stress changes on fault planes of small earthquakes. The seasonal loading analysis reveals earthquakes occurring more frequently during stress conditions that favor earthquake rupture. We infer that California seismicity rates are modestly modulated by natural hydrological loading cycles.

Changes in water storage and associated surface loads produce elastic deformation of Earth’s crust (13). In California, the accumulation of winter snowpack in the Sierra Nevada, surface water in lakes and reservoirs, and groundwater in sedimentary basins follow an annual wet and dry cycle. Seasonal vertical and horizontal surface displacements reflect the elastic response of the lithosphere under hydrospheric loads (4, 5) and have been linked to modulating large-scale regional seismic activity (69). Seismicity may also be correlated with annual variations in atmospheric pressure (10), surface temperature (11), snow accumulation (8), and crustal pore pressure (12). Seasonal earthquake modulation provides a natural experiment to investigate the earthquake nucleation process by probing the effects of low-amplitude stress perturbations on active faults (13).

The hydrosphere in California exhibits a strong annual cycle following the wet and dry seasons. Time series from the continuously operating Plate Boundary Observatory and Bay Area Regional Deformation global positioning system (GPS) networks provide dense spatiotemporal resolution of surface loading due to changing water storage (1, 3, 14). Seasonal vertical surface displacement peak-to-peak amplitudes in the Sierra Nevada and California Coast Ranges are 5 to 10 mm and peak late in the year following the dry summer months (2). The recent drought conditions in the western United States amplified the observed uplift signal as Earth’s crust responds to the loss of mass (3). In the Central Valley and other sedimentary basins, the observed seasonal vertical displacements are dominated by the poroelastic response to the head level in the local aquifer system and peak in the spring after the wet winter months (2). Stress changes in the crust from the seasonal water loads in the Sierra Nevada, Coast Ranges, and Central Valley may drive the periodic seismicity patterns, thus providing a physical connection to modulated earthquakes.

In central California, periodic signals in earthquake catalogs indicate that 12 months and the harmonics are statistically significant and document the modulation of earthquakes at natural loading cycles (15). Laboratory fault experiments simulating periodic loading conditions suggest that seismicity should more strongly correlate with the stressing rate if the stressing period exceeds the earthquake nucleation time (13). We expect a weaker correlation with the stress amplitude as the stressing period decreases, indicating a nucleation time exceeding the loading period (13). Numerical simulations predict that the strongest correlation should occur at periodic loading near the characteristic nucleation time (16). Earthquake triggering by the short-period daily tides, which are comparable in amplitude to the seasonal loads, suggests a very weak correlation of seismicity rates with tidal stresses (17), except in regions with unusually large tides (18). This supports the notion of earthquake nucleation times longer than the hourly tidal cycles, thereby obfuscating a detectable earthquake response to the short-period loading (13, 16). We might expect to find a correlation with the seasonal stressing rate if earthquake nucleation times for California faults are less than 1 year. We expect a somewhat weaker correlation with stress amplitude if annual periods are shorter than the nucleation time scale.

Determining the terrestrial water storage is an inverse problem that estimates the change in the surface mass load from the observed vertical deformation (1). To avoid large seasonal signals associated with deforming local aquifers, we examine the annual vertical displacement amplitudes and remove stations that exhibit uplift in the winter and subsidence in the summer. This includes many stations located in the Central Valley of California, where large volumes of groundwater are pumped seasonally for agricultural purposes and the GPS peak amplitude is anticorrelated with the Sierra Nevada and Coast Ranges modulation (2). The data from the remaining 661 GPS sites are averaged to monthly time series of vertical motions from January 2006 to December 2014. Changes in elevation are mapped to a surface mass load using Green’s functions for a spherical gravitating Earth model (19). The inversion method (14, 20) estimates surface mass loads on a 0.25° by 0.25° grid from 32°N to 42.5°N and 125°W to 115°W (movie S1). Knowing the time-dependent hydrological surface loads, we can calculate the monthly stress variations at seismogenic (8 km) depth in a stratified spherical elastic Earth (9). There are competing effects of other sources of seasonal forcing—i.e., atmosphere, temperature, and tides—that could further change the state of stress on the fault (10, 11, 21) but are generally smaller in amplitude during the annual water-loading cycle (20). Our focus is seasonal forcing, and the long-term trend is removed from the stress time series (20).

Our analysis attempts to provide the best hydrologically induced stress-change estimate by selecting the most plausible failure plane for tectonically driven earthquakes. Detailed investigations of focal mechanisms in California reveal a heterogeneous stress field and variations in faulting style (22, 23) that need to be considered when investigating the effect of applied stress transients on fault activity. Plate boundary deformation in central California involves right-lateral strike-slip faulting along the San Andreas Fault (SAF) system and distributed regional shortening in fold and thrust belts of the Coast Ranges (24, 25). To fully capture the complex fault structures in California, we examine a declustered (26) catalog of focal mechanisms [magnitude (M) ≥ 2.0, N = 3,689] in the region 34.5°N to 41.5°N and 124.5°W to 117°W between January 2006 and December 2014 (fig. S1) (20). The focal plane geometry is calculated using first-motion data recorded by the Northern California Earthquake Data Center (Fig. 1). We exclude events from two geothermal and volcanic centers where variation in the seismicity is mainly a response to large-volume fluid injection (27) or the upwelling of volcanic fluids (28). The focal mechanism nodal plane selection process is nontrivial; while shear-stress changes are identical on the two orthogonal planes, normal stresses generally are not. We address ambiguity in the focal mechanism geometry by selecting the preferred plane using a geometric approach. For the strike-slip mechanisms, we select the plane that is best aligned with the SAF (325°) for a geometry consistent with the tectonic environment. For the normal and reverse mechanisms, we assume Andersonian geometry and select the plane with a dip angle closer to 30° for reverse events and 60° for normal events. Using the preferred focal geometry, the shear, normal, and Coulomb (μ = 0.4) stress is calculated for each earthquake (Fig. 1).

Fig. 1 Northern California focal mechanism catalog.

The hydrological surface loading Coulomb stress change (μ = 0.4) is indicated by the color of the shaded focal mechanism compression quadrant for M ≥ 2.0 events from 2006 to 2014. Events with a stress change <0.05 kPa are not shown.

If slight changes in static stress influence the timing of earthquakes, either promoting or discouraging nucleation, then earthquakes will occur more often during slip-encouraging loading conditions. We quantify the percent excess seismicity (NEx) using 0.5 and 6.0 kPa/year stress intervals centered on zero for the resolved stress amplitude and stressing rate, respectively (20). The focal mechanisms are separated into populations of dip-slip events containing the reverse, normal, and oblique earthquakes (N = 2345) and strike-slip (N = 1344) events (Fig. 2 and figs. S2 and S3). Figure 2 shows the results for the shear-stress amplitude and rate for the dip-slip and strike-slip events. A positive trend of 22.5 ± 18.2 NEx/kPa (at 95% confidence) between –1.1 and 1.6 kPa is observed for the shear-stress amplitude of the strike-slip events, indicating a seismic response to the peak shear amplitude. The dip-slip events exhibit a response to the shear stressing rate with a positive trend of 12.5 ± 5.9 NEx/kPa/year, ranging from –11.8 to 26.7 kPa/year. Excess events are not observed for the strike-slip events from an increase in the shear stressing rate, and, similarly, the dip-slip events do not show excess seismicity as the shear stress amplitude increases. Interestingly, in all tests performed, the excess seismicity is ~0% for the near-zero stress interval. We repeat the analysis for a series of sensitivity tests (20) and assess our assumptions regarding magnitude cutoff threshold, focal plane ambiguity, and the error in the focal plane solutions (figs. S4 to S6). The population of M ≥ 2.3 (N = 1861) events indicate similar results for the strike-slip events (N = 762) but no significant trend for the population of dip-slip events (N = 1099) (fig. S4). The NEx values indicate that earthquakes are occurring at different times of the shear-stress cycle for different fault types, suggesting that the strike-slip events are mechanically different from the dip-slip population.

Fig. 2 Percent excess seismicity for shear-stress amplitude and shear stressing rate.

(A) Dip-slip events NEx values for shear-stress amplitude. (B) Strike-slip events NEx values for shear-stress amplitude. (C) Dip-slip events NEx values for shear stressing rate. (D) Strike-slip events NEx values for shear stressing rate. Bin interval is 0.5 kPa for shear-stress amplitude and 6.0 kPa/year for shear stressing rate. The error bars show the ±2 SD from bootstrapping. The minimum and maximum stress and stressing rate are shown numerically for the outer bins. The slope of the best-fit line is shown in the top left of each panel with the 2-SD error. The chi-square value is shown for the goodness of fit of the best-fit line.

The dominant faults in California belong to the northwest-striking SAF system consisting mostly of right-lateral strike-slip and off-fault dip-slip structures. Excess dip-slip events are resolvable during times of peak stressing rate, which agrees with laboratory and numerical modeling (13, 16) if the nucleation time is less than 1 year. This follows a Coulomb failure threshold–triggering model that predicts slip to occur with no lag time at the highest stressing rates. However, the population of strike-slip earthquakes indicates excess events when the peak stress amplitude is greatest (Fig. 2). This could either suggest that the nucleation time is longer for the strike-slip population of events or that there is an increased time lag from the peak stressing rate due to higher tectonic loading rates (16). The strike-slip faults associated with the central SAF system are accumulating strain at rates 8 to 10 times as high as the dip-slip faults (24), consistent with the time-lag scenario.

The results do not indicate a trend with changes in the normal stress, independent of fault type (figs. S2 and S3). The major SAF system of strike-slip faults are postulated to be frictionally weak (29), and we expected a reduced sensitivity to normal-stress changes. Laboratory measurements, heat flow, and regional principal stress-orientation data all support friction values between 0.1 and 0.2 on the major strike-slip faults (29, 30). Our results support the weak-fault hypothesis for faults in California. However, we cannot assume that all faults in California are inherently weak due to varying geology and the overall maturity of different fault strands. If we perform the analysis and select the fault plane with maximum tension, we do find a positive trend with the fault unclamping stress, but cannot say with confidence that the plane selected is the correct failure plane and that a true response to the fault unclamping is occurring.

We calculated seasonal stress time series using the Uniform California Earthquake Rupture Forecast model geometry (UCERF3) (31), a multiyear community effort that provides earthquake rupture probabilities on California faults. For each fault location, we fit a seasonal model to the stress time series and obtain model coefficients to calculate an average annual stress cycle. We assume that this represents the monthly stress change for the past 230 years. Figure 3 shows the annual peak-to-peak Coulomb stress amplitude on faults throughout northern California [see fig. S7 for the corresponding shear, normal, and Coulomb (friction values μ = 0.1 and μ = 0.7) stresses]. The change in Coulomb stress is <1 kPa on most strike-slip faults in the SAF system. The stress cycles on oblique and dip-slip faults in the eastern Coast Ranges and flanking the Sierra Nevada are larger. On many faults the peak slip-encouraging shear and normal stresses are out of phase by 30 to 180 days (fig. S8). We assess the timing of historic seismicity from 1781 to 2012 for M ≥ 5.5 earthquakes (N = 137) and only consider events not identified as foreshock or aftershocks (32). We assumed that the UCERF3 geometry of the nearest fault location is representative of the historic event to evaluate the stress for the month the earthquake occurred. Our simplified treatment of the historical seismicity records suggests that more events occur during times of slip-encouraging failure conditions with a trend of 14.5 ± 9.6 NEx/kPa for the Coulomb (μ = 0.4) stress change (Fig. 4). The more recent earthquakes occurring during the instrumental era (1932 to 2012) produce similar results, with more events during periods of inferred positive Coulomb seasonal stress changes (fig. S9). Although the geometry and hydrological stress conditions of large earthquakes are less well known, the historic catalogs suggest that large earthquakes might also be triggerable from low-amplitude seasonal stress changes. We infer that annual hydrospheric loading is a contributing factor in the modulation of microseismic earthquakes in California.

Fig. 3 Northern California faults in the UCERF3 fault model shown with the annual peak-to-peak Coulomb stress change from hydrological loading.

Fig. 4 Percent excess seismicity from 1781 to 2012 for M ≥ 5.5 events (N = 137), using the UCERF3 historic seismic records and the average monthly annual stress estimated from the closest fault–model geometry.

The error bars show the ±2 SD from bootstrapping. The minimum and maximum extent of stress values is shown numerically for the outer bins. The slope of the best-fit line is shown in the top left of each panel with the 2-SD error. The chi-square value is shown for the goodness of fit of the best-fit line.

Supplementary Materials

Materials and Methods

Figs. S1 to S9

Movie S1

References (3353)

References and Notes

  1. Materials and methods are available as supplementary materials.
Acknowledgments: Funding for this project was provided by the U.S. Geological Survey National Earthquake Hazards Reduction Program (G15AP00106) and the Southern California Earthquake Center. This material is based on work supported by the National Science Foundation Graduate Research Fellowship under grant number DGE1106400 for C.W.J. GPS data were collected by the Plate Boundary Observatory operated by UNAVCO for EarthScope ( and supported by National Science Foundation grant no. EAR-0350028 and EAR-0732947. The data were processed using the GIPSY OASIS II software provided by the Jet Propulsion Laboratory. Waveform data, metadata, or data products for this study were accessed through the Northern California Earthquake Data Center (NCEDC), doi:10.7932/NCEDC. The monthly water-loading estimates and earthquake catalog with the calculated stresses are all available in the supplementary materials.
View Abstract

Navigate This Article