Fine-scale spatiotemporal variation in subsidence across California's San Joaquin Valley explained by groundwater demand

We adopt an idealized conceptual model to motivate and structure our statistical analysis. The model integrates our understanding of the mechanics driving land surface displacement with the fundamentals of surface water hydrology (Supplementary Text S1 (stacks.iop.org/ERL/15/104083/mmedia)). We consider displacements dD/dt to be the sum of two responses: a local poroelastic displacement response dD_p/dt, which arises from local changes in pore pressure caused by change in groundwater storage [40] (stores of water in subsurface pore spaces); and a broad, potentially non-local elastic load response dD_e/dt, which arises from changes in the downward force of stored surface and subsurface water on the solid earth, reflecting motion opposite to that induced by poroelastic effects [41–44]. Thus, displacements are related to groundwater storage (equations (S1)–(S2)). Where groundwater is pumped from confined or semi-confined aquifers, as in the San Joaquin Valley [16, 24], the elastic load effect is overwhelmed by the poroelastic effect [42]. We therefore assume that observed displacement reflects poroelastic effects over regions of active groundwater pumping ($dD/dt = dD_p/dt$) and elastic load effects elsewhere ($dD/dt = dD_e/dt$). Artificial groundwater storage, or managed aquifer recharge, is limited in our study region and therefore not considered in our analysis.

We describe local change in groundwater storage using a standard water balance framework that relates local storage change dS/dt to the difference between surface water inputs (precipitation, irrigation) and outputs (evapotranspiration, discharge) (equation (S3)). We combine the water balance model with our geophysical understanding of the displacement response, and consider idealized conditions representative of drought conditions in the San Joaquin Valley [36, 37]. This yields a simple proportionality between displacement and vegetation water demand, or the difference between evapotranspiration ET and precipitation P: $dD/dt = -\alpha (ET - P)$ (equations (S4)–(S6)), where α is a proportionality constant. Positive water demand results in negative displacement, or subsidence. When surface water deliveries are present, such as in interim seasons or wet periods, we expect the opposite: either reduced rates of subsidence or, in some cases, uplift via recharge.

This conceptual model relates groundwater storage change to displacement using a proportionality constant (α), which is a function not only of the poroelastic and elastic responses, but also various aspects of the subsurface (e.g. geology, soil composition) [25] and land use practices (e.g. soil management, irrigation methods) (equations (S7) and (S8)). Due to this heterogeneity in surface and subsurface features across the study region, we relied on a large-sample statistical approach that was purposefully designed to average over this heterogeneity [45], and thereby illustrate key dynamics present across the study region. Thus, our empirical analysis is focused on exploiting variation in high-resolution observational data (section 2.2), through which we explore the capacity of a simple linear relationship (section 2.3) to reveal fundamental features of an otherwise complex connected surface and groundwater system.

Our primary data include vertical land surface displacements, land cover classifications, and vegetation water demand for irrigation ('water demand'). Our vertical land surface displacement data are observations at a 100 m x 100 m spatial resolution and 6 - 48 day (19 day average, 24 day mode) temporal resolution. Displacement data are a combination of two independent and complementary geodetic datasets: synthetic aperture radar (SAR) acquisitions from the European Space Agency's Sentinel-1 missions [28, 46–48], which are corrected using available continuous Global Positioning System (cGPS) stations [49] following the cGPS-enhanced InSAR method [35] (Supplementary Text S2.1). We formatted 50 spatiotemporal displacement estimates (Supplementary figure S1 and Supplementary table S1) between April 1, 2015 and October 23, 2017 as cumulative vertical land surface change relative to April 1, 2015 (mm, 'total' displacement), and rates of change (mean mm/day within 6-48 day measurement periods; displacement 'rate'). Our land cover data are classifications [50, 51] grouped into seven categories aligned with those used in land use surveys [52] and with irrigation practices [53] (Supplementary figure S2). We limit our analysis to areas whose land cover types remain unchanged over our study period, and consider areas with oil wells [54] separately. These land cover types and their percentage of the valley (1.2 million grid cells total), are as follows: fruit and nut crops, e.g. deciduous, citrus, and vineyard crops (33%); native vegetation (30%); urban land (13%); idle land (10%); field crops, e.g. corn, cotton, and soy (6%); pasture crops excluding pasture grasslands, e.g. alfalfa and switchgrass (4%); oil fields (3%); and truck crops, e.g. tomatoes, lettuce, and melons (1%). Because we lack fine-scale information on water supply, we generate an estimate of water demand as the difference between potential evapotranspiration and precipitation (PET − P). PET is the product of daily 2 km resolution reference evapotranspiration [55] (http://cimis.casil.ucdavis.edu/cimis/), and vegetation-specific coefficients [56] matched to the annual 30 m land cover data [50, 51]; from PET, we subtracted daily 4 km precipitation [57] (Supplementary Text S2.2). Additionally, we constructed a measure of proximity to artificial surface water supply infrastructure (e.g. canals) (Supplementary Text S2.3 and Supplementary figure S3) [58, 59]. Ancillary data include: soil infiltration characteristics [60]; elevation [61]; well water table depth measurements [62]; and annual groundwater use amounts reported by agricultural water districts [63, 64].

We performed analyses at the spatiotemporal resolution of the displacement data by using modal and median resampled values of land cover and water demand. Our analysis exploits variation across two dimensions. First, temporal differences between dry and wet time periods figure 2(a) represent extremes in surface water supply [36]. We defined the dry period to be April 1, 2015 - October 4, 2016 (22 observations) and the wet period to be October 5, 2016 - October 23, 2017 (28 observations). The dry period spans more than a single water year (October - September of either 2015-2016 or 2016-2017), due to our inclusion of additional April-September 2015 displacement observations to more closely balance the number of wet and dry period observations. Second, spatial differences between land cover types figure 2(b) represent differences in water demand. We split the study region into broad cultivated and uncultivated land regions using the extent of the cultivated alluvial valley [65] as a boundary (figures 1(a) and (b)), and 2(b). Where noted, we use either the full study region or the subset cultivated valley, as well as the land cover types therein.

Our analysis proceeds in three parts. First, we evaluate how displacement dynamics are differentiated by land cover. We provide summary statistics of displacement for different land cover types for the entire study period, and wet and dry year periods separately. We also calculate the correlation between displacement and water demand rates at the grid cell level across the full study region, which identifies unique temporal dynamics of displacement and water demand within areas with either strong positive or negative correlation, and identifies the land cover types that generate those dynamics. See Supplementary Text S3.1 for details.

Second, we explore the capacity of water demand to directly explain variation in displacement. We empirically model displacement as a linear function of water demand in the cultivated valley at the grid cell level, and for wet and dry time periods separately. This analysis yields estimates of grid-cell level mean displacement rates, mean change in displacement with a change in water demand, and the the fraction of variation in displacement explained by water demand. See Supplementary text S3.2 and equation (S9) for details.

Third, we evaluate how displacement response to water demand varies by land cover. We again model displacement as a linear function of water demand for separate wet and dry periods, but at the land cover type group level. This analysis yields estimates of land cover type- and time period-specific mean displacement rates, and mean change in displacement with a change in water demand. Lastly, we model displacement as a linear function of water demand interacted with land cover type and a measure of proximity to surface water supply infrastructure (Supplementary Text S2.3 and Supplementary figure S3). This allows simulation of displacements for each land cover type across a range of water demand rates, as well as for low and high levels of access to surface water supply infrastructure. See Supplementary Text S3.3 and equations (S10)–(S12).

Our attribution of heterogeneity in displacement to land cover and associated groundwater demand faces several possible sources of confounding. Therefore, we explore the capacity for variation in displacement with land cover to be attributable to variation in other features of the environment, scale, or model specification (Supplementary Text S3.4). Additionally, while there exist records of groundwater levels and withdrawals in the San Joaquin Valley, those that are publicly accessible are not temporally and/or spatially consistent across our large study region. Therefore, we only briefly explore the capacity of displacements to be used directly alongside a limited set of well water table depth measurements (Supplementary Text S3.4.5) and annual groundwater use amounts reported by agricultural water districts (Supplementary Text S3.4.6).

Statistical summaries over the cultivated valley reveal meaningful differences across land cover types, including different crops. Cumulative totals and average rates of displacement, expressed as the median of grid cell observations within each land cover type for the full time series, vary significantly across land cover types (figure 3 (a, b)). The largest cumulative magnitudes of subsidence occurred in the dry year and for field crops, followed by pasture crops, truck crops, and fruit and nut crops (figure 3(c)). For example, there was a median 272 mm of total cumulative subsidence for field crops over the study period, and a dry water year subsidence rate of 131 mm/year. For fruit and nut crops, there was a median 62 mm of total subsidence over the study period, and a dry water year subsidence rate of 31 mm/year. In contrast, for uncultivated lands (idle, native, urban, and oil), subsidence totals over the study period were between 4 and 38 mm, and dry water year rates of subsidence were between 5 and 14 mm/year. The lesser subsidence of fruit and nut crops relative to other crops is unexpected given high reported volumes of water consumption by fruit and nut crops in California. The area of land cultivating fruit and nut crops, which is large relative to field crops (33% vs. 6% of the cultivated valley), can nevertheless consume a greater total volume of water. For all cultivated lands, the subsidence rates experienced in the dry year were more than double that of the wet year, whereas dry and wet year subsidence rates for uncultivated lands differed only by up to 9%. Agricultural land types experienced an additional 20-88 mm of subsidence in the dry year relative to the wet year, indicating how much subsidence may be attributable to surface water deficits in a drought year.

Zoom In Zoom Out Reset image size

In general, land cover types with the greatest relative subsidence are those with higher positive water demand (figure 3(b)). Water demand is greatest for field crops, followed by truck, pasture, and then fruit and nut crops (figure 3(b), top). In general, we expect land cover types with greater relative water demand (e.g. field crops vs. fruit and nut crops) to have greater groundwater withdrawal and associated subsidence. However, differences in the amount of water applied as irrigation between agricultural land types are determined both by differences in physiological water demand and the efficiency of irrigation (section 2.1 and Supplementary Text S1). Field and pasture crops characteristically use irrigation methods that are less efficient and higher-volume (gravity-driven methods such as flood and furrow) than those used by fruit and nut crops (micro-sprinkler and drip irrigation) [66]. Thus, irrigation methods can increase or decrease the relative intensity of irrigation water use beyond what would be expected according to water demand alone. For example, pasture crops have only a slightly greater water demand than fruit and nut crops (figure 3(b), top), but are generally irrigated less efficiently [66]. This may explain the relatively greater subsidence during dry periods over pasture crops than water demand alone suggests (figure 3(b), bottom).

Correlation between time series of land surface displacement and water demand (figure 4 (a)) illustrates seasonal (phase) alignments between displacements and water demand that we expect based on the conceptual water balance model, and the presence or absence of poroelastic effects in agricultural and non-agricultural land, respectively (section 2.1 and Supplementary Text S1). Negative correlation indicates displacement that is out of phase with water demand (figure 4(b), blue lines), which is expected in areas where groundwater pumping results in subsidence during the dry summer season (during positive water demand). This pattern appears predominantly over cultivated land (figure 4(c), top). Positive correlation indicates that displacement is in phase with water demand (figure 4(b), orange lines), a pattern that occurs almost exclusively over rain-fed native vegetation (figure 4(c), bottom). Furthermore, the amplitude of the out-of-phase displacement signal, which corresponds to agricultural land cover, is greater than that of the in-phase non-agricultural displacement signal. Thus, cultivated areas are distinguishable from uncultivated areas according to the relationship between displacement and water demand. Additionally, the spatial extent of strong negative correlation increases in dry relative to wet periods (Supplementary figure S4), suggesting an observable increase in the spatial extent of agricultural reliance on groundwater during the dry period.

Zoom In Zoom Out Reset image size

We used a grid cell—level regression analysis (equation S9) to evaluate the extent to which spatial variation in water demand corresponds to spatial variation in displacement, across all land cover types within the cultivated valley. We assess spatial patterns in the direction and magnitude of mean displacement (figure 5(a), (d)), the marginal response of displacement to water demand figure 5(b), (e), and the degree to which water demand variation explains displacement variation figure 5(c), (f). Differences between dry (figure 5 top) and wet (figure 5 bottom) periods illustrate differences between cases of limited and abundant surface water supply in the growing season, respectively. Mean displacement, the marginal displacement response to water demand, and explained variance were all of greater magnitude in the dry period, suggesting greater relative reliance on groundwater throughout the valley. Meanwhile, there was substantial spatial variability in this water demand - displacement relationship. For example, negative displacement (subsidence) and displacement sensitivity to water demand were reduced but still present in the wet period in southwestern Tulare county (dark red and purple regions in figure 5(a) and (b), respectively), while there was uplift and greatly diminished displacement sensitivity in the Westlands region of southwestern Fresno county (dark gray and light purple regions in figure 5(d) and (e), respectively).

Zoom In Zoom Out Reset image size

In terms of the regression design, mean displacement (the intercept) and the marginal response of displacement to water demand (the slope) are independent from one another. Therefore, these results imply that even when local averages are accounted for, negative marginal responses of displacement to water demand (greater subsidence with greater water demand) have an extremely coherent spatial structure consistent with intensively farmed areas of the valley. Furthermore, the percent of variance in displacement that is explained by water demand (R-squared) is exceptionally high—up to 95%—in these same areas of the valley, even though displacement is not independent from other features of the environment, such as geology and soils, which are not accounted for in the empirical model. In the dry period, in areas with R-squared values ranging between 70%–95% (orange and red in figure 5(c)), water demand is the primary, or exclusive, driver of subsidence. In areas with lower R-squared values, including cultivated areas in the wet period (figure 5(f)), remaining unexplained spatial variability is likely due to surface water supply (specifically in the wet period), or other unobserved geophysical drivers of heterogeneity (e.g. subsurface geology). Discovery of these features is beyond the scope of the present study, but would be a fruitful area for future research, particularly with respect to understanding wet period dynamics.

First, we used a land cover type—level regression analysis (equation S10) to evaluate the extent to which land use—specific water demand explains land use - specific displacement in the valley, and within separate dry (figure 6(a)) and wet (figure 6(b)) periods. Average dry period rates of subsidence (figure 6(a), horizontal axis), and the average marginal response of displacement to change in water demand (figure 6(a), vertical axis), are both greater for those land cover types that also demonstrated greater cumulative subsidence and rates over time (figures 3(a) and (b)). For example, for field and pasture crops, average rates of displacement are not only more negative, but the marginal displacement response to increased water demand is also more negative, relative to other land cover types. Furthermore, the magnitudes of subsidence and marginal displacement response to water demand were reduced for all land cover types during the wet period (figure 6(b)), when surface water supplies were available.

Zoom In Zoom Out Reset image size

Second, we used fitted regression coefficients to simulate displacement across a realistic range of water demand (equation S11, Supplementary figure S5), and across both levels of water demand and levels (5th and 95th percentiles) of proximity to surface water supply infrastructure (equation S12), again within separate dry (figure 6(c)) and wet (figure 6(d)) periods. Without accounting for proximity to surface water infrastructure (equation S11), we observe that subsidence rates increase with water demand for agricultural land cover types in the dry period, as well as for select agricultural lands in the wet period, but to a lesser degree (Supplementary figure S5). Thus, agricultural lands in general are more sensitive to water demand in dry periods, but only select agricultural land cover types (field, pasture crops) remain sensitive in wet periods.

Accounting for access to surface water supply infrastructure (equation S12), which is done using a measure of the intensity of proximal artificial surface water features such as canals (Supplementary figure S3), reveals additional dynamics. The use of the surface water source intensity measure as a proxy for surface water supply was necessitated by the absence of spatially-explicit measurements of surface water deliveries at the resolution of our displacement or water demand datasets (Supplementary Text S2.3). The 5th percentile of water source intensities took a value of zero, indicating areas where no grid cells in the surrounding 5 km region contained artificial surface water supply infrastructure features; the 95th percentile took a value of 0.22, indicating areas where 22% of grid cells in the surrounding 5 km region contained an artificial surface water supply feature. In the dry period, subsidence increased with water demand for agricultural land cover types in areas with both low and high artificial water source intensities (figure 6(c)), when the measure of those intensities did not indicate access to surface water, as deliveries were limited or non-existent in that period. In contrast, in the wet period, when surface water deliveries returned, subsidence only continued to increase with water demand in areas with low artificial water source intensities—or in areas lacking surface water supply infrastructure. We note that figures 6 (c) and (d) present a subset of land cover types in order to highlight patterns specific to agricultural lands. Displacements for oil fields, a non-agricultural land type, should not be sensitive to the proximity of surface water supply networks, other than through the collocation of those oil fields with agricultural land that is sensitive to water supply. As anticipated, the slope for oil fields with low proximity to surface water infrastructure does not change between dry and wet periods as it does for agricultural land; the slope for oil fields with high proximity is consistent with collocated agricultural land.

This analysis demonstrates that the sensitivity of agricultural land subsidence to water demand is diminished with proximity to surface water supply networks, but only in wet periods when surface water is in fact accessible through those networks. This supports our understanding that water demand satisfied by pumping groundwater drives subsidence, and does so variably across different land cover types.

All robustness analyses generated results consistent with our primary findings. Alternative regression model specifications demonstrate that systematic differences in displacement and displacement sensitivity to water demand by land cover type are robust to model specification (Supplementary figures S6 and S7, Supplementary table S2, and Supplementary Text 3.4.1). For example, even when accounting for variation in proximity to artificial surface water supply infrastructure (Supplementary figures S3 and S8), soil infiltration characteristics [60], elevation [61], longitude and latitude, the field and pasture crop displacement responses in the dry period remain more negative than those of other land cover types (figure S6). Thus, key differences in displacement sensitivity to water demand by land cover type are not explained by those other features. Despite deviation in some estimates (e.g. we obtain variable results for truck crops and idle land, for which there are fewer samples and uncertain classification, respectively), it remains clear that subsidence rates are more sensitive to increased water demand, and subsidence is greater for agricultural land cover types, and in dry relative to wet periods.

While agricultural land cover types are relatively well-mixed spatially across the larger study region, there exists some degree of spatial segmentation by land cover type, which motivated analysis within a subset area of the valley. For the subset analysis, displacement remains differentiated by land cover type, with field and pasture crop subsidence totals and rates exceeding those of fruit and nut and other non-agricultural land cover types (Supplementary figure S9 (b, c)). In our analysis accounting for spatial correlation in displacement, we replicated our main descriptive analyses using displacement data adjusted for spatial correlation, and also replicated our land cover type-level regression analysis through bootstrap sub-sampling of the study region. In both analyses, the relative relationships between displacement for different agricultural land cover types, and specifically negative displacement (subsidence) during dry months, were retained (Supplementary figures S10 and S11). When incorporating land cover classification error into the main descriptive analyses, the magnitudes of displacement for land cover types varied. However, the relative rankings of agricultural land cover types remained the same (Supplementary figure S2 and Supplementary table S3).

Through comparison of displacement and monitoring well water table depth data (Supplementary Text S3.4.5), we found that displacement observations extracted at well locations were more correlated with each other than were the water table depth observations at those same well locations. However, temporal patterns in displacement were generally consistent with well water table depths (albeit lagged, as observed in previous research [18]). Both the well and displacement data similarly distinguished between agricultural and non-agricultural land types (Supplementary figure S13). It is nevertheless apparent that the information provided by point-located well data is fundamentally different in nature from that of gridded InSAR, despite general correspondence. In our analysis of annual groundwater use amounts reported by agricultural water districts that were incorporated into Groundwater Sustainability Agencies (GSA) following the passage of California's Sustainable Groundwater Management Act (SGMA) in 2014 [67], we found that volumes of displacement over GSAs were a fraction of the total groundwater withdrawal volumes, and that this fraction varied substantially between GSAs, and over time within individual GSAs (Supplementary Text S3.4.6). In 2015 and 2016, annual total volumes of displacement over GSAs ranged between 13% and 21% of reported groundwater use by volume, and displacement patterns varied substantially across adjacent agencies (Supplementary figure S14). Displacement is a function not only of groundwater withdrawals (themselves a function of irrigation demand and water use efficiency), but also long-term recharge rates, aquifer properties, and subsurface geology—all of which vary in space and time [12, 16, 18, 24, 25, 68]. Nevertheless, comparisons such as these demonstrate opportunity for evaluating subsidence data alongside other datasets and models of the aquifer system, which have the capacity to benefit local-level monitoring efforts, the assessment of permanently reduced storage capacity, and the identification of opportunities for managed aquifer recharge.