Hydrospheric mass changes exert direct influence over lithospheric deformation. Both seasonal and long-term changes to ice, snow or water loads may induce displacements of the Earth’s surface^{11, 12, 13, 14, 15} and can create stress perturbations that modulate activity on seismogenic faults^{16, 17, 18, 19}. A volume of groundwater approaching approximately 160 km³ in California’s Central Valley has been lost through pumping, irrigation and evapotranspiration over the past 150 years or so^{3, 4}. Historical and modern records demonstrate that groundwater depletion occurs primarily in the drier, hotter, southern portion of the basin^{1, 3} (the San Joaquin Valley), parallel to the central San Andreas Fault and adjacent to the high topography of the southern Sierra Nevada (Fig. 1). Previous studies demonstrating sensitivity to small-scale stress changes across this section of the San Andreas point to seasonal hydrologic⁶ or temperature²⁰ variations to explain observed changes in seismicity rates⁶ or strain in shallow boreholes²⁰.

Modest contemporary uplift rates of the southern Sierra Nevada observed using space geodesy^{7, 8, 9} have been attributed to various tectonic and geomorphic drivers, such as buoyant response to mantle delamination¹⁰, Basin and Range extensional faulting²¹, and erosional mass transfer from Pleistocene glaciation²². Inferred stress change from epeirogenic uplift in the southern Sierra Nevada is also potentially linked to the transition from locked to creeping sections of the San Andreas Fault²³. Here we seek to explore potential human impacts on contemporary deformation across this region through global positioning system (GPS) constraints on lithospheric flexure induced by anthropogenic groundwater level changes.

Continuous GPS networks spanning the southwestern USA provide a high-resolution framework for analysing crustal motion in three dimensions (Fig. 1). Longer available records and recently improved processing techniques enable determination of reliable GPS vertical velocities with a precision²⁴ of less than a millimetre per year, tied to the Earth system’s centre of mass to within 0.5 mm yr⁻¹ (refs 24 and 25). GPS station data were processed to create individual time series of vertical position and were then fitted with an empirical model including epoch position, velocity, and the amplitude and phase of annual and semiannual harmonic components (to model seasonal effects). This analysis includes only stations with at least two-and-a-half years of data and vertical rate uncertainties of ≤1.0 mm yr⁻¹ (566 stations). Details of the GPS analysis and processing techniques are included in the Methods.

Figure 1 shows the distribution of vertical GPS rates across central California and east of the Sierra Nevada into the Great Basin. Long-term rates exclude GPS stations located within the San Joaquin Valley groundwater basin that display comparatively large signals owing to their position on soft sediment. These stations are affected by compaction-driven subsidence through loss of pore water⁵, poro-elastic deformation related to groundwater level fluctuations²⁶, local irrigation and other processes unrelated to solid earth motion.

Vertical rates across a transect from the central Coast Ranges to the relatively stable Great Basin range up to about 3 mm yr⁻¹. The highest average velocities concentrate along the margins of the valley (Fig. 2), near the locus of the greatest historical changes in the deep confined aquifer¹ (Fig. 1). The uplift rates of stations west of the San Joaquin Valley are more variable than those at bedrock sites in the Sierra Nevada, reflecting a higher proportion of stations on or near agricultural basins or active faults in the central Coast Ranges. Scatter in the vertical rate transect probably reflects shifts in the locus of groundwater change in comparison with historical averages², as well as some aleatory variability. Nearby earthquakes, such as the San Simeon rupture in 2003 of moment magnitude 6.5, may also influence stations near the coast with anomalously high rates near 2 mm yr⁻¹ (Fig. 2). In the Sierra Nevada, the highest vertical rates occur along the western slope and steadily decay to the northeast in the Great Basin (Fig. 2). Previous studies show a similarly poor or inverse correlation between elevation and GPS uplift rates across the Sierra Nevada^{7, 8, 9}. Modelling of postseismic viscoelastic relaxation following nearby historical earthquakes suggests that such transient strain accounts for only a fraction of the total observed vertical signal in the southern Sierra Nevada⁹.

We compare vertical GPS rates with surface uplift predicted from an elastic half-space model simulating the response to load changes driven by variations in total water storage (Fig. 2) (see Methods). The model accounts for surface and subsurface deformation induced by line loads distributed across a 60-km-wide strip over the surface of an elastic half-space representing the San Joaquin Valley, the site of long-term, historical groundwater unloading. Using a range of elastic parameters, we fitted the GPS-derived vertical velocities with a rate of unloading of 8.8 (±1.3) × 10⁷ N m⁻¹ yr⁻¹ (Fig. 2). We compare this estimate with the current average unloading rate in the valley based on changes in total water storage, measured using satellite gravimetry from the Gravity Recovery and Climate Experiment (GRACE) (ref. 3; see Methods). This change, measured between October 2003 and March 2010, yields an average equivalent unloading rate of 7.2 (±2.1) × 10⁷ N m⁻¹ yr⁻¹, slightly lower but in good overall agreement with the GPS-derived estimate.

Assuming that changes in total water storage drive vertical motion surrounding the San Joaquin Valley, the GPS data reflect a combination of short-term, elastic response to ongoing groundwater depletion and longer-term viscous relaxation reflecting the history of hydrospheric mass changes. Slight overestimation of the unloading rate based on the GPS results leaves open the possibility of either a time-dependent (viscous) component of the uplift signal, or a (smaller) contribution of ongoing tectonic uplift. Both the overall match between the elastic model and observed GPS uplift (Fig. 2) and seasonal patterns inherent in the phase and amplitude of the GPS vertical time series, however, demonstrate the importance of the instantaneous (elastic) response to groundwater unloading.

Annual peak uplift for GPS stations in the Coast Ranges and Sierra Nevada occurs in the late summer and early autumn (Fig. 3), corresponding with diminished snow and surface water loads²⁷ and overlapping with the end of the summer growing season and peak groundwater pumping in the Central Valley (May to September; ref. 2). In contrast, recharge during the early spring drives larger-amplitude peak uplift within the valley through poro-elastic effects of aquifer water levels (Fig. 3 and ref. 2). Seasonal vertical displacements are more broadly distributed than average, long-term trends (Fig. 2), reflecting the variability in total water storage in upland catchments feeding the San Joaquin Valley. Notably, stations along the valley margins move upward annually in accordance with the Sierra Nevada and Coast Ranges during dry months and show a corresponding subsidence during the wet winter and spring (Fig. 3). These peripheral stations lead to the bimodal distribution for peak uplift times observed for the Central Valley (inset to Fig. 3), indicating that uplift from hydrospheric load changes dominates even the shortest-term signals in the vertical GPS measurements for stations unaffected by local irrigation or aquifer effects.

Flexure in the central Coast Ranges due to seasonal hydrospheric unloading in the San Joaquin Valley provides a viable mechanism to explain the annual modulation of seismicity on the San Andreas Fault. Both the locked and creeping fault sections at Parkfield (Fig. 1) exhibit an increase in the number of earthquakes greater than magnitude 1.25 during the late summer and autumn (Fig. 3), previously attributed to local changes in effective stress linked to the hydrologic cycle⁶. We explore the potential seasonal impacts of unloading and uplift of the Coast Ranges on short-term changes in the fault-normal stress resolved across the fault using the elastic half-space model (see Methods). Taking estimates of annual variations of 100–300 mm in equivalent water height in the San Joaquin Valley from GRACE data³ and a wider load encompassing the Coast Ranges and the western Sierra Nevada, the elastic model produces between 3 mm and 8 mm of maximum annual vertical ground motion, in good agreement with the peak-to-peak amplitudes of annual uplift measured by GPS (Fig. 2). The corresponding fault-normal stress variations on the San Andreas Fault are near 1 kPa at seismogenic depths (Extended Data Fig. 4), with peak unclamping during the dry summer and autumn. Although we cannot rule out the potential feedback of reduced effective normal stress due to diffusion of pore fluids into fault zone rocks during recharge⁶, our results suggest that unloading may contribute to seasonal modulation of seismic activity on the central San Andreas Fault. A similar mechanism was invoked to explain annual variations in seismicity in the Himalayas¹⁸. Stress changes due to groundwater unloading are somewhat higher for other historically active faults closer to the valley such as the Coalinga thrust system (Fig. 1), where peak-to-peak annual unloading gives rise to positive Coulomb stress changes of 1.0 kPa to 1.7 kPa (see Methods and Extended Data Fig. 4).

Given that relatively small annual stress perturbations (around 1 kPa) appear sufficient to modulate earthquakes along the central San Andreas Fault⁶, we use the elastic model to estimate the stress rate caused by ongoing human-induced groundwater removal (see Methods). Taking an estimate of the unloading rate due to loss in water storage in the San Joaquin Valley of 8.8 (±1.3) × 10⁷ N m⁻¹ yr⁻¹, we calculate a rate of unclamping of 0.07–0.18 kPa yr⁻¹ for the San Andreas Fault. Coulomb failure stress rates on the Coalinga thrusts are higher, averaging 0.13–0.57 kPa yr⁻¹. Estimates of total groundwater loss from historical and more recent water level records (around 160 km³; refs 3 and 4) correspond to a reduction in the effective normal stress of 2.7–9.5 kPa for the San Andreas Fault, and a Coulomb stress change of 10–15 kPa at Coalinga, since the beginning of groundwater extraction in 1860. These results suggest that human activity may give rise to a gradual increase in the rate of earthquake occurrence, as suggested by earthquake catalogues in central California⁶.

Future scenarios for groundwater in California suggest increasing demand for agricultural, urban and environmental use². Climate change will probably exacerbate the stress on this resource through altered precipitation patterns, more frequent droughts, earlier snowmelt, larger floods, and increasing temperatures and evapotranspiration^{28, 29}. We demonstrate how long-term and seasonal hydrospheric mass changes due to groundwater removal significantly affect regional crustal deformation, which can lead to annual variations in small-earthquake frequency and longer-term changes of stresses on the San Andreas Fault. Our model explains recent vertical uplift rates previously attributed to epeirogenic deformation in central California¹⁰, indicating that much of this signal results from anthropogenic groundwater removal. Taken together, our results highlight new and underappreciated links between human activity and solid earth processes driven by hydrospheric change.

GPS data processing

We use data from continuously recording GPS stations between −122° and −114° longitude and 34° to 38° latitude that had time series longer than 2.5 years that were >50% complete (Extended Data Fig. 1). This data set includes 566 continuous GPS stations within this geographic region. No episodic campaign or semi-continuous stations were included. Percentage completeness is defined as the number of daily solutions for the station divided by the number of potential daily solutions given the duration of the time series. Larger values for percentage completeness are preferred in order to best resolve seasonal oscillations of position.

The data were processed as a part of a >11,000 station mega-network analysis system that retrieves data daily and updates solutions weekly. We use the GIPSY/OASIS software provided by the Jet Propulsion Laboratory to estimate station coordinates every 24 h using the Precise Point Positioning method³¹. Ionosphere-free combinations of carrier phase and pseudorange were obtained every 5 min. Calibrations were applied for all antennas, ground receivers and satellite transmitters. To model tropospheric refractivity, the Global Mapping Function was applied³², with tropospheric wet zenith delay and horizontal gradients estimated as stochastic random-walk parameters every 5 min (ref. 33). The observable model includes ocean tidal loading (including companion tides) coefficients supplied by Chalmers University³⁴. Ambiguity resolution was applied to double differences of the estimated one-way bias parameters³⁵, using the wide lane and phase bias method, which phase-connects individual stations to IGS stations in common view³⁶. Satellite orbit and clock parameters were provided by the Jet Propulsion Laboratory, who determine these parameters in a global fiducial-free analysis using a subset of the available IGS core stations as tracking sites. Output station coordinates are initially in the loose frame of the Jet Propulsion Laboratory’s fiducial-free GPS orbits. These were transformed into reference frame IGS08 using daily seven-parameter transformations supplied by the Jet Propulsion Laboratory. IGS08 is derived from ITRF2008 (ref. 37) and consists of 232 globally distributed IGS stations³⁸. Finally, the solutions are aligned with our custom reference frame (NA12) that co-rotates with stable North America²⁴. This alignment is applied to each daily solution and provides a spatial filter to suppress errors correlated at the continental scale. The mean formal uncertainty in daily vertical coordinate is <3.0 mm. A summary reference table for GPS solution methodology is given at http://geodesy.unr.edu/gps/ngl.acn. The resulting time series are used to estimate rates of motion with respect to stable North America.

GPS time series analysis

An outlier exclusion criterion was applied to remove daily solutions that were more than 20 mm from the expected position based on a provisional time series model, or having position uncertainty greater than 20 mm. For each of the three components (east, north, up) of each GPS time series we solved for the intercept b, rate v, steps at the times t_eqp of N known equipment changes that can introduce discontinuities of size D. We solve for annual and semiannual seasonal terms by solving for the cosine (C_i) and sine (S_i) terms. For example, the east component is:

where is the Heaviside step function, where t is in years. We use a damped linear inversion with weights on the prior uncertainties of the model parameters scaled to allow the expected range of values for each parameter. However, because only continuously recording GPS stations were used, all parameters in the estimation (including seasonal terms) were well constrained, so the damping parameters had little influence on the solutions. The amplitude of the annual terms were calculated with:

and the day of year from:

Finally, stations were excluded if their rate uncertainty was greater than 1 mm yr⁻¹ or their rate was clearly anomalous (that is, greater than 100 mm yr⁻¹). In all, we estimated solutions for 566 stations, 96% of which have time series duration greater than 4.0 years, and 92% of which are 80% or more complete (Extended Data Fig. 2). Example vertical GPS time series are shown in Extended Data Fig. 3. GPS time series parameters, including vertical uplift rate and the amplitude and timing of peak uplift are included as source data for Fig. 1.

Elastic half-space model

We use an elastic half-space model for calculating the elastic response of the crust to changes in total water storage. The model is two-dimensional and accounts for surface and subsurface deformation induced by a line load over the surface of the half-space along a finite-width strip (half-width a, Extended Data Fig. 4a). Model derivation can be found in ref. 30, with a small correction from ref. 13. The normal displacement of the surface (z =  0) along a profile perpendicular to the load centred at x = x_d is given by

where x and z are the horizontal and vertical dimensions, respectively, N₀ is a line load, ν is Poisson’s ratio, G is the shear modulus, a is the strip half-width and K is a far-field rate offset correction. We fix the position of the line load to the centre of the San Joaquin Valley, and use a load width of 2a = 60 km corresponding to the average width of the valley floor. Given a range of elastic parameters (ν = 0.25, G = 35 ± 5 GPa), we can fit the vertical velocities derived from the GPS data by a line-load rate with = 8.8 (±1.3) × 10⁷ N m⁻¹ yr⁻¹ (see Fig. 2). Goodness of fit is estimated using a reduced chi-square criterion. Our best-fit model has a reduced chi-square of 1.9, indicating a reasonably good fit.

Unloading from estimated mass loss

The combined Sacramento–San Joaquin River basins lost an estimated volume of 30.9 (±2.6) km³ in total water storage over a 6.5-year period between October 2003 and March 2010 (ref. 3). Of this volume, 20.3 (±3.8) km³ is attributed to loss of groundwater in the Central Valley, 80% of which occurred in the San Joaquin Valley. Assuming that changes in total water storage not related to groundwater (for example, soil moisture, surface water and snow) are distributed equally between the Sacramento and San Joaquin basins, the total water storage change in the San Joaquin Valley over the 6.5-year period is 21.5 (±6.2) km³. This change is equivalent to a rate of mass loss of 3.3 (±1.0) × 10³ gigatons per year, and therefore a rate of unloading of 3.2 (±1.0) × 10¹³ N yr⁻¹ distributed over the 450 km length of the valley. Taking into account the range of uncertainties in groundwater and total water storage change³, the equivalent line-load rate is thus = 7.2 (±2.1) × 10⁷ N m⁻¹ yr⁻¹, in agreement with the estimate from the elastic model based on the vertical GPS velocities. If we add the estimated 1860–2003 depletion of 140 km³ (ref. 4) to the 2003–2010 water loss, the cumulative historic load decrease amounts to ~3.5 × 10⁹ N m⁻¹.

Stress modelling

The two-dimensional stress components at any point at depth caused by a distributed line load N₀ (in N m⁻¹, the negative value indicates unloading) at the surface of the elastic half-space are calculated as³⁰:

where θ₁ and θ₂ are the angles from both edges of the load measured clockwise from the positive x direction, a is the load half-width and z is positive downward (Extended Data Fig. 4a). The shear and normal stress components resolved on a fault plane dipping at an angle ϕ with a strike direction normal to the xz plane are given by:

Using Coulomb failure assumptions³⁰, the Coulomb failure stress is expressed as:

where μ is the coefficient of friction, p is pore-fluid pressure and S is cohesion. Assuming that p, S and μ are constant over time, the change in Coulomb stress is given by³⁹:

where Δτ_n is positive for tension^{30, 39}. Similarly, the Coulomb stressing rate is given by:

The first term on the right-hand side is valid for an isotropic failure plane. We consider the change in shear stress in the slip direction (or rake) and vary the sign of the change in shear stress accordingly. For the San Andreas Fault we consider a dip angle of 90° and the unclamping stress (or stress rate) is given simply by (or ). For the Coalinga blind thrust we use dip angles ranging between 15° and 30° (refs 40 and 41) and a negative change in shear stress. We also consider a 65° northeast-dipping fault with a positive slip direction, simulating stress conditions on the Nuñez thrust fault. For each fault we calculate the normal and shear stress components at depths of 5 km and 15 km, representative of the seismogenic depth in Central California.

Rate of unclamping of the San Andreas Fault

The rate of unclamping on a vertical fault (change in fault normal stress) can be calculated as³⁰:

where θ₁ and θ₂ are the angles from both edges of the loading strip to any point at depth, measured clockwise from the positive x direction. Taking an estimate of the unloading rate due to loss in water storage of 8.8 (± 1.3) × 10⁷ N m⁻¹ yr⁻¹, as well as average spatial parameters for the San Andreas Fault (load width of 2a = 60 km; distance to centre of the distributed line load of around 70 km; seismogenic depths of 5–15 km), we estimate a rate of unclamping of 0.07–0.18 kPa yr⁻¹ at the seismogenic depth along the San Andreas Fault due to the elastic rebound caused by unloading (Extended Data Fig. 4b and c).

For the Coalinga blind thrust, we calculate the normal (τ_n) and shear (τ_s) stress rates resolved on a range of fault planes and depths. The Coulomb failure stress rate (ignoring pore pressure), (ref. 30), is then calculated with μ = 0.7.

Seasonal vertical displacements and stress changes

The elastic model can be used to estimate total vertical displacements caused by variations in total water storage. The model considered is more broadly distributed over the Coast Range, the San Joaquin Valley and the western Sierra Nevada (width of 200 km; Fig. 2), corresponding to the approximate total width of the San Joaquin River drainage basin. Taking an average of 200 mm of change in total water storage, we estimate a load of ~4 × 10⁸ N m⁻¹, which produces annual vertical displacements of 3–8 mm within the load extent (Fig. 2). The associated stress variations on the San Andreas Fault are about 1 kPa and the Coulomb stress changes on the Coalinga thrust are 1.0–1.7 kPa (Extended Data Fig. 4d and e). Varying the width of the distributed line loads by ±20 km and the load centre by ±10 km has only a small impact (<0.5 kPa) on the resolved seasonal stress changes on the faults.

1. Williamson, A. K., Prudic, D. E. & Swain, L. A. Ground-water flow in the Central Valley, California. US Geol. Surv. Prof. Pap. 1401-D. (1989)
2. Faunt, C. C. (ed.) Groundwater availability of the Central Valley aquifer, California. US Geol. Surv. Prof. Pap. 1766. (2009)
3. Famiglietti, J. S. et al. Satellites measure recent rates of groundwater depletion in California's Central Valley. Geophys. Res. Lett. 38, L03403 (2011)
   • Article
4. Scanlon, B. R. et al. Groundwater depletion and sustainability of irrigation in the US High Plains and Central Valley. Proc. Natl Acad. Sci. USA 109, 9320–9325 (2012)
   • PubMed
   • Article
5. Poland, J. F., Lofgren, B. E., Ireland, R. L. & Pugh, R. G. Land subsidence in the San Joaquin Valley, California, as of 1972. US Geol. Surv. Prof. Pap. 437-H. (1975)
6. Christiansen, L. B., Hurwitz, S. & Ingebritsen, S. E. Annual modulation of seismicity along the San Andreas Fault near Parkfield, CA. Geophys. Res. Lett. 34, L04306 (2007)
   • Article
7. Fay, N. P., Bennett, R. A. & Hreinsdóttir, S. Contemporary vertical velocity of the central Basin and Range and uplift of the southern Sierra Nevada. Geophys. Res. Lett. 35, L20309 (2008)
   • Article
8. Bennett, R. A., Fay, N. P., Hreinsdóttir, S., Chase, C. & Zandt, G. Increasing long-wavelength relief across the southeastern flank of the Sierra Nevada, California. Earth Planet. Sci. Lett. 287, 255–264 (2009)
   • CAS
   • Article
9. Hammond, W. C., Blewitt, G., Li, Z., Plag, H. P. & Kreemer, C. Contemporary uplift of the Sierra Nevada, western United States, from GPS and InSAR measurements. Geology 40, 667–670 (2012)
   • Article
10. Saleeby, J., Saleeby, Z. & Le Pourhiet, L. Epeirogenic transients related to mantle lithosphere removal in the southern Sierra Nevada region, California: Part II. Implications of rock uplift and basin subsidence relations. Geosphere 9, 394–425 (2013)
    • Article
11. Holzer, T. L. Elastic expansion of the lithosphere caused by groundwater depletion. J. Geophys. Res. 84, 4689–4698 (1979)
    • Article
12. van Dam, T. et al. Crustal displacements due to continental water loading. Geophys. Res. Lett. 28, 651–654 (2001)
    • Article
13. Jiang, Y., Dixon, T. H. & Wdowinski, S. Accelerating uplift in the North Atlantic region as an indicator of ice loss. Nature Geosci. 3, 404–407 (2010)
    • CAS
    • Article
14. Nof, R. N. et al. Rising of the lowest place on Earth due to Dead Sea water-level drop: evidence from SAR interferometry and GPS. J. Geophys. Res. 117, B05412 (2012)
15. Fu, Y. & Freymueller, J. T. Seasonal and long-term vertical deformation in the Nepal Himalaya constrained by GPS and GRACE measurements. J. Geophys. Res. 117, B03407 (2012)
16. Heki, K. Snow load and seasonal variation of earthquake occurrence in Japan. Earth Planet. Sci. Lett. 207, 159–164 (2003)
    • CAS
    • Article
17. Luttrell, K., Sandwell, D., Smith-Konter, B., Bills, B. & Bock, Y. Modulation of the earthquake cycle at the southern San Andreas fault by lake loading. J. Geophys. Res. 112, B08411 (2007)
    • Article
18. Bettinelli, P. et al. Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology. Earth Planet. Sci. Lett. 266, 332–344 (2008)
    • CAS
    • Article
19. González, P. J., Tiampo, K. F., Palano, M., Cannavó, F. & Fernández, J. The 2011 Lorca earthquake slip distribution controlled by groundwater crustal unloading. Nature Geosci. 5, 821–825 (2012)
    • CAS
    • Article
20. Ben-Zion, Y. & Allam, A. A. Seasonal thermoelastic strain and postseismic effects in Parkfield borehole dilatometers. Earth Planet. Sci. Lett. 379, 120–126 (2013)
    • CAS
    • Article
21. Christensen, M. N. Late Cenozoic crustal movements in the Sierra Nevada of California. Geol. Soc. Am. Bull. 77, 163–182 (1966)
    • Article
22. Small, E. E. & Anderson, R. S. Geomorphically driven late Cenozoic rock uplift in the Sierra Nevada, California. Science 270, 277–281 (1995)
    • CAS
    • ISI
    • Article
23. Le Pourhiet, L. & Saleeby, J. Lithospheric convective instability could induce creep along part of the San Andreas fault. Geology 41, 999–1002 (2013)
    • Article
24. Blewitt, G., Kreemer, C., Hammond, W. C. & Goldfarb, J. M. Terrestrial reference frame NA12 for crustal deformation studies in North America. J. Geodyn. 72, 11–24 (2013)
    • Article
25. Altamimi, Z., Collilieux, X. & Métivier, L. ITRF2008: an improved solution of the international terrestrial reference frame. J. Geod. 85, 457–473 (2011)
    • Article
26. Schmidt, D. A. & Bürgmann, R. Time-dependent land uplift and subsidence in the Santa Clara valley, California, from a large interferometric synthetic aperture radar data set. J. Geophys. Res. 108, 2416 (2003)
    • Article
27. Argus, D., Fu, Y. & Landerer, F. Seasonal variation in total water storage in California inferred from GPS observations of vertical land motion. Geophys. Res. Lett. 41, 1971–1980 (2014)
    • Article
28. Dettinger, M. D. & Cayan, D. R. Large-scale atmospheric forcing of recent trends toward early snowmelt runoff in California. J. Clim. 8, 606–623 (1995)
    • Article
29. California Department of Water Resources. California water plan update 2005. Volume 1: Strategic Plan. Dep. Water Res. Bull. 160–05. (2005)
30. Jeager, J., Cook, N. & Zimmerman, R. Fundamentals of Rock Mechanics 4th edn (Blackwell, 2007)
31. Zumberge, J. F., Heflin, M. B., Jefferson, D. C., Watkins, M. M. & Webb, F. H. Precise point positioning for the efficient and robust analysis of GPS data from large networks. J. Geophys. Res. 102, 5005–5017 (1997)
    • Article
32. Boehm, J., Niell, A., Tregoning, P. & Schuh, H. Global mapping function (GMF): a new empirical mapping function based on numerical weather model data. Geophys. Res. Lett. 33, L07304 (2006)
    • Article
33. Bar-Sever, Y. E., Kroger, P. M. & Borjesson, J. A. Estimating horizontal gradients of tropospheric path delay with a single GPS receiver. J. Geophys. Res. 103, 5019–5035 (1998)
    • Article
34. Scherneck, H.-G. A parametrized solid earth tide model and ocean tide loading effects for global geodetic baseline measurements. Geophys. J. Int. 106, 677–694 (1991)
    • ISI
    • Article
35. Blewitt, G. Carrier phase ambiguity resolution for the Global Positioning System applied to geodetic baselines up to 2000 km. J. Geophys. Res. 94, 10187–10283 (1989)
    • Article
36. Bertiger, W. et al. Single receiver phase ambiguity resolution with GPS data. J. Geod. 84, 327–337 (2010)
    • Article
37. Altamimi, Z., Collilieux, X. & Métivier, L. ITRF2008: an improved solution of the international terrestrial reference frame. J. Geodesy 85, 457–473 (2011)
    • Article
38. Rebischung, P. et al. IGS08: The IGS realization of ITRF2008. GPS Solut. 16, 483–494 (2012)
    • Article
39. Harris, R. A. Stress triggers, stress shadows, and implications for seismic hazard. J. Geophys. Res. 103, 24347–24358 (1998)
    • ISI
    • Article
40. Stein, R. S. & Ekstrom, G. E. Seismicity and geometry of a 110-km long blind thrust fault. 2. Synthesis of the 1982-1985 California earthquake sequence. J. Geophys. Res. 97, 4865–4883 (1992)
    • Article
41. Toda, S. & Stein, R. S. Response of the San Andreas fault to the 1983 Coalinga-Nuñez earthquakes: an application of interaction-based probabilities. J. Geophys. Res. 107, 2126 (2002)
    • Article

Download references

Funding for this work comes from NSF EarthScope award number EAR-1252210 to G.B. and W.C.H. GPS data were collected using the EarthScope Plate Boundary Observatory, SCIGN, BARGEN, BARD, CORS and IGS networks. We are particularly grateful to UNAVCO for operating the vast majority of GPS stations used in this project. GPS data were processed using the GIPSY OASIS II software and data products from the Jet Propulsion Laboratory.