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

: 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 quantiﬁed ∆ 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 2 –3500 km 2 ) 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 ﬁdelity 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 ﬂow 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.


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 3 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: 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 2 -3500 km 2 ). 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.

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 2 . 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 2 . 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 3 , respectively, when compared to reservoirs in the northern Central Valley.

Figure 1.
Map showing study areas. The Kaweah-Tule Subbasin and the Butte Subbasin (black outlines), Subregions 18 and 5, respectively, within the C2VSim groundwater model, are outlined in black (other subregions are outlined in red). Locations of reservoirs (cyan circles), and stream gauges measuring inflow (blue circles) and outflow (purple circle) for each study region are depicted. The Feather River, the source of Lake Oroville, is highlighted in blue. The Friant-Kern Canal, which transports water from Millerton Lake to the Kaweah-Tule Subbasin, is highlighted in gray. Locations of reservoirs (cyan circles), and stream gauges measuring inflow (blue circles) and outflow (purple circle) for each study region are depicted. The Feather River, the source of Lake Oroville, is highlighted in blue. The Friant-Kern Canal, which transports water from Millerton Lake to the Kaweah-Tule Subbasin, is highlighted in gray.
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 3 ) within and delivered throughout the region.

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.

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 ∆S KT gw , as well as the Butte Subbasin, referred to as ∆S B gw . 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: 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: 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 managementrelevant 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].

Data Used in Water Balance
The sources of data used in this study are shown in Table 1, 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. 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

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.

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 groundbased 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).

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.

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.

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.

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: where x i and y i are the elements of each variable, and − x and − 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 (σ 2 ) associated with each dataset. For three area-averaged collocated estimates of the same variable, the random error variances (σ 2 ) are given by: 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 σ 2 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.

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.
(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.

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.

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)

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 3 /yr) since 1960, while Butte has experienced a minimal loss (~0 km 3 /yr) of groundwater [65].

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 regionspecific 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 2 ), which includes the Sierra Nevada mountains, while accurate estimates of ET were essential to resolve ∆S gw for the Central Valley (~55,000 km 2 ), 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

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.

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 Figures 2 and 3 and for Butte in Figures 4 and 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.

Precipitation (P)
Generally, the estimates of precipitation within the Kaweah-Tule Subbasin (Figures 2A and 3A) and Butte Subbasin (Figures 4A and 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 Figures 2A and 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.

Evapotranspiration (ET)
Among input water balance components, evapotranspiration contains the highest variance between datasets.  Figures 2 and 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 Figures 2A and 4A

Precipitation (P)
Generally, the estimates of precipitation within the Kaweah-Tule Subbasin (Figures  2A and 3A) and Butte Subbasin (Figures 4A and 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). 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.

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 Figures 2C and 3C and for the Butte Subbasin in Figures 4C and 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 Figures 2C and 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 Figures 2-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]).

Runoff (R), Discharge (Q in and Q out ), and Reservoir Storage (∆S R )
A comparison of Figures 2D and 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 3 , and therefore constitutes a critical component of regional surface water dynamics in and around the Butte Subbasin.

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

Comparison and Evaluation between Estimates of ∆Sgw Derived from Remote Sensing, Wells, and Groundwater Model
We compared monthly ∆Sgw estimates during 2002-2021 for the remote sensing scenarios (described in Section 3.4) against the independent ∆Sgw estimates derived for well data (Section 3.5.1) and C2VSim (Section 3.5.2), available from 2002-2019 and 2002-2015, respectively. The ∆Sgw 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

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) 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 3 /year; shown in Figures 4 and 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 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.

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  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.

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 Sections 4.2.1 and 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 .

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 2 to~3500 km 2 ) 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.  Data Availability Statement: All data used in this study is from publicly available sources. The code written to query the data and perform the analysis can be found at https://github.com/ kashingtonDC/kaweah_wb (accessed on 21 March 2023). and https://github.com/kashingtonDC/ butte_wb (accessed on 21 March 2023).