Modeling 65 Years of Land Subsidence in California’s San Joaquin Valley


 [copied directly from first paragraph of paper]

Land subsidence, caused by groundwater extraction and subsequent subsurface compaction, is an issue of global concern. Since the 1920s, there have been numerous periods of subsidence in California’s San Joaquin Valley leading to widespread sinking of the land surface which has locally exceeded 9 m. The most recent period of severe subsidence, which was triggered by the 2012-15 drought, is now causing damage which threatens the long-term viability of critical water distribution infrastructure in the Valley. However, there is neither a continuous monitoring record of the subsidence nor high-quality records of the hydrologic head changes in the subsurface which have caused the subsidence, making it impossible to understand, and thus mitigate, the subsidence. Here, we leverage subsidence and hydraulic head data from a variety of sources to create and validate a one-dimensional model of subsurface compaction and subsidence over the 65 years between 1952-2017. This model, which simulated up to 7.5 m of subsidence since 1952, provides a complete record of subsidence in our study region by filling crucial gaps in the observed record. Our model reveals the long-term processes causing subsidence, which operated over decades-to-centuries and caused exceptionally high rates of baseline subsidence in 2017, resulting in a critical risk of future subsidence. This risk is exacerbated as the Valley moves into drought conditions again in Spring 2021. We demonstrated an approach which provided the understanding of subsidence in the Valley needed to directly inform sustainable groundwater management, and which is applicable in subsiding regions around the World.


Main text
In California's San Joaquin Valley, land subsidence has been a recurring problem since the 1920s [9][10][11] . Similar to many sedimentary aquifer systems worldwide, the major source of subsidence in the San Joaquin Valley is the compaction of subsurface clays in response to groundwater extraction 12,13 . Between the mid-1920s and early-1970s, a period of extreme subsidence occurred in the Valley, resulting in a maximum observed subsidence exceeding 9 m 3 and millions of dollars of damage to wells, water-conveyance structures and other infrastructure 14,15 . During this period, a comprehensive monitoring program was established under the leadership of Joseph Poland of the U.S. Geological Survey, including regular levelling surveys and a network of extensometers which measured the compaction at different depths in the aquifer system 16 . By 1970, increased availability of surface water greatly reduced groundwater extraction, and subsidence rates diminished to near-zero by 1980 3 ; the monitoring program was largely abandoned in the decades that followed 17 . However, some limited subsidence was again observed during both the 1976-77 and 1987-92 droughts 3,17 . Far worse occurred during two periods of severe drought in 2007-10 and 2012-15, when subsidence returned at rates comparable to those observed in the 1920s-1970s 5 , triggering damage to canals, bridges and other infrastructure 18,19 . This recent subsidence, along with other negative environmental impacts of excessive groundwater extraction, led to the passage of a groundbreaking piece of legislation in 2014 -the Sustainable Groundwater Management Act (SGMA), which requires groundwater managers to take action to reduce future significant subsidence. The passage of SGMA, along with renewed drought conditions in the San Joaquin Valley in Spring 2021, makes understanding the controls on subsidence in the Valley a critical and urgent issue.
While our understanding of the recent subsidence has been hindered by the absence of a comprehensive monitoring network like that of the Poland era, we now have available geodetic data from Interferometric Synthetic Aperture Radar (InSAR) satellites. InSAR has been applied to create maps displaying subsidence in the San Joaquin Valley during several time windows between 2002 and present [20][21][22][23] ; the resulting images of the sinking of land in the Valley are a compelling reminder of the urgent need for sustainable groundwater management. However, InSAR observations have two significant limitations: 1) they provide measurements of total subsidence but contain no direct information about the specific depth intervals over which compaction is occurring and 2) they leave us without detailed information about subsidence outside of the relatively short periods for which high-quality data are available, which are the 2007-10 and 2015-present time periods 24 .
The first limitation is particularly important given the multi-layered aquifer system of the San Joaquin Valley. At most locations, it is comprised of an upper and lower aquifer separated by the regional confining unit known as the Corcoran Clay. While extensometer data gathered in the Poland-era showed that compaction primarily occurred in the lower aquifer 3,9 , the extensometer network was abandoned in the 1980s and there are almost no extensometers active in the areas affected by recent subsidence (Figure S 2). As a result, we have no direct monitoring of the depth at which compaction occurs, which is a key consideration in understanding the controlling mechanisms.
The second limitation has reduced our ability to assess the long-term features in the subsidence record and determine how groundwater extraction in the past has contributed to the recent subsidence. Modeling during the Poland era predicted decades to centuries of so-called deferred subsidence, which is subsidence that occurs due to past, but not current, groundwater extraction 3,25 . However, several recent InSAR-based studies have described deferred subsidence as lasting for only one-to-two years after the 2007-10 and 2012-15 droughts 24,26,27 . This discrepancy likely stems from the short monitoring period provided by the InSAR satellites, rather than a change in the nature of the compaction.
A powerful approach to address these gaps in monitoring and understanding is to use physical models to simulate the subsurface compaction which causes the subsidence. The aquitarddrainage conceptual model, developed during the Poland era 28 , remains our best theoretical description of how groundwater pumping triggers compaction and therefore subsidence. It describes how water pumped from the parts of the aquifer system dominated by coarse-grained materials (sand and gravel) causes a drop in the hydraulic head and corresponding increase in grain-to-grain stresses in those parts of the aquifer system. These increased stresses only cause limited compaction due to the low compressibility of sand and gravel. Over time, however, the hydraulic head drop spreads into adjacent clay interbeds and aquitards, triggering large compaction rates in those highly compressible parts of the aquifer system. This compaction is responsible for the majority of observed subsidence. The timescale over which a head drop propagates from the sands and gravels into a clay layer is related to the clay's hydraulic properties and thickness, and is typically large. This is the cause of deferred subsidence. In this study, our approach was to use a one-dimensional numerical implementation of the aquitard drainage model to simulate subsidence, in order to fill the observational gaps and improve our fundamental understanding of subsidence in the San Joaquin Valley.
Our modelling focused on a small study region just south of the town of Hanford (Figure S 1). The major breakthrough that launched our modelling effort was the creation of the 65-year head timeseries from 1952-2017 in the coarse-dominated material for both the upper and lower aquifers. To make these timeseries, we analyzed and combined head data from eight different wells within the study region, a procedure which is fully described in Methods 1. The remaining model inputs, and the model equations, are described in Methods 2-6. Where uncertainty was identified in input parameters, our approach was to use a grid-search through all possible inputs and validate the model output with observed subsidence data. We validated against observed subsidence data from two InSAR datasets (2007-10 and 2015-17) and from two levelling surveys (1966-70 and 1960s-2004), which provided excellent constraint on both short and long-term subsidence rates. Out of a grid search of over 6000 parameter combinations, 27 model runs agreed with the validation data and were retained; these runs provided our modelled subsidence record. Equipped with this record, we were then able to make important insights into the origin of subsidence in the Valley.

A 65-year record of subsidence history
Our model provides a complete history of subsidence rates over the past 65 years, without the substantial gaps of the observational record. We present the 65-year modelled subsidence record, together with the hydraulic head series used as model input, in Figure 1. The information in the middle panel, which shows the available subsidence data, reveals the temporal gaps which we have filled with our model. We are able to clearly discern a recurring pattern over the past 65 years: cycles of head decline in one or both aquifers, accompanied by elevated subsidence rates, followed by periods when the head was stable or undergoing partial recovery and the subsidence continued but at a decreased rate. These cycles can be linked to known hydrologic events that caused the head variations, shown by the shading in Figure 1.
The first cycle started with head decline caused by high groundwater extraction in the 1940s and 1950s 29,30 , which led to high subsidence rates between 10-30 cm/yr in the 1950s and early 1960s. As surface water supplies became available in the 1970s and 1980s 9 , head levels recovered and subsidence rates fell to 1 cm/yr by 1987. The second cycle began with the 1987-92 drought, where head declined again and subsidence rates increased to 15 cm/yr. After the drought, head recovered and subsidence rates fell again to 3 cm/yr by 1999.
There were three, rapid cycles in the 2000s. The first began in 2001-04, where we observe head declines in both aquifers. This period was not a hydrologic drought: the likely reason for the head declines in this period was a combination of moderately-dry conditions 31 and a reduction in deliveries of surface water to this area 32 . These head declines led to a spike in subsidence rates, reaching 15 cm/yr in 2004. Head then recovered until the next cycle which started with the 2007-09 drought, during which we simulated subsidence rates peaking at 17 cm/yr in 2009, and which finished with two years of stable head and reduced subsidence rates from 2010-2012. The final cycle began with the 2012-15 drought, which caused further, substantial head declines and the greatest simulated subsidence rate of the 65-year modelled period: 37 cm/yr in 2015. In the final two years of our modelled record, 2015-17, head recovered, primarily in the lower aquifer, and subsidence rates fell to 14 cm/yr in 2017. Note that since this 2017 subsidence rate was due entirely to past head drops, it was therefore deferred subsidence. Our model provides the first picture of the subsidence cycles over the past 65 years, including identification of high subsidence rates in the 1987-92 and 2001-04 periods of head decline for which there were no subsidence data in our study region. This 65-year record provides the historic context required to understand the controls on subsidence.  Poland et al. (1975) 9 and the Inter-agency Committee on Land Subsidence in the San Joaquin Valley (1958) 16 . Hydrologic events which caused the head variations are labelled at the top of the figure, and are marked with light orange (head decline) and cyan (head recovery) shading.

Compaction occurs over long timescales
The modelling that we have conducted provides valuable insight into the timescale over which compaction occurs following groundwater extraction. The concept of compaction extending over long timescales, as the head decline propagates into the clay interbeds and aquitards, is widely accepted. However, it is impossible to quantify those timescales without a sufficiently long temporal record of hydraulic head and subsidence. The aquitard-drainage model allows us to calculate an effective time constant ( ) which describes the characteristic timescale over which head propagation, and hence compaction, occurs (described in Methods 8). We tabulated these effective time constants for all 6000+ simulations in our grid search, which gave a range of values from 6 to 2500 years. However, for the 27 simulations which fitted the validation data, the range of time constants contracted to between 60 and 1300 years, which is comparable to those derived from the aquitard-drainage model in the Poland-era 25 . We note that the shortest value that we obtained, of 60 years, is far larger than the one-to-two year timescales mentioned in recent studies based on short periods of subsidence and head data 24,26,27,33 . We believe our longer timescales provide an accurate description of the duration of compaction, which was underestimated in recent studies due the limited length of the observational record. Since thick clay aquitards and interbeds are widespread in the San Joaquin Valley, we conclude that timescales on the order of decades to centuries are required to characterize compaction throughout the Valley.
An important implication of the long timescales of compaction is deferred subsidence. The deferred subsidence rate in 2017 was 14 cm/yr, which is as large as peak subsidence rates in the 1987-92 and 2001-04 periods of head decline. Our simulations allowed us to quantify the link between this large deferred subsidence rate and the multiple cycles of head decline and partial recovery since the 1940s in our study region. Specifically, we quantified the amount of subsidence in 2017 that was deferred from each cycle by fixing the head at the level observed at the end of the cycle and simulating the subsidence that would result, with no further changes in head, in 2017. The results, shown in Figure 2a, reveal that the deferred subsidence incrementally increased with each cycle. This cumulative effect is due to the long compaction timescales associated with the slow diffusion of head through thick clay interbeds and aquitards. For example, Figure 2b shows the 2017 head profile through a 13 m thick clay. The head levels in a 4 m thick zone in the center of the clay would still be greater than (so unequilibrated with) those in 1987, and for each cycle of head decline and recovery, the thickness of this zone increased. The cumulative effect of repeated cycles of head decline and partial recovery in a system governed by very long timescales has led to extremely high rates of deferred subsidence by 2017.
These high levels of deferred subsidence, ongoing today and previously unrecognized, represent a significant concern for groundwater managers for two main reasons. Firstly, they represent long-term subsidence rates which, in the absence of head recovery, will continue for decades into the future. Secondly, they are now established as the baseline so that any future head decline will add additional subsidence. Such high baseline rates mean the system is highly vulnerable to extremely large and damaging subsidence rates in the event of further head decline. This is of particular concern as the San Joaquin Valley enters another period of drought in Spring 2021. The high deferred subsidence rates in our study region are likely to be present widely across the Valley, where similar cycles of head decline and partial recovery have been observed for decades, and represent an acute challenge in the sustainable management of groundwater.

Depth of compaction
The final important insight provided by our modeling approach is the depth at which compaction is occurring. Our model divided the aquifer system into the three major hydrostratigraphic layers: the unconfined to semi-confined upper aquifer, the Corcoran Clay and the confined lower aquifer, and we simulated the compaction in each. The apportionment of compaction between these three layers is shown in Figure 3. These data fill an important data gap left by the sparse extensometer network. Figure 3 is separated into two time periods, 1955-1980 and 2000-2017. In both periods, compaction of the Corcoran Clay contributed less than 10% of the total subsidence, and the majority of the subsidence occurred in the lower aquifer. This is consistent with the extensometer measurements from the Poland era as well as previous numerical experiments 3,34,35 . However, an unexpected result was the relatively large amount of subsidence, up to 30%, which originated in the upper aquifer from 1955-1980. During this period, lower aquifer head levels were relatively constant whereas upper aquifer head levels showed a net decline. While it has often been assumed that subsidence is only a problem at locations where the Corcoran Clay is present, due to the confined conditions it creates in the underlying aquifer, our results show that under conditions of significant head decline, it is possible for substantial subsidence to originate within the unconfined-to-semi-confined upper aquifer. Our modelling approach could be generalized to other locations in the Valley, providing an alternative to the expense of installing extensometers and contributing critical information on the depth of compaction.

Conclusions
The sustainable management of groundwater can be guided by insights gained from observations of the past. We have used a one-dimensional numerical compaction model to obtain a 65-year record of head and simulated subsidence from 1952-2017 in a study region in the San Joaquin Valley. Our model provides insights into the long timescales over which subsidence occurs and highlights the presence of extremely high deferred subsidence rates. While we saw that the lower aquifer was the main source of subsidence, our model was able to highlight the capacity for both aquifers to contribute subsidence under conditions of sustained head decline. As the San Joaquin Valley enters another period of drought, and groundwater managers begin implementing their SGMA plans, our findings provide the critical quantitative insight into the causes of subsidence required to inform the sustainable management of groundwater.

Methods 1: Preparation of head inputs
In the San Joaquin Valley, it is generally difficult to construct long-term head time series corresponding to the upper and lower aquifers. The primary challenges are that well completion reports are seldom available, wells are often screened across both aquifers, and measurements are often either privately owned or contain unmarked errors 36 . We had access to a large amount of groundwater data and were therefore able to resolve these challenges in our study region. The available datasets were: the statewide well data from the California Statewide Groundwater Elevation Monitoring (CASGEM) database 37 ; a complementary local database of water level measurements in the Kaweah groundwater subbasin, compiled by GSI Water Solutions and provided with the permission of several local groundwater agencies; water-level measurements from Poland et al. (1975) 9 ; and water-level data from a dedicated monitoring well Well H, situated 2 km northeast of our study region, which had 4-hourly measurements of head in the upper and lower aquifers for the period 2015-present and were provided by Kings County Water District.
We categorized wells in our study region into four classes: gold standard upper aquifer wells, gold standard lower aquifer wells, supplementary upper aquifer wells and supplementary lower aquifer wells. Gold standard wells were those for which we had very high confidence of which aquifer the head measurements represented. To that end, gold standard upper aquifer wells were defined as those with available well completion information, well depth shallower than the Corcoran Clay, and seasonal amplitude variations similar to upper aquifer variations at Well H. Gold standard lower aquifer wells were those with well completion information available, perforated intervals exclusively below the Corcoran Clay, and seasonal amplitude variations similar to those in the lower aquifer at Well H. Corcoran Clay depths were determined from the borehole data as described in Methods 3. The two gold standard upper aquifer and two gold standard lower aquifer wells which we identified in our study region are highlighted in Table S 1. Unfortunately, using only these four wells would have left a head timeseries with gaps in the upper aquifer between 1978-2010 and gaps in the lower aquifer between 1981-2006. We identified supplementary wells to fill these gaps. Supplementary wells were additional wells which did not have completion information available but satisfied other criteria to give us confidence of the aquifer which they were recording. Supplementary wells were defined as those with head measurements which showed overlap with those from gold standard wells, and which showed seasonal variability similar to the corresponding aquifer at Well H. We identified two supplementary wells in each aquifer which, when combined with the gold standard wells, gave a continuous record of measurements from 1952-2017, shown in Figure  To convert the well head measurements into model input for each aquifer, we extracted a Fall minimum and Spring maximum for each year, then linearly interpolated between those measurements. We excluded Well UA2 from this process due to anomalously low head measurements in 2007-10, and likewise removed two measurements taken in 1995 and 1976 from UA3 which were also anomalously low (Figure S 4). For each year, we defined the Fall minimum as the minimum datapoint between the months of August and November (inclusive), and we defined the Spring maximum as the maximum datapoint between the months of January and May (inclusive); in the extremely rare case that no data were available within Spring or Fall, we used the measurement from the previous year (for instance, Upper Aquifer Fall 1977 and Fall 1976 are identical for this reason). We tested the validity of this approach by applying it to the continuously monitored data at Well H, and comparing subsidence simulated using our approach as input to that simulated using the "true", continuously-recorded measurement as input. The differences, which can be seen in Figure S 5, were relatively minor.
The final head timeseries used as model input are presented in Figure S 6.

Methods 2: Model equations and numerical solver
The conceptual model of the subsurface, which forms the basis for our numerical model, is shown in Figure S 3. Our numerical implementation of the aquitard-drainage model follows that of Helm (1975) 38 , and is summarized only briefly here. Clays appear in our model in three places: upper aquifer interbeds, lower aquifer interbeds and the Corcoran Clay aquitard. Within each clay layer, we solve the governing equation for effective stress: where ′ is the effective stress, Kv is the vertical hydraulic conductivity and Ss is the specific storage, related to skeletal specific storage by the relation = + , where Ssw is the specific storage due to the compressibility of water. A detailed explanation of the concepts of specific storage can be found in p. 4-5 of Sneed (2001) 39 . Equation M1 is solved using a finite difference scheme; our model code is available here: https://github.com/mlees0209/compaction_model_StaticPaperVersion. The initial conditions for each clay layer are discussed in Methods 5. The effective stresses of the two aquifers, which act as boundary conditions and drive the system's dynamics, are time-varying and are computed from the observed head in the aquifer layers using the relation ′ = + ℎ, in which is the overburden stress, is the density of water and is the gravitational acceleration. Overburden stress is calculated based on the upper aquifer head according to = ℎ where Sy is the specific yield; this relationship is concisely explained in p. 5-6 of Poland et al. To apply Equation M2, our model keeps track of effective stress history for each node and updates the value of Ssk accordingly. Helm (1975) 38 , which is reproduced here, modified slightly to match our notation:

Calculation of compaction for each clay layer follows Equation 8 in
where ( ) refers to the total compaction of the clay at timestep i, is the mesh spacing, and there are J nodes. The nodes used are the midpoints of the nodes used to solve Equation M1; values of effective stress at the midpoints are obtained through linear interpolation. Superscripts on effective stress refer to the value at a given timestep. Note that, in the condition that effective stress at the current timestep is equal to its maximum, the Sske term equates to zero and compaction is calculated purely as the inelastic sum.
In the coarse-dominated units, effective stress history was calculated as described previously. Compaction in these units was calculated using Equation 3 with Sskv = 0, an elastic value of Sske for coarse-dominated material and a single 'node' of thickness Δ equal to the aquifer layer thickness minus the total thickness of clay interbeds within the aquifer.
The values used for all hydrologic parameters are provided in see Methods 6. The hydrostratigraphy of the subsurface is described in Methods 3 and 4.

Methods 3: Estimating the number and thickness of subsurface clays and the top and bottom boundaries of the Corcoran Clay
Borehole data, in the form of a driller's log and a resistivity log, were both obtained at Well 19S22E08D002M. These logs were provided to us by Kaweah Delta Water Conservation District. The borehole was 210 m deep and both logs spanned its entire depth. The driller's log contained generalized lithologic descriptions of the subsurface at 5-foot intervals but did not describe individual clay interbeds; an illustrative extract from the driller's log is shown in panel a) of Figure S 7. The resistivity log was a 16-inch normal log which we interpreted presuming that clays are the least resistive layers.
To obtain the number and thickness of subsurface clays and the boundaries at the top and bottom of the Corcoran Clay, we applied resistivity cutoffs to the resistivity log. These resistivity cutoffs were values below which we assumed material to be clay. To choose appropriate cutoff values, we applied trial cutoffs and compared the result with our identification of the Corcoran Clay and thicker clay lenses from the driller's log. We determined that the smallest and largest cutoffs which gave good agreement were 30.5 and 38.5 ohm-m respectively, shown in Figure S 7 panels b) and c). There was no variation in the depths to top or bottom of the Corcoran Clay across this range, which were 108 m and 135 m respectively. However, there was variation in the number and thickness of clay interbeds in each aquifer. We retained these variations in our grid search, and included two intermediate cutoffs of 33 and 36 ohm-m. The four resulting clay distributions are shown in Figure S 7 panels d-g. One complication was that, for most simulations, the lower aquifer extended below the base of Well 19S22E08D002M (see Methods 4 for our approach to estimating the aquifer system base). In these cases, on the basis that there is no evidence for a change in depositional environment below the log, we scaled the number of clays from within the log to correspond to the lower aquifer thickness. For instance, the log covered a 50 m thick region of the lower aquifer; in a simulation with a lower aquifer thickness of 100 m, we doubled the number of clay interbeds identified in the log.

Methods 4: Assigning depth to aquifer system base and defining top of upper aquifer
We defined the lower aquifer as the region below the Corcoran Clay over which significant head changes occurred. Although the deepest pumping well in our study region is 322 m deep, we expect head changes to propagate for some depth below this. This propagation depth will depend on the hydraulic conductivity of the lower aquifer and the time scale being considered. From resistivity logs, an impermeable layer called the San Joaquin Formation has been identified in our study region at about 800 m below the surface (see, for instance, Figure 14b in the Mid Kings River GSA 32 ), which we consider to be the maximum possible depth to which head changes propagate. However, to accurately determine the propagation of head changes in the region between the base of wells (322 m) and the base of permeable units (800 m), we would require a detailed hydrogeologic model in a region for which we have neither head measurements for validation nor well data to describe the subsurface structure. Instead, we made the simple approximation that head changes in the lower aquifer propagated for some depth-ofinfluence di below the deepest pumping well, as annotated in Figure S 3. Since we could not reliably know di, we allowed it to be a variable in our grid search. The base of the lower aquifer in our model was therefore defined as zbase = zdeepestwell + di. Unlike in the previous hydromechanical modelling studies of Helm (1975) 38 and Smith and Knight (2019) 23 , this approach gave us a temporally-varying base of the lower aquifer, since zdeepestwell has increased as more wells were drilled over time, as shown in Figure S 8. Data on the depths of wells drilled over time came from the California Department for Water Resources' Online System for Well Completion Reports (OSWCR), and is available at: https://data.cnra.ca.gov/dataset/wellcompletion-reports.
In our grid search, we let di take the values [50, 100, 150, 200] m. This gave a maximum base of the lower aquifer as 522 m, which is similar to the range used by Smith and Knight (2019) 23 .
We defined the upper aquifer as the permanently saturated region throughout our modelled time period. We therefore defined the top of the upper aquifer to be the upper aquifer head in 2017, which was the lowest head during our modelled period, and was 41 m depth. This is illustrated in Figure S 3.

Methods 5: Assigning initial hydraulic head in clay layers
The initial hydraulic head condition within clays interbeds and aquitards cannot be measured, but it can be estimated based on the prior history of head in the aquifers 38 . Given that the Poland-era record reveals that subsidence rates only became large at our study region in 1953 9 , we assumed that the initial head profile in the clay interbeds (ℎ 0 and ℎ 0 ) were constant within each aquifer. To select the values for ℎ 0 and ℎ 0 we considered what was known about head levels pre-1952.
In the upper aquifer, we identified CASGEM well 33014 within the study region which contained head data from 1925-1963. This well showed head approximately constant at a value of 68 masl from 1925-1948. From 1948, it dropped steadily to the model start value in 1952 (55 masl). We therefore bounded ℎ 0 between 55 masl and 68 masl. For our grid search, we allowed values of 55, 61 and 68 masl.
In the lower aquifer, we were unable to identify any wells providing head information prior to 1952. Instead, we used the literature to give us a picture, albeit limited, of lower aquifer head levels. Mendenhall (1916) 41 documented artesian wells tapping the deep aquifer in the study region in 1905; hence we know that lower aquifer head was at least the level of the land surface, or 71 masl, in 1905. By 1925, Harding (1927) 42 reported that artesian wells in the valley had stopped flowing apart from in very wet winters. We do not, however, have any constraint on the deep aquifer head between 1925 and 1952. We therefore bounded ℎ 0 as being between the 1952 level (39 masl) and the land surface (71 masl). For our grid search, we allowed values of 39, 49, 54, 60, 65 and 71 masl.
For the Corcoran Clay aquitard, we set the initial head ℎ 0 − to be constant, and equal to the average of the head in the upper and lower aquifers.

Methods 6: Selection of input hydrologic parameters
In the model equations described in Methods 2, there are seven hydrologic parameters: vertical hydraulic conductivity of clays Kv, elastic skeletal specific storage of clays and of coarsedominated material (denoted here as and respectively), inelastic skeletal specific storage of clays Sskv, the specific storage due to the compressibility of water Ssw, the density of water and the acceleration due to gravity g. Three of these (Ssw, and g) are well-known physical constants so did not vary in our grid search. Additionally, since ≪ and ≪ , we did not expect or to have significant influence on our simulated subsidence, and we therefore did not vary them in our grid search. The values used for these five, nonvarying parameters are given in the top part of Table S 2.
The remaining two variables (Kv and Sskv) were important to our simulations and we varied them in our grid search. Our approach was to define the range of values reported in the literature and run simulations which covered a large portion of that range. We gave extra weighting to values from previous models, as opposed to values from laboratory experiments, as it is known that values do not agree at these two scales 39 . The literature ranges and the values used in our simulations are given in the bottom part of Table S 2.

Methods 7: Validation of model outputs
The four validation datasets, and the validation criteria, are described below.

2015-17 InSAR data
The 2015-17 InSAR data were acquired with the European Space Agency Sentinel satellite and processed by TRE Altamira under contract with the California Department of Water Resources (DWR), who provided them to us. A detailed description of how the data were processed can be found in TRE ALTAMIRA Inc. (2020) 43 , and a version of the data is publicly available through the DWR website: https://data.cnra.ca.gov/dataset/tre-altamira-insar-subsidence. The raw data provided measurements of vertical displacement every 7-14 days, starting in February 2015 and continuing until the end of our study period, at the pixel locations seen in Figure S 9a and b. Since our model input head was based on a Fall minimum and a Spring maximum for each year (i.e., temporal resolution of six months), we degraded the temporal resolution of the validation data from 7-14 days to annual by using October-to-October values. This allowed us to validate our simulated subsidence over two water years: Water Year 2016 and Water Year 2017. The validation criteria for these two years were: 1) Water Year 2016: simulated subsidence between 20 and 34 cm, which was the range of measured subsidence at the pixels in our study region (Figure S 9a). 2) Water Year 2017: simulated subsidence between 9 and 15 cm, which was the range of measured subsidence at the pixels in our study region ( Figure S 9b).

2007-10 InSAR data
The 2007-10 InSAR data were the ALOS data described in Smith and Knight (2019) 23 . Although the raw ALOS data contained 15 measurements from between June 1 st 2007 and July 2 nd 2010, we elected to use the mean rate of subsidence over that period for our validation because this gave a more reliable estimate of the actual subsidence than the individual measurements, which likely contained some noise. The validation criteria rate between these dates was a mean subsidence rate between 7.0 and 16 cm/yr, which was the range of measured subsidence in our study region. (This range is shown visually in Figure 44 . Considering that "the 1960s" could refer to anytime between 1960 and 1969, we converted 1.7 m into a range of deformation rates from 3.9-5.0 cm/yr, equivalent to 1.7 m over 44 or 34 years. The survey point is just north of our study area, as shown in Figure S 1. At the survey point, measurements from the two InSAR datasets record subsidence rates at 75% of the rates within our region. To account for the fact that our study area likely experienced slightly higher subsidence rates than the measurements at this survey point, we extended the bound slightly, and therefore applied the following validation criteria: the simulated subsidence between 1965 and 2004 was required to be between 90% and 160% of that recorded by the Highway 198 survey point; that is, it had to have a mean rate over that period between 3.9 and 6.9 cm/yr.

Poland-era survey data
The Poland-era survey data are a series of measurements made during the years 1933, 1947, 1953/4, 1957, 1958, 1962, 1966 and 1970 as reported in Poland et al. (1975) 9 and the Inter-Agency Committee on Land Subsidence in the San Joaquin Valley (1958) 16 . They show low rates of subsidence until 1953/4, and high rates of subsidence from then until the end of the surveyed period in 1970. Since the early part of our simulations is more strongly influenced by the initial conditions, we only used the 1966-1970 measurement as model validation. The survey recorded an average subsidence rate of 11 cm/yr during those four years. However, the survey location is slightly south of our study region (see location in Figure S 1). At the survey location, measurements from the two InSAR datasets recorded subsidence rates 125% of the rates measured within our region. To account for the fact that our study area likely experienced slightly lower subsidence rates than the measurements at this survey point, we extended this bound slightly, and applied the following validation criteria: the simulated subsidence rate between 1966 and 1970 was required to be between 65% and 100% of that recorded by the survey point; that is, it had to have a mean rate over that period between 7 and 11 cm/yr. where we use the values of Sskv and Kv for that particular simulation. To calculate the equivalent thickness beq, we follow the approach of Helm (1975) 38 , using a version of their equation 9, modified to fit our notation:

Methods 8: Calculating characteristic time constant of compaction for a simulation
where there are i unique interbed thicknesses ti in each model, and ni is the number of each thickness. To be consistent with the results of Helm (1978) 25 , we do not include the Corcoran Clay in this calculation, hence only consider interbeds; the time constants including the Corcoran Clay would be slightly larger. Figure S 1