Apportioning deformation among depth intervals in an aquifer system using InSAR and head data

Land surface subsidence due to excessive groundwater pumping is an increasing concern in California, USA. Interferometric Synthetic Aperture Radar (InSAR) is a remote sensing technique for measuring centimeter-to-millimeter surface deformation at 10–100 m spatial resolution. Here, a data-driven approach that attributes deformation to individual depth intervals within an aquifer system by integrating head data acquired from each of three screened intervals in a monitoring well with InSAR surface deformation measurements was developed. The study area was the Colusa Basin in northern Central Valley. To reconstruct the surface deformation history over the study area, 13 ALOS-PALSAR scenes acquired between 2006 and 2010 were processed. Up to ~3-cm year−1 long-term subsidence and up to ~6 cm seasonal subsidence were observed using the InSAR technique. The technique developed in this paper integrates the InSAR-observed seasonal deformation rate and the co-located head measurements in multiple depth intervals to estimate the elastic skeletal storage coefficient, the time delay between the head change and the observed deformation, and subsequently the deformation of each depth interval. This technique can be implemented when hydraulic head measurements within each depth interval are not correlated with each other. Using this approach, the depth interval that contributed the most to the total subsidence, as well as storage parameters for all intervals, are estimated. The technique can be used for identification of the depth interval within the aquifer system responsible for deformation.


Introduction
In the Central Valley of California, USA, groundwater pumping for irrigation has increased dramatically due to the recurrent droughts in the past decade. Much of the Central Valley's aquifer system is composed of fine-grained materials that are highly prone to compaction due to water extraction (Leake and Prudic 1991;Williamson et al. 1989). As a result, excessive pumping in some areas has led to aquifer compaction and associated subsidence that can be observed at the ground surface (Farr and Liu 2014;Faunt et al. 2009Faunt et al. , 2016Poland et al. 1975). In addition to the potential impact on infrastructure, aquifer compaction due to over-pumping can be permanent if the groundwater level declines beyond the maximum historical stress (Riley 1969). In one region of the San Joaquin Valley, the southern part of the Central Valley, there has been 25 cm year −1 of permanent compaction over the 2007-2010 and 20122015 drought periods, with significant delayed compaction for several years following the droughts (Chaussard and Farr 2019;Smith and Knight 2019;Smith et al. 2017).
While it had long been recognized that there was a link between groundwater pumping and subsidence of the ground surface in the Central Valley, the use of Interferometric Synthetic Aperture Radar (InSAR) to map surface deformation throughout the valley (e.g., Farr and Liu 2014;Smith et al. 2017) significantly raised the visibility of this issue. InSAR is a remote sensing technique for measuring surface deformation with 10s-100s of meters spatial resolution and millimeter-tocentimeter line-of-sight (LOS) measurement accuracy, which has been used on a regular basis to derive the spatial distribution of deformation throughout the valley. These high-quality surface deformation data, along with recent legislation requiring sustainable groundwater management in California (Senate Bill SB 1168 (Pavley), Assembly Bill AB 1739 (Dickinson), and Senate Bill SB 1319 (Pavley), are motivating local and state agencies to take on the task of proactive groundwater management to limit the extent of subsidence associated with the extraction of groundwater. The focus of this study was the development of a new method for the integration of InSAR data with head data that could be used to support decision-making related to sustainable groundwater management.
Previous studies have shown the use of InSAR for mapping long-term and seasonal surface deformation due to pumping and recharge of aquifers in both urban and agricultural regions (e.g., Galloway et al. 1998;Hoffmann et al. 2003;Hoffmann et al. 2001;Miller and Shirzaei 2015;Reeves et al. 2014a;Reeves et al. 2014b;Watson et al. 2002). InSAR surface deformation data, when integrated with in-situ hydraulic head data and geomechanical models, can be used to estimate the skeletal storage coefficients as well as the amount of permanent storage loss associated with the excessive pumping (e.g., Chaussard et al. 2017;Chen et al. 2017Chen et al. , 2016Smith et al. 2018Smith et al. , 2017. Because InSAR measures the total amount of deformation over the entire depth of the aquifer system, most previous studies have treated an aquifer system as a single deforming unit. In order to apportion the deformation at the ground surface among discrete subsurface layers, Smith and Knight (2019) developed a method to integrate InSAR data with geophysical data, acquired using an airborne platform and imaging sediment type to a depth of~400 m. This allowed them to model deformation as a function of the clay thickness and hydrologic properties at specific depths. While effective, this method requires high-quality geophysical data, which are not always available. The method also assumes that the aquifer system is interconnected, or that the change in head in high-permeability aquifer material is equivalent for all depths of the aquifer, which is not always the case.
In this study, a method was developed to determine the depth interval(s) over which aquifer deformation is occurring at specific well locations. The study site for this paper is a well in the Colusa Basin in the Sacramento Valley, in the northern part of the Central Valley. The Colusa Basin is the rice-growing capital of California, with groundwater supplementing surface water supplies, when needed, for irrigation. While subsidence throughout the Sacramento Valley has been much less than that observed in the San Joaquin Valley, due to more abundant surface water and less pumping of groundwater, there is still considerable interest among water agencies in assessing the vulnerability to subsidence under future scenarios of climate change and declining surface-water supplies. The depth interval responsible for observed seasonal deformation from 2006 to the end of 2010, and properties controlling deformation, were estimated by integrating InSAR data with head data from a well in the northern Central Valley. Application of this method could provide critical information for groundwater management, making it possible to strategically select the optimal locations and depth intervals for both the extraction and targeted recharge of groundwater.

Background: the deformation to head relationship
In an aquifer where the hydraulic head is fluctuating above the preconsolidation head, the aquifer deforms elastically, and the subsidence due to pumping is reversible when recharge occurs (Hoffmann et al. 2001;Riley 1969). Deformation can occur in both unconfined and confined aquifers (Leake and Galloway 2007;Poland and Davis 1969), but typically is greater in confined aquifers, which typically have higher drawdowns. When the head remains above the preconsolidation head (Riley 1969), the relationship between the change in head, Δh, in a confined aquifer and the associated elastic deformation, Δb e , can be written as: where S ke is the elastic skeletal storage coefficient. S ke can be approximated as S ske × b, where S ske is the specific elastic skeletal storage coefficient of the aquifer, and b is the initial thickness of the aquifer. In an unconfined aquifer, deformation is also affected by the change in overburden stress. The relation between deformation and change in head is given by Leake and Galloway (2007): where S y is the specific yield, or the ratio of the volume of water that can be drained by gravity to the total volume of the unconfined aquifer (Fetter 2001). Typical values for specific yield in the study area are on the order of 0.1 (Faunt et al. 2009;Terzaghi 1925). In this study, Eq. (1) was used for all aquifers, which results in a slight (~10%) underestimation of S ke in unconfined aquifers.
In the case that the hydraulic head declines below the preconsolidation head, inelastic deformation occurs, and the aquifer permanently loses some of its storage capability (Riley 1969;Terzaghi 1925). Under this scenario, the total amount of inelastic deformation Δb v can be written as: where S kv is the inelastic skeletal storage coefficient. S kv can be approximated as S skv × b, where S skv is the specific inelastic skeletal storage coefficient of the aquifer. S skv is typically 1-3 orders of magnitude higher than S ske (Sneed 2001).
In addition, an aquifer that contains materials of variable sediment type can exhibit a complex relationship between the changes in head measured in the aquifer and the timing of the associated deformation (Pavelko 2004). Within the interconnected sandy region of an aquifer, head tends to equilibrate rapidly, and the associated deformation occurs along with the head changes almost simultaneously. There can, however, be a significant delay between the change in head in the sandy regions and the change in head experienced within interspersed clay-rich regions, resulting in a time delay in the deformation of these regions. In the case of a decrease in hydraulic head of equal magnitude on both the top and bottom boundaries of a clay lens, the time it takes for the unit to equilibrate with the hydrologic head above and below it, τ, can be approximated by the following equation (Riley 1969): where K v is the vertical hydraulic conductivity, S s is the specific storage for the clay lenses, and b 0 is the thickness of the compacting clay.

Study site and available data
The location of the study site is in the Colusa Basin in the northern Sacramento Valley, shown in Fig. 1a. A geologic cross-section is also shown in Fig. 1c. The main formation of the upper 300 m in the study area, where most wells are drilled, is the Tehama Formation. This formation consists of Pleistocene-aged coarse-and fine-grained volcanic deposits (California Department of Water Resources 2014). ALOS PALSAR SAR data were used, acquired over the 2006-2010 period, as the basis for the InSAR analysis to quantify the total deformation of the ground surface. Head data were obtained through the California Statewide Groundwater Elevation Monitoring (CASGEM) program. Head data from 24 monitoring wells located within the SAR scene in the Colusa Basin were reviewed. The approach of this paper relies on attributing surface deformation to changing head at different depths, so head measurements at multiple depths were required. Only 14 of the 24 monitoring wells had a measurement history in three depth intervals that overlapped with the 4-year InSAR measurement period. In these 14 wells, multiple head measurements, between 20 and 31, were available during the 2006-2010 period, when the ALOS PALSAR data were acquired. These wells were all constructed as nested wells, which involved the installation of three casings, within one borehole, placed and containing a screened interval so as to measure head in three discrete depth intervals. Since the head measurements and SAR data were not acquired simultaneously, the measured data were first interpolated to a daily time series. Then a moving average, with a 15-day timeframe, was calculated to estimate head levels between the consecutive head measurements for comparison with the irregular SAR measurements.
The method to apportion deformation among various depth intervals required uncorrelated head data from multiple depth intervals. Correlated head data suggest that the corresponding depth intervals may be treated as a single aquifer unit. A review of the data yielded one monitoring well ('monitoring well' in Fig. 1) where the head data were uncorrelated. In all other wells, the measured heads in at least two depth intervals were correlated, with a correlation coefficient of 0.6 or greater between at least two of the three nested wells. In the selected well, uncorrelated head measurements over three distinct depth intervals were available, with 24 measurements at each depth in the time period December 2006 to November 2010. In addition, information was available from a detailed lithologic log- Fig. S1 of the electronic supplementary material (ESM). The top depth interval is screened over the depth range of 15-18 m, with the geological material described as poorly graded gravel and medium to coarse grained sand. The middle depth interval, screened over the depth range of 39-45 m and 51-54 m, is separated from the top interval by a 6-m-thick clay layer. The lithology over the screened interval is described as a mixture of poorly graded gravel and sand, with clayey sand and silt in between the two depth ranges. The bottom depth interval, screened over the depth range of 244-250 m, is separated from the middle interval by a 45m-thick clay layer, and is described as clay-rich, mixed with poorly graded sand and claystone.

Processing deformation time series with InSAR
The InSAR technique uses the phase difference of two SAR images to map surface deformation along the LOS of the satellite (Hanssen 2001). For this study, 13 L-band ALOS PALSAR scenes, acquired between December 2006 and November 2010 from the Alaska Satellite Facility, were downloaded. These scenes were used to generate 44 small baseline interferograms (less than 3,000-m spatial baseline and less than 3 years temporal baseline) using a motion compensation imaging radar processor developed by the Stanford radar interferometry group (Zebker et al. 2010). Topographic phases were removed using SRTM 3-arc sec digital elevation database. An adaptive interpolation filter was applied between Persistent Scatterer (PS) pixels to reduce the impact of vegetation decorrelation in the InSAR phase measurements (Chen et al. 2015a, b), and unwrapped the interferograms using SNAPHU and a minimum-cost flow algorithm (Chen and Zebker 2001). A long-wavelength linear ramp was removed from each interferogram, and InSAR pixels over the mountain areas with elevation greater than~250 m amsl were masked out. The permanent GPS station (P270), operated by California Special Reference Center (CSRC), near Corning, California, was used as the reference point ( Fig. 1), where no substantial surface deformation due to groundwater pumping and recharge was recorded during the 2006-2010 period to calculate the deformation. A linear deformation model was used to solve for the longterm LOS deformation rate v c at each pixel of interest using a Small Baseline Subset (SBAS) approach (Berardino et al. 2002): where δϕ = [δϕ 1 , δϕ 2 , …δϕ n ] is the n small baseline InSAR measurements derived from SAR acquisitions at the pixel of interest. Matrices B and P are defined as in Berardino et al. (2002).
To extract the seasonal deformation signal over the study area, the seasonal vertical displacement rate between each consecutive SAR acquisition v = [v 1 , v 2 , …v m−1 ] was calculated using the following equation: where δϕ is the InSAR data vector δϕ in Eq. (5) with the longterm subsidence removed. An additional integration of the velocity solution v over time yielded the seasonal deformation time series at the pixel of interest. Under the assumption that there is no significant variation in horizontal motion over the study area, the LOS displacement was converted to vertical displacement using the radar LOS direction unit vector e = [0.62 0.11 -0.78] at the midswath as d LOS = -0.78d vertical . The look angle variation over the study area is not substantial. Galloway and Burbey (2011) note that in areas with regional groundwater pumping as opposed to localized wells or wellfields, vertical subsidence is the dominant form of deformation. Since this study area has regional groundwater pumping, and no faults with significant slip rates (>1 mm) are present in the study area, it is considered that the assumption that all deformation is vertical is a reasonable one.

Joint analysis of InSAR and head data
A data driven method that integrates all head measurements at different depths with the InSAR deformation measurement to solve for the deformation, and values of S ke and τ for each of the depth intervals was developed. In this study, it was assumed that all deformation is elastic, including seasonal and long-term (multiyear) deformation. This assumption was  (2014). UPVF stands for Upper Princeton Valley Fill made because regional head data (within a 5 km radius of the monitoring well) show that groundwater levels over the study period are higher than the lowest levels previously experienced in the 1970s (Fig. 2). For this reason, it is unlikely that any substantial inelastic deformation occurred in the study period.
The 4-year trend was also removed from the head measurements, presuming that this trend in head would correspond to the long-term signal removed from the InSAR measurements. It was presumed that the deformation of the thick clay layer, between the middle and bottom depth intervals, was negligible relative to that of thinner clay lenses, which, in aggregate, in much of the Central Valley typically make up a larger portion of clay in an aquifer system (e.g., Faunt et al. 2009;Smith and Knight 2019), and deform much more rapidly in response to changes in head because they are not as thick.
The following equation shows the sum of calculated deformation D, in each depth interval a, as equivalent to the InSAR integrated vertical displacement d: Expanding D results in the following equation, with the two unknown parameters, S ke and τ, for each of the three depth intervals: where t denotes the InSAR acquisition time, and H is the head level at the InSAR acquisition date after accounting for the delay, τ a , between head change in the depth interval and thickness change of the clay lenses that are principally responsible for the deformation observed at the ground surface. The inverse problem was set up to solve for the parameter values that would minimize the difference between the InSAR-observed deformation, d, and ∑ 3 a D a . A grid search method was used to achieve a range of possible values of S ke and τ for each of the depth intervals. In the grid search method, the model is iteratively run for a range of parameters. Parameters that result in model output with lower error are considered more representative of the true parameter values. The model error in predicting d from ∑ 3 a D a calculated as the 2-norm of the residuals for each combination of parameter values.
The input grid of parameters needed to be sufficiently large so as to ensure that the correct values are captured. Assuming elastic deformation, S ke, a , is defined for each depth interval as: For S ske the extent of the grid search was set to fall within the range of the values reported from a literature review by Sneed (2001) of 2.0 × 10 −6 to 2.3 × 10 −5 m −1 . As the deformation is dominated by interspersed clay-rich regions, a range of 1-250 m was estimated for the total clay thickness b 0 in each depth interval; this range was taken from California's Groundwater Bulletin (2003) and represents the range in values of total observed clay thickness in the Sacramento Valley. These ranges were used with Eq. (9) to set the range of possible values for S ke . The study period is shown in the red ribbon. Note that the depth to groundwater in the 1970s exceeded the depths experienced during the study period Values ranging from 0 to 140 days were used for τ a . This is a conservative range, as it covers a much larger range of values of τ than was observed by Chen et al. (2016), whose maximum observed τ was 41 days in a similar study that evaluated time lags in elastically deforming clays. Elastically deforming clays have much shorter time lags than inelastically deforming clays (Smith et al. 2017) and thus are likely to equilibrate relatively quickly with aquifer head levels.

Results and discussion
Long-term and seasonal deformation results derived from InSAR Figure 3 shows the average long-term subsidence rate between December 2006 and November 2010 as inferred from ALOS PALSAR data over the northern Sacramento Valley. No substantial long-term subsidence was observed in the northern and southern portions of the Colusa Basin (Fig. 2), over the study period. A maximum long-term subsidence rate of about 3 cm year −1 occurred in the portion of the Colusa Basin between Orland and Willows cities, where the head also dropped~12 m between 2004 and 2015 (DWR, CA, Northern Region office, 2015).
The results of both long-term and seasonal subsidence suggested that there is intense groundwater pumping due to farming activity between the cities of Orland and Willows that caused the subsidence bowl in this area of the basin (Fig. 3). The higher amount of subsidence occurred closer to the foothills, which are distant from the Sacramento River that runs through the center of the Sacramento Valley. The portion of the aquifer adjacent to the river is directly recharged from the riverbed and bank, thus reducing the subsidence in response to groundwater pumping.
The InSAR results were validated using the GPS time series recorded at a permanent station P344 (as shown in Figs. 1  and 3) in the Colusa Basin. Figure 4 shows a comparison of InSAR and GPS surface deformation time series at the GPS station P344 (39.929°N, 122.028°E). At this location, no substantial long-term deformation was observed. The InSAR and GPS seasonal deformation estimates were similar in magnitude, with less than~1 cm difference between the two measurements. While the InSAR time series and GPS station data match within roughly 1 cm, there is more than 1 cm variation at some time periods, indicating some error in the InSAR data at this GPS station. Fig. 3 The long-term subsidence map over the northern Central Valley, California, between December 2006 and November 2010 as inferred from ALOS PALSAR data. The studied monitoring well is depicted on the map. The study area is Colusa Basin, which is located at the western portion of the map defined by the rectangle. Negative values in the scale bar (warm colors) denote subsidence. The maximum subsidence rate is 3 cm year −1 during a 4-year period The two major error sources in the ALOS Butte time series are vegetation decorrelation (and the associated phase unwrapping errors; Zebker and Villasenor 1992) and tropospheric noise (Zebker and Villasenor 1992). It is important to note that the uncertainty at different pixel locations often varies. The GPS validation was used as a rough estimate of error magnitude (~1 cm) to justify that the actual noise level at the well location is likely much smaller than the validation GPS location.
To mitigate the phase decorrelation noise, an adaptive interpolation between persistent scatterer (PS) pixels, those points at which surface scattering properties do not change much over time (e.g., road and man-made structure) was employed and thus decorrelation artifacts are minimal (Chen et al. 2015a, b). The accuracy of this PS-interpolation depends on the density of the PS pixels as well as the false-positive rate of the selected PS pixels (Agram and Zebker 2009). At the well of interest, very minimal PS-interpolation artifacts were observed and phase unwrapping errors in all interferograms used in the analysis, and therefore the uncertainty due to decorrelation noise at this location is not substantial.
Tropospheric noise consists of a stratified component that correlates with topography (Doin et al. 2009) and a turbulent component that is random at time scales longer than 1 day (Emardson et al. 2003). The stratified tropospheric noise is trivial given that the area of interest is relatively flat. Tropospheric turbulence noise increases with the distance from the InSAR reference point (Staniewicz et al. 2020). Because the monitoring well is much closer to the reference GPS station, the tropospheric noise at the well location is expected to be smaller than the centimeter-level noise observed at the GPS validation point. To be conservative, for this study, a 1-cm error for the InSAR time series data was estimated.

Estimates of aquifer properties and deformation
Using the method described previously, the head and InSAR data were analyzed to estimate the amount of deformation in each depth interval by solving for the skeletal storage and time lag values required to match the observed deformation. The estimated S ke and τ for each of the three depth intervals was assumed to be constant during the study period. Figure 5a shows the InSAR-observed deformation vs. the modeled deformations for the 2006-2010 period. Figure 5b depicts the modeled deformation (in meters) for the top, middle, and bottom depth intervals. Figure 5a represents the comparison between the deformation measured using the InSAR technique and the cumulative deformation of all three depth intervals calculated by the model. Figure 5a,b shows a range of possible modeled cumulative deformation, and deformation for each depth interval, that are considered plausible outcomes based on their model error. The threshold in choosing which parameters to include in this category was determined based on the root mean squared error (RMSE) between modeled deformation and InSAR-observed deformation. Models with RMSE less than the 10th percentile of all RMSE values from the model iterations in the grid search were considered to have reasonably good fit to observed deformation. The range of values from these models were then used to plot the ribbons (shaded areas) in Fig. 5a,b. In addition to ranges in modeled deformation, error bars are plotted on the InSAR data. A 1-cm error in the InSAR data was assumed based on the comparison with the GPS station (Fig. 4) and similar work from previous studies (i.e., Chen et al. 2017).
Not surprisingly, the top depth interval shows very little deformation due to the low magnitude of its head change relative to the head changes in the deeper depth intervals. While there are comparable head changes in the middle and bottom depth intervals (Fig. 2a of the ESM), it is the bottom depth interval that exhibits the largest deformation and is the principal source of the observed surface deformation. As can be seen in Fig. 5a, there is a good agreement between the total modelled deformation and integrated deformation measured by InSAR. Assuming a 1 cm uncertainty in the InSAR result, derived from the comparison between InSAR and GPS (Fig.  4), the total modeled deformation matches the InSAR measurement reasonably well within the expected error during the study period. The root mean square error (RMSE) and correlation coefficient were used to check the discrepancy and fit between the total modelled and observed deformation. The RMSE was 0.012 m (1.2 cm), and the correlation coefficient between InSAR and total modelled deformation was 0.65.
In Fig. 5c-e, the results of the grid search are shown as heat-maps to illustrate the accuracy of the estimates of S ke and τ. Cool colors (blue) indicate low error in the prediction of the InSAR-measured deformation, d. Parameter values corresponding to low error are more likely to be representative of the properties of the depth intervals. As can be seen in the results for the middle depth interval (Fig. 5d), low error was obtained with a wide range of values for S ke , and essentially all values of τ; i.e., the values of the parameters are not very well constrained by the measurements of InSAR deformation and head in this depth interval. The top interval (Fig. 5e) shows low error is obtained from essentially any combination of S ke and τ, meaning these parameters are very poorly constrained by measurements of InSAR and head, likely because there was very little head fluctuations at this interval, and little contribution from this interval to the InSAR deformation. The estimated values of S ke determined in the region of the plots with the lowest model error, are between 0.0 and 5.0 × 10 −3 for the top depth interval and between 0.0 and 2.5 × 10 −3 for the middle depth interval. Figure 5c shows the results for the bottom depth interval. This is the interval that is responsible for most of the observed deformation. In contrast to the wide region in the heat map with low error found for the other two intervals, in this interval low error is found in a relatively narrow range of S ke values between 2.0 × 10 −3 and 3.3 × 10 −3 , with the lowest model error found at a value of 3.3 × 10 −3 . The range of τ values is between 64 and 107 days, with no significant difference in model error of values within that range. Thus, the model is most sensitive to the S ke and τ of the bottom depth interval, which is contributing the most to total deformation.
The finding of this study that most deformation occurs at or near the depth of the deepest interval is consistent with independent geologic information from the driller's log for the monitoring well (see Fig. S1 of the ESM). Higher deformation is typically associated with a large amount of compressible sediments such as clays, that are experiencing a change in head. As noted previously, the top, middle, and bottom intervals are located at depths of 15-18 m, 39-54 m and 244-250 m, respectively. The top and middle intervals are located at fairly shallow depths, with relatively little clay in their vicinity (10% or less fine-grained material in the 30 m above and below each interval). The top interval has no recorded confining layers above it, and has low seasonal variations in head, consistent with an unconfined aquifer. The middle and bottom intervals both have clay layers above them that could be confining, and both also have significant seasonal variations in head, so are likely semiconfined or confined. While both the middle and bottom intervals have some confinement, only the bottom interval is surrounded with significant amounts of fine-grained material, with over 50% of the  Fig. 1). a Observed deformation by InSAR compared to the cumulative modeled deformation from all three depth intervals. b The contribution of each depth to the total deformation.
Ribbons indicate error bars (5th-95th percentiles). c-e Shows the error (RMSE, in m) for S ke and τ parameters for the lower, middle and top layers, respectively. Cool colors (blue) indicate low model error material in the 30 m above and below the interval composed of fine-grained sediment. For these reasons, based on the data from the geologic log alone, the bottom interval is the only interval that is likely to experience significant deformation. These geologic observations agree with the independent analysis of this study on depth-dependent deformation. Sneed (2001) compiled parameter estimates for S ske (S ke divided by the thickness of the compacting material) in the Central Valley, and found that the values ranged from 2.0 × 10 −6 m −1 to 2.3 × 10 −5 m −1 , with an average of 9.0 × 10 −6 m −1 for aggregate combinations of sand and clay. The implemented method does not allow a direct comparison with these values, since the total thickness of compacting material is unknown. However, if the best-constrained layer (the deepest layer) is considered and estimate the total thickness that is compacting to be less than the distance between this layer and the middle layer, since the two do not show significant correlation, it can be estimated that up to roughly 200 m of sediments could be experiencing a change in head and thus compacting. The amount of sediments compacting could be much lower than this, so this thickness is likely an upper bound. With this value, the S ske from was estimated to be 1.7 × 10 −9 m −1 or greater (if the thickness compacting is lower), slightly higher than the values estimated by Sneed (2001). While the estimate in this study is slightly higher than the values reported by Sneed (2001), those values were mainly collected in the southern portion of the Central Valley, and may not be representative of geologic materials in the north. Chen et al. (2017) reported S ke values in the San Luis Valley, Colorado, ranging from 2.9 × 10 −3 to 1.8 × 10 −2 , which are slightly higher than those estimated from this study. It is apparent from these variations that the geologic setting plays a key role in the geomechanical properties of sediments.

Two-layer deformation
It was assumed in section 'Estimates of aquifer properties and deformation' that each layer is capable of producing a significant amount of deformation. Since the upper layer at the study well appears to be unconfined and experiences little changes in head over time, with few clay layers, it is reasonable to assume that no significant deformation occurs in that layer. As noted previously, only one of the 14 available wells had head data that were not correlated (meaning the correlation coefficient was lower than 0.6) between each of the layers. However, if it is assumed that no deformation occurs in the upper layer, there is one additional well (well code 397325 N1221233, location shown as 'secondary monitoring well' in Fig. 1) that has heads in the lower two layers with a correlation coefficient below 0.6. The seasonal and long-term head variation in the top layer of this well is minor compared with the middle and lower layers (Fig. S2b of the ESM). This well was analyzed following the method above, with only the Fig. 6 Modeled deformation vs. observed InSAR deformation for the secondary well (location in Fig. 1). a Observed deformation by InSAR compared to the cumulative modeled deformation from all three depth intervals. b The contribution of each depth to the total deformation.
Ribbons indicate error bars (5th-95th percentiles). c-d Shows the error (RMSE, in m) for S ke and τ parameters for the lower and middle layers, respectively. Cool colors (blue) indicate low model error bottom two layers included. The results of this analysis are given in Fig. 6.
While the monitoring well presented in Fig. 5 had low sensitivity to the properties of the top and middle layer, this two-layer model has reasonably good sensitivity to parameter values of the middle and bottom layers. The deformation pattern is similar to Fig. 5a, which is expected as the wells are located~4 km away from each other. The parameter estimates for S ke and τ range from 2.6 × 10 −3 to 3.2 × 10 −3 and from 116 to 134 days, respectively, for the bottom layer, and from 2.1 × 10 −3 to 2.7 × 10 −3 and from 103 to 116 days for the middle layer. The range in modeled deformation values shown in Fig. 6a,b were determined as described in section 'Estimates of aquifer properties and deformation'. The middle and lower layers at this location contribute similar amounts of deformation, but the lower layer in most model runs contributes slightly more deformation.

Conclusion
As groundwater is a major water source for both drinking and irrigation supply, the demand for new technologies and approaches for facilitating long-term and reliable groundwater management has, recently, brought the use of satellite imagery into the field of groundwater hydrology.
The InSAR technique provides spatial information and maps the ground surface deformation over wide areas with very high spatial resolution. Using ALOS PALSAR SAR data, a maximum long-term subsidence rate of about 3 cm year −1 was observed in the region between Orland and Willows, northern California, for the period between December 2006 and November 2010. The seasonal subsidence rate for the same area reached about 6 cm.
In order to relate ground subsidence to groundwater drawdown and sediment compaction, most of the previous studies assumed that the entire aquifer system behaves as a single unit causing subsidence. This study presents a new way of connecting InSAR to head levels, highlighting the potential risks of assuming that the entire system is compacting. Using the developed approach, the depth interval that contributed the most to the total deformation was identified. The results were verified independently by comparing this analysis with existing information on clay content and interval thickness from the geologic log of the study site.
The enormous challenge of sustainable management of groundwater resources encourages us to explore new ways in which remote sensing data can be used to support both groundwater science and management. The application of the methodology developed in this study leads to an improved understanding of the relationship between the withdrawal and recharge of groundwater, and the ground deformation that can be observed with InSAR data. Understanding this relationship is essential for the proactive management of groundwater in areas vulnerable to subsidence.