Research ArticleGEOLOGY

Short-lived pause in Central California subsidence after heavy winter precipitation of 2017

See allHide authors and affiliations

Science Advances  29 Aug 2018:
Vol. 4, no. 8, eaar8144
DOI: 10.1126/sciadv.aar8144

Abstract

The Tulare Basin in Central California is a site of intensive agricultural activity and extraction of groundwater, with pronounced ground subsidence and degradation of water resources over the past century. Spatially extensive observations of ground displacements from satellite-based remote sensing allow us to infer the response of the aquifer system to changes in usage and to marked recharge events such as the heavy winter rainfall in 2017. Radar imagery from the Sentinel-1a/b satellites (November 2014 to October 2017) illuminates secular and seasonal trends modulated by changes in withdrawal rates and the magnitude of winter precipitation. Despite the increased precipitation in early 2017 that led to a marked decrease, or in some areas, reversal, of subsidence rates, subsidence returned to rates observed during the drought within a matter of months.

INTRODUCTION

The management of water resources is critical to global issues of food security, drought, and climate change (1, 2), particularly at a time when most of the large groundwater aquifers around the world are being used far more rapidly than they are naturally replenished (35). Adverse effects of groundwater overdraft have been observed for decades, including contamination and saltwater intrusion, the financial impacts of drilling new and deeper wells, damage to ecosystems (6), and large-scale, rapid, land subsidence and associated damage to infrastructure. Geodetic observations of land subsidence in the San Joaquin Valley reach back to the mid 1920s (7), with magnitudes of subsidence reaching approximately 9 m by the 1970s. Subsidence rates increased further, beginning in the winter of 2012, during California’s most extreme drought in recorded history. Between October 2012 and September 2016, approximately 40 km3 of water was lost from California’s Central Valley, as estimated using NASA’s Gravity Recovery and Climate Experiment (GRACE) (5, 8, 9). Global Positioning System (GPS) observations within the same region document regional crustal uplift in response to this mass deficit (1012). Space-based and airborne interferometric synthetic aperture radar (InSAR) allow characterization of land subsidence at scales that range from a single aquifer to a watershed, with high spatial resolution (meters to tens of meters) and temporal resolution of days to months (1316). Here, we examine C-band imagery from the Sentinel-1a/b satellites, part of the European Space Agency’s (ESA) Copernicus program, and describe the response of the Tulare Basin to the marked winter precipitation that occurred in 2016–2017.

California’s Central Valley is one of the world’s most productive agricultural regions and is highly dependent on groundwater resources. The Tulare Basin, which comprises the southern two-thirds of California’s Central Valley, is bounded by the Sierra Nevada mountain range to the east and Coast Ranges to the west (Fig. 1). The basin began to form during the Early Paleocene and is filled with sediments of marine, deltaic, and continental origin (17, 18). Roughly half of the sediments are fine-grained clays and silts susceptible to long-term compaction (19). The most areally extensive fine-grained deposit is the Corcoran clay layer, which has a thickness of up to 50 m (20) and acts as the confining layer for the confined aquifer utilized by the majority of wells in the western part of the basin (21). There is very little rainfall within the Tulare Basin, so natural groundwater recharge is dominated by percolation into the semiconfined and confined aquifers through coarse-grained alluvial deposits in the foothills of the Sierra Nevada via streams including the Tule River, Poso Creek, Deer Creek, and White River (22). This watershed once fed into the Tulare Lake, which was the largest freshwater lake in the western United States before tributaries were diverted at the end of the 19th century.

Fig. 1 Subsidence rates and seasonal amplitudes.

(A) Secular rates of displacement in the satellite’s line of sight (LOS) inferred from 40 Sentinel-1 acquisitions (November 2014 to November 2016), as well as the Corcoran clay extent (gray, white outline), rivers (black lines, dotted where inferred), major canals (black dashed lines), and faults (49) (red lines). Inset indicates location of (A) and (B) (red square), Tulare Basin (blue), and Corcoran clay extent (gray). (B) Amplitude of best-fit sinusoidal term. SAF, San Andreas Fault.

The immediate response of Earth’s shallow crust to the extraction of groundwater results in poroelastic compression of coarse-grained sands and gravel deposits in the aquifer, which results in a decrease in pore water pressure. This decrease in pore water pressure increases the effective lithostatic stress in the shallow crust and causes subsidence of the overlying surface. In many cases, the resulting deformation can be recovered with sufficient recharge to the aquifer (23, 24). Periodic subsidence and uplift cycles associated with the annual water cycle, which occur with annual periodicity, are commonly observed geodetically in many areas around the globe (25). However, in cases where extraction rates lower the hydraulic head below the preconsolidation level, there can be unrecoverable (inelastic) deformation due to compaction of the fine-grained sediments constituting the aquitard layers. This compaction can occur over the course of tens to hundreds of years—long after cessation of pumping as equilibrium is reestablished (2628)—and can cause prolonged subsidence on the order of tens of centimeters per year for large groundwater systems. These large magnitudes of subsidence can contribute to damage to infrastructure, flooding, faulting, decreased aquifer storage capacity, and degradation of wetland environments (15, 2931).

Growing demand from the Central Valley’s agriculture industry has increased groundwater dependence and pumping rates over the last century, leading to average rates of depletion with magnitudes as high as 1.85 km3/year since 1962 and more in recent years (table S1) (32). Water diversion systems, such as the Friant-Kern Canal, import surface water from the north and the valley margins (32), and can be impacted by subsidence and affect the distribution of water that would otherwise contribute to aquifer recharge. Temporary regulations were set during the drought to require usage reporting and reduced water use in urban areas. However, no restrictions or reporting requirements exist for agricultural usage, which constitutes approximately 80% of water consumption in California. In 2014, state laws were emplaced requiring agencies to develop plans for sustainable groundwater use by the year 2042 (Sustainable Groundwater Management Act).

MATERIALS AND METHODS

The end of the drought in California was declared in April 2017, after higher than normal rainfall during the winter months. To study the spatial and temporal evolution of ground deformation in response to the drought and its recovery in the Central Valley, we analyzed data from 69 SAR images from the C-band (5.6-cm wavelength) Sentinel-1a/b satellites (path 144; Fig. 2)—part of ESA’s Copernicus Program (33, 34). With the explosion of freely available SAR data acquired by Sentinel-1a/b, which image California every 6 to 12 days, we can now observe complex displacement histories in many areas previously limited by low coherence due to insufficient temporal coverage of SAR data relative to the time scale of surface change. We processed SAR imagery with the InSAR Scientific Computing Environment (ISCE) (35) through the stages of interferogram formation and removal of topographic effects at full resolution (20 m in azimuth and 5 m in range). We generated all sets of sequential interferograms between adjacent dates (68 interferograms) and used the Shuttle Radar Topography Mission (SRTM) digital elevation model to remove topographic effects (36).

Fig. 2 Temporal variability of ground deformation.

Displacement history at location of GPS station CRCN (Fig. 1) for Sentinel-1 (red dots) and GPS (black dots) projected into Sentinel-1 LOS (~38° from vertical), and best-fit model of secular + seasonal signal (green line), with daily stream discharge (light blue curve) at the U.S. Geological Survey (USGS) monitoring site at Deer Creek (https://waterdata.usgs.gov). The light gray–shaded area indicates the approximate time range of drought in California, with purple shading highlighting the rate changes that began in early 2017. We use GPS position solutions in the NA12 reference frame, available for download on the data products pages of the Nevada Geodetic Laboratory (http://geodesy.unr.edu).

We faced three primary challenges when constructing and interpreting reliable InSAR time series in the Tulare Basin: (i) constructing high-quality interferograms in the presence of rapid changes in ground surface properties due to agricultural activity, (ii) referencing displacement to a common datum in the presence of large magnitude signals that extend over a large portion of the SAR frame—a problem that can bias rates by up to several centimeters per year, and (iii) deconstructing the complicated displacement histories into a useful set of basis functions.

We addressed the first challenge as follows: Agricultural areas, such as the Central Valley, are often associated with low interferometric coherence because of the constantly changing vegetation and ground surface properties (37). Approaches that identify stable pixels and regions within a study area (38, 39) have been shown to improve the quality of the InSAR time series, as have methods that use the quality of individual, full-resolution pixels as weights during spatial averaging (40). We followed a variation of the latter approach, spatially averaging (six times in range and three times in azimuth) each interferogram and giving higher weight to the most stable pixels (fig. S1) (39). The most stable pixels are commonly located along roads and at locations of buildings and other more permanent structures (fig. S1). This weighted downsampling approach reduces the impact of decorrelation in regions where some pixels are relatively stable over time (for example, roads surrounded by agricultural fields).

However, while the weighted downsampling improves coherence along individual roads and corridors, the algorithm that we used for phase unwrapping (41) is a minimum cost flow approach that does not always perform well at unwrapping linear regions of good coherence surrounded by areas of poor coherence. Therefore, we also applied an additional masking and interpolation step to enhance connectivity between “good” regions during phase unwrapping (42). Our masking threshold involves removal of downsampled pixels where fewer than 6% of the original, full-resolution pixels that contributed to the downsampling were above a threshold of 0.4 for phase stability (fig. S1). We interpolated values within the masked regions using a Gaussian-weighted average of the remaining pixels within a five-pixel radius. We then performed phase unwrapping with SNAPHU (Statistical-cost, Network-flow Algorithm for Phase Unwrapping) (41). The areas that were interpolated were masked out after phase unwrapping, so they do not contribute to our final analysis except for their impact of resolving the phase unwrapping ambiguity for nearby points. Finally, we inverted our suite of interferograms for the temporal history of displacement at each pixel using the small baseline subset (SBAS) approach (43). Because we used all sequential pairs of dates as our set of interferograms, the inversion for the temporal history of displacement is trivial—a more complex set of interferograms could be used with the SBAS approach (43) or similar.

The second challenge is related to the definition of zero displacement for each time period (that is, each date relative to the first date). Interferograms are necessarily differences between the state of the surface and the propagation delays within the troposphere/ionosphere at two separate times, so there is, at minimum, an ambiguity in the overall shift of each interferogram relative to the rest that are being used in a time series analysis. Errors from propagation delays within the atmosphere and errors in the inferred position of the satellites can also introduce signals with large spatial scales. When these effects are neglected, the time series at all pixels will have an increased degree of scatter due to the spatially correlated offsets of the individual interferograms. Two primary approaches have been used in the past to resolve this ambiguity—if the deforming region is small relative to the size of the study area, then it is masked out, and a simple function (an offset, plane, or quadratic surface) is fit and removed from the remainder of each interferogram (44, 45). Alternately, if the deformation field is smooth and is sampled sufficiently by independent geodetic observations [for example, from GPS, then a prediction can be removed from each interferogram, followed by the same sort of removal of a simple function. In California’s Central Valley, the deformation signal is sufficiently large and spatially complex that it would be aliased by the sparse distribution of existing GPS sites.

We identified “nondeforming” regions by examining the frequency content of the SBAS time series at each pixel and then removed the mean value averaged over just those nondeforming areas. Signals with high power at low temporal frequencies (for example, periods longer than 180 days) are related to processes that occur smoothly over longer time periods such as subsidence and seasonal variations in water storage. Conversely, signals with low power in low-frequency bands are typically dominated by other processes related to noise such as atmospheric delays, which are approximately random in time. To identify areas with real deformation, we initially applied a low-pass filter to the time series at each pixel with no removal of a mean and then took the SD of each time series. Pixels that have SDs at low frequencies greater than 10 cm over the entire time range are considered pixels with a real deformation signal. We excluded that subset of pixels to calculate the mean and remove it from the SBAS time series results at each date (Fig. 3). This approach allows us to automatically reference the displacement in individual dates with respect to each other without introducing bias from the longer time scale signals of interest, and without making prior assumptions about where those signals are located. The success of this approach can be seen through the lower scatter that results in the nondeforming regions (Fig. 3, A and B).

Fig. 3 Displacement profiles and time series reference area.

(A and B) Profiles of cumulative displacements colored by acquisition time [profile location in (C)]. Every fifth date is plotted. (A) Cumulative displacements after removing a planar function and mean offset from each date. Note that the large region of subsidence has biased this estimate, resulting in an artificial signal of apparent uplift along the margins, particularly near point P. (B) Cumulative displacements after removing mean offset calculated on the basis of subset of pixels defined by SD criterion. (C) SD of time series at each pixel after applying a low-pass filter (allowing only periods greater than 180 days). The black contour line (10 cm) represents the threshold SD used to separate deforming areas from nondeforming areas in offset correction. Nondeforming areas were used as reference pixels in the displacement time series. Note that even some of the smaller subsidence features associated with hydrocarbon extraction in the southwest are automatically identified using our approach.

To address the third challenge (characterizing the displacement history), we solved for the secular rate and sinusoidal term (amplitude and phase) that best characterize the time series at each pixel during the time period of peak drought (before November 2016). We used a nonlinear least-squares solver based on the interior-reflective Newton method (39) to find the model parameters that, in a least squares sense, best fit the data points spanning November 2014 to November 2016 (Fig. 1). We also defined a value, α, that we refer to as “deviation from the multiyear trend,” which is the difference between any measured displacement and the model fit to only the time period of peak drought (Fig. 2). This calculation allows us to quantitatively assess how the aquifer system responds to the increased rainfall at the end of the drought as inferred from the complex displacement history. We formulated this problem asEmbedded Image(1)where αij is the deviation from the multiyear trend at the ith pixel at time tj, with time measured relative to the first date (8 November 2014), and uij is the displacement at the ith pixel at time tj, measured in the direction of the satellite’s LOS. Vi, Ai, and φi are the secular rate, amplitude, and phase shift of the seasonal term inferred from only data acquired during the period of drought.

RESULTS

The displacement time series over the Tulare Basin has a spatial resolution of ~60 m, temporal resolution of 6 to 12 days, and spans the time period November 2014 to September 2017. The time series is dominated by subsidence related to groundwater extraction in the northern half of our area of interest with much smaller, isolated areas of subsidence and uplift in the south near Wheeler Ridge/Lost Hills related to hydrocarbon production (Fig. 1) (46). Maximum rates reach over 45 cm/year within the satellite LOS (Fig. 1A). If we assume that all displacement is in the vertical direction (that is, subsidence), then this LOS value would represent ~55 cm/year of subsidence. We present our results as displacements projected along the near-vertical satellite LOS to avoid any potential bias associated with the assumption that all displacements are vertical. Records from continuous GPS stations within and on the margins of the Central Valley demonstrate that seasonal displacements in the east-west or north-south directions are on the order of a centimeter or more (fig. S2). Subsidence accelerates in the drier summer and fall months when less surface water is available for irrigation and farmers rely more heavily on groundwater extraction for agricultural production. Conversely, subsidence rates are lower in wetter winter and spring months due to more recharge to the aquifers and a lower need for pumping. We observed this seasonal pattern in the displacement time series at individual pixels (Fig. 2). The spatial distribution of the highest secular rates (V; Fig. 1A) overlaps broadly with the region that exhibits the largest seasonal amplitudes (Fig. 1B); however, the regions where these two metrics have their maximum values differ significantly. The peak subsidence rates occur further northwest compared to the peak seasonal amplitudes to the southeast. This pattern may provide insights into the properties of the aquifer and general patterns of pumping behavior, which will be discussed in more detail.

To better understand the aquifer system’s response to the changing recharge and extraction associated with the end of the period of drought, we explore the spatiotemporal characteristics of the displacement history in two ways: first by examining deviation from the multiyear trend (α) described above (Figs. 2 and 3 and fig. S5), and second by directly intercomparing the displacement histories of the three individual years of observations versus day of year. The magnitude of α evolves over time in two primary stages: the first stage in the spring season of 2016 and the second stage in the spring season of 2017 (Fig. 4). The initial stage is concentrated in the southeastern section of the aquifer system. It follows an approximately linear trend in space that coincides partly with the path of the Friant-Kern Canal and partly with the Deer Creek and White River flowing from Sierra Nevada tributaries to the east, where daily discharge is monitored by the USGS (Figs. 2 and 4). The second stage of increased deviation from the multiyear trend (α), which began in early 2017, continued in the southeast but expanded to the northwest section of the aquifer system, where the highest secular rates of subsidence during the drought were observed. It is correlated in time with the anomalously heavy winter precipitation. We observe time lags of 1 to 5 months between the onset of pronounced precipitation as measured at Deer Creek and the onset of pronounced increase in the deviation from the multiyear trend (α), which is roughly consistent with time lags observed between drought index changes derived from GRACE versus those derived from precipitation during the onset of the drought (9).

Fig. 4 Spatiotemporal progression of deviation from the multiyear trend (α).

Filled contours colored by date when LOS displacements deviated by more than 13 cm from model fit during peak drought. Time ranges are shown in inset plot (October 2016 to July 2017). Inset plot: α versus time at points 1 to 5. Light blue curve is a time series of daily stream discharge at the U.S. Geological Survey (USGS) monitoring site at Deer Creek (https://waterdata.usgs.gov) filtered with a 2-week moving average.

This analysis assumes that the behavior before 2017 closely follows a form similar to that in Eq. 1, without any variation in time of the quantities Vi, Ai, and φi. For an alternative approach, we directly compare the displacements each year from 2015 to 2017 versus day of year (Fig. 5). This approach eliminates potential bias and inherent assumptions associated with fitting a model (Eq. 1) to a subset of data as we do in our analysis of α. Both results are consistent and indicate marked deceleration of subsidence in 2017. We also see a two-stage recovery in the southeast area of the aquifer system that aligns with the Deer Creek and White River drainage (point 1; Figs. 4 and 5), and a single later stage of recovery in the remainder of the aquifer system (point 4; Figs. 4 and 5). Both approaches demonstrate that parts of the central area of the aquifer are associated with a magnitude of the deviation from the multiyear trend (α) of 20 cm or more (Figs. 2, 4, and 5).

Fig. 5 Year-to-year differences in subsidence behavior.

(A) Time series at two representative points as a function of the day of year for each of the three years, beginning on the date of the first available SAR acquisition (November 8; day of year, −53). Locations of examples 1 and 4 are the same as in Fig. 3 and are shown in map view in (B) and (C). Cumulative difference of first and second years (B) and second and third years (C) at each pixel, averaged over shaded time period shown in (A). Map symbols are the same as in Fig. 1. Differences less than 2 cm are transparent. This analysis assumes no model (for example Eq. 1) to characterize the data but still illustrates the spatial and temporal variability in subsidence, indicating that the effects that we observe are not just an artifact of the time period used in fitting Eq. 1.

DISCUSSION

The complex spatial and temporal characteristics of ground deformation observed during and after the drought require consideration of processes related to both natural and anthropogenic influences. The seasonal cycles of precipitation and the corresponding elastic deformation of the aquifer/aquitard system occur in conjunction with variations in anthropogenic water use during and after the drought—including the diversion of flowing water through canal systems, the storage of water in reservoirs, changes in agricultural activities (that is, crop rotation, fields left fallow, etc.), and the direct extraction of groundwater. The variations that we observe between the northwest and southeast regions of the aquifer system are likely related to differences in the pattern of water usage and recharge, as well as the properties of the aquifer itself. The northwest section of the aquifer is farther from the valley margins and is located over the confined section of the aquifer. The higher abundance of clay layers and observations of continued subsidence throughout 2016 may suggest that this part of the aquifer is more susceptible to compaction. Conversely, the southeastern section is closer to water sources coming out of the mountains, such as Deer Creek, which are associated with subsurface deposits of high permeability (22), and is located above a semiconfined part of the aquifer. The observed α magnitude and higher seasonal amplitudes in this area are likely related to a combination of higher availability of natural recharge, less stress from pumping, and lower amounts of compaction.

In the winter of 2016 precipitation returned to a normal level; however, only the areas shown in Fig. 5B, primarily in the southeastern section of the aquifer system, exhibited significantly less subsidence than in 2015. These areas may be where (i) sufficient water directly recharged the aquifer and/or (ii) there was sufficient renewal of surface water storage to require less groundwater pumping for agricultural and municipal use. Furthermore, recovery in the confined section of the aquifer would show a less distinct subsidence deceleration as the aquifers continue toward equilibrium following any compaction that occurred in previous years. Nevertheless, much heavier precipitation throughout the winter of 2017 increased the availability of surface water for irrigation and lowered groundwater dependence across the remainder of the watershed, resulting in lower subsidence rates relative to the time period of peak drought (Fig. 5C). However, while ~47% of wells in the Tulare Lake region reported groundwater levels more than 1.5 m higher in Spring 2017 than in Spring 2016, 86% were still more than 1.5 m lower than in Spring 2011(47). This suggests that recovery from the drought has not been complete and agrees with our results that show continued subsidence throughout the remainder of 2017 (fig. S6).

The growing demand for groundwater resources in California and across the globe increases at an even higher rate during periods of drought. The sustained subsidence observed in the geodetic time series suggests that unrecoverable inelastic compaction in aquitard layers likely occurred (24, 48), permanently lowering the storage capacity of the aquifer system. California, like many areas, has few restrictions on groundwater use, and reporting on pumping/usage is sparse and mostly unregulated. Efforts to assess and manage water scarcity challenges globally are supported by the existence of independent, free, and open data sets, particularly when there is inadequate reporting of groundwater usage. Space-based geodesy can provide insight into the response of the Central Valley groundwater system to variations in precipitation and usage, enabling evaluation of the efficacy of future policy choices and efforts to mitigate damage to this vital resource.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/8/eaar8144/DC1

Fig. S1. Data quality metrics.

Fig. S2. GPS time series.

Fig. S3. Time series at points.

Fig. S4. Example seasonal signal.

Fig. S5. Deviation from multiyear trend (α) averaged over March 2017.

Fig. S6. Continued subsidence, 2017.

Table S1. Previous observations of subsidence in the Tulare Basin.

Movie S1. Spatiotemporal progression of deviation of subsidence from the multiyear trend.

Data file S1. Text file containing a list of all dates of SAR acquisitions used in this study.

References (5053)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank T. Farr and an anonymous reviewer for their thorough and helpful comments on the manuscript. Funding: K.D.M. and R.B.L. were partially supported by NASA grants NNX16AL20G and NNX16AK57G. Author contributions: K.D.M. conceptualized the study, performed the research, and wrote the manuscript under the help and guidance of R.B.L. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Sentinel-1a/b data used in this study have been made freely available by ESA and were downloaded through an agreement between ESA, NASA, and the Alaska Satellite Facility. Stream flow observations are available from the USGS at https://waterdata.usgs.gov. GPS position solutions used here are available for download on the data products pages of the Nevada Geodetic Laboratory (http://geodesy.unr.edu). Most data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from authors.
View Abstract

Navigate This Article