Integrating shallow head measurements and InSAR data to quantify groundwater-storage change in San Joaquin Valley, California (USA)

Monitoring groundwater storage is essential for sustainable groundwater management. Storage can be quantified by considering the two main components through which storage change is expressed: saturation changes and deformation of aquifer materials. Here, these components were quantified using a selected area in California’s San Joaquin Valley (USA). First, this involved following existing observational approaches: quantifying the component expressed through saturation changes by identifying head measurements from shallow wells and scaling by specific yield. In the San Joaquin Valley, existing approaches to estimate the deformation component are to ignore it or approximate it with a simple linear relation to measured head. However, head and deformation measurements made at extensometers revealed that assuming a linear relationship between deformation and head might provide a poor estimate, particularly during periods in which measured head is rising. Instead, InSAR-derived surface deformation measurements were used to quantify the deformation component of storage changes. This showed that the two components—saturation and deformation—accounted for storage declines of equal magnitude over 2015–2021, suggesting that the deformation component should not be neglected when estimating storage changes in regions with subsidence. Summing the two calculated components gave a new estimate of the total storage change that captured the major trends seen in independent estimates, while better accounting for the deformation component. An additional benefit is that this method accounts for the deformation component in the unconfined aquifer. This method to quantify total storage change can be a practical and effective tool to support groundwater management.


Introduction
Projected population increases coupled with the impacts of climate change will drive stress on global groundwater supplies, which currently support ~20% of domestic and 40% of irrigation water needs (Siebert et al. 2010; UNESCO World Water Assessment Programme 2022).One key aspect of groundwater management is the monitoring of the volume of stored groundwater, commonly referred to as storage.Although there is a large array of monitoring methods ranging from direct measurements in wells to satellite remote sensing, it is generally challenging to monitor storage as there are seldom comprehensive in situ observations of key parameters.An improved ability to monitor groundwater storage would be a valuable contribution to the efforts of water managers and legislators to successfully implement sustainable groundwater management.
The San Joaquin Valley in California (USA) is home to over 4 million people and has an annual agricultural output valued over $35 billion (Fernandez-Bou et al. 2021).Having a semiarid climate and regular droughts, the region is dependent upon groundwater for a large portion of its water supply (California Department of Water Resources 2020a).There have been significant efforts in recent years to quantify changes in storage at the Valley scale, leading to the identification of major storage declines during the 2007-2009 and 2012-2016 droughts from which there was limited recovery during wetter years (Ahamed et al. 2022;Alam et al. 2021;Xiao et al. 2017).These studies highlight the vulnerability of the groundwater supply.Recognizing the importance of groundwater resources to state water use, California passed legislation in 2014 known as the Sustainable Groundwater Management Act (SGMA), which requires "significant and unreasonable declines" in groundwater storage to be identified and, by 2042, stopped.Fulfilling the requirements of the SGMA falls on groundwater sustainability agencies (GSAs), which are small, local groundwater management units.As a part of this legislation, annual reports containing estimates of storage changes must be provided by each GSA.There is thus a need for reliable and accurate estimates of groundwater storage change, both at a regional (Valley-wide) and local (GSA) scale.
A commonly used approach to quantify storage change in the San Joaquin Valley equates storage change with the product of shallow, assumed unconfined, head change and specific yield (Brewster et al. 2015;Harder 2022;North Kings Groundwater Sustainability Agency 2019;Scanlon et al. 2012).This approach effectively assumes that storage change is fully expressed through changes in saturation in the unconfined aquifer.
There is, however, a second component through which storage change can be expressed, which previous studies have shown can be important in the San Joaquin Valley (Lofgren 1978;Williamson et al. 1989).This is the deformation of aquifer materials, which occurs in response to a change in the effective stress within the aquifer.It is typically assumed that this deformation only occurs in confined aquifers, and that the magnitude of the associated storage change can be calculated as the product of head change, thickness and specific storage of the confined aquifer.
An approach to estimating storage that would appear to account for both the saturation and the deformation components is to sum the product of head and specific yield for the unconfined aquifer and the product of head, specific storage and aquifer thickness for the confined aquifer.Since the third storage component, the deformation of water itself, is typically very small in the San Joaquin Valley (Williamson et al. 1989), the approach would serve to quantify the total storage change.However, a lack of information about the depth from which head measurements are taken in the San Joaquin Valley makes it impossible to separate those made in the unconfined aquifers from the confined.In addition, abundant clay interbeds in the aquifer system create semiconfined conditions in many parts of the aquifer (Planert and Williams 1995), and there is no clear theoretical framework for using head to quantify storage change in semiconfined aquifers.
As a simplified approach to account for both storage components, a number of regional-scale studies have linearly related all head measurements (with no attempt to accurately determine the depth or level of confinement) to storage change by defining an effective storativity through calibration with GRACE satellite data (Alam et al. 2021;Rateb et al. 2020).Similarly, at the local scale, several GSAs have implicitly assumed that head measurements are linearly related to storage change by defining "significant and unreasonable declines" in storage as occurring when head measurements drop below specified thresholds at individual wells (Eastern Tule Groundwater Sustainability Agency et al. 2020;Greater Kaweah Groundwater Sustainability Agency et al. 2022).
The focus of this study was the development of a method for quantifying storage change with an improved level of accuracy.For the component expressed through saturation, the approach described above was followed, equating storage change with the product of shallow head change and specific yield, acknowledging that there are two primary sources of uncertainty-the lack of accurate information about the depth from which head measurements are taken, i.e. it cannot be said with certainty whether that the measurements are from an unconfined aquifer; and the value used for specific yield.
For the component expressed through deformation, it was first explored whether, in clay-rich confined aquifers, the measured head is a good predictor of deformation.There is a complicated relationship between measured headtypically representative of head in coarse-grained materials-and head in the clay interbeds, where the majority of the deformation occurs.In certain aquifer systems, where deformation is primarily elastic, this relationship can be approximated as linear (Chen et al. 2017;Reeves et al. 2011).However, in the San Joaquin Valley, the numerous, thick interbeds combined with the prolonged history of head decline mean there are time lags between the head change in coarse-grained materials and that in the clays.Here, without long records of head and detailed hydrostratigraphy that permit modelling of the head in the clays (as was done in the recent study by Lees et al. 2022), it is unclear whether head measurements can accurately quantify the deformation.To circumvent the difficulties of using head to quantify the deformation component, interferometric synthetic aperture radar (InSAR) data were used.These data measure surface deformation which, under certain circumstances, can be used as a measure of the component of storage change expressed through the deformation of aquifer materials.
An estimate of the total storage change over time was obtained by combining the two components-storage changes due to saturation changes estimated with shallow head measurements and storage changes due to deformation estimated using InSAR.This approach, which was recently described by several GSAs (Mid-Kings River Groundwater Sustainability Agency et al. 2022;Provost and Pritchard Consulting Group 2022), has the potential to provide a useful method for monitoring total storage change.In this study, the supporting theoretical framework that justifies taking this approach was developed.The estimates, derived using this approach, were then compared to two independent estimates, one obtained from the GSAs and the other from a groundwater model, to assess the extent to which the estimates differ and to identify possible benefits.
The focus of this study was the Tule and Kaweah subbasins of the San Joaquin Valley as a study area containing high-quality head and aquifer deformation data.In particular, data were used from three sites containing extensometers, down-well instruments that record aquifer deformation, with colocated head observations.Additionally, InSAR measurements of surface deformation were available across the study area from the Sentinel-1 satellites, spanning 2015-2022, and databases of contemporaneous head measurements, including those that are specifically from unconfined/shallow parts of the aquifer system, were also available.
Although this study is specific to the Tule and Kaweah subbasins of the San Joaquin Valley, the difficulties faced using head measurements to monitor storage changes are common in many clay-rich aquifer systems.The findings of this study are expected to have transferable applications, as well as immediate groundwater management implications in the San Joaquin Valley.
The organization of this paper is as follows.First, general information about the study area and the theory of storage change is provided.Then, three main sections follow, each containing data, methods, results and discussion; this organization was chosen because the second and third sections, 'The deformation component of storage' and 'Estimating the total storage change', rely upon results from the preceding section(s).In the first of these three sections, the component of storage change expressed through saturation change is estimated.In the second, the extensometer sites are used to assess whether head is an appropriate estimator of the deformation component and InSAR is used to provide an improved estimate.In the section 'Estimating the saturation component of storage change', the storage components from the first two parts are summed and the total is compared to independent storage estimates.Through these three sections, an approach is demonstrated that can provide improved estimates of storage change, from the local to regional scale, to support the implementation of the SGMA.

Study area
The study area consists of the Kaweah and Tule subbasins, located within the south-eastern part of the San Joaquin Valley, as shown in Fig. 1a.Typical of the Valley, these subbasins host a multibillion dollar (USD) agricultural industry (California Department of Food and Agriculture 2022), receive relatively low annual precipitation of between 20-25 cm (California Department of Water Resources 2003), and are highly dependent on surfacewater imports and groundwater (providing ~30 and 70% of water use respectively) (California Department of Water Resources 2020a).Severe hydrologic conditions occurred during 2015-2021, with an extreme drought from 2015-2016, wet-to-average years from 2017-2019, and a return to severe drought conditions from 2020-2021 (California Department of Water Resources 2020b).These recent droughts have led to lowering of the water table and major land subsidence in the study area, which follows a previous, historic period of subsidence between 1920-1970(Poland et al. 1975)).The spatial extent of the recent and historic subsidence, which reaches beyond the study area and across much of the valley, is shown in Fig. 1b.
The study area is organized into 10 GSAs, with three in the Kaweah subbasin and seven in the Tule subbasin (Fig. 1c,d).Although each GSA is required individually to describe how storage will be measured and how "significant and unreasonable" declines will be defined, there is substantial coordination between GSAs within each subbasin.The Kaweah subbasin has a coordination agreement between all three GSAs, which sets out a shared approach defining an "undesirable result" as a situation where at least 33% of monitoring wells across the three GSAs measure heads below a given threshold.In the Tule subbasin, there is also a coordination agreement between the seven GSAs which sets out a similar, but slightly more lenient, criterion.In their agreement, they define an "undesirable result" to have occurred if 50% of monitoring wells across the GSAs measure heads below a specified threshold for two consecutive years.
The key hydrogeologic features of the study area are shown in Fig. 2, and are described below.The aquifer system, which ranges from ~10-400 m in thickness, consists primarily of unconsolidated alluvial and fluvial sediments (Petersen and Nicely 2020), with a high proportion of clay, estimated by Faunt et al. (2010) to be ~40% by volume.This clay is mostly concentrated in clay interbeds: these are lenticular regions of clay which, while sometimes 10s of meters thick, do not occur as laterally continuous units.The eastern side of the study area abuts the granitic bedrock of the Sierra Nevada; along the western side lie the Kings and Tulare Lake subbasins.In the western part of the study area, there is a laterally continuous clay, the Corcoran Clay, which acts as a confining layer.Underneath the Corcoran Clay, the aquifer system is described as confined, and the terminology "upper aquifer" is used to describe the region above the Corcoran Clay.In the upper aquifer, as well as in the aquifer system east of the lateral extent of the Corcoran Clay, the aquifer system is typically described as an unconfined layer that grades to semiconfined with increasing depth, although there is no consensus on how to define the depths at which this gradation occurs (Petersen and Nicely 2020;Thomas Harder and Co. 2020).Groundwater is pumped from unconfined, semiconfined and confined layers across the study area, apart from a portion in the south-eastern corner where the unconfined aquifer is considered dry and is not used (Harder and Lewis 2020).
Due to the heterogeneity of the vertical hydrogeological structure of the San Joaquin Valley, the depths of the screened interval for wells are crucial for understanding which part of the aquifer system is being sampled by head measurements.While there are many thousands of groundwater wells in the study area, only a small number report the screened interval.A histogram of the reported screened intervals is provided in Fig. 3.While some wells are screened over small <20-m-long intervals, others are screened over 100 m or longer intervals, and some, screened over intervals of 300 m or greater, tend to be screened across the entire thickness of the aquifer system.
One of the advantages of working in the Kaweah and Tule subbasins is the ability to build on previous investigations (Kang et al. 2022;Knight et al. 2018;Lees et al. 2022;Smith and Knight. 2019) and continue collaboration with the water agencies which can provide valuable knowledge of local conditions and access to data.Although there are differences, the generalized hydrogeology of the study area has many commonalities with the rest of the San Joaquin Valley, including the presence of unconfined, semiconfined and confined layers.Therefore, it is expected that this study will have findings which are transferable throughout the Valley, in addition to hydrogeologically similar aquifer systems globally.

Storage change equations
At any given scale, the total storage change ΔS tot in an aquifer system is the sum of the storage changes expressed through three components: saturation (ΔS sat ) deformation of aquifer materials (ΔS def ) and deformation of water itself (ΔS w ; Weight 2019) and can be written To quantify these three terms, it is common to relate them to hydraulic head using storage coefficients, which takes a different form for unconfined and confined aquifers.In unconfined aquifers, in most practical situations, the storage change expressed through saturation change is far larger than the change expressed through the other two components (Heath 1983).The storage change in an unconfined aquifer is therefore written as: where ΔS′ unc is the storage change per unit area of the unconfined aquifer, S y is the specific yield and Δh meas is the measured head value (Fetter 2007).This relationship has   three implicit simplifications.First, it assumes that negligible deformation of water or aquifer materials occurs.In the San Joaquin Valley, it is generally thought that the deformation of water is a very minor part of storage changes, and although the deformation of aquifer materials can account for a significant portion of total storage changes (Williamson et al. 1989), this deformation primarily occurs in deeper, confined parts of the aquifer system (Ireland et al. 1984).
Second, it assumes that the measured head in the unconfined aquifer is representative of the water-table elevation, an assumption which can be difficult to validate given the poor quality of depth information associated with head measurements in the San Joaquin Valley.Third, it assumes that an appropriate value of specific yield can be assigned; this is usually done through either field measurements or using the results of calibrated numerical models.
It is common to treat confined aquifers as remaining fully saturated at all times.Hence, confined aquifer storage change is the sum of the aquifer material deformation component and the deformation of water component (Lohman 1972).Under the assumption of a homogeneous aquifer of thickness b, these components are typically expressed in terms of head as where ΔS′ conf is the storage change per unit area of confined aquifer, S sk is the skeletal specific storage coefficient and S w is the specific storage of water.The high clay content in the San Joaquin Valley aquifer system means that the component of storage change expressed through the deformation of aquifer materials is far larger than the component expressed through the deformation of water (S sk b ≫ S w b) (Sneed 2001).This is particularly the case when the head level declines below the historic minimum, as observed in many parts of the Valley following recent droughts, and deformation occurs inelastically, which is characterized by large compaction and permanent deformation (Faunt et al. 2016).Thus, the component of storage changes expressed through the deformation of water can be neglected and storage change in a confined aquifer can be simply written as Contained within Eq. ( 4) is the prediction that measured head change is linearly proportional to the component of storage change expressed through the deformation of aquifer materials, which is explicitly written in Eq. ( 5): Equation ( 5), as originally derived by Jacob (1940), only applied to situations where "there are a sufficient number of clay-laminae interbedded with the sand so that the release (5) ΔS � def = S sk bΔh meas .
of stored water from the clay is virtually instantaneous".If there are thicker clay interbeds, this leads to time lags between measured head and the storage change expressed through deformation, and Jacob indicated that Eq. ( 5) should not apply in such circumstances.More recently, this condition has often been omitted; for instance, Freeze and Cherry (1979) state that the equation can be "specified for aquitards as well as aquifers", but do not comment on the application to an aquifer with interbedded clays, while the textbook of Fetter ( 2007) makes no mention of interbedded clays when describing Eq. ( 5).In the San Joaquin Valley, Eq. ( 5) has been applied, despite evidence for clay interbeds and long time lags, estimated in the study area to be decades-to-centuries by Ireland et al. (1984) and Lees et al. (2022).For groundwater management, Eq. ( 5) has been used directly by GSAs to quantify changes in confined aquifer storage using Eq. ( 5) (e.g.(2022).In the published scientific literature, it has been used implicitly in those regional-scale studies that assume a linear coefficient can be used to relate all head measurements and total storage (Ahamed et al. 2022;Alam et al. 2021).A goal of this study was to test the accuracy of using Eq. ( 5) with head measurements in the San Joaquin Valley to quantify storage change.
In Eqs. ( 4) and ( 5), it is common to separate the skeletal specific storage coefficient into its elastic and inelastic components, S ske and S skv (Ojha et al. 2019;Smith et al. 2017).The elastic component refers to the portion of deformation which is reversible, which is characterized by a relatively small coefficient S ske .The inelastic component, which is assumed to only apply when head levels fall below historic minima, refers to the permanent, irreversible component of deformation and the coefficient is much larger (Sneed 2001).It can be written that: In a semiconfined aquifer, there is no established mathematical representation linking a head change to a storage change.It is expected that storage change in the semiconfined aquifers of the San Joaquin Valley will be expressed partially through saturation changes-from exchange with the water table-and partially through deformation of aquifer materials.The partitioning between the two will depend on (1) the degree of confinement and (2) the timescale of interest; however, no detailed study has explicitly investigated this.On this basis, and in response to the situation where head measurements cannot be separated into unconfined/semiconfined/confined aquifers, it has been proposed (6) S sk = S ske + S skv to use the following approximation for total storage change (Alam et al. 2021;Rateb et al. 2020): where ΔS′ tot is the total storage change per unit area across the entire aquifer system depth (including unconfined, semiconfined and confined layers), S e is so-called effective storativity, assumed to take a value intermediate to that of S y and bS s , while Δh meas is a spatial/temporal average of head measurements in the area where measurements generally include those from any/all of confined, unconfined and semiconfined aquifers.To use Eq. ( 7), S e is determined by calibration with a storage change calculated using independent methods, most commonly GRACE satellites (Alam et al. 2021;Rateb et al. 2020).Equation ( 6) has been applied in the San Joaquin Valley, and its use contains within it the implicit assumption that the deformation component of storage has a linear correlation with head measurements (i.e., the relationship given in Eq. 5).

Theory of head measurements and interbedded clays
In the San Joaquin Valley, wells are commonly screened across multiple layers each of which may have different head values.In particular, the heads in the clay interbeds and the Corcoran Clay are presumed to be substantially different from those in the coarse-grained material due to long timescales of equilibration between the two (Ireland et al. 1984;Lees et al. 2022).Three head values of importance are defined: head in the coarse-grained layers, head in clay layers (interbeds or the Corcoran Clay), and the head value measured in wells.The head in the clay layers is closely related to the deformation component of storage change since the majority of deformation occurs in clays, whereas the head measured in wells is the value observed and used to estimate storage changes.It is therefore important to consider how the head measured in wells relates to the head in the clay layers.
If the head measured in a well is allowed to equilibrate with the heads in the layers across which it is screened, so it no longer varies with time, then this 'steady-state' head measured in the well is the transmissivity-weighted average of the head in each of the layers (Sokol 1963).It follows that the measured head is insensitive to the head in the clay layers, as their transmissivity is orders of magnitude lower than that of the coarse-grained materials.In the San Joaquin Valley, there is frequent pumping which disturbs the establishment of steady state.In this situation, the head measured in the well will change over time, but it is still expected that measurements made will be more representative of the heads in the coarse-grained materials than in the clays, since (7) ΔS � tot ≈ S e Δh meas the rate of flow from clays into the well bore will be slow.In summary, head measurements made in wells in the San Joaquin Valley do not provide a good representation of head conditions in clay layers, and are far more indicative of the head in coarse-grained units.

Direct measurement to quantify the deformation component of storage
Although it is most common to use head measurements and a storage coefficient to quantify the component of storage changes expressed through deformation, multiple studies have used measurements of deformation within an aquifer system to represent the deformation component directly (Jiang et al. 2018;Konikow and Neuzil 2007;Rezaei et al. 2020).These direct measurements offer estimates of this component of storage change independent of head measurements.Mathematically, the relationship between the component of a storage change expressed through deformation and the deformation (the change in thickness, Δb) across some interval z 1 to z 2 can be written as where Δϕ(z) is the change in porosity at a given depth.Note that, here, recent convention is followed-e.g.Reeves et al. (2011) and Smith et al. (2017)-in defining Δb such that a positive Δb implies expansion (increasing aquifer thickness) and a negative Δb implies compaction (decreasing aquifer thickness).In using Eq. ( 8), it is assumed that all storage changes occur in the vertical dimension.

Measurements of aquifer deformation
In this study, deformation measurements provided by extensometers and InSAR are used.A brief introduction to both is provided in the following.Extensometers are instruments installed in wells or boreholes to measure compaction over an interval of the subsurface.The extensometers used in this study measure the changes in tension on a cable running down a well and anchored in the aquifer materials at the base, allowing measurement of deformation with sub-cm or even sub-mm accuracy (Gambolati and Teatini 2021).Detailed information on the design and operation of San Joaquin Valley extensometers is provided in Poland et al. (1984) and Riley (1986).Extensometers are particularly useful for this study as, where head is also measured within the borehole, they can provide an observational record of colocated head and deformation across a known depth interval.
Since 2014, many studies have reported the use of InSAR to study subsidence in the San Joaquin Valley (Chen et al. 2015;Farr and Liu 2014;Ojha et al. 2018Ojha et al. , 2019;;Smith et al. The accuracy of measurements from the first studies, using data from older satellites, was impacted by inconsistent satellite geometry, signal decorrelation, data gaps and a lack of independent deformation estimates for calibration and/or validation.However, the launch of the Sentinel-1A and 1B satellites in 2014 has provided imagery of the valley with high temporal resolution and a reliable imaging geometry (Jones et al. 2021), which has been combined with a growing network of global navigation satellite system (GNSS) sensors in the Valley allowing for superior calibration and validation of InSAR data products (e.g., Neely et al. 2020).InSAR measures the surface deformation, which is the sum of all deformation, including that due to storage changes, tectonics, and other causes.For this study, InSAR is most beneficial due to its near-continuous spatial coverage, which other methods cannot provide.

Saturation component: methods and data
To estimate the storage change expressed through saturation, a dataset of shallow head measurements was assembled by combining data from two sources.The first source, the Department of Water Resources (DWR)'s Seasonal Groundwater Reports Well Elevations dataset (California Department of Water Resources 2023a), is a subset of the DWR's Periodic Groundwater Levels dataset and is produced by filtering the latter to only include measurements representative of the unconfined to shallow semiconfined aquifer system.The filtering process is partially documented in Brewster et al. (2015) and via the SGMA Data Viewer Docs (California Department of Water Resources 2023c).In summary, the data are first separated into spring and fall measurements, defined as January 1st to May 31st and September 1st to November 30th, respectively.Then, two filtering steps aim to extract only shallow measurements.In the first step, wells with a total depth less than 99 ft (~30 m) are removed, to avoid unreliable data from perched water tables, as well as those with a top perforation deeper than 400 ft (~120 m), which are considered too deep to be unconfined.Where well depths and/or perforations are not known, the wells are retained.In the second step, a triangulated irregular network surface is constructed for each season and anomalous measurements are individually identified through inspection and removed (B.Brewster, personal communication, 2022).The dataset is complete up to 2021.For the years and location of the study, 2015-2021, this dataset contains 385 wells with a total of 1,693 measurements.The well locations are shown in Fig. 4a.
These measurements were supplemented with a second head dataset which was constructed from shallow measurements made by local agencies.This second dataset was added since some agency measurements are not included in the Seasonal Groundwater Reports Well Elevations dataset, and improved spatial and temporal coverage was desired.The GSAs in the study area operate monitoring networks, and these networks include well measurements identified as representing shallow conditions, labelled "Upper Aquifer System" or "Assumed Unconfined" in the Kaweah subbasin and "Upper Aquifer" in the Tule subbasin.These data are available through annual reports in the two subbasins (Harder 2020(Harder , 2021(Harder , 2022;;Provost andPritchard Consulting Group 2020, 2022) as well as through the SGMA monitoring sites portal (California Department of Water Resources 2023d).The shallow measurements were manually extracted from both sources and, after removal of measurements which were already included in the first dataset, an additional 26 wells with 71 measurements were found in the Tule subbasin and 66 wells with 270 measurements were found in the Kaweah subbasin within the time period of 2015-2021.The locations of these wells are also shown in Fig. 4b.
The measurements from the new, combined dataset of shallow head measurements were interpolated onto uniform grids with 1-km spacing in latitude and longitude.This was done for spring and fall of each year from 2015-2021, the period for which the head datasets were complete, using the definitions of spring and fall as in the Seasonal Groundwater Reports Well Elevations dataset.To interpolate, the arithmetic mean was taken of measurements within a certain search radius, using the smallest search radius that still resulted in 90% of the grid cells averaging two or more measurements.When interpolating using a search radius, any grid cells that did not have a measurement within the search radius were assigned the head value of the nearest grid cell.
To estimate the saturation storage changes, each interpolated grid of head was differenced from the interpolated grid of spring 2015 (the first measurement date) and the resulting head changes were scaled by a spatially constant specific yield value of 9% from Bertoldi et al. (1991).Although specific yield has been estimated in multiple ways, the estimation described in Bertoldi et al. (1991) was selected as it was based on a compilation of laboratory measurements made on samples from the Central Valley, rather than other estimates that may rely on modelling.The next step was integration across the study area to obtain the time series of storage change.The integration was done using the "grdvolume" function from the Generic Mapping Toolbox (Wessel et al. 2019).To estimate the size of uncertainty introduced by error in the measurements or biases in spatial coverage, a bootstrapping approach was used, in which the storage calculation was repeated 1,000 times, withholding between 1-10% of 1 3 the head measurements each time.The interquartile range of the bootstrapping results was used as an estimate of the error.

Results and discussion on the saturation component
The resulting estimate of the component of storage expressed through saturation changes is shown in Fig. 5.The 14 interpolated grids of shallow head measurements, one each for spring and fall from 2015-2021, used to estimate the saturation component are shown in Figure S1 in the electronic supplementary material (ESM1) and the search radii are reported in Table S1 in ESM1.Over the entire period 2015-2021, the saturation component contributed a net, long-term storage decline of ~2-3 km 3 , with large seasonal changes between spring and fall of each year.

The deformation component of storage
The analysis of the deformation component of storage change contains two parts.The first part uses colocated head and deformation measurements at extensometer sites in confined and semiconfined aquifers to explore whether head is a good estimator of the deformation component of storage.The second part demonstrates the method of using InSAR to quantify the component.

Methods and data for colocated head and deformation measurements
The deformation measurements used were made by extensometers installed in boreholes where head measurements  1998).Of these 11 instruments, those where a long record of deformation was accompanied by water level measurements were desired.An additional requirement was for both the deformation and head measurements to be deep, in parts of the aquifer system which were confined or very semiconfined, since head in such aquifers is related to the storage change expressed through deformation (Eqs.4 and 5) as opposed to that in unconfined aquifers where head is primarily related to the saturation component of storage change (Eq.2).Four extensometers covering three sites satisfied these criteria, and information about these extensometers' period of record, depth and the screened intervals for head measurements are detailed in Table 1.Following Ireland et al. (1984), the Pixley Upper and Lower extensometer data were combined, by differencing the two, to obtain the resulting deformation over the narrow, confined interval 130-230 m.The locations of the extensometers are shown in Fig. 6 and the data sources for each are provided in Table 1.
For the three extensometer sites, the head measurements were made across a subsection of the interval for which deformation was known, but they were always from depths where the aquifer system is semiconfined or confined.These head and deformation measurements were used to test, at the local scale, the applicability of Eq. ( 5), which describes a positive linear relationship between head measurements and the storage component expressed through deformation.Equation ( 8) states that a measured/known deformation across an interval is, by definition, equal to the component of a storage change across that interval expressed through deformation.Deformation measurements were therefore  treated as a direct measure of the deformation component of storage changes.
It is noted that the head and deformation measurements cover different (albeit similar) depth intervals.This means that the analysis is not a perfect test of Eq. ( 5).However, the analysis mirrors the practice in which measurements of head made in a well (with some associated screened interval) are used to infer information about the storage change due to deformation.Hence, the analysis is applicable to practical situations in groundwater management where head measurements are used to make inferences about storage changes.

Results and discussion on colocated head and deformation measurements
The records of head and deformation are shown in Fig. 7, which contains three plots showing cumulative deformation on the y-axis and measured head on the x-axis.Figure 7 also shows the depth intervals over which deformation and head were measured.To match the annual storage-reporting requirement of SGMA, a 365-day moving average was applied to both the head and deformation data prior to plotting.Note that an increase in cumulative deformation (i.e., In each panel, time is not explicitly plotted, but always proceeds from top to bottom, apart from where indicated by arrows (c).In all panels, certain years are labelled to aid interpretation moving upwards in Fig. 7) corresponds to expansion of the aquifer materials and a decrease in cumulative deformation (moving downwards on Fig. 7) corresponds to compaction.Note also that time, which is not plotted directly, proceeds from top to bottom in Fig. 7, except where noted in Fig. 7c.
In each plot, time periods were identified by eye within which two types of relationship between head and cumulative deformation were observed: 1 Type 1 is behavior that is in agreement with Eq. ( 5) and the general understanding of the relationship between head and the deformation component of storage.Two styles of type 1 behavior were seen: a. Type 1a: This is when cumulative deformation is decreasing (i.e.compaction is occurring) and head is decreasing.Type 1a behavior was very common, and some examples include 1960-1962and 1974-1977at Pixley, 2019-2022at Deer Creek or 1959-1962 at Friant Kern 1. b. Type 1b: This is when cumulative deformation is increasing (i.e., expansion is occurring) and head is rising.Type 1b behavior was much less common than type 1a, and was only seen at Friant Kern 1for instance, during 1977-1980.
To highlight the correspondence with Eq. ( 5), linear gradients for each example of type 1 behavior have been approximated-these are the dashed lines in Fig. 7, where the colors distinguish type 1a from type 1b.These lines highlight how, as well as the general relationship in Eq. ( 5) holding (i.e., deformation and head change in the same direction), the specific relationship-that of linearity between the two-also holds, at least over time periods of several years. 2 Type 2 is behavior that is not in agreement with Eq. ( 5), running contrary to the general understanding of the relationship between head and the deformation component of storage.Only one form of type 2 behavior was seen: it was always characterized by increasing head and cumulative deformation decrease (i.e., compaction).It was common at Pixley (e.g. 1962Pixley (e.g. -1964) ) but also occurred at Deer Creek (e.g. 1977Creek (e.g. -1980) ) and Friant Kern 1 (e.g. three episodes between 1962-1965).In Fig. 7, type 2 behavior is indicated with red shading.Note that the type 2 behavior typically has a shallower slope than type 1a, which represents less compaction per change in head than the type 2 behavior.
First, a physical explanation for type 2 behavior is considered: why does Eq.( 5) frequently not hold?This is explained in the context of a phenomenon widely reported in the literature: residual compaction.Residual compaction occurs when head decline in the clays lags behind that in the coarsegrained materials, which means drainage and compaction of clays can continue after a fall in measured head and even during rising head as sampled through measurements (e.g., Galloway et al. 1999;Lees et al. 2022).It is suggested that during periods of rising measured head, although the head is rising and expansion is occurring in the coarse-grained materials, the compaction caused by delayed head decline in clays dominates so that the storage component expressed through the deformation of aquifer materials is one of storage loss.This can explain why Eq. ( 5) does not hold in the case of rising head, which was the time type 2 behavior was identified.It also explains why the slope of cumulative deformation during type 2 behavior is shallower than during type 1a, because compaction of clays during type 2 behavior is occurring alongside expansion in coarse-grained materials, which lessens the total compaction.It is noted that, were head measurements available within the clays, it would be expected that Eq. ( 5) would hold; it is only due to the fact that available head measurements sample the coarse-grained materials, and not the clays, that times are identified during which Eq. ( 5) does not hold.
Next, the practical implications of the observed type 1 (a or b) and type 2 behavior for monitoring storage are considered.To that aim, two scenarios are considered: periods where measured head was decreasing, and periods where it was increasing.In the former case, only type 1a behavior was seen to be consistent with Eq. ( 5); hence, it is concluded that a reasonable estimate of the storage change component expressed through deformation can be made by assuming a linear relationship to measured head during measured head declines.However, in the latter case (measured head increasing), two different behaviors were seen: type 1b, which was consistent with Eq. ( 5), and type 2, which was not.Type 2 was the predominant behavior, which can be expected when residual compaction rates are large, and it is suggested that residual compaction is therefore generally high in the study area.The type 1b behavior at Friant Kern 1 may have been the result of the large and sustained measured head rise between 1962-1980 halting the residual compaction at that location.The predominance of type 2 behavior during rises in measured head is significant for storage estimation.If storage estimations were made using a linear relationship during such periods, the estimated storage change would be an increase whereas the true storage change would be a decrease.This suggests that storage estimations made using a linear approximation in periods where measured head is increasing are prone to underestimate storage losses or even mischaracterize them as storage recovery.Hence, a reliance on head measurements and Eq. ( 5) is considered to be an unreliable approach to estimate the component of storage expressed through deformation in the study area.This conclusion is expected to apply to other areas where there is substantial water coming from drainage of clay interbeds with associated large time lags, thus making invalid the use of Eq. ( 5), as derived by Jacob (1940).This idea is revisited when comparing the new storage estimates made here to existing estimates in section 'Estimating the total storage change'.

Methods and data for calculating the deformation component
The results of section 'Analysis of colocated head and deformation measurements' indicate that head measurements may not be a good predictor of the storage change component expressed through deformation.Instead of head measurements, surface deformation measurements derived from InSAR data were used.Starting with Eq. ( 8), reproduced below as Eq. ( 9) for ease of viewing, and setting the depth interval to be from the water table z wt to the base of the aquifer system z base , puts the equation into a form describing the deformation storage change for the entire aquifer system (Eq.10).It is assumed that this quantity could be approximated by the InSAR measurements of surface deformation Δs InSAR , which is also written in Eq. ( 9).
The assumption that the InSAR measurements correspond to ΔS def is dependent on three conditions: (1) that the deformation within the aquifer system propagates to the surface without attenuation, (2) that there are no tectonic or other nonhydrological causes of deformation, and (3) that there is no deformation occurring in the unsaturated zone (since that zone is not part of the aquifer system).In the study area, these three conditions are satisfied.The aquifer system is no more than 400 m thick, and in such shallow systems it is common to assume there is no attenuation of the deformation signal during propagation to the surface (Galloway et al. 1998;Reeves et al. 2011;Smith et al. 2021).Although there is tectonic subsidence in the San Joaquin Valley, due to deep lithospheric detachment, its rate is a few mm per year, far smaller than the deformation due to compaction of aquifer materials (Saleeby et al. 2013).Finally, while a form of shallow compaction known as hydrocompaction has been reported in the San Joaquin Valley (Lofgren 1969), there are no reports of shallow deformation in the study area, and the (9) unsaturated zone is thin compared to the thickness of the total aquifer system, so deformation in the unsaturated zone will not be significant.
It is noted that this method applies where it can be assumed that the majority of deformation is vertical, and storage changes are not expressed in the horizontal direction.It is believed that this is a valid assumption, as previous studies have noted that the predominant deformation component in the San Joaquin Valley is vertical (Carlson et al. 2020;Lees et al. 2022;Smith et al. 2017).
To obtain Δs InSAR , InSAR surface deformation measurements collected by the European Space Agency Sentinel-1 satellites and processed by TRE Altamira under contract to the DWR were used.The dataset provides vertical deformation every 5-7 days across California's groundwater basins spanning the period February 2015 to October 2021 and is available at the California Natural Resources Agency data repository (California Natural Resources Agency 2022).The methods for generation and validation of the vertical deformation dataset are provided in TRE Altamira (2021) and Towill Inc. ( 2021).The deformation product was used as provided, without altering the workflows described in the two reports.A summary of pertinent information from TRE Altamira (2021) and Towill Inc. ( 2021) is provided here.To obtain vertical deformation, TRE Altamira took several steps.First, the raw phase data were processed in both ascending and descending tracks using the SqueeSAR algorithm, obtaining line of sight (LOS) deformation estimates.For pixels where both LOS directions were available, the two LOS measurements were combined to obtain vertical displacement vectors.At points where only a single imaging geometry was available, or the gap between the ascending and descending measurement was too long, the LOS measurement was directly projected into vertical.The resulting vertical displacement dataset was calibrated against a subset of ~40 of the Valley's GNSS stations by planar removal and absolute calibration.The dataset was validated by Towill Inc. ( 2021), who used a further subset of GNSS stations including ~20 not used in the calibration procedure.The result of the validation exercise was an estimate of the InSAR error within a 95% confidence interval, based on a California-wide analysis, of 18 mm.This dataset was used as provided as it is freely distributed by the DWR and was produced using modern, rigorous processing techniques and subjected to an independent validation analysis using GNSS.
The SqueeSAR algorithm means that certain pixels do not have deformation measurements for every time channel; instead, the presence of a measurement depends on the phase stability of the pixel.In the dataset, this can be seen by the ~40,000 InSAR pixels with values reset to zero in 2020, due to there being a large gap in temporal coverage.For simplicity, all these pixels were removed for the analysis in this paper, leaving 86,000 pixels.The progression of subsidence recorded by these 86,000 pixels-in other words, the maps of measured deformation-is shown in the video in ESM2.
The InSAR measurements of vertical deformation are provided at a point/pixel-scale.To convert these into a study area-wide estimate of storage change, the deformation at each pixel was first averaged to obtain one deformation value for each spring and fall, consistent with the shallow head datasets, by taking the mean of the surface deformation measurements which fell within each season.Any pixels with fewer than five measurements for any given season were removed.Next, the InSAR was interpolated onto uniform grids with 1-km spacing in latitude and longitude for each season.To interpolate, the arithmetic mean of measurements within a certain search radius was taken, using the search radii as used for the head data, described in section 'Saturation component: methods and data'.These steps put the InSAR data into a form with the same spatial and temporal sampling as the head data.Lastly, the grdvolume function from the Generic Mapping Toolbox (Wessel et al. 2019) was used to calculate the volume of subsidence for each grid across the study area.This resulted in estimates of the storage change expressed through deformation for each spring and fall season between 2015-2021.The uncertainty of 18 mm that TRE Altamira determined from GPS calibration was used for the entire study area to estimate the error for this component.
An estimate of the elastic and inelastic parts of the deformation component was also made.This is of interest to water managers because inelastic deformation is nonrecoverable and represents a permanent loss in storage capacity.Various studies in the past have used simplified physical modelling (Smith et al. 2017), principal component analysis (Chaussard et al. 2014), wavelet transforms (Miller et al. 2017) and spectral techniques (Ojha et al. 2018(Ojha et al. , 2019) ) to separate the elastic and inelastic components.Here, a simple decomposition was used, similar to Chaussard et al. (2014), where an exponential decay was fit to the total InSAR-measured storage change to estimate the inelastic component, and the remainder was assumed to be the elastic component.The full description of the method used to separate the elastic and inelastic deformation is provided in the Appendix.

Results and discussion on calculating the deformation component
The calculated component of storage expressed through deformation, and the estimates of the elastic and inelastic parts of this component, are shown in Fig. 8.The 14 interpolated grids of InSAR deformation used to estimate the deformation component, one each for spring and fall from 2015-2021, are shown in Figure S3 in ESM1.Over the entire period 2015-2021, the deformation component exhibited a monotonic decline, contributing a net storage decrease of 2.3 km 3 .This decline was primarily through the long-term, inelastic part, with the elastic part contributing to seasonal oscillations of magnitude <0.2 km 3 year -1 .
In section 'Estimating the saturation component of storage change', it is seen that the saturation component of storage contributed a net storage decline of between 2-3 km 3 .The deformation and saturation components, therefore, contributed a similar decline.
The relative magnitude of the two components is consistent with the estimate made by Williamson et al. (1989) in a study of the region.They estimated that between 1961-1978, approximately the period of maximum historical subsidence, between 20-70% of storage change was expressed through deformation.However, while authors at the time suggested that any future large storage losses would unlikely be expressed through deformation, since inelastic compaction is an irreversible process (e.g.Poland et al. 1975;Lofgren 1978), the findings here show that the expression of storage declines through deformation of aquifer materials continues today.This is most likely due to the shifting location of the subsidence bowl (Fig. 1b), which indicates that today's compaction is occurring inelastically from sediments which have not experienced previous major stress.Another contributing factor is the continued drive to ever-deeper pumping, which has also exposed new clays to inelastic deformation.
That the size of the decline was similar for the two components indicates that neglecting the component of storage change expressed through the deformation in this area is a poor assumption.It is anticipated that the deformation component of storage change will be important in areas of the Valley beyond the study area, especially in those areas with recent deep well installations pumping and those areas where historical subsidence was not large.It is suggested that the component of storage change expressed through the deformation of aquifer materials should not be neglected when monitoring storage in the San Joaquin Valley, and that estimates which do so (e.g.Brewster et al. 2015;Scanlon et al. 2012) may substantially underestimate storage declines.
This finding may have some global implications.While the land subsidence rates in the San Joaquin Valley are high compared to those of many aquifer systems worldwide, there are other notable examples with subsidence in the tens of centimeters per year, including Mexico City (Chaussard et al. 2021) and locations in Iran (Haghshenas Haghighi and Motagh 2019), which are similar hydrogeologic settings to the San Joaquin Valley.Further work must be done to determine whether the current findings generalize to these aquifer systems, but it is expected that the storage component expressed through the deformation of aquifer materials may be important in any rapidly deforming aquifer system.
The monotonic decline of the deformation component contrasts with the large variations of the saturation component.The absence of recovery in the deformation component is consistent with the study of Smith et al. (2017), which focused on a similar area to this study and showed that the compaction of aquifer materials was primarily permanent (inelastic) between 2007-2010.In contrast, uplift has been observed in the Westlands Water District (Neely et al. 2021), indicating that storage recovery can be expressed through the deformation of aquifer materials.It is suggested that recharge in the study area has primarily been occurring in the shallow aquifer zone, where it can be expressed as a change in saturation.The absence of any recovery in the deformation component shows that, at the scale of the entire study area, the deeper zones are not seeing a storage recovery, even in wetter years.This is despite some observations of rises in head.However, it is not expected that this finding is applicable across the Valley, and the possibility of future uplift in the Tule and Kaweah subbasins cannot be ruled out.
The estimates of the elastic and inelastic parts of the deformation component indicate that almost all the storage losses expressed through deformation occurred inelastically.The mean annual rate of inelastic deformation was 0.35 km 3 year -1 .In a single year, this was at least twice the size of the elastic deformation (magnitude <0.2 km 3 year -1 ), and throughout the entire study period, the total inelastic storage loss was over 20× larger than the elastic seasonal deformation.Although the method to separate these components was simplistic, the finding is similar to other published studies which have found the majority of the subsidence to be inelastic (Chaussard and Farr 2019;Ojha et al. 2019;Smith et al. 2017).As previous authors have noted, this finding implies that the subsidence represents a permanent loss in storage capacity, i.e., the pore space lost through subsidence would not be recovered if hydraulic head rises."

Total storage change: methods and data
An estimate of total storage ΔS tot was obtained by combining the estimates of the saturation and deformation components of storage change calculated in sections 'Estimating the saturation component of storage change' and 'The deformation component of storage', according to the following equation: This sum was calculated and compared with storage estimates from two, independent sources.The first independent estimate came from the California Central Valley Groundwater-Surface Water Simulation Model (C2VSim; Brush et al. 2013).The available version of C2VSIM gives storage estimates for the Tule and Kaweah subbasins up to 2015 but a version produced by Alam et al. (2021), which extended the record to 2019, was used here.The second came from the GSAs in the study area, which have made estimates of storage over the entire 2015-2021 period as required by the SGMA.The GSA estimate of storage change was produced using an amalgam of methods.For the period 2015-2017, the GSAs in both subbasins used a groundwater budget approach (Harder and Lewis 2020; Petersen and Nicely 2020), but for 2018-2021, the methods differed by subbasin.For 2018-2021, the Kaweah subbasin summed the product of shallow head and specific yield and InSARderived subsidence (Provost and Pritchard Consulting Group 2020, 2021, 2022).In the Tule subbasin, 2018-2021 estimates were made using only shallow head measurements and specific yield, omitting the storage changes expressed through deformation (Harder 2020(Harder , 2021(Harder , 2022)).

Total storage change: results and discussion
The estimates of total storage change are shown in Fig. 9.The new estimate is shown as the thick black line and the components (deformation and saturation) are additionally shown as thinner black lines.The main trends, which are (11) ΔS tot = ΔS sat + ΔS def captured by all estimates, are a decline in storage during 2015-2017, a recovery during 2017-2019 and, with the exception of the estimate from C2VSim, which does not extend beyond 2019, further decline in 2019-2021.The new method estimates the total storage loss over the entire period to be 5.0 km 3 .
The calculated estimate of storage change, obtained by summing the two components estimated using head and InSAR, matches the general trends seen by the C2VSim model and the GSA estimate in the study area.However, it is expected that the approach taken here gives a more accurate result due to its treatment of the deformation component.Although the C2VSim model does simulate the deformation of aquifer materials, the code used (IWFM v3.02) does not allow for delayed drainage of interbeds and assumes the deformation of aquifer materials is synchronous with head changes in the coarse-grained materials (Central Valley Modeling Unit Modeling Support Branch Bay-Delta Office 2018).This means the C2VSim estimate is prone to the same overprediction of recovery during periods of rising head as can be introduced by assuming the linear head-deformation relationship of Eq. ( 5).In the GSA estimate, which uses an amalgam of approaches, the storage change expressed through aquifer deformation is neglected in the Tule subbasin during the period 2019-2021.As shown previously, that component accounts for ~50% of storage losses, and the GSA approach likely underestimates storage decline from 2019-2021 for this reason.The GSA and C2VSim estimates are most similar to the saturation component in Fig. 8, whereas the new estimate containing both components diverges from this, particularly during the storage recovery period 2017-2019.
Beyond the study area, storage estimates are needed at the scale of the entire San Joaquin Valley.Both InSAR data and shallow head measurements are available at this scale, so the method demonstrated here would readily upscale.It is suggested that the current approach would provide improved storage estimates compared to the existing regional estimates which use head and effective storativity (e.g.Alam et al. 2021or Ahamed et al. 2022) because it would avoid the assumption that the deformation component has a linear relationship to measured head.In particular, based on the findings from section 'The deformation component of storage', it is thought that these estimates are most prone to overestimate storage recoveries during periods of rising head measurements.It is suggested that future work should apply the approach to update existing regional estimates.
Because InSAR provides an integrated measure of aquifer deformation, another benefit of the method is that it accounts for any deformation occurring in the unconfined aquifer.This term tends to be ignored in other estimations, but data and modelling both show it to be nonzero (Lees et al. 2022).Although evidence still shows upper aquifer deformation is a relatively minor component of the total deformation, its inclusion provides a more complete picture of storage change than would be obtained by neglecting it.
For the purposes of SGMA, which requires annual reporting of storage estimates, the methodology applied here has practical appeal.The input data required-shallow groundwater measurements and InSAR data-are available with high latency; although one area in which data availability could be improved would be the existence of a single, centralized repository for shallow groundwater measurements including all data gathered by local agencies.The methodology is also relatively simple to apply.Flux-based approaches, such as the ones used by the Kaweah and Tule GSAs to estimate 2015-2017 storage change, require a greater number of data types and are typically complex to compute.The same is true for groundwater models like C2VSim, which additionally require substantial computational resources.It is expected that the demonstrated methodology to estimate storage changes can be beneficial for GSAs in developing Fig. 9 The temporal evolution of total storage estimated using the new approach, along with the independent estimates from the GSAs and the C2VSim groundwater model is shown.The major periods of decline, recovery and decline are marked, as are the start and finish of major droughts plans for, and monitoring compliance with, SGMA, and it is anticipated that the current study can speed adoption of this approach.
For all its advantages, the new methodology in this study contains uncertainty from several sources.First is the fact that shallow head measurements are imperfect, as is knowledge of specific yield.Some error is also introduced by the assumption that the InSAR-measured surface deformation is exactly equal to the storage change expressed through deformation; this could be false in the presence of tectonic deformation, attenuation of aquifer system compaction, or compaction in the unsaturated zone, although it was previously described how these factors are unlikely to be significant in the study area.It is suggested that improved practices for gathering shallow head measurements, and more detailed studies quantifying local-to-regional scale values for specific yield, would be important contributions to improve the robustness of storage estimates made using this new methodology.

Conclusions
The road to sustainable groundwater management in California has immense challenges, many of which must be met at the local level as GSAs decide upon and implement their sustainable groundwater management plans.Here, it has been demonstrated how existing methods to quantify storage change using hydraulic head measurements can give misleading results due to difficulties in capturing the component of storage change expressed through the deformation of aquifer materials.A method whereby the remote-sensing technology of InSAR is combined with shallow head measurements to obtain improved storage estimates has been proposed, and its application demonstrated through a detailed data-intensive study.The authors hope to see rapid adoption of this approach in the San Joaquin Valley and anticipate that the findings of this study can inform estimating groundwater storage in the many groundwater basins worldwide experiencing subsidence associated with unsustainable groundwater use.

Additional methods for separating elastic and inelastic parts of deformation component
To estimate the inelastic and elastic components of the storage component estimated through deformation, InSAR measurements for each time channel were taken (the time channels in the raw data are at 5-7 days, from March 2015 to October 2021) and gridded by averaging measurements within a search radius of 8 km.(This was selected as it was the mean search radius used for interpolating the head measurements.)Then, this surface was integrated under using the "grdvolume" function from the Generic Mapping Toolbox (Wessel et al. 2019), which resulted in a time series of total deformation at 5-7 days temporal resolution.To obtain the inelastic part, an exponential decay function of the following form was fitted to the time series: which is a similar formulation to the one given in Chaussard et al. (2014).To fit the function, the Trust Reflective algorithm implementation from the scipy Python package was used.n = 4 was used as the lowest value which gave a good fit.
It was assumed that the residual was an estimate of the elastic part.The regular periodicity of that part supports this assumption (Fig. 8).

Fig. 1 a
Fig. 1 a Map illustrating the location of the Kaweah and Tule subbasins, which together comprise the study area.The inset shows the location of the San Joaquin Valley within California.b The same map, this time showing the location of the subsidence which occurred between 1925 and 1970 and between 2015 and 2020.Note that the orange areas, i.e. the overlap of both the 1925-1970 and 2015-2020 ΔS tot = ΔS sat + ΔS def + ΔS w (2) ΔS � unc = S y Δh meas

Fig. 2
Fig. 2 Schematic illustration of the physical geography and hydrogeology of the study area.Figure modified from Faunt (2009)

Fig. 3
Fig. 2 Schematic illustration of the physical geography and hydrogeology of the study area.Figure modified from Faunt (2009)

Fig. 4 a
Fig. 4 a The locations of wells used in the new dataset of shallow measurements.Background color indicates topography.b An illustration of the interpolated shallow head map for spring 2018.The full

Fig. 5
Fig. 5 The calculated saturation component of storage change, with error bars representing the interquartile range obtained from the bootstrapping analysis

Fig. 6 Fig. 7
Fig. 6 The locations of extensometer sites used in this study.Background color indicates topography

Fig. 8
Fig. 8 The calculated evolution of the deformation component of storage change in the study area.Shown in the two grey lines are the estimated inelastic part (dashed) and the estimated elastic part (dot-dash)

Table 1
Details of the four extensometers used in this study