Groundwater Loss and Aquifer System Compaction in San Joaquin Valley During 2012-2015 Drought.

California's millennium drought of 2012-2015 severely impacted the Central Valley aquifer system and caused permanent loss of groundwater and aquifer storage capacity. To quantify these impacts within the southern San Joaquin Valley, we analyze various complementary measurements, including gravity changes from Gravity Recovery and Climate Experiment (GRACE) satellites; vertical land motion from Global Positioning System, interferometric synthetic aperture radar, and extensometer; and groundwater level records. The interferometric data set acquired by the Sentinel-1 satellite only spans the period January 2015 and October 2017, while the other data sets span the entire drought period. Using GRACE observations, we find an average groundwater loss of 6.1 ± 2.3 km3/year as a lower bound estimate for the San Joaquin Valley, amounting to a total volume of 24.2 ± 9.3 km3 lost during the period October 2011 to September 2015. This is consistent with the total volume of 29.25 ± 8.7 km3, estimated using only Global Positioning System deformation data. Our results highlight the advantage of using vertical land motion data to evaluate groundwater loss and thus fill the gaps between GRACE and GRACE-Follow-On missions and complement their estimates. We further determine that 0.4-3.25% of the aquifer system storage capacity is permanently lost during this drought period. Comparing groundwater level and vertical land motion data following September 2015, we determine an equilibration time of 0.5-1.5 years for groundwater levels within aquitard and aquifer units, during which residual compaction of aquitard and land subsidence continues beyond the drought period. We suggest that such studies can advance the knowledge of evolving groundwater resources, enabling managers and decision makers to better assess water demand and supply during and in-between drought periods.

Introduction

In the southwest United States, notably California, severe and frequent drought periods have become a norm that local communities try to cope with (Woodhouse et al., 2010). However, the recent drought during 2012–2015 is the most severe in a millennium due to reduced precipitations and record high temperatures; the latter is linked to anthropogenic warming (Diffenbaugh et al., 2015; Griffin & Anchukaitis, 2014). Unusually dry and warm weather persisted for about 4 years, and the supply of surface water to rivers, lakes, and reservoirs dramatically reduced (Xiao et al., 2017). These factors significantly impacted ecosystems relying on rivers and led to deficient lake levels and intensified wildfires (Diffenbaugh et al., 2015; Williams et al., 2015). In the year 2014, a drought state of emergency was declared in the state of California, which was triggered by a decline of 8.14 km3 in surface water available to agriculture, and this was compensated by increased groundwater pumping of 6.3 km3 (Howitt et al., 2014). As a result of this drought, the loss for California's economy in its vast agricultural sector is estimated to be ~$2.7 billion and 21,000 jobs (Howitt et al., 2014).

The higher extraction rates during periods of severe droughts cause the decline of groundwater levels (Alley et al., 1999; Faunt et al., 2009, 2015; Galloway & Riley, 1999; Russo & Lall, 2017; Wada et al., 2010). Measurements of gravity fields from the Gravity Recovery and Climate Experiment (GRACE) satellites allow the quantification of changes in the water stocks at a regional scale (Famiglietti, 2014). For the drought period, 2007–2009, Famiglietti et al. (2011) and Scanlon et al. (2012) estimated a groundwater loss of ~24 and ~31 km3 within the river sheds that encompass California's Central Valley aquifer (Figure 1a). During 2012–2016, which includes the period of recent drought, Xiao et al. (2017) found a groundwater loss of ~40 km3 as an agreement between estimates from GRACE observations and that obtained from a water balance model for the Central Valley. Investigating the correlation between various drought indicators and GRACE observations of total water storage (TWS) changes, including changes in groundwater, soil, snow, and reservoirs, Wang et al. (2016) found a longer delay for recovery of groundwater stocks during the recent drought.

Figure 1

(a) Overview of the study area within the State of California (bold black outline), the San Joaquin Valley (filled brown area) and the Gravity Recovery and Climate Experiment (GRACE) analysis region (bold red outline). (b) Footprint of the Sentinel‐1A/Sentinel‐1B synthetic aperture radar scene (path P137, white rectangle). Locations of wells providing records of groundwater level for the period January 2008 through May 2017 (blue squares). Locations of continuous Global Positioning System (GPS) stations (red triangles) and extensometers (orange squares). Light cyan polygons show major subbasins of the Central Valley. The background displays an optical image distributed together with Sentinel‐1 products.

Groundwater depletion is often accompanied by compaction of aquifer systems and land subsidence. A recent study using high‐resolution measurements of vertical land motion provided by interferometric synthetic aperture radar (InSAR) suggests a subsidence rate of ~300 mm/year in the southern Central Valley during the previous drought period of 2007–2009 (Ojha et al., 2018; Smith et al., 2017). Note that once the effective stress (normal stress minus pore fluid pressure) increases beyond a preconsolidation stress level, inelastic deformation occurs and aquifer system storage capacity is permanently lost (Burbey, 2001; Jacob, 1940; Terzaghi, 1925). Measurements of land subsidence in conjunction with groundwater level data can be used to quantify mechanical properties of aquifer system (Miller & Shirzaei, 2015) and to assess the success of water conservation plans (Miller et al., 2017). Recent works combine GRACE measurements with InSAR estimates of land subsidence to characterize groundwater dynamics in Brazil, Central Mexico, and California as well as to discuss their potentials and their limitations (Argus et al., 2017; Castellazzi, Martel, Galloway, et al., 2016; Castellazzi, Martel, Rivera, et al., 2016; Hu et al., 2017).

In California, drought and land subsidence are typically most severe in the southern part of the state. The San Joaquin Valley aquifer system covers 26,000 km2, about half that of the entire Central Valley. This region includes the counties of San Joaquin, Merced, Madera, Fresno, Kings, Tule, Tulare, Kaweah, and Kern, which are home to a population of ~4 million, as of 2009 (Howitt et al., 2014), and which contribute one fifth of the national agricultural food production. The arid to semiarid climate receives on average 127–406 mm of rainfall annually (Galloway & Riley, 1999). Most of the streamflow in the valley enters from the east side draining the western Sierra Nevada, where much of the precipitation occurs as snow during winter (Galloway & Riley, 1999; Knowles et al., 2006; Lundquist et al., 2008; Pandey et al., 1999). Precipitation and streamflow in the valley vary significantly from year to year, and these variations substantially affect changes in groundwater stocks in the region (Dettinger et al., 2004; Kattelmann, 1997; Pandey et al., 1999).

Quantifying the response of aquifer system to drought and understanding its link to land subsidence in San Joaquin Valley is of great importance to water managers, policymakers, and the scientific community. Here, we investigate meteorological, hydrological, and geodetic data sets to quantify impacts of the 2012–2015 drought on groundwater resources within San Joaquin Valley. We begin with analyzing precipitation records and drought indicators together with GRACE‐based estimates of changes in groundwater stock (sections 2.1, 2.2, and 3.1). Next, we investigate the compaction of the aquifer system during and after the drought period manifested in the measurements of vertical land motion obtained from Global Positioning System (GPS), extensometer, and InSAR observations (sections 2.4). In combination with measurements of groundwater levels across the valley (sections 2.3), we solve for mechanical properties of the aquifer system (sections 2.5, 3.2, and 3.3). Further, we apply a 1‐D poroelastic model to the GPS‐derived vertical displacements to determine changes in groundwater volume (sections 2.5 and 3.4).

Data and Methods

Precipitation and Drought Indicators

Three‐hourly Multi‐satellite Precipitation data at 0.25 × 0.25° grid resolution provided by the Tropical Rainfall Measuring Mission are used to investigate the spatiotemporal variation of precipitation in San Joaquin Valley (Tropical Rainfall Measuring Mission, 2011). The daily precipitation product (data identifier: 3B42) is aggregated to obtain a monthly estimate of precipitation. Then, monthly precipitation anomalies are calculated as the monthly estimate minus the monthly decadal mean. Figure 2a illustrates the time series of precipitation anomalies for the period 2004 to 2016 over the area covered by two tiles (southern and northern) as well as only the southern tile of the GRACE analysis region, defined in Figure 1a and Figure S1 in the supporting information.

Figure 2

(a) Monthly precipitation anomalies based on Tropical Rainfall Measuring Mission data in the GRACE analysis region shown in Figure 1a (solid black line) and its southern half (red dashed line). (b) Temporal variation of the drought state. Light orange to dark red colors indicates area fraction of the state of California characterized by USDM drought severity states D0 to D4 (see Table S1 for descriptions). Graphs show average SPEI drought index within the GRACE analysis region (solid black line) and its southern half (red dashed line). GRACE = Gravity Recovery and Climate Experiment; USDM = U.S. Drought Monitor; SPEI = Standardized Precipitation‐Evapotranspiration Index.

A drought severity classification for the state of California is provided by the U.S. Drought Monitor (USDM, 2018). It combines various metrics, including the Palmer drought severity index, standardized precipitation index, and measurements of soil moisture and streamflow states to describe the dryness and state of drought in a given area. The drought severity is classified into five categories from abnormally dry to exceptional drought (see Figure 2b and Table S1) and provided along with the fraction of the region covered by each class.

The Standardized Precipitation‐Evapotranspiration Index (SPEI; Vicente‐Serrano et al., 2010) is also presented as an alternative drought indicator that includes a temperature component. Record high temperatures, such as those of the 2012–2015 drought, intensify drought impact by increasing a region's water loss through evapotranspiration. SPEI intensity is calculated on a range from −3 to +3 from precipitation and temperature records on multiple time scales from 1 to 48 months. Negative values correspond to a drought period. For this study, the 12‐month SPEI is acquired until the end of 2015 from the global SPEI database on a 0.5° resolution (SPEI12, 2017). The SPEI index is aggregated for the entire GRACE analysis region (shown by a red polygon in Figures 1a and S1) as well as for the southern part of that region comprising the San Joaquin Valley.

GRACE Observations of Water Storage Changes

Measurements of changes in the Earth's gravity field obtained from GRACE satellites are used to estimate the long‐term loss of groundwater storage over California's Central Valley. To this end relevant storage variations in soil (), snow () and surface water () are averaged for GRACE analysis region and then their effect is removed from TWS variations () obtained by GRACE observations. Here we do not explicitly account for the possible groundwater loss in the crystalline basement in the Sierra Nevada, suggested by Argus et al. (2017). Accordingly, the annual component and a long‐term trend ( can be calculated for each time series. Therefore, the long‐term changes in groundwater ( ) for a region is given by the following:

All data sets are either acquired or interpolated on a grid of 0.5° resolutions at the monthly temporal resolution and then averaged for the analysis region. Soil moistures variations are retrieved from the latest Noah model of the Global Land Data Assimilation System Version 2 (GLDAS‐2, 0.25° resolution; Beaudoing & Rodell, 2016; Rodell et al., 2004). Snow storage variations are acquired from the Snow Data Assimilation System (National Operational Hydrologic Remote Sensing Center, 2004). Surface water storage is accumulated from storage changes in the largest reservoirs located within the GRACE analysis region. The reservoir records are retrieved from the California Department of Water Resources (California Department of Water Resources, 2017). However, our data set does not include the reservoir and lakes with a volume of less than 0.9 km3.

GRACE estimates of TWS variations are obtained from Jet Propulsion Laboratory (JPL) providing level‐3 data based on a mass concentrations approach with 3° spherical caps (mascons, version 2, coastal resolution filter applied; Watkins et al., 2015). This GRACE solution relies on geophysical constraints to suppress noise in the data, making a posteriori filtering unnecessary. Previous studies indicate that this GRACE data set exhibits less signal leakage and bias on a scale above 3°, when compared to GRACE solutions based on spherical harmonics, specifically, if the coastal resolution filter is applied to compensate for leakage from land to ocean (Scanlon et al., 2016). The long‐term changes of TWS are estimated as a spatially averaged signal within the study region (Figures 1a and S1).

To correct for the amplitude bias caused by applying an inherent smoothing function to the 3° × 3° mascons, JPL provides also a 0.5° gridded scaling factors, which are calculated using hydrological models (Figure S1; Wiese et al., 2016). Rescaling of the mascon data adjusts for the signal variation within the original 3° mascon caps. These factors were not specifically designed for refining the long‐term trend and are only suitable for adjusting the TWS time series (Watkins et al., 2015). The location of the Central Valleys aquifer (~52,000 km2) does not match placements of the individual spherical mass caps because its area falls well below the area of either a spherical mass cap (~110,000 km2) or the usually referenced GRACE footprint of ~150,000–200,000 km2 (Longuevergne et al., 2010; Rowlands et al., 2005). However, there is no other aquifer with a similar depletion rate in the vicinity of Central Valley and the effects of other hydrological compartments can be calculated and removed. Therefore, for calculating in equation (1) we do not cut through the original mascon tiles and perform our analysis on the larger region to capture the potentially missing groundwater signal distributed inside each mascon tile. We do not apply any scaling factors because we process GRACE data within two complete 3° JPL mascon tiles intersecting with the Central Valley aquifer system (see Figures 1a and S1). However, we consider a 10% scaling uncertainty () to account for the error in the signal amplitude due to signal leakage across mascon tiles. Note that the estimate of TWS and associated errors are correlated in space due to the limited spatial resolution of the GRACE data. Signal leakage across mascon tiles is minimized but not zero, because signal within a tile is not completely independent from its surroundings (Watkins et al., 2015). Finally, to estimate the uncertainty of the trend of basin‐averaged groundwater change ( ), we begin with estimating uncertainty ( ) of TWS trend () through bootstrapping (Werth et al., 2017). We then replace in equation (1) by the GRACE trend, introduce the 10% leakage error, and apply the concept of error propagation (Mikhail & Ackermann, 1976):

The application of the bootstrapping and the introduction of leakage error across mascon tiles provide a more conservative error estimate for GRACE‐based TWS than that of Wiese et al. (2016).

During the analysis, we accounted for data gaps in the GRACE time series. A conservative uncertainty of 30% is also applied for trends in snow, surface, and soil water storage because independent measurements to validate snow and reservoir data are not available. Moreover, the soil moisture estimates vary among different models, and the associated uncertainties are generally unknown.

Due to the limited spatial resolution of GRACE observations, the signal in the southern tile correlates to that of northern tile. Thus it is not possible to retrieve a sufficiently accurate estimate for the San Joaquin Valley. Although, we mainly discuss the GRACE‐based estimate of groundwater change for the entire Central Valley. We also inspect individual estimates for the southern and northern mascon tiles covering the study region (approximately covering the south and north valley, respectively). These estimates provide upper and lower limits on the expected water mass losses in the southern region (see discussion in section 3.1). In the following, all basin‐averaged time series and estimated trends for water budget components are referred to the GRACE analysis region, unless otherwise is stated. As seen, all components of water storage experience significant loss during the 2012–2015 drought period (Figures 3a and 3b and Tables 1 and 2).

Figure 3

Temporal variations in components of the water cycle within the Gravity Recovery and Climate Experiment (GRACE) analysis region (Figure 1a). (a) Changes in water storage in soil (brown, Global Land Data Assimilation System Version 2, GLDAS‐2), surface reservoirs (green, California Department of Water Resources) and snowpack (yellow, Snow Data Assimilation System), and total water storage (TWS) variations from GRACE (black) after applying a 12‐month moving average window (dashed black). (b) Groundwater storage variations without (solid lines) and with (dashed lines) a 12‐month moving average window applied and calculated for the GRACE analysis region (black), its northern (green), and southern tiles (red) separately. Gray areas in the background indicate periods characterized by one of following U.S. Drought Monitor (USDM) drought states (see also Figure 2): Light gray refers to more than 30% of California's area being in moderate drought, medium, gray means more than 30% being in exceptional drought and dark gray refers to more than 50% being in exceptional drought. These simplified drought periods are very similar in the San Joaquin County, except that the two latter classes are merged (medium and dark gray indicate D4 > 50% in San Joaquin).

Table 1

Estimates of Change in Various Components of the Water Budget as well as of TWS for the Period October 2011 to September 2015

Storage compartment	mm w. eq. ann. loss (mm/year)	Vol. ann. loss (km3/year)	Volume trend (km3)
Total	−87.1 ± 12.3	−18.9 ± 2.7	−75.5 ± 10.7
Soil water	−7.9 ± 2.4	−1.7 ± 0.5	−6.9 ± 2.1
Snow water	−1.4 ± 0.4	−0.3 ± 0.1	−1.2 ± 0.4
Reservoir	−13.2 ± 3.9	−2.9 ± 0.9	−11.4 ± 3.4
Groundwater (GRACE analysis region)	−64.6 ± 13.2	−14.0 ± 2.9	−56.0 ± 11.4
Groundwater (GRACE southern tile)	−64.3 ± 24.8	−6.1 ± 2.3	−24.2 ± 9.3

Note. Total water storage is from Gravity Recovery and Climate Experiment (GRACE), soil water storage is from model simulations (GLDAS‐2), snow water storage is from remote sensing products (SNODAS), and surface water storage in the form of reservoirs is from terrestrial monitoring products (see section 2.2). Estimates for groundwater storage are provided for the GRACE study region and its southern tile (Figure 1a). The second column shows millimeter water equivalent for annual loss in the water budget, the third column presents annual volume loss of each component, and the fourth column contains the trends of loss for each component. Also, 1‐sigma standard errors are shown for each estimate.

Table 2

Comparing GRACE‐Based Estimates of Groundwater Changes Within the GRACE Study Region (Shown in Figure 1a) for Various Periods

Period	Duration (year)	G.W. Sea. Var. (km3)	G.W. Vol. loss (km3)
January 2004 to June 2017	13.5	9.0	−52.5 ± 11.4
January 2007 to December 2009	3	n/a	−27.5 ± 5.5
January 2012 to December 2016	5	n/a	−43.0 ± 12.5
October 2011 to September 2015	4	n/a	−56.0 ± 11.4
February 2015 to June 2017	2.4	n/a	21.3 ± 12.0

Note. The periods include the following; January 2004 to June 2017: almost entire duration of Gravity Recovery and Climate Experiment (GRACE) observations, January 2007 to December 2009: a period of previous severe drought, January 2012 to December 2016: a period of recent severe drought and October 2011 to September 2015: a period of fastest depletion rate during the recent drought. Note that no significant seasonal signal is present in the GRACE‐based estimate of groundwater during very dry months. The first column shows the investigated period, the second column provides the duration that is examined, the third column contains the estimated volume of seasonal groundwater storage change, and the fourth column is the estimated total groundwater loss as well as 1‐sigma standard error.

Groundwater Level Data

From the database of more than 500 groundwater pumping wells across the San Joaquin Valley (Department of Water Resources, 2017), we select 175 wells with continuous measurements during our observation period (Figure 1b). Following Ojha et al. (2018) and Smith et al. (2017), we assume shallow and deep aquifers are hydraulically connected, and vertically separated units are part of one larger aquifer system. Figure 4a shows the spatial distribution of wells across the San Joaquin Valley color coded with the measured total groundwater level drop during October 2011 to October 2015. The maximum groundwater level drops are measured within the western and central part of the valley. Tulare Lake, Tule, Kaweah, Kings, and Westside subbasins experience most substantial groundwater level changes of >30 m. Examples of the time series of groundwater level change are shown in Figures 4b and 5. Note that between 2016 and 2017, despite the heavy precipitation near the end of 2015, groundwater levels remain stable and only show seasonal variations.

Figure 4

(a) Spatial distribution of wells used for monitoring groundwater levels color coded by the cumulative groundwater level drop during the period October 2011 to October 2015. Black polygons show major subbasins. Arrows indicate the approximate location of three Global Positioning System (GPS) stations, whose vertical displacement time series is shown in panels (b)–(d) by red lines. (b–d) Examples of time series of groundwater level change (blue lines) measured at the wells adjacent to the GPS stations.

Figure 5

(a–e) Time series of compaction measured at locations of extensometers (orange squares; Figure 1b) versus time series of groundwater level change (blue squares) at the nearby well.

GPS, Extensometer, and InSAR Data

To measure land subsidence associated with the compaction of the San Joaquin Valley aquifer system, we investigate data from extensometers and GPS stations and synthetic aperture radar (SAR) interferometry (Figure 1). Although sparsely distributed, GPS and extensometers provide temporally dense measurements for the entire drought period of 2012–2015 and beyond. The publicly available SAR data, acquired by Sentinel‐1A/Sentinel‐1B satellites over the San Joaquin Valley, span the period 31 January 2015 to 29 October 2017. Thus, although InSAR measurements provide high‐resolution estimates of land subsidence, they cover the last year of the drought period and the following ~2 years. Therefore, a comprehensive understanding of land subsidence in space and time is obtained through joint analysis of these data sets. The continuous observation of land subsidence at GPS and extensometer stations in combination with measurements of the groundwater level allows us to determine aquifer mechanical properties and estimate groundwater loss during the drought period. While the spatially dense measurements of subsidence during the postdrought period obtained by SAR interferometry enables to characterize the response of aquitard layers to the groundwater overdraft from the aquifer unit.

The 3‐D displacement field is measured at the location of 19 GPS stations distributed across the San Joaquin Valley (see Figure 1b). Similarly, Argus et al. (2017, Figure S2 therein) provide a displacement time series for 16 GPS sites within the Central Valley. The Nevada Geodetic Laboratory provides the daily solutions of GPS locations in the North American fixed reference frame 2012 (Blewitt et al., 2013). Figures 4b, 7a–7c, and S2 show examples of GPS vertical displacement time series. For the period 2012–2016 up to 1,000 mm of land subsidence is measured at some stations.

Figure 7

(a–e) Comparison of vertical land motion measured by GPS (red) and InSAR (black triangle) at various stations. The InSAR pixels within 200‐m radius distance of each station are averaged. Also, GPS station P544 (see Figure 1b) is used as the reference (i.e., zero displacement) for both InSAR and GPS measurements. See Figure S2 for time series of the remaining GPS stations. (f–j) Comparison of compaction measured by extensometers (orange) and InSAR vertical land motion (black triangle). STD is the standard deviation of the difference between InSAR and GPS or extensometer. InSAR = interferometric synthetic aperture radar; GPS = Global Positioning System.

In Figures 4b–4d, we also provide a comparison between the time series of GPS vertical displacement and groundwater level measured at the nearby well. There is a good overall correspondence between two data sets, as they both simultaneously show rise and fall superimposed on a long‐term declining trend. Exceptions appear following 2016 when groundwater level become stable and only show seasonal variations. We notice a lag of ~10–20 months between recovery of groundwater levels and land subsidence reversal. This observation provides an estimate of equilibration delay between fine grain aquitard head levels to neighboring aquifer levels.

We also acquire data from five extensometers, whose locations are shown in Figure 1b. An extensometer measures the vertical compaction and expansion of the aquifer system within certain depth. The extensometer depths range between ~100 and ~400 m, and they are located within the Westside subbasin, where depletion rate is slower. The associated time series with up to 200‐mm compaction are shown in Figures 5 and 7. Figure 5 also provides a comparison with the time series of groundwater level change at the nearby well. Similar to above, there is an overall correspondence between the temporal behavior of both data sets, except following 2016 when groundwater levels cease to decline, while extensometers continue to measure compaction.

We further use 68 SAR images obtained from the ascending track path 137 (heading angle ~347° and incidence angle ~33–43°; Figure 1b). Using a spatial and temporal baseline of 200 m and 100 days, we generated 228 interferograms. Applying a multilook factor of 30 and 5 pixels in range and azimuth direction, respectively, we obtain an average spatial resolution of 70 m × 70 m on the ground. A Shuttle Radar Topography Mission Digital Elevation Model of 1‐arc sec (~30 m) spatial resolution and precise satellite orbital information were used to calculate and remove the effect of topographic phase and flat earth correction (Franceschetti & Lanari, 1999). We then apply an advanced multitemporal SAR interferometric algorithm to reduce impacts of environmental artifacts and to constrain a temporally variable deformation field (Shirzaei, 2013; Shirzaei et al., 2017). To further reduce the effect of the topography correlated component of atmospheric phase delay and the spatially uncorrelated Digital Elevation Model error, a suit of wavelet‐based filters is applied (Shirzaei, 2013; Shirzaei & Bürgmann, 2012). Using a reweighted least squares estimation, we invert the set of corrected interferograms to generate deformation time series and velocities along the line‐of‐sight direction. Using the horizontal component of GPS stations, we estimate and remove horizontal tectonic signals due to plate movements (Shirzaei & Bürgmann, 2018). The remaining signal is then projected onto the vertical direction using the satellite unit vectors (Hanssen, 2001). Figure 6a shows the spatial distribution of long‐term subsidence velocity across the San Joaquin Valley during January 2015 to April 2017. The southern San Joaquin Valley is characterized by a faster subsidence rate than the northern part. In particular, the northern Tulare lake, northwest of Tule, southwest Kaweah, and central Chowchilla subbasins are subsiding rapidly at a rate of >150 mm/year. The estimated rate for some parts of Tule and Kaweah subbasins is as high as 250 mm/year. Further north, the Chowchilla subbasin is also subsiding at a rate of ~150 mm/year.

Figure 6

(a) Rate of vertical land motion across the San Joaquin Valley during January 2015 to April 2017 measured by multitemporal processing of Sentinel‐1A/Sentinel‐1B synthetic aperture radar acquisitions. Negative values mean subsidence. (b) Deformation time series at five locations with fastest rates, whose position is marked in panel (a).

Figure 6b shows the time series of subsidence at five selected sites with high rates. We find that an even faster subsidence rate of >300 mm/year occurs between January 2015 and September 2016. Then the subsidence nearly ceases, and moderate uplift of about 10–20 mm occurs at some locations from September 2016 to October 2017.

Figures 7 and S2 compare InSAR vertical displacement time series against that measured by GPS stations and extensometers. We improve the signal‐to‐noise ratio by averaging the value of InSAR pixels within 200 m of each GPS and extensometer stations. In contrast to GPS observations, InSAR measurements are typically provided in a local reference frame. Namely, they refer to a point with an assumed zero displacement within the SAR frame. Such observations are suitable for investigating local deformation due to compaction of aquifers, as is performed here. In contrast, elastic load modeling exercises that solve for variation in TWS (Argus et al., 2017; Johnson et al., 2017) require measurements of surface deformation within a global reference frame (i.e., with respect to the Earth's center of mass). Here, for the sake of comparison, we refer both InSAR and GPS vertical displacements to the GPS station P544. Thus, this station acts as a reference with zero displacements for both data sets. Figures 7a–7e compare InSAR and GPS time series obtained at the location of each GPS stations after removing that measured at site P544. We find a good agreement between InSAR measurements and the other data sets. The average standard deviation of the difference between InSAR and GPS measurements is 21.9 mm, and between InSAR and extensometer measurements it is 14.5 mm. The agreement with extensometer measurements indicates that the compaction below ~400 m is negligible.

Mechanical Properties of Aquifer System

Estimating aquifer properties and characterizing the impact of drought on the groundwater storage capacity during the period October 2011 to September 2015, we jointly analyze groundwater levels and vertical land motion data measured at the location of GPS stations and extensometers data following Miller et al. (2017). The skeletal storage coefficient () can be expressed as a combination of elastic () and inelastic () storage coefficients as (Hoffmann et al., 2003) where and are the seasonal elastic components of vertical land motion and groundwater level, respectively, while and are the long‐term irreversible compaction and water level change. The is a threshold of preconsolidation stress level, and is the compaction time constant. To extract seasonal elastic components of the signal the least squares spectral analysis is implemented (Vaníček, 1971). To estimate the aquifer system storage loss (ASL) in percent, we consider inelastic storativity for the thickness of fine grain materials () undergoing permanent compaction: where is groundwater level change below the preconsolidation level. Here, the preconsolidation groundwater level is established as the lowest level of groundwater during the previous drought period 2007–2009. Since accurate estimates of are not available, we consider a range of 50–400 m for the thickness of fine grain sediments and evaluate ASL as a function of this thickness.

Given the one‐dimensional poroelastic relation between infinitesimal volumetric strain and effective stress presented by Terzaghi et al. (1996), the loss of groundwater storage () can be expressed as follows: where E is elastic modulus, is the effective stress (total stress minus pore fluid pressure), is Biot‐Wills coefficient (0 <  < 1) here equal to 0.9, is the area subject to subsidence, is a volumetric strain, and and are the density of water and acceleration due to gravity, respectively. To estimate areas subject to subsidence, the GPS vertical components are interpolated over a grid of 200‐m × 200‐m resolution using an inverse distance algorithm. We consider an elastic modulus of 300 MPa, assuming coarse grain materials dominate the shallow aquifer system in San Joaquin Valley (Faunt et al., 2009).

The distribution of GPS stations (Figure 1) compared with the area of extreme subsidence, detected by InSAR (Figure 6) suggests that GPS observations may not capture the fastest subsidence rates, although they cover the entire drought period. On the other hand, the publicly available SAR data sets from Sentinel‐1A/Sentinel‐1B satellites only cover the last year of the drought period. Thus, reliance on GPS observations is unavoidable to solve for the change in groundwater stock during the period October 2011 to September 2015. This also means that the GPS‐derived rate of land subsidence for this period is likely underestimated. Therefore, the results obtained following equations (6) and (7) provide a lower bound on the estimated volume of lost groundwater and on the percentile of permanently reduced aquifer storage capacity.

Moreover, the sparse distribution of GPS stations does not provide sufficient resolution to solve for the spatial distribution of groundwater volume change. Thus, we only report a spatially aggregated estimates for the entire San Joaquin Valley. However, benefitting from continuous measurements of GPS stations, we can obtain a time series of monthly change and loss in groundwater across the valley. Figure 8 presents a comparison of this time series against that of previous studies, which is further discussed in section 3.4.

Figure 8

Time series of monthly groundwater loss estimated in this study using GPS (green squares) for the San Joaquin Valley as well as from GRACE data for the southern mascon tile (red line) and both tiles of the GRACE analysis region (black line; see Figure 1a). Bars indicate a 1‐ error range for the GPS estimates of monthly groundwater loss. Our estimates are compared to that of earlier studies summarized in Figure 1b of Argus et al. (2017) and digitized here for the sake of visualization. These studies include Xiao et al. (2017; blue line for Central Valley), Faunt et al. (2015; purple line for Central Valley), and Argus et al. (2017; purple dashed line for the combination of estimates by Xiao et al., 2017, and Faunt et al., 2015, and light gray line for a GRACE estimate of the entire Sacramento/San Joaquin River basins). GRACE = Gravity Recovery and Climate Experiment; GPS = Global Positioning System.

Discussion and Conclusions

Drought Indicators and Groundwater Loss in Central Valley

California's most severe drought on record begins around October 2011; at the end of 2011, the seasonal rainfall dropped below its average level indicated by precipitation anomalies (Figure 2a). According to the USDM, which incorporates water availability and vegetation health, almost the entire state of California is classified as under drought during the period 2012–2016 (Figure 2b). The drought reached its peak during mid‐2014 after precipitation dropped to its lowest. In the second half of 2014 more than half of the state is characterized by exceptional drought (Figure 2b). Meanwhile, the drought indicator that only accounts for TWS variations reaches its peak in September 2015 (Argus et al., 2017). At the beginning of 2016 precipitation returns to an above average level and most of the state is relieved from intense drought. Nonetheless, until the beginning of 2017, at least one third of the state remains under extremely dry conditions, and a larger part of it is still characterized as under severe drought status. The SPEI drought index predicts a very similar timing for the onset of the drought, but it determines worst drought conditions for the early 2014 and a slow recovery that begins near the end of 2014. Note that it is expected that various drought indicators provide different estimates for drought duration and intensity, owing to the data sets that are used for creating the indicator to serve a particular purpose or community. For example, an agricultural drought estimated by USDM indicator may persist longer than that of the meteorological or hydrological indicators. This is because the response of vegetation health to changes in water availability is neither immediate nor linear because several other factors are also involved. In the following, we focus on the hydrological drought, which reflects water availability on and below the surface.

The drought had a profound impact on water quantities stored in soil and reservoirs (Figure 3a). These resources promptly respond to lack of water inflow as well as increased losses and undergo a sharp decline in 2012. Snowpack was low during 2012 and 2013 and effectively zero during 2014 and 2015, consistent with the record low estimate of TWS obtained from GPS and GRACE (Argus et al., 2017). Such snowpack decline is not observed during the previous extreme drought of 2007–2009. Following a positive precipitation anomaly, the recovery process began during the wet season of 2015 for soil moisture and of 2016 for snowpack and reservoirs. Note that soil moisture trend is minimal during the 2012–2015 drought period (Table 1). But despite its early recovery, a substantial interannual response to drought is evident (Figure 3a).

The aggregated impact of drought conditions on components of the water cycle is reflected in variations of TWS obtained from gravity measurements by GRACE satellites. TWS decline begins at the end of 2011 and lasts through the end of 2015 (Figure 3b) with a total loss of ~76 km3 for the entire GRACE analysis region (e.g., two mascon tiles) shown in Figure 1a (see also Table 1).

Given the time series of soil moisture, snowpack, and reservoir volume change, monthly estimates of groundwater variations and its long‐term trend can be obtained (Figure 3b and Tables 1 and 2). Groundwater inferred from GRACE shows minimal seasonal changes during the drought periods, compared to 9 km3 during 2004–2017 (Table 2). This indicates that there has been very little net recharge to the aquifer system during the drought periods. Interestingly, the maximum rate of groundwater loss for the entire GRACE analysis region is 14 km3/year with a total of 56 km3, which is measured for the period October 2011 to September 2015 (see Tables 1 and 2). We consider this period as the duration of most intense groundwater drought, which is slightly different than that suggested by drought indicators for the period 2012–2016 (Table 2). The estimated groundwater loss for the period October 2011 to September 2015 is near twice the volume loss during the drought period of 2007–2009, being ~28 km3.

Due to leakage errors, the GRACE observations do not allow distinguishing between the signals from the southern and northern part of the Central Valley and some of the signal associated with the south mascon tile likely smears into the surrounding ones. GRACE‐based estimates indicate a water volume loss in the northern part comparable to that of the southern part of the study region (Table 2, rows 5 and 6). However, high‐resolution InSAR measurements of vertical land motion by Ojha et al. (2018) show that most of the aquifer compaction and groundwater depletion during the 2007–2009 drought occurs in the southern Central Valley. Therefore, it is reasonable to assume that the leakage bias in the GRACE signal for the southern tile (encompassing the San Joaquin Valley) is negative. This means the GRACE‐based amount of groundwater loss for the San Joaquin Valley is underestimated, while an amount obtained for the entire Central Valley is likely overestimating the loss in the San Joaquin Valley. Therefore, we consider the GRACE‐based estimate of ~24 km3 (~56 km3) for groundwater loss in the southern tile (south and northern tiles) as a lower (upper) bound on the amount of loss in the San Joaquin Valley (Table 2).

According to the GRACE observations, groundwater depletion flattens in 2015 when positive anomalies characterize precipitation, and less than 30% of the state is within the USDM classification of exceptional drought (Figure 3b). It, however, does not rise until almost a year later at the beginning of 2016, when precipitation anomalies reach ~20 mm above their average level. Some of the lost groundwater volume (~23 km3) is later recovered during the following period of 2015–2017. Given that GRACE records are discontinued during the recovery phase, it is unclear how much groundwater is replenished following mid‐2017.

The study by Argus et al. (2017) infers a substantial long‐term groundwater loss in the Sierra Nevada and Klamath Mountains. They also suggest that hydrological models underestimate the loss of groundwater in the mountains. Thus, any estimate of groundwater loss using GRACE measurements may include both contributions from depletion of Central Valley aquifers and groundwater loss in river alluvium and the crystalline basement beneath the Sierra Nevada. This conclusion, however, cannot be verified by using satellite gravimetric data, due to the limited spatial resolution of GRACE measurements. The valley‐wide study by Ojha et al. (2018), which explores large‐scale high‐resolution vertical land motion data through a 1‐D poroelastic model, estimates a groundwater loss of 21.3 ± 7.2 km3 for the entire valley from 2007 to 2009, comparable to the lost volume of 27.5 ± 5.5 km3 obtained from GRACE and hydrological data (see Table 2, row 2). This agreement may mean that contributions from other sources such as the groundwater loss in the Sierra Nevada and Klamath Mountains are negligible, given the model and data uncertainties, and thus cannot be resolved through our approach. However, we consider both results from Argus et al. (2017) and our study as complementary insights into hydrological changes in the entire California.

Investigating the temporal lag between the TWS and GW change and precipitation anomalies, we first remove the seasonal oscillations in the three data sets by applying a boxcar filter of 12‐month length. We then identify the lag at which the smoothed TWS and GW change have maximum cross correlation with smoothed precipitation anomalies. The delay was also calculated for smoothed TWS and GW change versus the 12‐month SPEI index. We find that for the entire GRACE observation period (January 2004 to June 2017) the precipitation anomaly preceded the TWS and GW change (as well as SPEI) by 4 and 9 months (as well as 3 and 12 months), respectively. Splitting the observation period into shorter periods of January 2004 to December 2010 and January 2011 to June 2017, we find that the time lags between TWS and GW change and precipitation anomaly (and SPEI) significantly increase following the recent drought. For the period before the recent drought, the precipitation anomaly preceded the TWS and groundwater change by 3 and 7 months (and 1 and 6 months for SPEI), respectively. While the corresponding values for the following period are 22 and 30 months (11 and 14 months for SPEI). These observations are consistent with the finding of Wang et al. (2016) and Lin et al. (2015) suggesting that following severe droughts and excess overdraft recovery of the aquifer system may be much slower than usual.

Aquifer System Mechanical Properties

To characterize the response of the aquifer system to the drought period of October 2011 to October 2015, we jointly analyze the measurements of vertical land motion obtained at locations of GPS and extensometer stations with groundwater levels at nearby wells. We use only those GPS stations that have a well with frequent measurements within a radius of 20 km distance. This criterion leaves us with 18 GPS stations that are adjacent to a well (see Table 3). Following the procedure discussed in section 2.5, we calculate the coefficients of elastic skeletal storage () and inelastic skeletal storage () for each stations (Table 3) and provide mean values for the San Joaquin Valley (Table 4).

Table 3

Estimated Mechanical Properties of Aquifer Systems During the Observation Period of October 2011 to September 2015, for Each GPS and Extensometer Station Shown in Figure 1b

	Elastic skeletal storage Ske (×10−4)	Inelastic skeletal storage Skv (×10−2)
GPS:
CAFR	9.93	3.51
CALB	0.59	5.37
CAM5	0.10	4.83
CHOW	5.66	1.89
CRCN	0.11	2.45
LEMA	0.11	0.77
MULN	0.11	1.14
P056	0.11	0.71
P303	0.11	1.87
P304	0.10	1.38
P307	0.10	1.14
P544	0.28	0.19
P545	1.32	0.33
P564	8.70	1.72
P566	0.63	0.76
P809	9.98	10.00
RBRU	0.10	4.93
TRAN	0.11	1.69
Extensometer:
12S/12E‐16H2	0.10	1.05
14S/13E‐11D6	0.11	0.51
13S/15E‐31J17	0.53	0.50
18S/16E‐33A1	0.10	0.56
20S/18E‐06D1	0.11	0.52

Table 4

Mean Estimate of Aquifer Mechanical Properties Estimated From Values Provided in Table 3, As Well As the Aquifer System Storage Loss Assuming a Range of 50–400 m for Aquifer Thickness and the Groundwater Storage Loss Estimated Using GPS Measurements for the Period of October 2011 to September 2015

Aquifer property	Mean
Elastic skeletal storage coefficient Ske	1.7 × 10−4
Inelastic skeletal storage coefficient Skv	2.00 × 10−2
Total aquifer storage loss (%)	0.4–3.25
Total groundwater storage loss (km3)	−29.25 ± 8.7

Estimates of elastic skeletal storage vary by an order of magnitude, but the range of 0.1 × 10–1 × 10 and the mean value of 1.7 × 10 are consistent with those of earlier works (Sneed, 2001; Poland et al., 1975; Helm, 1978; Riley, 1969; Ojha et al., 2018; Smith et al., 2017). For instance, assuming a 50‐m‐thick aquifer unit, Sneed (2001) determined a range of 3 × 10–13 × 10 for the specific elastic skeletal storage within the San Joaquin Valley. Poland et al. (1975), Helm (1978), and Riley (1969) reported values of specific elastic skeletal storage for Central Valley that, respectively, correspond with 0.75 × 10, 2.3 × 10 and 0.55 × 10, given 50‐m‐thick sediments. Our estimates for the inelastic skeletal storage coefficients () have a mean value of 2.00 × 10, and their order of magnitude corresponds to the range of 2.3 × 10–11.0 × 10 reported by Smith et al. (2017), assuming a 50‐m‐thick aquifer unit.

We note that our estimates are systematically near the lower end of these suggested ranges. Preisig et al. (2014) show that the storage coefficients are stress dependent, and thus, their values may decrease during severe droughts when the effective stress significantly increases. However, the stress‐dependent reduction seems to be only a few percent and may not entirely explain the tendency of our estimates toward the lower bound of the suggested range for the storage coefficients. The other possible explanations include the following: (1) GPS and extensometer stations are not collocated with wells and are often away from the zone of maximum subsidence rate, which are presumably active aquifer units, and (2) aquifer mechanical properties are a function of local geology, which is heterogeneous across the valley.

Residual Compaction

As seen in Figure 6b, the InSAR time series of the vertical land motion indicate that fast subsidence rates continue until September 2016. Compared with Figure 3b, where GRACE‐based lower and upper estimates of groundwater change in the San Joaquin Valley are shown, we find that groundwater recovery started sometime in 2015. A similar conclusion can be made by investigating groundwater level measures at various wells shown in Figures 4b and 5 (as discussed in section 2.4). Aquifer systems with thick aquitard lenses may undergo lagged equilibration as a result of past overdrafts (Miller et al., 2017). Therefore, we attribute this lag of 0.5–1.5 years to the delayed equilibration of hydraulic head level differences between fine grain aquitards and adjacent aquifers. Such a lag is equivalent to the compaction time constant () present in equation (5). The estimate is consistent with a delay of ~2 years suggested by Ojha et al. (2018) and Smith et al. (2017) for the drought period 2007–2009. The San Joaquin Valley aquifer system comprises many clay layers (Faunt et al., 2009), which compose the system's total thickness. Assuming a comparable thickness for vertically separated units, the estimated can be related to the thickness of a single equilibrated hydrologic unit (), given the hydraulic diffusivity () as following (Smith et al., 2017):

Smith et al. (2017) reports average values of in case of elastic and inelastic deformation of clays to be 7.60 × 10 and 1.32 × 10, respectively. Thus, is estimated to be in the range of 7.4–12.8 and 1.0–1.7 m, during elastic and inelastic clay deformation, respectively.

Aquifer Storage and Groundwater Loss

Considering the lowest level of groundwater during the previous drought period of 2007–2010 as the preconsolidation groundwater level for the 2012–2015 drought and using a range of 50–400 m for sediment thickness, we apply equation (6) to determine the amount of aquifer storage loss (Table 4). We find that during the period of October 2011 to October 2015, 0.4–3.25% of aquifer storage capacity is permanently lost within San Joaquin Valley. Ojha et al. (2018) suggest an aquifer storage loss of 0.5–2% during the 2007–2009 drought period. Thus, an accumulation of the effects from both drought periods results in a 0.6–5.25% loss of storage capacity. These repeated observations are alarming because more severe and frequent droughts are becoming a new norm (Woodhouse et al., 2010). The storage capacity reduction may increase the replenishment time for aquifers after each drought period. Such time may not be sufficiently available in the future, due to more frequent droughts. This can cause positive feedbacks, exacerbating degradation of the aquifer system and posing further challenges to water management efforts.

Considering equation (7) and using the long‐term component of vertical land motion obtained from GPS and extensometers, we estimate a cumulative groundwater loss of −29.25 ± 8.7 km3 for the period of October 2011 to September 2015 (Table 4). This estimate is consistent with the lower bound obtained from GRACE measurements for its southern tile covering the San Joaquin Valley (Table 2). In Figure 8 we provide monthly estimates of groundwater loss based on GPS and GRACE data, as obtained here. We also compare these results with those from previous studies relying on a water balance model (Xiao et al., 2017), the Central Valley Hydrologic Model (Faunt et al., 2015), a combination of both (composite model of Argus et al., 2017), and an alternative GRACE estimate for the combined basins of the Sacramento and San Joaquin rivers (Argus et al., 2017). Our GPS‐based time series of groundwater volume change for the San Joaquin Valley indicates a fast rate of depletion during the drought period. Our estimate agrees well with those of other studies, given the range of errors and uncertainties. Noteworthy, the trend of our estimate of groundwater depletion is very similar to those obtained from hydrological models (Faunt et al., 2015; Xiao et al., 2017).

Interestingly, during the winter seasons, when pumping is expected to be the lowest, the GPS‐based time series shows a slower depletion rate, which is also evident in the results from the water balance study (Xiao et al., 2017). Although the GRACE‐based time series is affected by more significant errors, the associated long‐term trend is comparable with that of the GPS‐based estimate. The long‐term trend for estimates of groundwater in the southern GRACE mascon tile is slightly lower than our estimate (Table 1, row 6 vs. Table 4, row 4), but given the range of errors, the time series (red curve in Figure 8) agrees well with our GPS‐based estimate. As we discussed above (sections 2.5 and 3.1), both amounts for groundwater loss obtained based on observations of GPS and GRACE southern tile are considered the lower bound estimates for the San Joaquin Valley.

Although the groundwater loss estimates based on hydrological models (from Faunt et al., 2015, and Xiao et al., 2017) pertain to the entire Central Valley, they are similar to our amount estimated based on GPS data. This agreement is due to a wide uncertainty range associated with various estimates, and it is also due to the fact that the amount of groundwater depletion in the entire Central Valley is dominating that of the San Joaquin Valley (Ojha et al., 2018). The GRACE‐based estimate for the whole Central Valley (black curve in Figure 8), however, indicates a faster rate of groundwater depletion compared with other results (see also Table 1, row 5 vs. Table 4, row 4). This difference is possibly due to uncertainties in the GRACE data and, to a lesser degree, caused by contributions from the northern tile to groundwater loss. Therefore, we consider the amount of groundwater loss for the entire Valley obtained based on GRACE observations the upper bound estimate for that of the San Joaquin Valley.

Also in Figure 8, GRACE‐based groundwater changes provided by Argus et al. (2017) do agree well with our GPS‐based groundwater depletion; however, their calculations are not fully comparable with ours, because their averaging region reaches farther north outside the Central Valley, where climate is wetter and groundwater depletion is much less severe.

In summary, we find an overall agreement among various estimates of groundwater loss within the Central and San Joaquin valleys. However, there are also discrepancies between different data sets, which need to be investigated to identify those that are statistically significant. A better understanding of the uncertainty ranges can further improve predictions of regional groundwater changes. After a lifetime of more than 15 years, the GRACE mission ended in October 2017. A follow‐on mission has been launched successfully in May 2018 (Gravity Recovery and Climate Experiment Follow‐On, 2018), and the signal gap has grown to over a year (from summer 2017 to spring of 2019). Our results suggest that deformation data from GPS stations and, if available, InSAR have the potential to fill the gaps and complement GRACE estimates of groundwater change.

Supporting Information S1

Click here for additional data file.