Sustained Groundwater Loss in California's Central Valley Exacerbated by Intense Drought Periods.

The accelerated rate of decline in groundwater levels across California's Central Valley results from overdrafting and low rates of natural recharge and is exacerbated by droughts. The lack of observations with an adequate spatiotemporal resolution to constrain the evolution of groundwater resources poses severe challenges to water management efforts. Here we present SAR interferometric measurements of high-resolution vertical land motion across the valley, revealing multiscale patterns of aquifer hydrogeological properties and groundwater storage change. Investigating the depletion and degradation of the aquifer-system during 2007-2010, when the entire valley experienced a severe drought, we find that ~2% of total aquifer-system storage was permanently lost, owing to irreversible compaction of the system. Over this period, the seasonal groundwater storage change amplitude of 10.11 ± 2.5 km3 modulates a long-term groundwater storage decline of 21.32 ± 7.2 km3. Estimates for subbasins show more complex patterns, most likely associated with local hydrogeology, recharge, demand, and underground flow. Presented measurements of aquifer-system compaction provide a more complete understanding of groundwater dynamics and can potentially be used to improve water security.

Introduction

Groundwater resources comprise about 20% of the global freshwater supply, and their availability and sustainable use are vital to the basic needs of more than 1.5 billion people in the world (Alley et al., 2002; Famiglietti et al., 2011; Wada et al., 2010). Groundwater delivers a worldwide average of 42% of the water needed for irrigation (Döll et al., 2012), and it benefits aquatic ecosystems (Taylor et al., 2013). In recent decades, the widespread decline of groundwater levels is associated with higher extraction rates (Russo & Lall, 2017) and periods of severe drought (Diffenbaugh et al., 2015). Thus, knowledge of large‐scale distribution of aquifer properties and the impact of droughts on storage capacity of a stressed aquifer‐system is of great importance for water management and for forecasting the availability of water resources during future drought periods (Farr et al., 2015; Faunt et al., 2016; Miller et al., 2017). Furthermore, groundwater overextraction causes land subsidence (Galloway & Burbey, 2011), which is detectable from ground and space (Amelung et al., 1999; Galloway et al., 1998; Miller & Shirzaei, 2015). Land subsidence not only poses severe threats to infrastructure but also is a potential indicator of permanent compaction of an aquifer‐system (Poland & Ireland, 1988), which may imply a reduction in aquifer‐system storage capacity and groundwater stocks. In situ measurements of groundwater levels at sparse observational wells together with low‐resolution estimates of regional groundwater storage variation and loss through analysis of Earth's gravity field, obtained from the Gravity Recovery And Climate Experiment (GRACE) satellites, are widely used to characterize groundwater budgets (Famiglietti et al., 2011; Richey et al., 2015; Scanlon, Longuevergne, et al., 2012a). These measurements, however, lack the spatial resolution needed for management plans. Thus, the extent to which aquifers respond to overextraction through compaction is generally unclear, owing to a lack of dense and large‐scale monitoring networks. Moreover, after a lifetime of more than 15 years, the GRACE mission ended in October 2017. Although its follow‐on mission GRACE‐FO is due to lunch in spring 2018 (Flechtner et al., 2014), a signal gap of at least 6 months is unavoidable. Thus, filling the data gaps and validating large‐scale estimates of the groundwater budget is of great importance.

Central Valley in California, with an area of ~52,000 km2, includes one of the world's largest and yet most stressed aquifer‐systems (Famiglietti, 2014; Richey et al., 2015). The cultivation of more than 250 different types of crops in the Central Valley claims ~75% of California's groundwater supply (Lassiter, 2015; Maupin et al., 2014). The remaining 25% of groundwater sustains wetlands and supplies urban needs (Lassiter, 2015). The significant demand for groundwater, accelerated by population growth (US‐Census‐Bureau, 2014) and extreme droughts, has resulted in depletion of groundwater resources since decades (Scanlon, Faunt, et al., 2012b). The occurrence of extreme and prolonged drought events have not only impacted the ecosystems but also severely affected the agricultural outcome of the region (Howitt et al., 2014). Dealing with a rising water crisis, the California Sustainable Groundwater Management Act passed in 2014. It is designed to steer the state toward sustainable groundwater use by 2042, targeting a balance among critical societal values including public health and safety, healthy economy, ecosystem vitality, and opportunities for enriching experiences (CNRA, 2016). Land subsidence and aquifer compaction due to groundwater depletion affect the sustainability of resources and undermine these societal values. Knowledge of aquifer‐system response to drought and overdraft will provide valuable information to water managers and policy makers for planning actions. While the occurrence of drought or the natural impact of drought on the groundwater resources is out of direct human control, their negative impact on the aquifer‐system can be minimized by guiding groundwater use in a way that prevents permanent compaction of underground storage capacity and thereby sustaining groundwater availability for future droughts and generations to come.

In this study, we use interferometric analysis of large sets of synthetic aperture radar (SAR) images acquired by Advanced Land Observing Satellite (ALOS) L‐band satellite in combination with observations of Global Positioning System (GPS) and groundwater level data to constrain the valley‐wide high‐resolution vertical land motion associated with aquifer‐system depletion and degradation as well as their mechanical properties. The notion of “high resolution” here pertains to the dimension of ground resolution cells at which actual observations of land motion or gravity are obtained, which is ~100 m for InSAR versus ~50 km for GPS and ~300 km for GRACE.

Background

The Central Valley's aquifer‐system comprises layers of confined and semiconfined sedimentary units (Belitz & Heimes, 1990), and drillers' logs indicate lateral and vertical variations in sediment texture (Faunt, 2009). Based on the geomorphological distribution of the Corcoran clay, gravel, silt, and sand, the thickness of deposits holding groundwater is estimated to vary from 300 to 400 m across the valley (Faunt, 2009). Nonetheless, similar patterns observed in head levels associated with the shallow and deep aquifers suggest that these units are hydraulically connected (Smith et al., 2017). Groundwater levels generally show a long‐term decline in the northern and southern regions of the valley (Figures 1b–1d). During 2007–2010 drought period, the valley experienced accelerated groundwater decline, triggering an escalation of permanent compaction of underground layers. However, the decline rate is fastest in the southern San Joaquin Valley (Figure 1d). Previous studies based on GRACE measurements estimated the total groundwater loss in the valley to be 6.4–10.4 km3/yr during the drought period of 2007–2010 (see Table S1 in the supporting information and Famiglietti et al., 2011; Scanlon, Longuevergne, et al., 2012a). Moreover, significant seasonal variations are associated with the annual recharge and discharge of the aquifer‐system. Models of elastic loading constrained by GPS measurements indicate a seasonal oscillation of ~8 km3 (Argus et al., 2014), consistent with GRACE‐based estimates of 9.2 km3 in a dry year (Scanlon, Longuevergne, et al., 2012a).

Figure 1

(a) One thousand six hundred four observational wells with observations during the study period of 24 December 2006 and 1 January 2010, illustrating overall hydraulic head level drop. Kern, Tule, Tulare Lake, and Kaweah subbasins in southern Central Valley show larger head level drops (>30 m), and similarly, the Colusa subbasin to the north shows maximum drops of ~15 m. The black polygons are the 39 subbasins named in Figure 4a, forming the Central Valley hydrological system. Examples of groundwater level measurements at selected wells in the (b) northern (ID‐388965N1219139W002), (c) central (ID‐390214N1214057W001), and (d) southern valley (ID‐359125N1194468W001). (e) Highlights the major distinct hydrological units of Central Valley consisting of Redding, Sacramento, and San Joaquin Valley.

Figure 4

(a) Interferometric synthetic aperture radar (InSAR) vertical velocity map over the ~52,000 km2 area of Central Valley. The 39 subbasins are also shown. (b) Bivariate plot showing the relationship between the rate of vertical land motion measured by InSAR and Global Positioning System. The standard deviation of differences is ~8 mm/yr. (c) Close‐up view of the vertical land motion over the northern part of Sacramento Valley.

Methods

This study based on a careful analysis of Interferometric SAR data sets in conjunction with observations of GPS stations and groundwater levels across the Central Valley, California. A large set of SAR images were acquired by the ALOS‐1 L‐band satellite, one of the worlds' largest Earth observing satellites launched in January 2006 by Japan Aerospace Exploration Agency. Earlier works describe the methods used for multitemporal InSAR analysis (Shirzaei, 2013), combination of InSAR and GPS to obtain vertical land motion data (Shirzaei & Bürgmann, 2018) and integration of land subsidence and groundwater level data to estimate mechanical properties of aquifer‐systems (Miller et al., 2017; Miller & Shirzaei, 2015). The chart in Figure 2 summarizes major steps of the data analysis presented in this study.

Figure 2

Schematic diagram of the overall processing chain combining interferometric synthetic aperture radar (InSAR), Global Positioning System (GPS), and groundwater levels to characterize mechanical properties of aquifer‐system in California's Central Valley. The red color highlights the InSAR processing step, where M is the number of interferograms generated, , and are the wrapped and unwrapped interferometric phase, respectively, and are displacements in LOS and vertical directions, and is the LOS velocity. The blue color represents the GPS processing step, where and show horizontal (easting and northing) and vertical velocities. The black color represents the process of characterizing mechanical properties of aquifer‐system, where is the groundwater changes; and are elastic and inelastic components of groundwater level change obtained from PCA analysis; and are the elastic and inelastic components of vertical surface deformation obtained through PCA analysis; and  are the elastic and inelastic skeletal storage coefficients; and , and represent the aquifer storage loss percentage, seasonal groundwater change, and cumulative groundwater storage loss, respectively.

Multitemporal InSAR Measurements of Land Subsidence

To measure surface deformation across the Central Valley, we applied an advanced multitemporal Wavelet‐Based InSAR algorithm (Shirzaei, 2013; Shirzaei & Bürgmann, 2013). We processed 420 L‐band SAR images acquired on ascending orbit tracks (heading angle ~240° and incidence angle ~34.5°) of the ALOS‐1 satellite during the period 24 December 2006 to 1 January 2010 (Figure 3 and Table S2). The SAR images are acquired in 36 frames (i.e., F690–F800) and six partially overlapping paths (i.e., P217–P222). We generated more than 1,600 interferograms with maximum perpendicular and temporal baselines of 2,000 m and 1,500 days, respectively. The pixel dimension of ALOS‐1 data at full resolution is ~12 m × ~4 m in range and azimuth direction, respectively. To enhance the interferometric phase signal‐to‐noise ratio we applied a multilook factor of 8 × 25 in range and azimuth direction, which yields a pixel of dimension ~100 m ×~100 m on the ground. The Shuttle Radar Topography Mission digital elevation model with 1‐arcsec (~30 m) spatial resolution was used to calculate and remove the effect of topographic phase and flat Earth (Franceschetti & Lanari, 1999). A wavelet‐based analysis was used to examine the complex interferometric phase noises and identify elite (i.e., less noisy) pixels (Shirzaei, 2013). To reduce the effect of a topographically correlated component of atmospheric phase delay and spatially uncorrelated digital elevation model error, a suite of wavelet‐based filters was applied (Shirzaei, 2013; Shirzaei & Bürgmann, 2012). Using a reweighted least squares estimation, the set of interferograms was inverted to generate deformation time series and velocities along the line‐of‐sight (LOS) direction. Displacement maps were initially obtained for each individual SAR frame. We then used the pixels within overlapping areas of adjacent frames as tie points and calculated a constant shift, accounting for the difference between the local reference frame of a given frame and that of the reference frame (i.e., F690). This procedure allowed aligning all frames within a given path. Once the frames belonging to a given path were aligned, we used the pixels within overlapping areas of adjacent paths to further align paths with respect to the reference path (i.e., P217). We found that applying a constant shift in all cases was enough to yield a seamless displacement signal across frames and paths.

Figure 3

The footprint of Advanced Land Observing Satellite 1 synthetic aperture radar (SAR) images (cyan rectangles) and Global Positioning System stations (triangles) within the Central Valley. The SAR data were acquired in ascending orbit tracks (heading angle = ~350°, incidence angle = ~34.5°).

To estimate the vertical components of land motion, we used a combination of InSAR LOS displacement and horizontal (easting and northing) component of displacements measured at the GPS stations (Shirzaei & Bürgmann, 2018). We used only GPS stations far from zones subject to rapid land subsidence. Note that GPS stations used here belong to Plate Boundary Observatory network and were essentially located far from zones affected by none‐tectonic signals. Using radar satellite unit vectors (Hanssen, 2001), the 2D components of the GPS displacement were projected onto the LOS direction, linearly interpolated at the location of elite pixels, and then subtracted from the measured LOS displacements. This operation yielded a LOS displacement that was solely due to vertical land motion, which can readily be back‐projected on the vertical direction using the satellite unit vectors. Figure 3 shows the spatial distribution of the interferometric data set, including six satellite paths. The baseline plots for each path are shown in Figure S1 in the supporting information.

Aquifer Mechanical Properties

Elastic and Inelastic Storage Coefficients

Specific storage of a confined aquifer, , is given by Jacob (1940): which relates the amount of water produced as groundwater level decreases and aquifer‐system compresses and water expands. is the density of water, is gravitational acceleration, is aquifer compressibility, is water compressibility, and is porosity. Here we assume water is incompressible. Assuming a constant overburden load, Δ, where Δ is the change in groundwater level (Poland & Davis, 1969), aquifer compressibility is given by where Δ is the compaction and is the initial thickness of the compacting layer (Jacob, 1940). Considering a thickness of 50 to 800 m and maximum subsidence rates of 25 and 5 cm/yr for the San Joaquin and Sacramento Valleys, the aquifer compressibility in the south and north are 3 × 10−9 Pa−1 − 5 × 10−9 Pa−1 and 5 × 10−9 Pa−1 − 8 × 10−9 Pa−1, respectively, which correspond with bulk moduli of 20–300 MPa and 15–200 MPa. Substituting equation (2) in 1 yields the dimensionless skeletal storage coefficient (:

Depending on whether effective stress is greater than preconsolidation stress ( threshold, the skeletal storage can be separated into elastic () and inelastic () skeletal storage coefficients (Hoffmann et al., 2003): where Δ and Δ are the elastic, seasonal components of the vertical displacement and water level, respectively. Also, Δ and Δ are the long‐term compaction and the water level change, respectively, and is the compaction time constant. Assuming porosity of 0.3 for a typical aquifer‐system and compressibility of ~10−10 Pa−1 for water, the assumption of incompressible water may cause an error of the order 10−6 m−1 in estimation of elastic skeletal specific storage coefficient, which is often negligible. To estimate the elastic and inelastic skeletal storage coefficients, we follow the procedure detailed in Miller et al. (2017). They used wavelet analysis to extract seasonal components of the signal, but here we applied a spectral principal component decomposition approach that uses trigonometric base functions to extract annual and semiannual components of the land subsidence and groundwater‐level time series. Since the subsidence and groundwater level times series are sampled irregularly, we implemented the so‐called least squares spectral analysis (Vaníček, 1971), which allows robust estimation of the spectrum of unevenly sampled data sets. Figure S3 shows the amplitude of long‐term and seasonal components of the vertical land motion averaged over subbasins.

Aquifer‐System Storage Loss

Permanent aquifer‐system storage loss indicates the volume of one‐time used water, mostly extracted from aquitards in the area of overdraft. The percentage of permanent aquifer‐system storage loss () is defined as the ratio between the volume of the aquifer‐system that is compacted permanently and the total aquifer volume and given by (Miller et al., 2017): where is the inelastic skeletal storage coefficient, is the total thickness of the fine‐grained materials undergoing permanent compaction, and Δ is the amount the groundwater level change below the preconsolidation level. To estimate Δ, preconsolidation head level was defined within the 10‐year period of 1996–2006. Since the thickness of fine‐grained materials varies across the entire valley and the accurate thickness of the compacted layer is poorly understood (Faunt et al., 2010; Smith et al., 2017), we considered a range of thickness including 50, 100, 200, and 400 m and accordingly estimated the aquifer‐system permanent storage loss (Figure 5c). The thinness range of 50–400 m, though locally variable, is in agreement with the general range of aquifer‐system thickness across the valley.

Figure 5

Mechanical properties of aquifer‐system across Central Valley, for the study period of 24 December 2006 to 1 January 2010. (a) Elastic skeletal storage coefficient averaged over subbasin areas. (b) Inelastic skeletal storage coefficient averaged over the 15 subbasins which are characterized by land subsidence >3 cm. (c) Aquifer storage loss percentile per subbasins for different scenarios of fine‐grained layer thicknesses ranging from 50 to 400 m. (d) Seasonal groundwater volume amplitude, summed over subbasin areas. (e) Cumulative groundwater loss over 3 years, summed over subbasin areas.

Seasonal Change and Permanent Loss of Groundwater Storage

To estimate the seasonal changes and permanent loss in groundwater storage, we use a one‐dimensional relation between infinitesimal volumetric strain () and effective stress (dσ) (Terzaghi et al., 1996): where is elastic modulus and is bulk volume. Considering effective stress , where  is the total stress, is the effective pore pressure, is the Biot‐Wills coefficient ( 0 <  , and an area of A for each InSAR pixel, the associated change in volume of water () is

Summing the estimated values of for subsiding pixels within each subbasin yields the corresponding groundwater seasonal change and permanent loss, depending on whether is the seasonal or long‐term volumetric strain, respectively (Figures 5d and 5e). To evaluate equation (7), we interpolated our deformation data on a grid of 200 × 200 m (i.e., A = 40,000 m2) using an algorithm with weights proportional to the inverse of distance and assumed . The ratio between fine and the coarse‐grained material is variable across the valley, and the south is more dominated by coarse‐grained layers, while in the north, fine‐grained layers are more abundant (Faunt, 2009). Thus, we consider an elastic modulus of 200 and 300 MPa for the Sacramento Valley and San Joaquin Valley, respectively, under fully saturated conditions (Wang, 2017) consistent with the upper limit of the above estimated bulk modulus.

Results and Discussion

Observation of Land Subsidence and Permanent Compaction

Figure 4a shows the rate of vertical land motion across the Central Valley, consisting of more than ~23 million elite pixels, with a spatial resolution of ~100 × 100 m on the ground. Figure S3a also shows rates averaged for subbasins. Although the map appears continuous, due to phase decorrelation processes such as land cover change associated with agriculture and precipitation, there were some gaps in each frame where interferometric phase was corrupted, and thus, associated pixels were discarded. To validate the InSAR‐derived vertical land motions, we used the vertical components of ~300 continuous GPS stations, continuously monitored during our study period (Figure 3a). Given that GPS stations and InSAR pixels are not necessarily collocated, we defined an area of 200 × 200 m surrounding each GPS station and averaged the estimated values of land motion at pixels located within this box. This averaging may further improve the signal‐to‐noise ratio. Figure 4b shows a bivariate plot visualizing the degree and pattern of the relation between InSAR and GPS vertical land motions. We found a standard deviation of ~8 mm/yr for the differences, indicating good agreement between the two data sets.

During the observation period, the southern part of the valley, near Tulare Lake subbasin (latitude of 36.14°) experiences the fastest subsidence rate of ~25 cm/yr (see also Figures S2 and 4a). In the northern Sacramento Valley (latitude of 39.62°) the velocity pattern shows a lower subsidence rate of ~5 cm/yr. Central sections of the valley exhibit a negligible subsidence rate. For the southern part of Central Valley, the San Joaquin Valley, the cross section that passes through the southeast part of Delta‐Mendota subbasin displays a maximum subsidence trend of 15–20 cm/yr. Similarly, the cross section across Tulare Lake subbasin (latitude of 36.14°), north of Tule, and south of Kaweah subbasins shows a displacement rate of 20–25 cm/yr (see Figure S2). The subsidence rate further in the south of Kern county drops to about 5 cm/yr. There is also a good correspondence between zones of rapid subsidence with those showing a rapid groundwater level drop (Figures 1a and S2). Although our search did not turn up any study of valley‐wide land subsidence using InSAR, our results are in a good agreement with earlier smaller‐scale studies in the San Joaquin Valley (Farr et al., 2015; Smith et al., 2017) and Sacramento Valley (Farr et al., 2015).

Aquifer Response to Drought and Overdrafting

Combination of vertical land motion and groundwater level change allows the estimation of elastic () and inelastic () skeletal storage coefficients at location of each well (see section 2.2 and Miller et al., 2017; Miller & Shirzaei, 2015). Note that although InSAR subsidence maps are of high resolution, the aquifer mechanical properties can only be estimated at sparse well locations. Thus, to improve the signal‐to‐noise ratio and obtain relevant values for each subbasin, we averaged the sparse measurements associated with wells within each subbasin. These subbasins are closely related to the water balance regions defined by California Department of Water Resources used for modeling efforts (Faunt et al., 2010). Thus, providing estimates of mechanical properties at subbasin scales is also useful for modeling purposes. We found that larger values of correlate with subbasins having larger percentages of coarse‐grained sediments in the middle and southernmost parts of the Central Valley (see Figure 5a and Table 1). Table 1 provides the values for different mechanical properties per subbasin shown in Figure 5. Moreover, Figure S4 provides a visual comparison between our estimates and those obtained by earlier works. The average estimated for the valley is 2.84 × 10, consistent with earlier estimates of 4.6 × 10−4 by Helm (1978), 1.1 × 10−4 by Riley (1969), and 1.3 × 10−4 by Poland et al. (1975) assuming ~100‐m thickness of sediments. The Delta‐Mendota subbasin has the largest elastic skeletal storage coefficient of 7.1 × 10 in the San Joaquin Valley, while in the southern‐most San Joaquin Valley, the Kern County subbasin has a high value of 5.56 × 10. The Kaweah subbasin has the lowest elastic skeletal storage coefficient of 1.69 × 10 in the San Joaquin Valley. In the Sacramento Valley, the largest and lowest estimates of are 3.99 × 10 and 4.2 × 10, which are obtained in the Sutter and Red Bluff subbasins, respectively. These results indicate that the aquifer‐system in the southern part of the Central Valley has a higher yield and produce more freshwater for a given decline in groundwater level without undergoing irreversible deformation. In the following, Figure 5b shows the average coefficient of inelastic skeletal storage coefficient () for the 15 subbasins (see also Table 1) that experience cumulative subsidence of >3 cm and a long‐term drop in groundwater levels, which are in the same range as suggested by an earlier works (e.g., Sneed, 2001). The average is 4.08 × 10 for the entire valley. We find in the San Joaquin Valley to be larger than that in the Sacramento Valley. The Kern, Tule, Tulare Lake, Kaweah, and Merced subbasins in the San Joaquin Valley have the highest values of 5.0 − 6.5 × 10. The Sutter subbasin has a large of 5.0 × 10 compared to the rest of the Sacramento Valley. Such valley‐wide estimates of allow the identification of the aquifers susceptible to permanent storage loss. Note that we were unable to reliably constrain the compaction time constants. However, the observed nearly simultaneous changes in groundwater level and nearby GPS vertical components suggest a time constant on order of days to months, consistent with pumping mostly from thinner aquitards and fine‐grained interbeds in the aquifer‐system.

Table 1

Numerical Estimates of Aquifer Mechanical Properties Shown in Figure 5

Subbasins	Valley	Elastic skeletal storage coefficient Ske (×10−4)	Inelastic skeletal storage coefficient Skv (×10−3)	Aquifer‐system storage loss (%)	Seasonal groundwater amplitude (km3)	Cum. groundwater loss (km3)
Kern County	San Joaquin Valley	5.56	65.20	0.23–1.84	3.06 ± 0.59	−9.59 ± 1.99
Tule	San Joaquin Valley	3.15	56.60	0.15–1.19	1.04 ± 0.57	−3.33 ± 1.87
Pleasant Valley	San Joaquin Valley	4.94	55.40	0.21–1.71	0.054 ± 015	−0.18 ± 0.08
Tulare Lake	San Joaquin Valley	3.20	51.00	0.19–1.58	0.50 ± 0.15	−1.43 ± 0.71
Kaweah	San Joaquin Valley	1.69	53.50	0.02–0.16	0.71 ± 0.30	−2.06 ± 0.98
Kings	San Joaquin Valley	2.09	50.00	0.08–0.70	0.45 ± 0.10	−1.23 ± 0.43
Merced	San Joaquin Valley	2.09	60.90	0.09–0.71	0.51 ± 0.12	−0.84 ± 0.32
Delta‐Mendota	San Joaquin Valley	7.10	50.00	0.005–0.04	0.25 ± 0.06	−0.36 ± 0.15
Modesto	San Joaquin Valley	3.83	10.00	0.01–0.13	0.05 ± 0.02	−0.04 ± 0.03
Suisun‐Fairfield Valley	Suisun‐Fairfield Valley	1.43	41.60	0.003–0.02	0.06 ± 0.04	−0.13 ± 0.06
Solano	Sacramento Valley	1.23	21.70	0.02–0.21	0.26 ± 0.04	−0.41 ± 0.12
Yolo	Sacramento Valley	0.82	20.00	0.02–0.15	0.11 ± 0.02	−0.26 ± 0.07
Sutter	Sacramento Valley	3.99	49.90	0.02–0.21	0.26 ± 0.04	−0.22 ± 0.05
Colusa	Sacramento Valley	1.17	16.90	0.05–0.45	0.41 ± 0.07	−0.90 ± 0.21
Red Bluff	Sacramento Valley	0.42	10.70	0.02–0.20	0.10 ± 0.02	−0.29 ± 0.08
					Total SGA = 10.11 ± 2.5	Total GWL = −21.32 ± 7.2

Note. Provided values are associated with 15 subbasins across the Central Valley, which are subsiding >3 cm during the observation period

Permanent aquifer‐system storage loss is a direct consequence of inelastic deformation and can be estimated knowing the inelastic skeletal storage coefficient (Figure 5b) and thickness of fine‐grained materials. The thickness of fine‐grained materials is variable across the Central Valley, and little is known from drillers' logs (Faunt, 2009). We therefore consider thicknesses of 50, 100, 200, and 400 m to estimate a range for storage loss percentile (see section 2.2.2). Given the linear relationship between storage loss percentile and the thickness of fine‐grained layers, this range provides a lower and upper bound on the expected aquifer‐system storage loss across the Central Valley (Figure 5c). The aquifer‐system in the southern part of the San Joaquin Valley, associated with faster depletion rates, demonstrates more substantial storage loss compared with the Sacramento Valley. The estimated value is higher in Kern County subbasin with ~0.23%–1.84% of aquifer capacity being lost during the study period, while a slightly lower percentage loss is measured in the adjacent subbasins like Tule and Tulare Lake (~0.15%–1.19% and ~0.19%–1.58%; Table 1). In the Central Valley, a total of ~0.16–1.3% of aquifer‐system storage volume has been lost permanently during the study period. However, hydrological models (Scanlon, Faunt, et al., 2012b) and InSAR studies in the San Joaquin Valley (Smith et al., 2017) suggest larger values of up to 8%. We speculate that the difference pertains to the inelastic storativity coefficient used. Here we estimated values that are relevant to this study period rather than using approximate values from earlier studies.

To estimate the seasonal groundwater volume oscillation, we extracted the amplitude of the annual vertical land motion (Figure S3b) and applied a first‐order 1‐D poroelastic calculation (see section 2 and Figure 5d). The estimated seasonal oscillation for the entire valley is 10.11 ± 2.5 km3, consistent with earlier studies by (Scanlon, Longuevergne, et al., 2012a) given the error ranges. The spatial pattern of the seasonal groundwater oscillation is, however, heterogeneous. The largest seasonal oscillation of 3.06 ± 0.59 km3 is observed in the southernmost San Joaquin Valley (in Kern County subbasin). The rest of the valley has a seasonal groundwater oscillation of less than 0.5 km3. Although the source of groundwater recharge in the south is agricultural runoff and is principally different from that in north driven by natural processes such as precipitation, the recharge volume in San Joaquin Valley is thought to be comparable to that of Sacramento Valley (Faunt, 2009). Therefore, the results shown in Figure 5d likely represent the annual recharge affected by groundwater and surface water flows. We next use the first‐order 1‐D poroelastic calculation to estimate the total volume of depleted groundwater from the measured long‐term land subsidence (Figure 5e). A total of 21.32 ± 7.2 km3 of groundwater storage loss has been estimated, amounting to a depletion rate of ~7 km3/yr for the period December 2006 to January 2010. However, total groundwater storage loss is more significant in the southern part of San Joaquin Valley, in particular in Kern, Tule, Tulare Lake, and Kaweah subbasins with an amount of about 16.41 ± 5.5 km3 of groundwater loss within the observation period. Subbasins in the Sacramento Valley have a lower groundwater loss of ~2 km3, where the maximum rate of depletion has been noticed in the Colusa subbasin in the Sacramento Valley, with ~0.3 km3/yr loss. The InSAR‐based estimate of total groundwater volume loss is in good agreement with that obtained from GRACE observations. Although our period is slightly different from that covered in earlier GRACE studies, it is still within the error range so that all results agree well (see Table S1 and Famiglietti et al., 2011; Scanlon, Longuevergne, et al., 2012a; Xiao et al., 2017). Figure S5 also shows the spatial distribution of the rate of groundwater storage loss estimated by GRACE and InSAR, in terms of equivalent water thickness. The GRACE estimates are obtained by applying a method similar to that of (Werth et al., 2017) and based on JPL mass concentration products oversampled on a regular grid of 0.5° resolution (Watkins et al., 2015) combining with soil moisture variations from the Noah model of the Global Land Data Assimilation System (Rodell & Beaudoing, 2016), snow storage variations from the Snow Data Assimilation System (NOHRSC, 2004), and surface water storage monitored by the California Department of Water Resources (CDWR, 2017). Note that although the total volumes obtained from the two data sets are comparable, their spatial distributions are significantly different, likely due to GRACE's poor resolution and leakage errors.

Conclusions

As a result of excess pumping rates (Russo & Lall, 2017) and low precipitation due to frequent and intensified droughts (Diffenbaugh et al., 2015), groundwater resources are stressed (Famiglietti, 2014). Accurate knowledge of aquifer response is limited to few sites. This study used novel radar remote sensing data and highlights that the Central Valley's aquifer‐system was gradually losing their storage capacity permanently. We presented the first valley‐wide observations of vertical land motion, and our unique observational evidence enabled us to address several outstanding questions regarding the sustainability of groundwater resources and response of the aquifer‐system to overdraft and drought periods. We estimated the extent to which the aquifer‐system storage capacity was diminished due to groundwater levels dropping below historic low levels, a key factor currently considered infrequently by regulators for development and management plans. Our observations further confirmed regional scale estimates of groundwater decline provided by GRACE data sets. Thus, as an alternative to GRACE estimates, our InSAR measurements provided detailed data on the groundwater storage conditions that were an essential prerequisite for contemporary water management. Given the unprecedented decline in groundwater in the U.S. and globally, the information provided here has tremendous value to water managers, policy‐makers, and the many communities living in water‐stressed regions.

Figure S1

Click here for additional data file.

Figure S2

Click here for additional data file.

Table S1

Click here for additional data file.