Contemporaneous Subsidence and Levee Overtopping Potential, Sacramento–San Joaquin Delta, California

provides open access, scholarly publishing services to the University of California and delivers a dynamic research platform to scholars worldwide. for this work the Commission's the Climate Scripps Additionally we Greg Smith for providing ground-water well records; Mark Merrifield and Jon Stock for discussion; and Cecily Wolfe, Steve Ingebritsen, and Devin Galloway for critical pre-submission reviews that significantly improved the paper. Abstract: The levee system in California’s Sacramento-San Joaquin Delta helps protect freshwater quality in a critical estuarine ecosystem that hosts substantial agricultural infrastructure and a large human eScholarship provides open access, scholarly publishing services to the University of California and delivers a dynamic research platform to scholars worldwide. ABSTRACT The levee system in California’s Sacramento–San Joaquin Delta helps protect freshwater quality in a critical estuarine ecosystem that hosts substantial agricultural infrastructure and a large human population. We use space-based synthetic aperture radar interferometry (InSAR) to provide synoptic vertical land motion measurements of the Delta and levee system from 1995 to 2000. We find that Delta ground motion reflects seasonal hydrologic signals superimposed on average subsidence trends of 3 to 20 mm yr -1 . Because the measurements are insensitive to subsidence associated with peat thickness variations over Delta-island length scales, it is most likely that InSAR rates reflect underlying Quaternary sedimentary column compaction. We combine InSAR rates with sea-level rise scenarios to quantify 21st century levee overtopping potential. If left unmitigat-ed, it is likely that 50 to 100 years from now much of the levee system will subside below design thresholds.


INTRODUCTION
The largest estuary in the Western United States, the Sacramento-San Joaquin Delta provides fresh water to ~1 million cultivated hectares and more than twothirds of California's human population (Figure 1). The Delta, also an important migratory pathway for birds and endangered fish species (Lund and others 2007), has been a tidal freshwater marsh for over 6,000 years since decreased submergence rates outpaced sea-level rise following the last glaciation (Atwater and others 1979), and it is underlain by the Great Valley sedimentary column that is as much as 9 km thick (Dickinson and Seely 1979;Atwater 1982). Delta islands, particularly the four westernmost, help maintain freshwater quality by acting as a barrier to seasonally-and tidally-driven saline incursions from the San Francisco Bay (Galloway and others 1999;Lund and others 2010).
The Delta islands are protected by more than 1,700 km of earthen levees whose structural stability is variable (over 100 levee failures have occurred since the 1890s) and poorly known (Mount and Twiss 2005;Suddeth and others 2010). If a levee fails, rapid island-infilling can draw brackish water back into the Delta and seriously degrade fresh water quality and supply. While structural stability-the resistance to seismic-shaking and strong differential land subsidence-is most often the principal factor determining levee vulnerability, susceptibility to high water overtopping is also considered a threat (CDWR 2009). For instance, during Hurricane Katrina in 2005, overtopping exacerbated by regional subsidence was one of the principal failure modes, and many overtopped New Orleans levees spatially coincided with high modern subsidence rates (Dixon and others 2006).
Land subsidence since the late 1800s, thought to occur primarily by oxidation of peat soils drained for agriculture purposes (Rojstaczer and others 1991;Deverel 1993, 1995;Deverel and Rojstaczer 1996;Deverel and others 1998;Mount and Twiss 2005;Deverel and Leighton 2010), has left Delta island elevations as much as ~ 8 m below sea level (Coons and others 2008) (Figure 1). Local geodetic studies demonstrate that some islands have subsided at rates of ~ 4 to 70 mm yr -1 during the past century, and that island interiors have subsided by greater amounts than their margins Deverel and Rojstaczer 1996;Deverel and others 1998;Deverel and Leighton 2010). Peatrelated subsidence in the Delta is caused by oxidation of soil carbon resulting from the drainage of agricultural lands that started in the late 19th century (Deverel and Rojstaczer 1996). Subsidence in some locations has slowed since the 1950s because of the decreasing organic content of island soils and the implementation of improved land-use techniques (Deverel and Rojstaczer 1996;Deverel and others 1998;Mount and Twiss 2005;Deverel and Leighton, 2010).
Given that maximum measured subsidence rates in the Delta are of the same order of magnitude or up to 70 times the currently observed rate of global sealevel rise (~1 to 3 mm yr -1 ) (Church and White 2006; Woppelmann and others 2009), quantifying Deltawide subsidence is essential to best inform mitigation and planning decisions that attempt to balance ecosystem health and anthropogenic needs. The practical difficulties in making ground-based measurements over large areas, however, precludes their use in Delta-wide investigation of the spatio-temporally varying subsidence. Furthermore, subsidence of the levees themselves has remained largely unquantified.
In this paper, we use a space-based geodetic method, synthetic aperture radar interferometry (InSAR), to measure Delta subsidence monthly from 1995 to 2000. The synoptic nature of the measurements and the myriad returns from the levees allow us to couple the InSAR measurements with sea-level rise projections, and to estimate the overtopping potential of the Delta-wide levee system in the 21st century.

INSAR ANALYSIS
We used space-based InSAR (Burgmann and others 2000) to measure Delta ground motion with sub-centimeter motion resolution from 38 descending pass ERS-2 image pairs that span 1995 to 2000 ( Figure 2, Table 1, Appendix A). Traditional InSAR performance is compromised when interferometric coherence is degraded from either the phase instability of specific portions of the reflective surface or satellite orbital baselines exceeding a critical distance threshold.
Here, to increase temporal and spatial resolution arising from poor interferometric coherence in the highly cultivated Delta, we employed the recently developed 'persistent scatterer' interferometry (PInSAR) technique (Ferretti and others 2001;Hooper and others 2004). PInSAR uses the fact that minimal baselinerelated de-correlation occurs for stable, point-like reflectors, and so interferometric phase may be interpreted even for scene pairs with long perpendicular baselines that may exceed the critical baseline. The method differs from traditional InSAR processing primarily by: (1) the selection of stable scatterers; (2) elevation model correction; and (3) time-series analysis that is focused on isolating deformation signals from other error sources that are dominated by atmospheric-phase delays caused by tropospheric water vapor content. For a general discussion of the technique, we refer the reader to a recent textbook (Kampes 2006). In the appendix, we provide more details specific to our analysis in the Delta.
We used the technique of Werner (Werner and others 2003;Brooks and others 2007) and determined point targets using joint measures of backscatter temporal amplitude variability and spatial spectral diversity of candidate points (Ferretti and others 2001). The method yields 117,235 stable radar-scattering targets ('PS targets'). Of these, ~31,00 are within the Delta, and they are derived primarily from buildings, electrical towers, and road guardrails: ~40% of the levee system yields measurements ( Figure 2).
We present the results as: (a) a map of each scatterer's average annual vertical displacement rate ( Figure 2); (b) representative and median vertical displacement-rate time series for localized regions ( Figure 3); and, (c) to characterize longer spatial Sets"). Moreover, if motion of the levees towards islands interiors resulting from differential subsidence caused significant horizontal motion, then, for a NW-oriented look direction, we would expect preferentially positive LOS range change values (motion towards the satellite) on the western sides of the islands. The lack of any such pattern suggests that horizontal motion from levee instability does not significantly affect the InSAR results.
wavelength ground motion signal, the median annual rate of all scatterers in 1 ×1-km-grid cells throughout the study area ( Figure 4). We attribute all line-ofsight (LOS) motion to vertical displacement, finding no evidence to suggest that InSAR-measured motion is significantly affected by differential horizontal motions associated with the greater San Andreas fault system, by hydrocarbon production, or by related injection (see also "Comparison With Other Data  The grid cell analysis highlights the general downward motion of the Delta with respect to its margins ( Figure 4), a pattern consistent with an independently derived InSAR deformation study of a similar time period over the same area (Burgmann 2009). Although individual scatterers and local regions may subside by up to 2 cm yr -1 (Figure 2), generally, Delta interior regions demonstrate subsidence on the order of 3 to 5 mm yr -1 (Figure 4). The time series also show the relative stability of the margin sites ( Figures 3A, 3B) and the downward trend of sites within the Delta ( Figure 3C, 3D). Indeed, the median time series of all scatterers within the Delta margins clearly shows a downward trend ( Figure 3E). The time-series data also indicate seasonal vertical oscillations (rising in winter months, falling in sum-mer months) upon which the trends are superimposed ( Figure 3). Below, we explore the correlation of these seasonal oscillations with concomitant hydrological data from the region.

Continuous GPS Data
For the 1995 to 2000 time period, contemporaneously collected ground-control data do not exist, precluding direct assessment of the stability of our solution in an absolute reference frame. However, four plate boundary observatory ( to place constraints on our InSAR results (see also Appendix A).
Because only one satellite orbital path and look direction are available for this study, the InSAR solution yields only one displacement component, and so differential horizontal motions from local tectonics (d'Alessio and others 2005)-particularly east-west directed (sub-parallel to the radar's look direction)could be incorporated into the range-change solution.
Because they are spatially distributed throughout the Delta, we use the PBO CGPS average velocities to assess differential horizontal motion. For each pair of  Figure 4 Colored boxes, median InSAR-derived vertical displacement rate from all scatterers in 1x1-km-grid cells. Colored bold circles, vertical displacement rate from Plate Boundary Observatory (PBO) continuous GPS stations, station names indicated. Green circles, re-located earthquake epicenters from Waldhauser and Schaff (2008). Red, white, blue circles, groundwater well-level annual trend (1995)(1996)(1997)(1998)(1999)(2000) from DWR groundwater records. InSAR and water trend data share the same color scheme though with different max/min values.
PBO sites, we calculated average horizontal velocity difference, and projected this onto LOS. The maximum range-change value has a magnitude of less than 0.5 mm yr -1 and so we consider the effect of differential horizontal motion to be negligible.
To compare CGPS displacement rate estimates with InSAR, we projected the 3-component CGPS velocity (Table 2) onto LOS, and then projected this onto vertical. Because CGPS site P257 exhibits vertical motion of only ~ -0.2 mm yr -1 we chose this as a reference point for the InSAR data and explored potential reference-frame biases by imposing a range of biases (± 5 mm yr -1 ) and calculating RMS error with the CGPS sites ( Figure 5). The results demonstrate that the choice of P257 as a reference (0 bias case) is good to better than a 2 mm yr -1 RMS ( Figure 5), which is particularly encouraging given that the measurement periods do not overlap. If we allow an upward bias of 1 mm yr -1 , the RMS error decreases slightly ( Figure 5), although for internal consistency in this analysis, we prefer to leave our InSAR reference point coincident with the P257 CGPS site. The results in Figure 5 demonstrate that we should not interpret rate differences at less than the 1 mm yr -1 level, and below we further choose to impose a ± 5 mm yr -1 threshold for projection purposes.

Hydrocarbon-related Data
It is well known that anthropogenic activities, including hydrocarbon extraction and injection as well as groundwater pumping can produce substantial ground-motion signals (Bawden and others 2001;Borchers 1998). Because there are a number of gas and oil fields in the Delta (Figure 6), it is critical to assess whether the InSAR ground-motion data are contaminated by the anthropogenic signals. We compare records of hydrocarbon monthly production and injection volumes (http://www.conservation.ca.gov) and displacement time series, both visually and quantitatively by evenly sampling and removing a linear trend from each time series, and computing the cross-correlation coefficient for a range of temporal lags (Figure 7, Appendix A). A causative relationship would be identified by a negative correlation coefficient (the surface subsides when hydrocarbon production increases) and a positive time lag (the forcing must precede the response). The analysis shows that there are no obvious longer-term displacement rate changes associated with either increases or decreases in production and/or injection at any of the sites from 1995 to 2000. Thus, we can find no evidence that the signal is contaminated by hydrocarbon extraction/injection activity, despite its widespread nature across the Delta, a conclusion reached independently by others (Rojstaczer and others 1991). We suggest that this is likely because of the combined facts that overall volumes are small compared to places where large subsidence signals have been seen (Gambolati and Teatini 1998), and that time periods are long (Hettema and others 2002) between produc-  tion peaks and our measurement period. This analysis leaves as unexplained the uplift signal on the western limit of the Delta near the Brentwood gas field and the town of Antioch (Figures 4, 6, A6, A22). Because this signal from the Delta's margin is secondary to the main subsidence signal we leave further exploration of its cause for future work, although we note that it is not likely tectonic-related due to the high rate and lack of observed local tectonic activ-

Hydrological Data
To explore the possible correlation of InSAR-derived ground-motion signals with hydrologic phenomena, we obtained groundwater-level time series from the  Map showing Delta gas fields with production activity (green polygons). Cyan circles represent a 5-km radius about the centroid of each gas field. Magenta circles are InSAR scatterers falling within each circle. The time-series from each of these scatterers is used to construct the median time-series, an example of which is displayed in Figure 7.
California Department of Water Resources database (http://wdl.water.ca.gov/gw/) for 258 wells that lie within our study area (Figure 4, Table A1). Generally, California's Central Valley aquifer system is divided into a number of basins (Williamson and others 1989) and administrative sub-basins (CDWR 2003), and we refer the reader to the previous references for a thorough overview and for specific information on lithology, water-bearing units, and flow trends.
The records are from maximum depths of ~40 m, and have variable sampling intervals ranging from monthly to quarterly. Additionally, we obtained river flow data for the Sacramento and San Joaquin rivers and the streams on the east side of the Delta (CDWR 'DayFlow' algorithm; http://www.iep.ca.gov/dayflow).
We used the InSAR time series data to investigate possible seasonal correlation with groundwater lev-Jan80 Jan85 Jan90 Jan95 Jan00 Jan05  Gas production records (red) and median InSAR vertical displacement time series (blue) for Lathrop field indicated in Figure 6. Gas volumes are reported in industry standard MCF (thousand cubic feet). The lower plot shows results of cross-correlation analysis of the temporally overlapping gas and individual InSAR time-series. In this case, causative relation would be indicated by negative correlation coefficient and positive time lag. Color-coding is number of records with a given peak correlation coefficient and time lag. No causative relationship is apparent. Both production and injection records from the other fields indicated in Figure 6 show results similar to this figure (see Figures A3-A22). els and surface flow of the principal rivers feeding the Delta. We plotted median time series from the entire InSAR and ground-water well data sets and compared these with the total river flow from the sum of gauged inflows from the Sacramento and San Joaquin rivers and the streams on the east side of the Delta (Figure 8). Because ground-water well levels were not recorded at the same time for each site, we created the median time series for the complete set of wells by sampling the data quarterly ( Figure 8B). Additionally, to better illustrate the seasonal ground-water fluctuations, we de-trended the median ground-water time series by estimating and removing a linear rate from the data ( Figure 8C). The median displacement time series comprises seasonal fluctuations with peak amplitudes of 30 to 40 mm ( Figure 8A) superimposed on an average subsidence trend of ~ 3 to 20 mm yr -1 . The seasonal signal is strongly correlated with both the ground-water well levels ( Figure 8B, 8C) and the Delta's total river flow ( Figure 8D); the CGPS time series from 2006-2009 confirms that the correlation is persistent over the decadal scale ( Figure 8A). Near the start of each year median displacement rates reach a maximum, and towards the end of each year they reach a minimum in close correspondence with maxima and minima in the ground-water ( Figure 8B, 8C) and river flow ( Figure 8D) time series. Taken together, these data suggest a simple conceptual poro-elastic model for seasonal Delta ground-motion seen in other similar environments (Schmidt and Burgmann 2003): the river flow that inundates the Delta during peak spring run-off quickly diffuses through shallow organic and mineral deposits in the Delta and on its margins, causing the land surface to rise elastically. During low-flow months, the excess groundwater moves either laterally or to lower levels in the Great Valley aquifer system (Williamson and others 1989), and the ground surface relaxes.
Analysis of the longer-term signal, however, demonstrates that annual groundwater-level change rates, sampled from 1995-2000, are predominantly positive, in the opposite sense to the InSAR-derived vertical displacement rates (Figure 4). Thus, long-term Delta subsidence is not controlled by hydrological phenomena, and we seek an explanation related to geological process.

Delta-wide Subsidence
The elevated average subsidence rates occur within a broad area defined by the Delta islands (Figure 4). The lack of sharp velocity gradients in our result, combined with the generally low rates of active tectonics in the region (d'Alessio and others 2005), leads us to attribute the first-order signal to a Deltaspecific sedimentary process rather than to a regional seismo-tectonic process (Burgmann and others 2006;Waldhauser and Schaff 2008). The InSAR-derived rates from 1995-2000 in the Delta are significantly less than the ~ 40 to 70 mm yr -1 measured in island interiors (Deverel and Rojstaczer 1996;Deverel and others 1998;Rojstaczer and Deverel 1995) and estimated by assuming that the islands were at sea level at the start of the 19th century (Mount and Twiss 2005) (dashed line, see Figure 9A). Moreover, there is no general correlation between the InSAR-measured rates and Delta elevations for the individual scattering targets ( Figure 9A). This is not consistent with the inverse relationship (higher subsidence rates associated with lower elevations) that would be expected if all InSAR-measured subsidence were caused by the same organic soil-oxidation process instigated by water pumping for agricultural purposes towards the end of the 19th century (Rojstaczer and others 1991). Additionally, the lack of correlation between subsidence rate and elevation holds, even when the scattering target population is subdivided into groups from    within 50 m of levee crests and from within island interiors ( Figure 9A). Therefore, the InSAR rates are not biased by the levee rates, and the InSARmeasured subsidence process that affect both levees and island interiors occurs over length scales greater than island dimensions. In contrast, previous studies provide strong evidence that significant lateral subsidence gradients are associated with peat thickening toward island interiors .
The InSAR-measured rates also are insensitive to subsidence controlled by peat layer thickness variations; where collocated observations exist, there is no correlation between peat thickness and InSAR-measured subsidence rates ( Figure 9B). Thus, instead of shallow peat compaction, we suggest the InSAR-derived rates more likely record the ongoing, slower compaction of the entire underlying sedimentary column deposited throughout the later Quaternary Delta environment (Shlemon and Begg 1975). As a comparison, recent sub-surface measurements (Dokka 2006) and numerical modeling estimates (Meckel and others 2007) for Mississippi Delta sediments with similar mechanical properties, column thickness, and age of deposition as those in the Delta, show that non-peat-controlled compaction rates would be the same order of magnitude as our InSAR-measured rates. Because the radar returns are preferentially from anthropogenic structures that rest on the levees as well as the island interiors, we suggest that InSAR's insensitivity to peat thickness is likely related to conditions associated with the construction of these structures. These conditions are capping by asphalt or packed gravel, nearly complete consolidation of the upper peat layer from the structure's weight, or the structure being anchored below the peat layer. This does not discount visual observations-common within the Delta engineering community-of coupled peat mobilization and levee deformation; rather, it is consistent with levee-related deformation being highly localized. Future work and ground truthing is needed to detect and understand the signals of levee-related deformation in the InSAR data.

Sea-level Rise and Levee Overtopping Potential
Because, as described above, the data most likely measure regional subsidence associated with compaction of the Quaternary sedimentary column that affects island interiors and margins alike, the InSAR results should not be used in analyses of levee stability based on differential Delta-island subsidence (e.g. Mount and Twiss 2005). Rather, it is most appropriate to use the InSAR results to evaluate levee-related subsidence directly. The spatial extent of our measurements allows us to project future levee subsidence and potential overtopping throughout much of the Delta, given sea-level rise scenarios, and assuming that future subsidence rates follow the InSAR-measured ones; we recognize, for instance, that in the Santa Clara Valley, the basin has been uplifting in the 1990s (Schmidt and Burgmann 2003). Levee design standards in the Delta are variable (for a review, see Suddeth and others 2010). For the purposes of a self-consistent analysis, however, we assign to all locations the PL 84-99 standard wherein levees are designed to have ~0.5 m (1.5 ft) freeboard above the 100-year flood stage (Betchart 2008). PL 84-99 is the federal target standard for two-thirds of Delta levees, although few meet this criteria: more conform to the Hazard Mitigation Plan (HMP) criteria of 1.0 ft of freeboard above the 100-year flood stage (Suddeth and others 2010).
We assume that sea level rise in the Delta will not significantly differ from global estimates, and follow Rahmstorf's methodology (Rahmstorf 2007) to project 21st century sea-level rise for low, medium, and high rates ( Figure 10A). These curves, particularly in the latter half of the 21st century, predict significantly higher sea-level rise than the 2 mm yr -1 rate assumed by most previous analyses in the Delta (e.g. Mount and Twiss 2005). Further, for each sealevel rise scenario we present estimates considering 0 and 5 mm yr -1 absolute geodetic reference-frame biases, which are conservative estimates given the good agreement between the InSAR and CGPS data ( Figures 10B-10J). By 2025, for all sea-level rise scenarios and reference-frame biases, the projections indicate that only small, isolated regions in the Delta will have subsided more than 0.5 m below current levels ( Figure 10B, 10E, 10H). By 2050, however, we project widespread subsidence below the 0.5-m threshold level for all but the most extreme downward reference frame bias ( Figure 10C, 10F, 10I). For the 0 mm yr -1 reference-frame bias case, at least 10%, 31%, and 38% of the Delta's entire levee system is projected to be below the 0.5 m level for low, medium, and high sea level rise rates, respectively. These percentages are minima, because of the lack of InSAR returns from some portions of the Delta. Although we cannot use topography to predict levee subsidence rate for these non-reflective regions, because of the lack of correlation between elevation and subsidence rate (Figure 9), we know of no reason to suspect that the Delta-wide process that controls the InSAR-measured subsidence will affect the un-measured portions of the levees differently from the measured ones. Thus, at any given time, the percentage of levees falling below design threshold from which there are InSAR returns (10%, 31%, and 38% in Figure 10F, for example) can be used to estimate the percentage of all levees projected to fall below the design threshold (26%, 78%, and 97%). By 2100, all scenarios-except the lowest sea-level rise rate combined with the lowest reference-frame biasproject that at least ~38% and likely closer to ~97% of all levees will subside below the 0.5-m threshold. These projections are conservative in that they assume historical 100-year flood levels hold into the future. In fact, peak floods are very likely to increase because of a shift from snowfall to rainfall over the Delta's watershed in response to a warming climate (Dettinger and others 2008).
Previous concern in the Delta has understandably focused on the susceptibility of levees to catastrophic events, such as structural failure from high peatrelated subsidence gradients or seismic shaking. Our work shows that overtopping resulting from continued levee subsidence from Deltaic sediment compaction and future sea-level rise also poses a substantial and increasing threat to levee integrity and, ultimately, the Delta's fresh water supply. Our analysis has focused only on levee overtopping potential, though we recognize that the effect of increasing sea-level rise amplified by regional subsidence will imply further negative effects on overall Delta levee integrity that remain to be studied. For instance, the predicted increased frequency and magnitude of 21st century extreme tidal events (Cayan and others 2008) will be superimposed on continuously increasing overtopping potential. Additionally, subsidence not only increases the likelihood of overtopping, but also increases the potential for levees to fail because of underseepage and piping from increases in static hydraulic head (Mount and Twiss 2005). Because the increase in overtopping potential is continuous and can be monitored synoptically via InSAR, however, scientists, engineers and policy-makers can use these results to help devise a solution that will most effectively consider competing ecosystem interests and economic constraints.