Remote Sensing-Based Estimates of Changes in Stored Groundwater at Local Scales: Case Study for Two Groundwater Subbasins in California’s Central Valley

Abstract

Sustainable groundwater management requires high-quality and low-latency estimates of changes in groundwater storage (∆S_gw). However, estimates of ∆S_gw produced using traditional methods, including groundwater models and well-based measurements, typically lag years behind the present because collecting the required on-the-ground data is a time consuming, expensive, and labor-intensive task. Satellite remote sensing measurements provide potential pathways to overcome these limitations by quantifying ∆S_gw through closing the water balance. However, the range of spatial scales over which ∆S_gw can be accurately estimated using remote sensing products remains unclear. To bridge this knowledge gap, this study quantified ∆S_gw for the period of 2002 through to 2021 using the water balance method and multiple remote sensing products in two subbasins (~2700 km²–3500 km²) within California’s Central Valley: (1) the Kaweah–Tule Subbasin, a region where the pumping of groundwater to support agriculture has resulted in decades of decline in head levels, resulting in land subsidence, damage to infrastructure, and contamination of drinking water and (2) the Butte Subbasin, which receives considerably more rainfall and surface water and has not experienced precipitous drops in groundwater. The remote sensing datasets which we utilized included multiple sources for key hydrologic components in the study area: precipitation, evapotranspiration, and soil moisture. To assess the fidelity of the remote sensing-based model, we compared estimates of ∆S_gw to alternative estimates of ∆S_gw derived from independent sources of data: groundwater wells as well as a widely used groundwater flow model. The results showed strong agreement in the Kaweah–Tule Subbasin in long-term ∆S_gw trends and shorter-term trends during droughts, and modest agreement in the Butte Subbasin with remote sensing datasets suggesting more seasonal variability than validation datasets. Importantly, our analysis shows that the timely availability of remote sensing data can potentially enable ∆S_gw estimates at sub-annual latencies, which is timelier than estimates derived through alternate methods, and thus can support adaptive management and decision making. The models developed herein can aid in assessing aquifer dynamics, and can guide the development of sustainable groundwater management practices at spatial scales relevant for management and decision making.

1. Introduction

Groundwater depletion is a major problem in many areas around the world, including India, China, Iran, and the United States. Hotspots of groundwater depletion include northwestern India, the Malwa region in Punjab, the Central Valley of California, and the Shabestar Plain in Iran. Studies of these regions have shown the rapid depletion of groundwater storage in northwestern India, accompanied by spatial heterogeneity in groundwater levels and storage changes [1]. In the Malwa region of Punjab, more than 30% of tubewells showed an average depletion rate of about 40 cm/year [2]. At least 50 km³ of groundwater has been abstracted from the Central Valley of California since 1950 [3,4,5,6]. A study in the Shabestar Plain, Iran evaluated the impact of climate change on groundwater levels and found that future climate change scenarios will lead to a decline in groundwater levels [7]. These studies highlight the urgent need for better groundwater management and sustainable practices to ensure the future availability of groundwater resources.

Remote sensing has recently emerged as a powerful tool to advance groundwater science and management. A rich array of recent studies has demonstrated that remote sensing measurements from multiple sources of spaceborne data can be used to draw novel insights about groundwater systems and thus support management and decision making [4,6,8,9]. For instance, Argus et al. (2017) showed that Global Positioning System (GPS) data can be used with a hydrological model to obtain spatially distributed estimates of changes in groundwater storage for large regions. Vasco et al. (2022) showed how measurements of ground deformation from Interferometric Synthetic Aperture Radar (InSAR) and estimates of mass changes from the Gravity Recovery and Climate Experiment (GRACE) can resolve changes in confined and unconfined aquifer storage. Kim et al. (2021) highlighted the utility of remotely sensed data to support timely management and decision-making in the context of the Sustainable Groundwater Management Act, which was enacted in California in 2014. Ahamed et al. (2021) showed that remote sensing data can be used in a water balance approach to derive low-latency estimates of groundwater storage changes at coarse spatial scales.

Despite these advancements, remote sensing approaches applied at spatial scales and with latencies relevant for adaptive groundwater management and decision-making remain limited in practice. The de facto technique through which groundwater managers and decision-makers assess the status of groundwater systems continues to be through use of observations made in monitoring wells. Monitoring well networks are costly to maintain, and the measurements of groundwater levels can be sparse in both space and time due to the time and labor-intensive data collection process [10,11]. Timely information on groundwater resources is critical, especially given that managers are typically balancing the competing groundwater needs of agriculture, industry, ecosystems, and public supply at any given time. Remote sensing methods and data can alleviate some of the dependence on networks of monitoring wells, and serve as supplementary tools to groundwater managers and decision makers.

The growing availability of remote sensing data describing the components of the hydrologic cycle opens the possibility to estimate area-integrated changes in groundwater storage (∆S_gw) through a water balance approach:

∆S_gw = P  − ET + Q_in − Q_out − Q_{gw net} + R −∆S_SM  −  ∆S_SWE  − ∆S_R,

(1)

where P is precipitation; ET is evapotranspiration; Q_in and Q_out are surface water flow in and out of the area; Q_{gw net} is net groundwater flow in and out of the area; R is net runoff; ∆S_SM is the change in soil moisture, defined as moisture contained within the soil and the top few meters of the vadose zone; ∆S_SWE is the change in snow water equivalent; and ∆S_R is the change in surface water storage. Recent studies [6,12,13] have highlighted the potential of this approach and demonstrated that low-latency estimates of ∆S_gw at coarse spatial scales are achievable if remote sensing data is used to solve a form of the hydrologic water balance shown in Equation (1).

The viability of the remote sensing estimates of ∆S_gw, however, depends on many factors, including the ability of the available datasets for each variable to capture the key hydrologic conditions specific to a given area. Factors that complicate the use of certain remote sensing products can include irrigation (some remotely sensed ET datasets struggle to capture irrigated conditions), elevation (remotely sensed P products can underestimate magnitudes, especially in mountainous regions), surface water infrastructure (imported and exported water must be accounted for), and groundwater inflow and outflow (must be a relatively small component of overall water balance; groundwater flux is challenging to reliably estimate by any method). Moreover, remote sensing datasets derived from a variety of sources can have different spatial and temporal sampling frequencies which might not be suited to adequately capture dynamic conditions specific to smaller regions.

Here, we go beyond the remote sensing water balance framework described in Ahamed et al. (2021) by applying the method to local scales and assessing the accuracy, uncertainty, and applicability of remote sensing data to estimate changes in groundwater storage (∆S_gw) at spatial scales practical for management and decision making (~2700 km²–3500 km²). We select two regions within California’s Central Valley, which are characterized by contrasting surface hydrologic and hydrogeologic conditions, to assess the range of suitable conditions over which the remote sensing water balance may be applied.

2. Study Area: Kaweah–Tule Subbasin and Butte Subbasin

We selected two subbasins defined in the California Central Valley Groundwater-Surface Water Simulation Model (C2VSim) [14], shown in Figure 1, to test the approach of solving the water balance through the use of remotely sensed data. The first subbasin (C2VSIM Subregion 18) is located within the southern part of the Central Valley and constitutes an area of 3537 km². It falls within Tulare County and coincides roughly with the Kaweah and Tule Subbasins defined by the Department of Water Resources [15]. Thus, we will refer to this region as the Kaweah–Tule Subbasin. The aquifer system, which underlies the Kaweah–Tule Subbasin, is comprised of two principle layers: an unconfined aquifer composed primarily of unconsolidated alluvium which extends from the surface to the Corcoran Clay lithologic unit (located at approximately 100–300 m depth). The Corcoran Clay acts as the confining unit for a deeper confined layer composed mostly of clay [3]. There are 2003 available wells with periodic groundwater level measurements in the well database maintained by DWR [16]. This region is characterized by severe groundwater overdraft [14], which has led to land fallowing [17], significant land subsidence [18] resulting in damage to the infrastructure, and contamination of the groundwater by chemicals such as arsenic [19].

The second subbasin (C2VSim Subregion 5) is located within the northern part of the Central Valley and encompasses an area of 2794 km². The entire northern half of the area is located within Butte County, whereas the southern portion is split between Sutter County to the west and Yuba County to the east. Thus, we refer to this region as the Butte Subbasin. In contrast with the Kaweah–Tule Subbasin, the Butte Subbasin is characterized by a shallow semiconfined aquifer [3] comprised mostly of unconsolidated alluvium, and has not experienced dramatic declines in groundwater levels. For the Butte Subbasin, the California Department of Water Resources well database [16] contains 587 wells with periodic groundwater level measurements. Most of these wells are located in the northern and eastern portions of the region, with far fewer wells in the central and western portions of the study area.

The two regions have a contrasting hydrography, which makes them suitable test sites to evaluate the range of conditions over which the remote sensing water balance approach may be applied. The Kaweah–Tule Subbasin receives little natural surface water inflow; its primary source of surface water inflow is the Friant Kern Canal, shown in Figure 1 in gray, which transports water from Millerton Lake along the eastern margin of the Central Valley. As an internally drained basin, the Kaweah–Tule Subbasin does not have a natural surface water outlet. Principal reservoirs in the region include Lake Kaweah and Lake Success (Figure 1, cyan circles), which have relatively small storage capacities of 0.23 and 0.1 km³, respectively, when compared to reservoirs in the northern Central Valley.

In contrast, the Butte Subbasin receives abundant surface water from rivers originating in the Cascade and Klamath Mountains to the north and the Sierra Nevada to the east, including the Upper Feather River (shown in blue in Figure 1). The flow in some of these rivers is measured by stream gauges operated by the United States Geological Survey (Figure 1, blue circles). The Feather River and most other rivers flowing out of the region eventually drain to the Sacramento River, which lies west of the region. The headwater portion of the Feather River supplies water to Lake Oroville, the second largest reservoir in California. Lake Oroville lies on the eastern margin of the Butte Subbasin and constitutes the majority of the total surface water stored (total capacity of 4.3 km³) within and delivered throughout the region.

3. Methods and Data

We adopted the methodological approach described in Ahamed et al. (2021), with the minor modification of using exclusively remotely sensed data in the water balance. First, for each water balance component shown in Equation (1), we compared summary statistics (annual and monthly mean) and the full time series for all available data sources, which included ground-based measurements, remotely sensed datasets, land surface models, and reanalysis data. A key finding of previous studies [6,20] was that an evaluation of multiple datasets describing each variable can help with uncertainty reduction, and can be used to identify inaccurate datasets as well as guide the selection of appropriate datasets for a given region. For each water balance component, we used all available datasets to qualitatively assess seasonality and magnitude, and quantified the uncertainty and correlation among datasets describing each variable.

Because the focus of this study was the use of remotely sensed data, we exclusively selected remotely sensed data, where available, to construct ensemble water balance estimates of ∆S_gw. This contrasts with the approach taken in Ahamed et al. (2021), where the other available data derived from ground-based measurements and land surface models were also used to obtain estimates of each water balance component. ∆S_gw was estimated by using different combinations of available remotely sensed data. After estimates of ∆S_gw were computed, the net changes, trends during droughts, and seasonality was assessed. The fidelity of the remotely sensed ∆S_gw estimates were evaluated by comparing trends, magnitudes, and seasonality with independent estimates from well data as well as outputs of the C2VSim groundwater model. In subsequent sections, we formulate the water balance method, report the data sources used to estimate each variable, and describe the comparison data used in the context of the above procedure.

3.1. Method to Compute Changes in Groundwater Storage (∆S_gw) from Remote Sensing Mass Balance

The hydrologic water balance method was applied by solving Equation (1) to estimate ∆S_gw in the Kaweah–Tule Subbasin, referred to as ${\Delta S}_{gw}^{KT}$, as well as the Butte Subbasin, referred to as ${\Delta S}_{gw}^B$. For each subbasin, minor modifications to Equation (1) were required in order to accurately reflect their condition. Neither region receives a significant amount of snow, so we removed changes in the snow water equivalent (∆S_SWE). Moreover, we did not account for net groundwater flow (Q_{gw net}) for either subbasin because it is thought to be relatively small for regions within the Central Valley [14].

The Kaweah–Tule (KT) Subbasin is internally drained, having no natural outlet; hence, we also removed the Q_out component from Equation (1) resulting in the following equation:

$${\Delta S}_{gw}^{KT}=P - ET + Q_{in} + R - \Delta S_{SM} - \Delta S_R.$$

(2)

We note that the Q_in term encompasses water deliveries from the Friant-Kern Canal, shown in gray in Figure 1 on the eastern margin of the Kaweah–Tule Subbasin. Though the source of the water, Millerton Lake, is relatively far outside the boundary of the Kaweah–Tule Subbasin, the Friant-Kern Canal traverses much of the eastern portion of the basin, supplying significant volumes of water for agriculture occurring within the basin; therefore, it was factored into the estimates of inflow.

In the case of the Butte Subbasin, we found that we had inadequate measurements of surface water inflows and outflows to reliably quantify Q for every river entering and exiting the study area. To overcome this data limitation, we included the changes in the volume of water stored in Lake Oroville in the ∆S_R component, instead of considering inflows and outflows from the Feather River, the primary source of water feeding and draining Lake Oroville. Though Lake Oroville falls slightly (<4 km) outside the study region, the water stored in the reservoir is driven by a flow originating in the mountains to the north and east of the region, capturing the inflow component, while releases from the reservoir subsequently flow downstream and dictate the streamflow exiting the region, capturing the outflow component. This gives an equation to estimate ∆S_gw within the Butte Subbasin:

$${\Delta S}_{gw}^B=P - ET + Q_s + R - \Delta S_{SM} - \Delta S_R.$$

(3)

For each of the two subbasins, we evaluated ∆S_gw estimates produced by Equations (2) and (3) in order to better understand the behavior of ∆S_gw for management-relevant local areas with significant differences in hydrologic conditions, and the availability of validation datasets for each study area. For both the Kaweah–Tule and Butte Subbasins, the performance of remotely sensed estimates of ∆S_gw was evaluated through a comparison to independent ∆S_gw estimates derived from water level measurements made in groundwater wells [21] and outputs of the C2VSim groundwater model [14].

3.2. Data Used in Water Balance

The sources of data used in this study are shown in

Table 1. Data used in the water balance models developed for this study. Acronyms in the Data Type column refer to: Ground-based Measurements (GBM), Remote Sensing (RS), Reanalysis data (RE), and Land Surface Models (LSM). The optimal data selected for the Preferred Remote Sensing (PRS) scheme, which is described in Section 3.4, are denoted with a *. For more information about each source of data, including acronyms, retrieval methods, forcing data, model outputs, resolution, and availability, see the references in the last column, and the full table in Ahamed et al. (2021).

Precipitation (P)
Source	Data Type	Spatial Resolution	Availability	Reference
PRISM AN81	GBM	800 m or 4 km	1895–Present	[22]
Daymet v4	GBM	1 km	1980–Present	[23]
GPM v6 *	RS	10 km	2001–Present	[24]
CHIRPS	RS and GBM	5 km	1981–Present	[25]
PERSIANN-CDR	RS and GM	25 km	1983–2021	[26]
Actual Evapotranspiration (AET)
Source	Data Type	Spatial Resolution	Availability	Reference
MODIS (MOD16) *	RS	1 km	2001–Present	[27]
SSEBop (MODIS) *	RS	1 km	2002–2020	[28]
Terraclimate	LSM	4 km	1958–Present	[1]
FLDAS v1	LSM	10 km	1980–Present	[29]
GLDAS v2.1	LSM	25 km	1980–Present	[30]
Potential Evapotranspiration (PET)
Source	Data Type	Spatial Resolution	Availability	Reference
MODIS (MOD16) *	RS	1 km	2001–Present	[27]
Runoff (R)
Source	Data Type	Spatial Resolution	Availability	Reference
Terraclimate *	LSM	4 km	1958–2021	[31]
ECMWF ERA5 *	RE	10 km	1981–Present	[32]
GLDAS v2.1 *	LSM	25 km	1980–Present	[30]
FLDAS v1 *	LSM	10 km	1980–Present	[29]
Soil Moisture (SM)
Source	Data Type	Spatial Resolution/Number of Layers/Maximum Depth	Availability	Reference
SMAP *	RS	25 km/-/15–30 cm	2015–Present	[33]
SMOS *	RS	25 km/-/21 cm	2010–2019	[34]
GLDAS v2.1	LSM	25 km/4/2 m	1980–Present	[30]
Terraclimate *	LSM	4 km/1/6 m	1958–Present	[31]
Reservoir Storage (SR)
Source	Data Type	Temporal Resolution	Availability	Reference
California Data Exchange Center (CDEC) Monthly Reservoir Storage *	GBM	Monthly	Varies	[35]
Landcover and Crop Cover
Source	Data Type	Temporal Resolution	Availability	Reference
USGS National Landcover Database (NLCD) *	RS	Varies; typically once per 3 years	1992–Present	[36]
USDA Cropland Data Layer (CDL) *	RS	Annual	2008–Present	[37]
Discharge (Q)
Source	Data Type	Spatial/Temporal Resolution	Availability	Reference
USGS Stream Gauges *	GBM	Varies, typically daily or 15-min	Varies	[38]

, and can be grouped into the following categories: (1) remote sensing data and models used in the water balance; (2) ground-based measurements including those made in groundwater wells as well as rain, stream, and reservoir staff gauges; (3) land surface models (LSMs); and (4) reanalysis datasets. Each of these categories is denoted in the “Data Type” column shown in Table 1. All of the data sources shown in Table 1 used to estimate the terms in Equations (2) and (3) are publicly available. The following sections provide further details on the sources used to estimate each of the variables used in the analysis.

3.2.1. Precipitation (P)

At local scales, both the footprint of remotely sensed measurements and the number and spatial distribution of available rain gauges in the region will have an important impact on the quality of satellite-based and gauge-based precipitation estimates, respectively. Here, we included precipitation datasets based on remote sensing data as well as products based on the interpolation of ground-based rain gauge measurements. Ground-based products included the Daymet meteorological model [23] as well as the Precipitation Regression on Independent Slopes Model (PRISM) [22]. Remotely sensed precipitation estimates included the Global Precipitation Measurement Mission (GPM) [24]. Precipitation products that combine estimates from ground-based and remote sensing measurements included: (1) the Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks (PERSIANN) [26] model, and (2) the Climate Hazards Group InfraRed Precipitation with Station Data (CHIRPS) [25].

The PRISM dataset utilizes an interpolation scheme which incorporates measurements from five rain gauges falling within the Kaweah–Tule Subbasin and four additional stations within 10 km of the Kaweah Subbasin. There are no rain gauges located within the southern half of the Kaweah Subbasin, which could lead to uncertainties in this region. PRISM also contains measurements from six rain gauges falling within the Butte Subbasin and an additional fifteen rain gauges that fall within 10 km of the Butte Subbasin. The GPM satellite footprint is approximately 10 km, meaning that there are about 35 pixels which intersect or fall entirely within the Kaweah–Tule Subbasin, while there are approximately 27 pixels which fall entirely within or intersect the Butte Subbasin. Fractional pixels which intersect with each study area are scaled proportionally to remove the fraction which falls outside the region. For instance, if one-third of a pixel falls within the study area, the rainfall amount is scaled by 0.33 to reflect an accurate sumtotal and remove the precipitation which falls within the pixel but outside of the region from the estimates of P.

3.2.2. Actual Evapotranspiration (AET)

Actual Evapotranspiration (AET) is defined as the amount of water transferred to the atmosphere via evaporation from the land surface and shallow subsurface and via transpiration from plants. Regardless of the measurement or modeling approach, AET can be challenging to reliably estimate [39,40,41,42] due to the landscape as well as the assumptions made throughout the modeling process. Ground-based measurements, remote sensing data, and LSMs can all be used to estimate AET. One commonly used source of ground-based data to evaluate AET at local scales is AmeriFlux [43]. However, the flux tower data from AmeriFlux are sparsely available, both spatially and temporally—there are three available AmeriFlux stations within the Kaweah–Tule Subbasin, and a single AmeriFlux station within the Butte Subbasin, none of which contain long-period records of AET. Meteorological data can also be used to estimate PET, which can then be used with crop coefficients to estimate AET. The California Irrigation Management Information System (CIMIS) [44] is a network of meteorological stations widely used in California for irrigation management and scheduling. CIMIS contains five weather stations within Kaweah–Tule, and two proximal weather stations which could be useful in an interpolation. The Butte Subbasin contains only two CIMIS stations, and four additional proximal stations which may be used in an interpolation scheme.

Remotely sensed and LSM-based estimates of AET have the advantage of broader spatial coverage relative to in situ approaches. LSMs estimate AET by simulating the physics of mass and energy transfer. LSMs are typically forced using weather and radiation data and calibrated against ground-based measurements through data assimilation. However, as many studies have noted (e.g., [5,6,45]), LSM-based estimates of AET tend to perform poorly in irrigated and highly engineered areas, because irrigated conditions and water deliveries are typically not incorporated into the models. Hence, they are included in this study as a basis for comparison, but are not used to produce estimates of ∆S_gw through the water balance approach.

For remotely sensed estimates of AET, we consider data from the Operational Simplified Surface Energy Balance (SSEBop) model [28], as well as NASA’s MOD16 ET product [27]. Both products have demonstrated favorable agreement when compared against ground-based measurements of AET made in flux towers. Moreover, both products have global coverage and may be applied to local areas due to relatively high (1 km) spatial resolution. However, these products have also been shown to underestimate AET in irrigated agricultural settings [46,47] like much of the study area. To overcome this, we applied a third remote sensing-based approach to estimate AET, which estimates AET through scaling the potential evapotranspiration (PET) by a dimensionless crop coefficient [48] for irrigated regions, as determined from the National Landcover Database [36], and uses MOD16 actual evapotranspiration data to estimate AET for non-irrigated regions. The approach is described in greater specificity in Ahamed et al. (2021).

3.2.3. Stream Discharge (Q) and Surface Runoff (R)

Stream discharge is comprised of the inflows and outflows to and from the study areas resulting from natural rivers and streams as well as manmade aqueducts and canals. For both the Kaweah–Tule Subbasin and the Butte Subbasin, discharge was quantified using streamflow data from the USGS [38] (gauges shown in Figure 1). In total, we identified six stream gauges measuring inflows to the Kaweah–Tule Subbasin, which drains internally and hence does not have any outflows, as well as three stream gauges measuring inflows and one stream gauge measuring the outflow from the Butte Subbasin.

Despite the availability of data, an estimation of the true inflow and outflow for the Butte Subbasin is complicated, due to boundary conditions and ungauged canals and diversions. Regarding the former, both the western and southern margins are proximal to the Sacramento River and the Bear River, respectively, the source areas of which differ geographically from the Butte Subbasin, which is not a closed watershed. Regarding the latter, there are numerous diversions from the Sacramento River and the Feather River which transport water in and out of the study region; these include the Sutter Bypass, Sutter Extension Canal, Butte Slough, the West Borrow Ditch, and East Borrow Ditch, for which there are limited available data to rigorously determine inflow and outflow. Rather than explicitly account for the complex engineered network of diversions, canals, and aqueducts, we elected instead to use estimates of changes in reservoir storage for Lake Oroville as the estimate of inflow (Q_in) and outflow (Q_out) for the Feather River. This reservoir is 4 km outside the study region, but acts as the primary regional driver in terms of regulating flows to and from the Butte Subbasin.

The runoff accounts for a relatively small percentage of inflow relative to Q [3,49]. However, it can still contribute meaningful volumes of water to a given area during periods of intense precipitation. Some smaller watersheds cross the subbasin boundaries and thus contribute to diffuse surface water runoff to the basins, as opposed to streamflow at a single point [14]. Thus, we include runoff entering the Kaweah–Tule and Butte subbasins by including the small watersheds defined in C2VSIM adjacent to both subbasins, and the collected runoff data available from Terraclimate, GLDAS, and FLDAS LSMs for these areas. We also included runoff from the European Center for Medium-Range Weather Forecasting (ECMWF) produced as part of the ECMWF Re Analysis version 5 (ERA5) [32], which uses historical observations and numerical weather models to estimate runoff.

3.2.4. Reservoir Storage Changes (∆S_R)

Reservoirs are an important component of changes in surface water storage [50], especially within the Central Valley [6]. For the reservoirs falling within each study region, monthly data describing volumetric changes in reservoir storage were accessed from the California Data Exchange Center [35]. Within the Kaweah–Tule Subbasin, there are two reservoirs: Lake Kaweah and Lake Success. For the Butte Subbasin, CDEC (2022) contains data for five reservoirs, including Lake Oroville, Thermalito, Thermalito Afterbay, Thermalito Forebay, and Thermalito Divers Pool. Lake Oroville constitutes >95% of the total ∆S_R for the Butte Subbasin, and falls 4 km outside of the Butte Subbasin, but was included within our analysis to account for the inflows and outflows from the Feather River, which were difficult to quantify using gauge-based data (as described in the previous section). To determine the monthly fluctuations in water stored within the Kaweah–Tule Subbasin and the Butte Subbasin, the total storage volume for all reservoirs falling within each region was compiled and summed. To determine storage changes, which serve as inputs to the water balance model, summed reservoir storage volumes were differenced from the starting value at the beginning of the study period (1 October 2001). Because water storage in rivers, lakes, aqueducts, and streams is a negligible portion of the overall water budget, we did not include estimates of volumetric changes in these surface water bodies in our calculations.

3.2.5. Soil Moisture Water Storage Changes (∆S_SM)

Soil moisture measurements acquired from in situ sensors are sparsely distributed—the Soil Climate Analysis Network (SCAN) [51], which is a network of monitoring sites maintained and operated by the Natural Resource Conservation Service, contains just 15 stations within California; most have data for the last 8–12 years. We therefore considered alternative data sources, namely remote sensing and LSM-based datasets, in order to estimate ∆S_SM. In terms of remote sensing, reliable estimates of soil moisture can be obtained using passive microwave-based radiometers. These include the Soil Moisture Ocean Salinity (SMOS) [34] satellite, which contains data for 2010 to 2019; as well as the Soil Moisture Active-Passive satellite (SMAP) [33,34] which contains data for 2015 to present. Land Surface Models which simulate ∆S_SM dynamics are available for the entirety of our study period between 2002–2021; these include GLDAS [30] and Terraclimate [31].

Within land surface models, soil moisture is a state variable, and is often partitioned between depth intervals. The total soil moisture can be obtained by summing the soil moisture across all depth intervals. Hence, for each LSM, we summed the soil moisture for each time step across all of the layers in order to determine the total soil moisture. Depth intervals are model specific: GLDAS contains four depth intervals: 0–10 cm, 10–40 cm, 40–100 cm, and 100–200 cm. In contrast, Terraclimate uses spatially distributed single-layer depth intervals, as described in Wang-Erlandsson et al. (2016), which, within the study regions, range from 1 m to 6 m in depth.

3.3. Estimate of Correlation among Water Balance Components and Uncertainty within Components

We quantified the correlation and uncertainty associated with each water balance component as described in Ahamed et al. (2021), first determining the correlation between the datasets describing each water balance variable, then assessing the uncertainty in the water balance components. To estimate correlation, we used the Pearson’s correlation coefficient [52] between two variables x and y:

$$P_{xy}=\frac{\sum _{i=1}^n{(x}_i-\stackrel{-}{x}){(y}_i-\stackrel{-}{\mathrm{y}})}{\sqrt{{\sum _{i=1}^n{(x}_i-\stackrel{-}{x})}^2}\sqrt{{\sum _{i=1}^n{(y}_i-\stackrel{-}{\mathrm{y}})}^2}}$$

(4)

where x_i and y_i are the elements of each variable, and $\stackrel{-}{x}$ and $\stackrel{-}{\mathrm{y}}$ are the mean of each variable.

To estimate uncertainty, we used Triple Collocation (e.g., [53,54,55,56,57]), which has been shown to be a suitable uncertainty and error estimation approach for data describing hydrologic and meteorological measurements. Given at least three independent and collocated datasets describing the same variable, Triple Collocation approximates the random error variances (𝜎²) associated with each dataset. For three area-averaged collocated estimates of the same variable, the random error variances (𝜎²) are given by:

𝜎²_𝑥 = ⟨(𝑥−𝑦_s)(𝑥−𝑧_s)⟩

𝜎²_𝑦 = ⟨(𝑦−𝑥_s)(𝑦−𝑧_s)⟩

𝜎²_𝑧 = ⟨(𝑧−𝑥_s)(𝑧−𝑦_s)⟩,

(5)

where brackets indicate the temporal mean, x, y, and z are the three datasets, and the subscript s denotes the scaling of the x, y, or z dataset to the mean and standard deviation of the reference dataset [54]. If more than three datasets are available for a single variable, we considered each unique combination of triplets, estimated 𝜎² as in Equation (5), and reported the mean 𝜎 for all triplets containing each specific dataset. For example, if six datasets of the same variable were available, the number of unique triplets is 6!/(6 − 3)! = 120. Of the 120 unique triplets, each individual dataset appears (3/6) × 6!/(6 − 3)! = 60 times. For the 60 triplets where a dataset appears, the random error variance (𝜎) is computed. We then report the mean and standard deviation of 𝜎 (in mm) across the 60 triplets as the uncertainty associated with the dataset in question.

3.4. Estimation of Changes in Groundwater Storage

Estimates of ∆S_gw derived from solving Equations (2) and (3) can vary depending on the data sources used. What is not well known is to what extent the use of various input datasets can affect the ∆S_gw in terms of magnitude, seasonality, and trends during extremes. Moreover, it is not well understood how remotely sensed estimates of ∆S_gw, derived through the water balance, compare with estimates of ∆S_gw derived from other methods (e.g., wells and groundwater flow models). Ahamed et al. (2021) compared estimates of ∆S_gw derived from LSMs, ground-based measurements, and remote sensing data with independent estimates from the Gravity Recovery and Climate Experiment, as well as wells and C2VSim, finding that remote sensing-based estimates of ∆S_gw were able to effectively capture seasonal fluctuations, net changes, and trends during droughts when compared to validation datasets. Therefore, in this study, we focus on the application of remote sensing data using the water balance approach to recover ∆S_gw at local scales for the Kaweah–Tule and Butte Subbasins. This will allow us to better understand the applicability limits of remotely sensed data and analyses, which typically are applied to larger spatial scales than the Kaweah–Tule Subbasin and Butte Subbasin analyzed in this study. We considered two scenarios to estimate ∆S_gw:

(1)

Ensemble Remote Sensing scenario (ERS): Estimates of ∆S_gw were computed using an ensemble approach which considers the inputs of remotely sensed data for all hydrologic variables when available, those being three for P, three for AET, and three for SM. Because remotely sensed estimates of SM are available only from 2010 onwards, we used Terraclimate from 2002–2021, SMOS from 2010–2019, and SMAP from 2015–2021 for the ERS scenario ensembles. Lastly, R was derived from the three LSMs and the ERA5 Reanalysis dataset, resulting in 108 distinct scenario combinations for both subbasins (N_P × N_ET × N_SM × N_R × N_Q × N_SR = 3 × 3 × 3 × 4 × 1 × 1 = 108).

(2)

Preferred Remote Sensing (PRS) scenario: Estimates of ∆S_gw were computed using a subset of all remote sensing datasets used in the ERS scenario. We selected the remote sensing data for different hydrologic components that are generally found to exhibit a lower error than others in the selected study regions. The basis of the selection of the preferred datasets includes two factors: (1) a review of the literature for each data source, in order to determine the range of applicable conditions, and (2) an assessment of the seasonality, magnitudes, correlation, and uncertainty present for the estimates of each variable shown in Table 1 and described in the previous section. For precipitation, we used the state-of-the-art and newly released GPM Integrated Multisatellite Retrievals (IMERG) as our remote sensing estimate because of higher accuracies relative to legacy products [58]. For AET, we used the temporal mean of SSEBop and MODkc estimates, because of the sensitivity to ET dynamics in irrigated regions [28], which is a key feature of the study regions. The MOD16 product, in contrast with SSEBop and MODkc, is known to underestimate ET in croplands [45] For the case of SM, there is no available remotely sensed estimate that covers the entire study period. To approximate a continuous timeseries of SM, which is required in order to complete our analysis over the study period, we computed the temporal average across all data sources in the ERS scenario (Terraclimate, SMOS, and SMAP), and used that as our final SM estimate. Similarly, for the case of R, we used three LSMs and ECMWF ERA5 reanalysis datasets to compute the temporal mean, which was used in the water balance model. The dataset used in the PRS scenario is marked with a * in Table 1.

Both the trend and seasonality of ∆S_gw estimates generated from the remote sensing and water balance approach were used to assess performance. To better identify these dynamics, we differenced each scenario’s time series estimate of ∆S_gw from the mean of the first five years. This procedure effectively decreases the variability among estimates while still capturing the prevailing long-term trends. The 5-year mean differencing procedure is identical to the 5-year calibration period used in Ahamed et al. (2021), where it was found to provide a good balance between characterizing the seasonal cycle and recovering long-term trends. In addition, the 5-year calibration period is commonly applied to time-varying gravity data from GRACE in order to estimate the changes in stored groundwater for larger regions. For the ERS scenario specifically, we computed the mean ∆S_gw across all ensembles, and used the mean ∆S_gw as the final estimate.

3.5. Validation of Remotely Sensed ∆S_gw Estimates

In order to be deemed effective for quantifying ∆S_gw at any given spatial and temporal scales, a method should ideally characterize the long-term changes, seasonal fluctuations, and capture extreme events, whether those are droughts or wet periods. Validating the estimates of ∆S_gw remains challenging, since well-based estimates can have spatial and temporal data gaps, biases, and uncertainties (Section 2) [10,11], and estimates derived from groundwater models are calibrated on available well data. Despite these factors, wells and groundwater models are the most widely used sources of data currently used to measure and monitor changes in groundwater storage. Thus, in order to determine the fidelity of ∆S_gw estimates determined using the remote sensing water balance approach, we compared each scenario to estimates of ∆S_gw derived independently, which we describe here. The independent sources of data used to estimate ∆S_gw are based on: (1) groundwater levels measured in wells over the period 2002–2019 [21], and (2) outputs of the C2VSim groundwater flow model [14], which are available during the time period from 2002–2015. In the next section, we describe each of the comparison sources in more detail, and detail the metrics which are used for validation of the remote sensing approach.

3.5.1. Changes in Groundwater Storage Estimated from Wells

We used data obtained from wells describing groundwater levels for the time period 2002–2019 [21] in order to calculate ∆S_gw for the both the Kaweah–Tule and Butte Subbasins. Past studies (e.g., [6,21,59]) have also used groundwater level measurements to estimate ∆S_gw for the entire Central Valley. The data describing groundwater levels were obtained from the California Department of Water Resources (DWR) California Statewide Groundwater Elevation Monitoring (CASGEM) database [16]. This database contains measurements from wells used for a variety of purposes, including groundwater monitoring, irrigation, and domestic uses. The Kaweah–Tule Subbasin contains measurements from 2003 wells, while the Butte Subbasin contains measurements from 587 wells. In order to estimate ∆S_gw from water level data, we multiplied water level changes with the area of each subbasin and a storage coefficient, which captures the spatial variability in an aquifer’s volumetric yield for a unit decrease in water level. A three-step procedure (described in detail in [6]) was employed in order to estimate ∆S_gw from the well data: In the first step, well-based measurements of groundwater levels were spatially interpolated using inverse distance weighting at grid cell size of 10 km. It should be noted that this interpolation scheme permits the inclusion of wells which fall outside the subbasin boundaries in order to estimate groundwater level changes near the borders. In the second step, the water level data were used to estimate monthly water level changes for each grid cell. In the final step, monthly water levels were multiplied by spatially distinct storage coefficients weighted based on (1) the magnitude of water level fluctuation at a given grid cell and (2) the specific storage and specific yield values from C2VSim and the Central Valley Hydrologic Model (CVHM) [3]. The methodological details of steps 1–3 are given in Alam et al. (2021).

3.5.2. Model-Based Estimates of Changes in Groundwater Storage

The Central Valley Groundwater-Surface Water Simulation Model (C2VSim) [14] is a widely-used model created by California’s Department of Water Resources (DWR). C2VSim simulates water fluxes through the land surface, root zone, and groundwater system on a finite-element grid. A fine scale version of C2VSim (C2VSimFG), calibrated on observed surface water flows, groundwater heads, groundwater head differences between well pairs, and stream–groundwater flows for the period between September 1975 and October 2003 [60] has a spatial resolution of ~1.5 km, covering the 1975–2015 hydrologic years. The model’s elements are grouped into 21 subregions, shown in Figure 1. We used the outputs for subregions 18 and 5, referred to as the Kaweah–Tule and Butte Subbasins (Section 3; Figure 1), to estimate model-based ∆S_gw in our study regions. Since it was first calibrated in 2013 [14], many studies have used C2VSim to analyze surface and groundwater dynamics within the Central Valley (e.g., [61,62,63,64]. Of all the subbasins within the Central Valley, C2VSimFG suggests that the Kaweah–Tule Subbasin has experienced the most groundwater loss (−0.75 km³/yr) since 1960, while Butte has experienced a minimal loss (~0 km³/yr) of groundwater [65].

3.5.3. Evaluation of Remotely Sensed ∆S_gw Estimates

A comparison between various approaches to estimate ∆S_gw is essential to determine the range of suitable estimation approaches, which in turn depends on the size of the area, current and historic instrumentation, boundary conditions, and a host of region-specific factors (e.g., rain vs. snow dominated, irrigated vs. natural vegetation, natural vs. engineered surface water bodies). An important finding of Ahamed et al. (2021) was that in order to produce reliable estimates of ∆S_gw through the remote sensing approach, there must be accurate data describing the water balance components with the greatest contribution, in terms of magnitude, to the overall water budget in an area. For example, accurate estimates of the snow water equivalent were essential to resolve ∆S_gw for the entire Central Valley Watershed (~150,000 km²), which includes the Sierra Nevada mountains, while accurate estimates of ET were essential to resolve ∆S_gw for the Central Valley (~55,000 km²), a highly irrigated agricultural region.

There are a number of criteria which could qualify a high-quality estimate of ∆S_gw. These include the ability to recover: (1) seasonality, (2) long term changes, and (3) dynamics during extreme events (e.g., multiyear droughts). A comparison of multiple approaches to estimate ∆S_gw for a given area can reveal the strengths and weaknesses of various measurement approaches (e.g., wells, models, water balance, gravity, etc.) and enhance our understanding of which are best suited to capture the dynamics of hydrologic systems, especially across spatial scales and domains. Each single source of data describing ∆S_gw has inherent limitations, but earlier studies have successfully compared ∆S_gw from multiple approaches in order to understand ∆S_gw behavior in the Central Valley region [3,4,5,6,14,59,66,67,68]. Thus, in this study, we included multiple ∆S_gw estimates from independent sources, and adopt the following metrics to evaluate the performance of ∆S_gw estimates: (1) the net change in ∆S_gw for three time windows (long term: 2002–2021; multi-year droughts: October 2006–September 2009 and October 2011–September 2015), (2) the rate of change in ∆S_gw (km³/year) during each drought period (2006–2009 and 2011–2015), and (3) the seasonality ∆S_gw estimates derived from each water balance scenario and the two sources of validation data.

4. Results and Discussion

The remote sensing water balance approach was used to estimate ∆S_gw for the Kaweah–Tule and Butte Subbasins for the hydrologic years 2002–2021. First, we compare the magnitudes of water balance components for each of the subbasins and assess the correlations and uncertainties among water balance variables. Next, we discuss the ability of the remote sensing approach to produce high-quality estimates of ∆S_gw by comparing water balance results with those derived wells and C2VSim. Lastly, we discuss the limitations associated with the approach, and highlight future research opportunities and challenges.

4.1. Correlations and Uncertainty for Water Balance Components

The datasets shown in Table 1 were collected and analyzed for the Kaweah–Tule and Butte subbasins. For each variable, the magnitudes, seasonality, uncertainty, and time series for each of the data sources is shown for Kaweah–Tule in Figure 2 and Figure 3 and for Butte in Figure 4 and Figure 5. In general, we observe that the Kaweah–Tule Subbasin has less annual precipitation, runoff, and changes in reservoir storage relative to the Butte Subbasin, consistent with the regional hydroclimatic differences described in Section 3. The following sections highlight the key differences, including correlations and uncertainties, among datasets describing each water balance component in the study areas.

4.1.1. Precipitation (P)

Generally, the estimates of precipitation within the Kaweah–Tule Subbasin (Figure 2A and Figure 3A) and Butte Subbasin (Figure 4A and Figure 5A) agree between remotely sensed and ground-based datasets in terms of seasonality and magnitudes. This is consistent with the findings of Ahamed et al. (2021), which suggested strong overall agreement between ground-based and remotely sensed precipitation datasets for the Central Valley and Central Valley Watershed. However, there are region-specific differences between precipitation datasets. In the case of the Kaweah–Tule Subbasin, we find that the gauge-based estimates of P derived from PRISM and Daymet are marginally higher (~100 mm/year) than the remotely sensed estimates derived from GPM, PERSIANN, and CHIRPS (Figure 2A). This may be due to the lack of rain gauge coverage in the southern half of the Kaweah–Tule Subbasin, described in Section 4.2.1, which could lead to higher uncertainties associated with gauge-based products in this region. In contrast, we find that for the Butte Subbasin, the remotely sensed estimates of P derived from GPM, PERSIANN, and CHIRPS tend to be higher than the ground-based estimates of P derived from PRISM and Daymet (Figure 4A). We attribute this difference to two factors: (1) relatively sparse rain gauge coverage within both study regions, which may not capture localized storms, and (2) uncertainties associated with remotely sensed precipitation data tend to decrease over larger spatial domains, which can persist at smaller scales [69] such as the Butte Subbasin. For both the Kaweah–Tule and Butte Subbasins, PERSIANN had the highest Triple Collocation errors relative to other sources of data, shown in the fourth column of Figure 2A and Figure 4A respectively. PERSIANN also has the highest uncertainties across the Central Valley and Central Valley Watershed, which encompass both the Kaweah–Tule and Butte Subbasins [6], suggesting that other products may be better suited to capture regional precipitation patterns.

4.1.2. Evapotranspiration (ET)

Among input water balance components, evapotranspiration contains the highest variance between datasets. Ahamed et al. (2021) also highlighted ET as the most variable and uncertain water balance component for the Central Valley and Central Valley Watershed areas. As noted in Section 4.2.2, LSM-based estimates of ET struggle to accurately simulate the highly irrigated and engineered conditions specific to the Central Valley, including the Kaweah–Tule and Butte Subbasins. Figure 2 and Figure 4 show that LSM-based estimates of ET, which include GLDAS, FLDAS, and Terraclimate, suggest a maximum ET in April and May, which is inconsistent with the timing of irrigation in the region. This problem is not exclusive to LSM-based data; the MOD16 ET product [27], known to underestimate ET in croplands [45], also suggests peak ET during the spring season, both for the Kaweah–Tule Subbasin (Figure 2B) and Butte Subbasin (Figure 4B). In contrast, remotely sensed estimates of ET, which are constructed to be more sensitive to irrigation, the SSEBop [28] and MODkc approaches register a more realistic ET peak for irrigated conditions occurring during the summer months of June–August for both Kaweah–Tule (Figure 2B) and Butte (Figure 4B). The correlation matrix in the 4th column of Figure 2A and Figure 4A also show that there is a high correlation between ET datasets which are sensitive to irrigation (SSEBop and MODkc), and a high correlation between ET datasets which are insensitive (GLDAS, FLDAS, Terraclimate, MOD16). Triple Collocation suggests a mean error of approximately 20–40 mm across all datasets for both the Butte Subbasin and Kaweah–Tule Subbasin.

Our findings support the emerging critical need for high-quality datasets describing ET over longer timescales, which do not require landcover maps or crop coefficients to produce accurate estimates across a range of conditions. One such project that has the potential to meet this need, especially for irrigated regions, is Open ET [70], which simplifies the dissemination of high-resolution remotely sensed ET datasets. The numerous ET models within OpenET (i.e., [71,72,73,74,75,76,77]) improve the spatial resolution and ability to capture irrigation relative to LSMs and legacy ET products.

4.1.3. Soil Moisture Changes (∆S_SM)

Remotely sensed estimates of soil moisture [33,34] are a relatively recent advancement in remote sensing science, with ∆S_SM estimates from the Soil Moisture Ocean Salinity (SMOS) satellite being available from 2009–2019 and estimates from the Soil Moisture Active Passive (SMAP) satellite available from 2015 until present. Along with the estimates of ∆S_SM derived from LSMs, satellite-based estimates are shown for the Kaweah–Tule Subbasin in Figure 2C and Figure 3C and for the Butte Subbasin in Figure 4C and Figure 5C.

Passive microwave-based remotely sensed data are directly sensitive to moisture content, allowing for high quality estimates of soil moisture, albeit with a coarse resolution of approximately 25 km. In contrast, LSMs represent SM as a state variable, allowing for surface and subsurface mass exchanges, but typically contain a highly simplified version of the subsurface [78,79]. Moreover, the coarse resolution of most LSMs makes it challenging to capture surface–subsurface interactions [80] that are driven by topographic and geologic features at much finer scales. The estimates of soil moisture derived from remotely sensed data and LSMs are consistent in terms of seasonality and differ in terms of magnitude—LSMs tend to have a greater ∆S_SM magnitude than remotely sensed datasets. This could be due to (1) a relatively shallow penetration depth of satellite-based estimates and/or (2) simplified LSM subsurface interactions—some of the soil moisture fluctuation may instead be attributable to groundwater, as opposed to soil moisture. One LSM with lower magnitude fluctuations in ∆S_SM is Terraclimate, which has a high correlation (0.75–0.9) with both passive microwave remote sensing data sources, for both the Butte and Kaweah–Tule subbasins (shown in Figure 2C and Figure 4C, respectively). In contrast with other LSMs, Terraclimate is a relatively high-resolution model (4 km) and contains a spatially varying depth component [81], which may be able to capture some of the local scale subsurface interactions and therefore improve upon other LSM-based estimates of ∆S_SM.

Agreement between Terraclimate with SMOS and SMAP, suggested by a strong correlation, consistent seasonality, and comparable magnitudes, facilitated the temporal extrapolation performed in the PRS scenario, where Terraclimate LSM estimates of ∆S_SM were for the time period which predates the availability of passive microwave remote sensing products (2002 to 2010). As shown in Figure 2, Figure 3 and Figure 4, estimates of ∆S_SM have a higher variance than estimates of P, but a lower variance than estimates of ET—this may simply reflect the maturity of remote sensing science and algorithm development for each of these hydrologic variables. P is the most studied due to its importance in weather forecasting; ∆S_SM is the next most studied, due to its importance in climate modeling; and ET was the least studied until recent times, when a growing focus on water management has led to an increase in ET research and product development (e.g., [70]).

4.1.4. Runoff (R), Discharge (Q_in and Q_out), and Reservoir Storage (∆S_R)

A comparison of Figure 2D and Figure 4D shows that the Butte Subbasin receives considerably more runoff than the Kaweah Subbasin, which is consistent with regional hydroclimatic conditions. For both the Butte Subbasin and the Kaweah–Tule Subbasin, estimates of runoff agree in terms of magnitude and seasonality, with the majority of runoff occurring during the spring period from approximately February to May. Runoff is the smallest component in the water balance for both the subbasins, and has an overall low triple collocation error of approximately 10–20 mm across datasets.

Estimates of inflow and outflow can be challenging to compile for highly engineered regions where structures divert and channel water. These systems often do not have reliable flow records (described in Section 3.2.3). However, these data are important components of the overall water balance, and reasonable estimates are needed in order to reliably estimate the inflow and outflow traveling to and from a region. For the case of the Kaweah–Tule Subbasin, we included the Friant-Kern Canal as the inflow, because much of this water travels from Millerton Lake to the Kaweah–Tule Subbasin. For the Butte Subbasin, we could not find adequate gauges in order to accurately estimate the inflows and outflows specific to the Feather River. Noting these limitations, and recognizing the highly engineered nature of the system, we elected to incorporate changes in reservoir storage (∆S_R) for Lake Oroville as a proxy for the inflow and outflow of the Feather River. Though Lake Oroville falls slightly outside of the study area, the reservoir is the second largest in California with a storage capacity of 4.3 km³, and therefore constitutes a critical component of regional surface water dynamics in and around the Butte Subbasin.

4.2. Comparison and Evaluation between Estimates of ∆S_gw Derived from Remote Sensing, Wells, and Groundwater Model

We compared monthly ∆S_gw estimates during 2002–2021 for the remote sensing scenarios (described in Section 3.4) against the independent ∆S_gw estimates derived for well data (Section 3.5.1) and C2VSim (Section 3.5.2), available from 2002–2019 and 2002–2015, respectively. The ∆S_gw results are shown in Figure 6; the ERS scenario ensembles are shown in light red, with the ensemble mean as a red dashed line, the PRS scenario is shown as a blue line, C2VSim estimates are shown as a green line, and well-based estimates are shown as an orange line. Figure 6A shows the results for the Kaweah–Tule Subbasin, and Figure 6B shows the results for the Butte Subbasin. Trendlines for each estimate of ∆S_gw, along with annual rates of depletion, are shown with the same color scheme during the 2006–2009 and 2011–2015 droughts.

4.2.1. Evaluation of Remote Sensing Results for Kaweah–Tule Subbasin

Remote sensing results show a strong agreement in terms of seasonality, net change, and drought dynamics when compared to the validation estimates of ∆S_gw determined from wells and the C2VSim groundwater model (Figure 6A;

Table 2. Summary of net ∆S_gw (km³) changes derived from each dataset within the Kaweah–Tule and Butte Subbasins for the comparison period (second column), 2006–2009 drought (third column), and the 2011–2015 drought (right column).

Kaweah–Tule
Estimate	Net ∆Sgw (km3)	2006–2009 Net ∆Sgw (km3)	2011–2015 Net ∆Sgw (km3)
ERS (2002–2021)	−8.4 (−0.44/yr)	−2.7 (−0.9/yr)	−5.6 (−1.4/yr)
PRS (2002–2021)	−7.3 (−0.38/yr)	−2.7 (−0.9/yr)	−5.6 (−1.4/yr)
C2VSIM (2002–2015)	−8.1 (−0.62/yr)	−3.6 (−1.2/yr)	−5.2 (−1.3/yr)
Wells (2002–2019)	−6.5 (−0.39/yr)	−2.0 (−0.7/yr)	−6.0 (−1.5/yr)
Butte
Estimate	Net ∆Sgw (km3)	2006–2009 Net ∆Sgw (km3)	2011–2015 Net ∆Sgw (km3)
ERS (2002–2021)	−1.3 (−0.07/yr)	−2.1 (−0.7/yr)	−0.8 (−0.2/yr)
PRS (2002–2021)	−6.6 (−0.35/yr)	−2.7 (−0.9/yr)	−2.4 (−0.6/yr)
C2VSIM (2002–2015)	−0.6 (−0.04/yr)	−0.6 (−0.2/yr)	−0.8 (−0.2/yr)
Wells (2002–2019)	−0.7 (−0.06/yr)	−0.6 (−0.2/yr)	0 (0/yr)

) within the Kaweah–Tule Subbasin. The PRS-based estimate has high Pearson’s correlations (Section 3.3) of 0.89 when compared to the well-based estimates and 0.93 when compared to the C2VSim-based estimates of ∆S_gw, while the ERS-based ensemble mean has a Pearson’s correlation of 0.91 when compared to the well-based estimates and 0.93 when compared to C2VSim-based estimates. The strong agreement between remotely sensed estimates of ∆S_gw with wells and C2VSim, defined in terms of consistent seasonality, correlation, and comparable net change, is similar to the findings for the remote sensing approach applied to the entire Central Valley, described in Ahamed et al. (2021). These findings suggest that the remote sensing approach is able to recover ∆S_gw very similar to existing methods, but with the benefits of lower latencies and fewer on-the-ground data requirements: estimates of ∆S_gw shown in Figure 6 are available until the hydrologic year 2022, whereas well-based and C2VSim-based estimates are years out of date.

The Kaweah–Tule Subbasin receives far less surface water from rivers, reservoirs, canals, and aqueducts; these sources collectively constitute a far lower magnitude volumetric surface water flux (<0.5 km³/year; shown in Figure 4 and Figure 5) and a far smaller proportion of the overall water balance. Similarly, inflows, outflows, and changes in the reservoir storage represent a relatively small magnitude compared to other water balance components at the larger scale of the Central Valley (55,000 km²).

4.2.2. Evaluation of Remote Sensing Results for Butte Subbasin

Within the Butte Subbasin study area, there are significant variations between remotely sensed and validation datasets describing ∆S_gw, both in terms of seasonal fluctuations and longer-term trends (Figure 6B; Table 2). Remote sensing estimates of ∆S_gw, shown in Figure 6B as red (ERS) and blue (PRS) lines, suggest larger annual fluctuations, and greater long-term declines in groundwater than estimates based on wells and C2VSim, shown in green and orange, respectively, which show almost no seasonal variation or longer-term trend. Lack of monitoring wells within the western portion of the Butte Subbasin may result in an inability to capture the true magnitude of seasonal fluctuations. Unrealistically high magnitude fluctuations produced by the remotely sensed ∆S_gw estimates are likely due to the inclusion of changes in reservoir storage within the water balance model. As discussed in Section 4.1, Lake Oroville has a profound impact on surface water storage, as well as the inflows and outflows to and from the Butte Subbasin. The impact of storage changes at Lake Oroville on the overall water balance is especially visible during 2016–2017 in the remote sensing estimates of ∆S_gw for the Butte Subbasin, shown in Figure 6B. The correlation between remotely sensed and validation data is considerably lower for the Butte Subbasin: the PRS-based estimate has modest Pearson’s correlations (Section 3.3) of 0.65 with a well-based estimate and 0.59 with a C2VSim-based estimates of ∆S_gw. The ERS-based ensemble mean has a slightly lower Pearson’s correlation of 0.56 with a well-based estimate and 0.53 with the C2VSim-based estimates, but matches more closely in terms of net ∆S_gw over the length of the study period (Table 2).

The difference in the seasonal magnitude of ∆S_gw estimates derived from the remote sensing estimates and validation estimates of ∆S_gw highlights an important limitation of the water balance method: accurate estimates of inflow and outflow are essential in order to close the water balance. The inability of the remote sensing method to recover reasonable seasonal changes comparable to validation datasets within the Butte Subbasin demonstrates the importance of boundary conditions on estimates of ∆S_gw. Ideally, the remote sensing water balance may be applied with confidence to a watershed with a single outlet, for which there are continuous, long-term, and high-quality estimates of discharge available. Confidence in the water balance method decreases with an increasing number of engineered or natural sources of inflow or outflow for which measurements are not available.

4.2.3. Evaluation of Net Changes and Trends during Droughts for Study Areas

To assess the quality of the remote sensing approach, ∆S_gw estimates for each of the remote sensing and validation approaches were evaluated over three timescales: (1) the entire period of availability for each dataset, and each of the drought periods occurring (2) between 2006–2009 and (3) from 2011–2015. Table 2 shows the ∆S_gw results during each time interval for both the Butte and Kaweah Subbasins. Every approach suggests that the Kaweah–Tule Subbasin experienced greater groundwater loss compared to the Butte Subbasin during the entire period between 2002–2021, as well as during each of the droughts occurring from 2006–2009 and 2011–2015.

The Kaweah–Tule Subbasin experienced rapid declines in groundwater storage during both the 2006–2009 drought and 2011–2015, shown in Figure 6 and Table 2. All remotely sensed and validation estimates of ∆S_gw suggest that annual rates of groundwater decline during the second drought were 50–100% greater than the rates of groundwater decline during the first drought. Our analysis suggests that ∆S_gw in Kaweah–Tule accounts for a significant proportion of the total decrease in groundwater storage for the Central Valley over the same time period. This finding is consistent with C2VSim results [14] and other studies [6,65]. The finding of more severe groundwater depletion for the 2011–2015 drought in comparison to the 2006–2009 drought is also consistent with many previous studies (e.g., [4,5,6]) which quantified drought dynamics for the Central Valley region, including both the Butte and Kaweah–Tule subbasins, for the same time periods.

In contrast to the drought dynamics in the Kaweah–Tule Subbasin, data describing ∆S_gw for the Butte Subbasin suggest that the 2006–2009 drought caused more groundwater loss than the 2011–2015 drought, shown in Table 2. This is corroborated by the results from the ERS scenario, PRS scenario, and well data, while C2VSim suggests comparable amounts of groundwater loss during each drought, with the 2011–2015 drought causing slightly more groundwater storage loss. The widely used Palmer-Drought-Severity-Index (PDSI) [82], suggests that the Butte Subbasin was in a drought for 67% of the 2006–2009 period, while the region was in a drought for only 55% of the 2011–2015 period, indicating that the 2006–2009 drought may have been slightly more severe than the 2011–2015 period within the Butte region. It is possible that the limited number of wells, especially in the western portion, reduces the quality of well-based and model-based estimates of ∆S_gw, since dense coverage is needed to compute robust well-based estimates of ∆S_gw and model-based estimates are calibrated on the measurements made in wells.

We observed differences in ∆S_gw between the two study domains, both in terms of trend and net change. Such differences can be explained by the difference in seasonal water demand coupled with the magnitude and timing of the available water supply. For instance, water demand is significantly high in late spring and summer in both regions, but due to greater agricultural production in the south and lower precipitation, more groundwater is extracted in Kaweah compared to Butte [14,21]. Human intervention in these regions impacts groundwater use and supply in important ways—reservoirs store winter flows and release them to the Central Valley during late spring and summer when crop water demand is at its highest. The amount of water availability is thus highly dependent on the water stored in the reservoirs, and reservoir releases have a large impact on groundwater recharge through (1) the provision of irrigation water, some of which is not consumed by crops and (2) recharge occurring from losing stream channels and leaking aqueducts, which also recharges groundwater.

4.3. Limitations and Outlook for Remote Sensing Water Balance

This study showed that local-scale remote sensing water balance estimates of ∆S_gw can be comparable in trend, magnitude, and seasonality to independent estimates of ∆S_gw obtained through well-based measurements, and through groundwater flow models. Despite advantages of the remote sensing approach, including rapid computation, low latency relative to current methods, and a growing number of datasets available to assess uncertainty in the estimates, there are four key limitations associated with the approach that must be considered when evaluating the potential application to alternative areas.

Firstly, the method produces area-integrated estimates, while well-based estimates and groundwater models can produce a higher spatial resolution and gridded estimates of ∆S_gw, as long as dense well coverage is available. Secondly, as shown in this study and as discussed in Section 4.2.1 and Section 4.2.2, boundary conditions describing surface water inflow and outflow to a region can be a critical component of the overall water balance. In regions without such data, the water balance technique may be subject to large errors, or may prove impossible altogether. This data requirement may be alleviated by the Surface Water Ocean Topography (SWOT) satellite mission [83]. Launched in December 2022, SWOT will measure terrestrial surface elevations via altimetry, thereby providing the supplementary data needed to accurately estimate Q and S_R, critical components required to accurately estimating ∆S_gw through the water balance approach.

Thirdly, errors associated with input datasets can propagate through to estimates of changes in groundwater storage. For example, some remotely sensed estimates of ET struggle to capture irrigated dynamics, which skew ∆S_gw estimates computed using these data sources. This problem can be alleviated through ensemble modeling, and through the selection of the most accurate available data to adequately reflect the most critical water balance components specific to a given region (e.g., irrigation, snow, surface water infrastructure). The satellite footprint and resolution of various remotely sensed datasets is also a limitation–errors and uncertainties in remotely sensed data can be averaged out over larger regions, which may explain the strong remote sensing results at the Central Valley Scale described in Ahamed et al. (2021), but uncertainties and errors tend to persist at smaller scales (e.g., [69]), which may explain the poor results for the Butte Subbasin. Lastly, the application of the remote sensing water balance approach is not advised for regions where net groundwater inflow/outflow is a relatively large component of the overall water balance compared to other components. Regions with significant groundwater fluxes likely require advanced modeling approaches as well as additional measurements in order to accurately quantify groundwater flux and estimate ∆S_gw.

5. Conclusions

This study assessed the potential of a remote sensing water balance approach to accurately recover groundwater storage changes, ∆S_gw, at local scales (~2700 km² to ~3500 km²) highly relevant for the management of groundwater resources. The method was applied to the Butte and Kaweah–Tule Subbasins, which both fall within California’s Central Valley, but have substantially different hydrologic conditions in terms of surface water dynamics, available reservoir storage, and groundwater loss over the last 20 years. The remote sensing method was compared with estimates of ∆S_gw independently derived from groundwater wells and a widely-used groundwater flow model. The results showed a strong agreement between remotely sensed ∆S_gw estimates in the Kaweah–Tule region, but a limited agreement and higher seasonal variability within the Butte Subbasin. Part of the explanation for unrealistically high seasonal variability is the requirement of accurate inflow and outflow estimates for a study region, a data gap which SWOT may help to address. Future work will apply the analysis to basins with natural conditions, which will allow for an assessment of the influence of water management practices such as irrigation and water storage and distribution in reservoirs and canals. While more work is needed to understand the limitations and potential application of remotely sensed datasets to recover the changes in groundwater storage at local scales, the results suggest that the remote sensing approach can provide accurate estimates in some regions with the added benefits of low latency, rapid computation, and scale independence. These positive attributes can improve our understanding of aquifer dynamics, and facilitate up-to-date monitoring and adaptive management of groundwater resources.