Impounded Marshes on Subsided Islands: Simulated Vertical Accretion, Processes, and Effects, Sacramento–San Joaquin Delta, CA USA

There is substantial interest in stopping and reversing the effects of subsidence in the Sacramento–San Joaquin Delta (Delta) where organic soils predomi-nate. Also, the passage of California Assembly Bill 32 in 2006 created the potential to trade credits for carbon sequestered in wetlands on subsided Delta islands. The primary purpose of the work described here was to estimate future vertical accretion and understand processes that affect vertical accretion and carbon sequestration in impounded marshes on subsided Delta islands. Using a cohort-accounting model, we simulated vertical accretion from 4,700 calibrated years before present (BP) at a wetland area located within Franks Tract State Recreation Area (lat 38.059, long −121.611, hereafter, “Franks Wetland”)—a small, relatively undisturbed marsh island—and at the Twitchell Island subsidence-reversal demonstration project since 1997. We used physical and chemical data collected during the study as well as literature values for model inputs. Model results compared favorably with measured rates of vertical accretion, mass of carbon sequestered, bulk density and organic matter content. From 4,700 to model-estimated 350 years BP, the simulated rate of vertical accretion at Franks Wetland averaged about 0.12 cm yr -1 , which is within the range of rates in tidal wetlands worldwide. Our model results indicate that large sediment inputs during the last 150 to 200 years resulted in a higher accretion rate of 0.3 cm yr -1 . On Twitchell Island, greater organic inputs resulted in average vertical accretion rates as high as 9.2 cm yr -1 . Future simulations indicate that the managed impounded marsh will accrete highly organic material at rates of about 3 cm yr -1 . Model results coupled with GIS analysis indicate that large areas of the periphery of the Delta, if impounded and converted to freshwater marsh, could be restored

From 4,700 to model-estimated 350 years BP, the simulated rate of vertical accretion at Franks Wetland averaged about 0.12 cm yr -1 , which is within the range of rates in tidal wetlands worldwide. Our model results indicate that large sediment inputs during the last 150 to 200 years resulted in a higher accretion rate of 0.3 cm yr -1 . On Twitchell Island, greater organic inputs resulted in average vertical accretion rates as high as 9.2 cm yr -1 . Future simulations indicate that the managed impounded marsh will accrete highly organic material at rates of about 3 cm yr -1 . Model results coupled with GIS analysis indicate that large areas of the periphery of the Delta, if impounded and converted to freshwater marsh, could be restored to tidal elevations within 50 to 100 years. Most of the central Delta would require 50 to 250 years to be restored to projected mean sea level. A large portion of the western Delta could be restored to mean sea level within 50 to 150 years (large areas on Sherman, Jersey, and Bethel islands, and small areas on Bradford, Twitchell, and Brannan islands, and Webb Tract). We estimated that long-term carbon sequestration rates for impounded marshes such as the Twitchell Island demonstration ponds will range from 12 to 15 metric tons carbon ha -1 yr -1 . Creation of impounded marshes on Delta islands can substantially benefit levee stability as demonstrated by cumulative force and hydraulic gradient calculations.

INTRODUCTION
The 300,000-ha Sacramento-San Joaquin Delta is a critical natural resource, important agricultural region and the hub for California's water supply.
Within an area of about 81,000 ha in the central Delta, an estimated 4.5 billion m 3 of soils with a high organic matter content (peat) accreted during the last 6,700 years as sea level rose (Shlemon and Begg 1975;Drexler et al. 2009aDrexler et al. , 2009bMount and Twiss 2005). In tidal marshes, peat vertically accretes concomitantly with sea level rise (Jelgersma et al. 1993;Mitsch and Gosselink 2000) and when organic matter accumulates at a faster rate than it decomposes (Boelter and Verry 1977;Mitsch and Gosselink 2000). Climate, plant community dynamics, tectonics, and hydrologic processes influence long-term peat accumulation (Siegel 1983;Jelgersma et al. 1993;Mitsch and Gosselink 2000). Using physical and chemical data for cores collected throughout the Delta, Drexler et al. (2009aDrexler et al. ( , 2009b provided insight into vertical accretion of peat and demonstrated fundamental differences among islands where mineral input was more dominant (e.g., Browns Island) because of sediment input from the Sacramento River and those where organic accretion predominated (e.g., at a wetland area located within Franks Tract State Recreation Area [lat 38.059, long −121.611], hereafter, "Franks Wetland") near the San Joaquin River (Figure 1).
Subsidence of peat began during the late1800s and early 1900s when Delta islands were levied and drained for agriculture (Thompson 1958). The primary cause is microbial oxidation of peat (Deverel and Rojstaczer 1996) and present-day rates range from less than 5 to over 30 mm yr -1 (Deverel and Leighton 2010). Consequences include increased seepage and hydraulic forces on levees, decreased levee stability, increased subsurface drainage loads of dissolved organic carbon, methyl mercury and salinity, and higher drainage costs due to increasing pumping lifts. Also, farming is becoming more difficult and some subsided islands contain areas that are no longer arable because of excessive wetness.
Since the late 1980s, there has been substantial interest in stopping and reversing the effects of subsidence by creating impounded marshes. Under the hypothesis that construction of permanently flooded impounded marshes would stop subsidence and carbons loss, experiments were conducted in 1,000-m 2 enclosures on Twitchell Island beginning in 1993 under Steven Deverel's direction at the U.S. Geological Survey (USGS), in cooperation with the California Department of Water Resources. Deverel et al. (1998) reported a net carbon gain in permanently flooded impounded marshes in the enclosures and, thus, demonstrated their ability to stop and reverse the effects of subsidence. In contrast, they reported a net carbon loss for seasonally flooded wetlands drained from early fall through winter and drained during spring and summer thus indicating that this wetland water management regime will result in continued subsidence. These results and those of Miller et al. (2000) led to the conversion of 6 ha of agricultural land to the impounded marsh demonstration project on Twitchell Island  in 1997 by personnel of the California Department of Water Resources, HydroFocus, Inc., Reclamation District 1601, and U.S. Geological Survey California Water Science Center (USGS-CWSC). The 6-ha plot was split into two similarly sized areas, the East Pond and the West Pond, and the ponds were flooded to depths of 25 cm and 55 cm, respectively in the fall of 1997. Vertical accretion in the Twitchell Island marsh areavaried spatially, and depended on water depth, plant community composition and colonization, degree of marsh maturity, and water residence time (Miller et al. 2008). From 1997 to 2008, wetland-surface elevations increased by an average of about 3 cm yr -1 , but rates varied spatially from -0.5 to + 9.2 cm yr -1 Miller and Fujii, 2010). The largest rates occurred within dense stands of Schoenoplectus acutus (hardstem bulrush) and Typha (cattail) species.
Since the passage of California Assembly Bill 32 in 2006 with its requirements to reduce greenhouse gas emissions, interest has increased in carbon dioxide J e r s e y I s l a n d  (Deverel et al. 1998;Miller et al. 2000Miller et al. , 2008 and there is a large accommodation space below sea level on Delta subsided islands of over 2.3 billion m 3 (Deverel and Leighton 2010). Improved understanding of peat accumulation rates and processes will aid in the development of impounded marshes. Moreover, there is an unfulfilled need to integrate, synthesize, and effectively use data sets for tidal wetlands and the impounded marsh on Twitchell Island to answer questions about future accretion and levee stability benefits. Therefore, the primary objectives of the study described here were to: (1) simulate accretion rates and better understand processes affecting vertical accretion and carbon accumulation in impounded marshes and; (2) quantify the potential benefits to levee stability. The work and data described here were part of the larger study designed to evaluate marsh accretion in the Delta (Drexler et al. 2009a(Drexler et al. , 2009bDrexler 2011).
Vertical accretion rates of admixtures of inorganic and organic material which form peat determine the elevation of a marsh relative to the tidal frame (e.g., Redfield 1972;Stevenson et al. 1986). Based on data in Drexler et al. (2009b) and preliminary analyses, we surmised that data for long-term accretion during the last 4,700 years on Franks Wetland, a remnant channel tidal wetland (Figure 1), would provide insight to simulate future accretion in impounded marshes on subsided islands such as Twitchell Island. This is because accretion is dominated by organic deposition evidenced by highly organic deposits (generally over 60% organic matter content). Recently accumulated sediments in the impounded marsh on the Twitchell Island demonstration project ) also are over 75% organic matter. We used data described in Drexler et al. (2009aDrexler et al. ( , 2009b to simulate accretion on Franks Wetland during the last 4,700 years. We also analyzed data collected during the first 12 years in the Twitchell Island demonstration project and simulated vertical accretion and carbon accumulation rates. Much of the previous modeling of marsh accretion processes in wetlands in the San Francisco Estuary focused primarily on inorganic processes (Krone 1985) and used a constant rate of organic accretion (Orr et al. 2003;French 1993). In contrast, Callaway et al. (1996) used a cohort accounting approach to more explicitly account for organic accretion processes and successfully simulate approximately 300 years of tidal-marsh accretion in the southeastern United States and England. The model processes discrete annual wetland volumes or cohorts and simulates their vertically downward movement in the marsh. Rybczyk et al. (1998) extended the Callaway et al. (1996) approach to simulate elevation changes as affected by sea level rise in a forested wetland in the southeastern United States. Rybczyk and Cahoon (2002) and Day et al. (1999) used a similar approach to simulate effects of sea level rise on marsh accretion in Louisiana and the Venice Lagoon. Most recently, Swanson et al. (2013) used a modified version of the Callaway et al. (1996) model to evaluate changes in marsh surface elevation and the impact of these elevation changes on marsh habitat for specific species of concern. We modified the Callaway et al. (1996) model to integrate data and simulate Delta marsh vertical accretion.

METHODS, DATA SOURCES AND EXPLANATION OF MODEL INPUTS
The Callaway et al. (1996) model simulates surface organic matter and mineral deposition, subsurface organic matter decomposition, below-ground organic matter production and consolidation, and movement of mineral and organic matter accumulated at the surface to older age classes ( Figure 2). The model also calculates changes in cohort composition, mass, and porosity, and tracks the depth and elevation of the cohort in the accumulating sediment. Model inputs include surface organic matter and mineral inputs to the marsh, rates of organic matter decomposition, below-ground organic matter production and consolidation, initial and limiting porosities, sea level rise and organic and mineral particle densities. We modified the Callaway et al. (1996) model to simulate accretion in Franks Wetland and the Twitchell Island wetland demonstration project ( Figure 1). The overall approach was to develop simple models with input parameter values based on available data. We calibrated as few input variables as possible. We used results described in Drexler et al. (2009aDrexler et al. ( , 2009b and Miller et al. (2008) for model calibration and validation (Tables 1 and 2). Through extensive analysis of peat cores collected throughout the Delta, Drexler (2009aDrexler ( , 2009b addressed questions about the processes that affect rates of Delta marsh formation.

Franks Wetland Simulation
We estimated mineral deposition rates for input to the model using geochemical data from cores collected during the study. For model organic matter deposition and decomposition rates, we used values from the literature to initially estimate and place bounds on model input values. Data for age, porosity, bulk density and organic-matter content described in Drexler et al. (2009aDrexler et al. ( , 2009b were used to estimate inputs and for model calibration. We also used lead-210 ( 210 Pb) dating to identify the peat horizons deposited in 1963 and 1953, and we used heightened concentrations of the heavy metals Pb and Hg at the initiation of the California Gold Rush to estimate the depth of the 1850 horizon (Alpers et al. 2008). For simulation of Franks Wetland accretion from 4,700 years (calibrated years before present, BP) to within the last 350 years, we modified the Callaway et al. (1996) model Fortran code to accept varying annual inputs for mineral and organic matter.
For organic matter input for the 4,700-year simulation, we used a calibrated value of 0.025 g cm -2 yr -1 (Table 2), which was constant throughout the simulation period until 350 years BP. Core data (Drexler et al. 2009b) demonstrated significant accretion rate and peat compositional changes within the last 350 years. We, therefore, modified the model to accept time-varying inputs for first-year decomposition rates, organic and mineral inputs, initial and final porosity, and subsidence rates. We used measured porosity and bulk density data to estimate values for initial porosities (h2oin) and limiting porosities (porelim) ( Table 1). The following sections describe data and model inputs.

Age Dating
We compared data shown in Appendix A (Table A-2) with model results for Franks Wetland to calibrate the model. Carbon-14 ( 14 C) data were from Drexler et al. (2009b). In the San Francisco Estuary, lead contamination from hydraulic mining occurred from 1852 to 1884 (Hornberger et al. 1999;Domagalski 2001;Heim et al. 2007;Bouse et al. 2010). Alpers et al. (2008) established the approximate elevation of the 1850 horizon by determining the earliest occurrence of lead contamination (normalized to titanium to account for variable sediment input through time) at -29 cm MSL ( ± 4 cm) in the Franks Wetland core. Lead-210 was used to determine the elevation of sediments deposited in about 1963 and 1953 by applying the constant rate of supply model in cores described in Drexler et al. (2013).   (Miller and Fujii 2010) Surface organic matter deposition (g cm -2 ) Orgin Calibrated for Franks Wetland using data reported in Drexler et al. (2009b). Values for Twitchell Island were from Miller and Fujii (2010) Below-ground organic matter production (g cm -2 ) undpro Data for Twitchell Island impounded marsh. For Franks Wetland, we assumed undpro equaled orgin. For Twitchell, we varied undpro over time from 50% to equal to orgin as the marsh matured based in data provided by the USGS.
Decomposition rate constants for year 1, 2, 3, and later (yr -1 ) k decomp ,1,2,3 For Franks Wetland we specified the value of 1.0 based on the literature review. For Twitchell, we used spatially and temporally variable values determined from decomposition bag experiments reported by Miller and Fujii (2010). For modeling during the last 4,700 years on Franks Wetland, decomposition declined to close to zero within 100 years after introduction to the marsh.
Initial pore space (fraction) h2oin For Franks Wetland, we used the largest reasonable porosity values from data reported in Drexler et al. (2009a). For Twitchell we estimated values from bulk density data and literature values for organic and mineral particle density.
Final pore space (fraction) porelim Core data from Franks Wetland and Twitchell Island marsh. For recent accretion simulation for Franks Wetland, we varied porelim with time to accurately simulate measured bulk density. Input values were constrained by data for Franks Wetland.
Consolidation constant (unitless) k cons A final calibrated value of 10 was used for all simulations based on model calibration and simulation of mesocosm accretion described in Appendix A.
Sea level rise (cm yr -1 ) slr Specified as 0.17 cm yr -1 for all Franks Wetland simulations based on Atwater et al. (1977). No sea level rise was included for modeling at the Twitchell site because the site is impounded.
Subsidence (cm yr -1 ) from gas withdrawal subsid Calculated from 210 Pb data, Rojstaczer et al. (1991) and information about gas production and pressures. See Appendix A for explanation of development of time-variable subsidence rates.

Surface Organic Inputs and Sediment Deposition
For the Franks Wetland simulation, we reviewed the available data and literature for surface organic inputs. This included data from the Twitchell Island demonstration project (Miller and Fujii 2010) and other areas with similar climate and vegetation (Table A-1). Average organic matter deposition in the Typha (cattail)-dominated Twitchell Island demonstration project for East and West ponds ranged from 0.05 to 0.6 g cm -2 yr -1 from 1998 to 2006. During a 2-year study, Reed (2002) reported rates that ranged from 0.65 to 2.83 g cm -2 yr -1 for tidal freshwater marshes in the western, northern and central Delta. Reed (2002) reported consistent vegetative growth among study areas, dominated by Schoenoplectus (bulrush), Typha, and Phragmites (common reed) species. All other literature values for tidal marshes were substantially lower than those reported by Reed (2002) and Miller and Fujii (2010) (Table A-1). For tidal marshes, the average of literature values was 0.22 g organic matter cm -2 yr -1 . Excluding the Reed (2002) data resulted in average value of 0.05 g organic matter cm -2 yr -1 . Reed (2002) reported inorganic deposition rates for tidal marshes in the Delta ranging from 1.8 to 22.6 g cm -2 yr -1 . To provide a quasi-independent estimate of model inorganic sediment input, we used titanium data from Alpers et al. (2008) for Franks Wetland. Titanium is strongly correlated with the mineral, non-organic fraction of the peat matrix (Alpers et al. 2008) and Kabata-Pendias and Pendias (1992) reported that it does not accumulate in plant tissue; therefore, we used it to estimate the mass of mineral sediment in cores (see Appendix A).

Surface Organic Matter Decomposition and Underground Production
The model simulates organic decomposition as an exponential decay function (Equation 1).

Fraction of initial mass remaining
The value of k decomp is typically determined from measurements of decomposition bags that are harvested at varying times after initial placement. We surveyed the literature for values for k decomp for similar climatic conditions and plant species; primarily Typha, Schoenoplectus, and Phragmites species (see Table A For the Franks Wetland simulation, we assumed underground production rates equal to above ground organic production. This is based primarily on the data for the Twitchell Island demonstration project which indicated rates of underground production similar to above ground production in recent years (Miller and Fujii 2010).

Sea Level Rise
Estimates of sea level rise during the Holocene varied from 1.33 to 1.7 mm yr -1 . In their study of cores collected in San Francisco Bay, Atwater et al. (1977) showed a rate of sea level rise of about 1.7 mm yr -1 from 6,000 years to the present. Using a compilation of radiocarbon dates from the Delta and other locations in central California, Rosenthal and Meyer (2004) reported 1.33 to 1.54 mm yr -1 during 7,000 years BP. Gornitz et al. (1997) estimated relative rates of sea level rise during the 20th century on the West Coast of the United States to be from 1.4 to 1.5 mm yr -1 . For future simulations, we used global sea level rise projections from Meehl et al. (2007).

Deep Regional Subsidence
Since the mid-1900s, western Delta tidal marsh accretion has been affected by deep regional subsidence from natural gas withdrawal in the Rio Vista Gas Field, which underlies Twitchell, Bradford, Ryer, Staten, Jersey and Bethel islands ( Figure 1) (Rojstaczer et al. 1991) and by possible consolidation of Quaternary sediments (Brooks et al. 2012). The Rio Vista Gas Field is the largest gas field in California and has been under production since 1936. Production peaked in 1946 and then declined (Cummings 1999). (See data in Appendix A).
Although gas withdrawal-induced subsidence is relatively rare, significant compaction can occur in relatively shallow, thin, and wide reservoirs (Martin and Serdengecti 1984) such as those in the Rio Vista Gas Field (Rojstaczer et al. 1991;Oldenburg et al. 2001).
Delta marshes have remained at sea level so vertical accretion has offset sea level rise and deep regional subsidence. Using cesium-137 ( 137 Cs) (Delauney et al. 1978) to date sediment cores at undisturbed sites on Delta channel marsh islands, Rojstaczer et al.
(1991) estimated regional subsidence rates of 0.35 to 0.54 cm yr -1 from gas and groundwater withdrawal from 1963 to 1989. For cores collected in Franks Wetland, Rojstaczer et al. (1991) showed a subsidence rate of 0.4 cm yr -1 from 1963 to 1989. They attributed this subsidence to groundwater or gas withdrawal. However, Kerr and Leighton (1999) showed about 0.05 cm yr -1 of subsidence from groundwater withdrawal at extensometers on Bethel and Bacon islands. We assumed that gas withdrawal was the primary cause of subsidence and used 210 Pb data (see Table  A-2) and Rio Vista gas field production to estimate temporally variable subsidence rates since 1936 (see Appendix A for details).
where c is the model-calculated density of the material above the cohort and k cons (unitless) controls the rate of volume decrease. For marshes dominated by organic matter deposition, k cons , is greater than 1 (Callaway et al. 1996). We obtained a calibrated value of k cons of 10 from bulk density data described in HydroFocus, Inc. (2007) and Appendix A.

Effect of Accretion on Hydraulic Force on Delta Levees and Seepage
Accretion on subsided Delta islands will reduce the force upon-and seepage through and under-levees.
To approximate the potential benefit of vertical accretion to Delta levees, we first estimated the change in the cumulative force (CF in Newtons) on levees (Mount and Twiss 2005): where ρ is the density of water (1,000 kg m -3 ), g is gravitational acceleration (9.8 m s -2 ), H is the difference between the average channel water surface elevation (including sea level rise) and the average elevation of the island and L is the levee length for the island. Each levee is considered as a wall in which H determines the magnitude of hydrostatic force.
We estimated the under-seepage hydraulic gradient in drainage ditches nearest levees on Twitchell Island. As land-surface elevations decrease and sea level rises, hydraulic forces-which cause seepage under levees through underlying mineral sediments to drainage ditches-increase. Exit gradients in drainage ditches are an indicator of the potential for underseepage that can result in levee foundational piping: where WLE td and WLE ua are the water level elevations in the toe drain and underlying aquifer, respectively, and D is the vertical distance from the drainage ditch bottom to the bottom elevation of the organic soil.

Twitchell Island Simulation
Data collected from 1997 to 2008 provided temporally variable model input values for surface organic and inorganic marsh inputs, decomposition, underground production, and consolidation constants. Data used for model calibration and validation included sedimentation-erosion-table-measured elevation changes (Boumans and Day 1993), inorganic and organic inputs, decomposition rates and bulk density and organic matter content of accreting materials. Using weirs, water depths in the west and east ponds in the 6-ha Twitchell Island demonstration project were maintained at approximately 25 and 55 cm . The marsh received a constant supply of San Joaquin River water. Using data from Gamble et al. (2003), we calculated an average sediment deposition rate of 0.03 g cm -2 yr -1 . As the marsh vertically accreted, weir crest heights increased to maintain approximately constant water depths.
We used the Twitchell vertical accretion model and a Geographic Information System (ArcGIS) to predict future accretion for the entire Delta. We used land-surface elevation data collected by the California Department of Water Resources using LiDAR (Light Detection and Ranging) technology during January and February 2007. We used the average of the two intermediate sea level rise projections of 0.395 cm yr -1 from Meehl et al. (2007) for global sea level rise to the end of the 21st century. To estimate time to reach tidal elevations, we assumed that this rate of sea level rise rate would continue beyond the 21st century for estimating time to reach tidal elevations.

Consolidation
In the Callaway et al. (1996) model, consolidation is equal to For exit gradient calculations, we used groundwater level data reported in Deverel et al. (2007) and, based on field observations, approximated the elevation of the water in the drain ditch at the base of the levee (toe drain) at 1.3 m below land surface. We also assumed that the drain-water surface elevation will decrease with time because drainage ditches are typically deepened through time to compensate for on-going subsidence. To estimate future subsidence, we used the SUBCALC model (Deverel and Leighton 2010) and initial soil organic matter content reported for the soil nearest the levee at locations shown in Figure 1.

Franks Wetland Accretion
We simulated vertical accretion starting 4,700 calibrated years BP (Figures 3 to 5). The average modeled rate of vertical accretion during this period was 0.12 cm yr -1 . The measured bulk density and organic matter content were relatively consistent: between -80 to -500 cm elevation mean sea level (MSL). Consistent with relatively quiescent accretion, bulk density varied from 0.04 to 0.115 g cm -3 in this part of the peat profile ( Figure 4). Larger values below elevation -500 cm were from the influence of allochthonous sediments deposited from the greater watershed ( Figure 4), (Drexler et al. 2009a). The mean simulated bulk density value for the -80 to -500 cm elevation (MSL)-interval was 0.08 g cm -3 . For the same depth interval, the measured organic matter fraction varied from 0.62 to 0.93 ( Figure 5). The mean simulated value was 0.83.
Between elevation -80 cm (MSL) and land surface elevation (27 cm MSL), marsh sediment input, bulk density, and the average accretion rate were larger than in the underlying materials (Figures 3, 4, and 5). Also, subsidence from gas withdrawal began within the last 80 years. The estimated 1850 elevation (Figure 3) of -29 cm MSL is slightly shallower than the -39 cm MSL depth where bulk density initially increased ( Figure 4) and organic matter content decreased ( Figure 5). The elevation of the maximum bulk density (0.46 g cm 3 ) and minimum organic matter content (19.6%) was 0.7 cm (Figures 4 and 5). Simulated sediment input increased starting in the 1870s and peaked in the 1930s. Sea level rise was simulated as 0.17 cm yr -1 (Atwater et al. 1977). Model organic inputs varied from 0.001 to 0.21 g cm -2 yr -1 (Table 2).
Maximum bulk density and minimum organic matter content (Figures 4 and 5) were measured at elevations consistent with the estimated 1850 horizon. Within  Figure 6 shows the simulation results and measurements for elevation change. The 12 years of simulation and measured elevation data demonstrate temporally variable conditions associated with marsh establishment and subsequent plant colonization and maturation. Vertical accretion rates were greater in the West Pond than in the East Pond before 2003. Accretion progressed slowly initially in the East Pond at about 1.7 cm yr -1 . Rates increased substantially to over 7 cm yr -1 between 2003 and 2005 and 5.7 cm yr -1 between 2005 and 2008. Accretion in the West Pond was about 3 cm yr -1 during the entire period of record.

Twitchell Island Impounded Marsh
We achieved good agreement between simulated and measured elevation changes, bulk density, and organic matter content (Figures 6 and 7 and Table 3). Medians of measured bulk density values were larger in samples collected from the 15-to 25-cm interval relative to the 0-to 15-cm interval. This effect was more pronounced in the West Pond samples. Overall, lower bulk density values were measured in the East Pond, where water depths were deeper. Core samples from the West Pond had lower organic matter contents than those from the East Pond and organic matter content decreased with depth. In the East Pond, the median value for organic matter fraction was slightly lower in the 15-to 25-cm interval samples relative to the 0-to 15-cm interval samples.
We used different input variables for the East Pond and West Pond simulations (Table 2). For the West Pond, because of the shallow water levels and apparent higher sediment input, initial and final porosities were lower and mineral input was larger. Organicmatter decomposition rates were also greater for the West Pond, according to data provided by the USGS (2008 data file from Robin Miller, USGS-CWSC, to Steven Deverel, unreferenced, see "Notes"; Miller and Fujii 2010).
Using flow and suspended solids data collected by the USGS-CWSC (Gamble et al. 2003;Miller et al. 2008) for 1997 to 2001 in the Twitchell Island demonstration project, we varied the model sediment input to fit the bulk density and organic matter content of the cores because two additional factors affected sediment input and bulk density. Using the range of flow and suspended sediment data from Gamble et al. (2003), we estimated that sediment input likely varied from 0.01 to 0.07 g cm -2 yr -1 . Second, as the marsh developed, there was perturbation and uplifting by plants of the underlying, original peat. This underlying soil has lower organic matter content (18% to 24%) (Fleck et al. 2004) than the accreting material in the marsh (Table 3). This caused higher bulk density in deeper core samples especially in the West Pond (Table 3). Vegetation was denser in the West Pond, which probably resulted in more uplifting of the underlying organic soil during early years causing a greater bulk density change relative to the East Pond.
To account for these effects, we therefore adjusted the sediment input to as high as 0.2 g cm -2 yr -1 during the early simulation years.
We used USGS-CWSC plant biomass measurements Miller and Fujii 2010) in fully vegetated areas for marsh organic matter input (orgin). There were lower organic matter inputs to the East Pond even though more accretion was observed there, because of lower decomposition rates and higher porosities. We used measured decomposition rates as model inputs. Higher decomposition rates were measured earlier in the marsh history.
Using core data from Miller et al. (2008)   In light of generally good agreement between simulated and measured results for Twitchell Island and Franks Wetland, we estimated future elevation changes in impounded marshes on Delta subsided islands. We assumed that impounded marshes for subsidence mitigation and carbon sequestration will be managed similarly to the Twitchell Island demonstration project and will contain similar wetland vegetation.

Future Simulation of Impounded Marshes on Subsided Islands
We estimated 300 and 200 years for the West and East ponds to reach current sea level on the center of Twitchell Island which is about -630 cm. For these predictions, we used the model inputs for the most recent 4 years. These likely are more representative of long-term decadal accretion because decomposition rates slowed and accretion rates increased relative to initial marsh development. Because of faster measured accretion for the East Pond, we performed additional future simulations using the input parameters for the East Pond. The simulated carbon longterm accumulation rate was about 13 t C ha -1 yr -1 .
The extent to which accreting sediment will consolidate as it transforms from loose and highly porous material collected from the Twitchell Island impounded marsh to consolidated peat soil (see Figures A-3 and A-4 in Appendix A) is a key uncertainty in estimating future elevation change. Consolidation occurs as organic material decomposes, particles rearrange and pore space decreases. The key model variables that affect the rate and extent of consolidation are k cons (Equation 2) which controls the rate at which the deposited material consolidates and porelim, the limiting porosity.
The low bulk density, low mineral-and high organic matter-content peat that accumulated at Franks Wetland before at least 350 years BP is probably physically similar to materials that will eventually accumulate in an impounded marsh. Therefore, Franks Wetland core data and model simulations provide insight into probable values for model inputs that affect consolidation (initial and limiting porosities, h2oin and porelim, and k cons ) for future impounded marsh simulations. We used Franks Wetland bulk density and organic matter content (OMC) data to estimate particle density and porosity inputs as follows. The mean bulk density of 0.08 g cm -3 measured in the Franks Wetland core (Drexler et al. 2009a) results in a porosity of 0.94 using Equation 5 (Hillel 1980). Porosity = 1 − bulk density/particle density (5) We estimated the material particle density as equal to OMC × 1.14 g cm -3 + (1-OMC) × 2.61 g cm -3 where 1.14 and 2.61 g cm -3 are the particle densities for organic and mineral fractions, respectively (DeLaune et al. 1983). Using Equation 5, the average bulk density of 0.03 g cm -3 measured for recently accreted material in the Twitchell Island  east pond and mean OMC of 92% results in a porosity of 97.6% which was equal to the initial input for porosity (h2oin) for the Twitchell Island East Pond simulation. Therefore, we set the value of the initial porosity (h2oin) and limiting porosity (porelim) to 97% and 94%, respectively, for the future simulation of accretion on impounded marshes. Consistent with the Franks Wetland simulation and analysis of consolidation in the Twitchell Island mesocosms, the consolidation constant, k cons was 10 (see Appendix A).
Modeled accretion is sensitive to input porosities and k cons (Callaway et al. 1996). To estimate uncertainty in accretion estimates we performed multiple simulations to assess how changing the values of these variables affected consolidation (Figure 8 and Table 4). The uncertainty for time required to reach current sea level on Twitchell Island varied from 106 to 257 years, based on results that used a reason-able range in porelim and k cons values (Table 4 and Figure 8). Callaway (1996) stated that k cons values greater than unity are applicable for organically dominated sediments and they used a maximum calibrated value of 250. Therefore, we assumed that these values set the range for our sensitivity analysis.
Model results (Figure 8) indicate varying degrees of reductions in future accretion rates. Data provided by the USGS (2008 data file from Robin Miller, USGS-CWSC, to Steven Deverel, unreferenced, see "Notes") support this; the mean rate of Surface Elevation Table  ( For the 202-year estimate (run 11) to reach tidal elevations, simulation results indicate that rates will decline and stabilize at about 2.9 cm yr -1 within about 40 years. Using larger values for k cons resulted in a slower accretion rate decline. For k cons equal to 250 (runs 15 and 16), modeled rates declined to about 5 cm yr -1 within 100 years. For the lower porelim value of 92.5%, rates declined to 3.3 and 2.2 cm within 100 and 40 years for k cons values of 250 and 10, respectively. For all simulations, we assumed 0.03 g cm -2 yr -1 of mineral input consistent with measured values for the Twitchell Island impounded marsh. We predict that accreted material will be over 95% organic matter and have a bulk density of less than 0.1 g cm -3 .
Model results indicate that large areas in the periphery of the Delta could be restored to projected sea level within 50 to 100 years ( Figure 9). This includes Ryer, Grand, Tyler and Staten islands and in the southwestern, northern, eastern and southeastern Delta. Most of the central Delta would require 150 to 250 years.
In the western Delta where there is heightened interest in subsidence mitigation because of water supply vulnerability, model results indicated that a large portion of Sherman Island could be restored to projected sea level within 50 to 150 years. Deeply subsided portions on the southwestern and southeastern parts of Sherman would require 150 to 200 years to reach sea level. Large portions of Jersey and Bethel islands could be restored to sea level within 51 to 100 years. Small areas on Bradford, Twitchell, and Brannan islands, and Webb Tract would accrete to sea level in 51 to 100 years, but most of the area on these islands is deeply subsided and would require 150 to 250 years.  Table 4.

Hydraulic Force and Exit Gradients
We used Equation 3 to estimate the change in the cumulative force on levees (CF) from subsidence and impounded marsh accretion on Twitchell Island. Using 3.8 m as the initial average value for H and 19,000 m for L for Twitchell Island, we estimated the percent reduction in CF with time as an impounded marsh accretes similar to the East Pond during 50 years. The estimated percent reduction in CF with time after flooding of an impounded marsh using the results for simulation 11 (Table 4) is 13% of the current value from 2007 to 2057 ( Figure 10). In contrast, with continuing business as usual (i.e., subsidence), we estimated the CF will increase by 44%.
Two areas on Twitchell (Figure 1) are vulnerable to under-seepage because peat soils are: (1) high in organic matter content; (2) subsiding between 1 and 2 cm yr -1 ; (3) near the levee; and (4) relatively thin (less than 4 m). We estimated that exit gradients will approach critical levels for under-seepage sediment movement during the next several decades ( Figure 11). The Delta Risk Management Strategy levee vulnerability analysis reported (URS Corp. and J. R. Benjamin & Assoc. 2008) increased risk of levee instability as exit gradients approach -1.0. Exit gradients more negative than -1.0 can be associated with sand boils, which result from movement of levee foundation materials towards drainage ditches. In the east area, we estimated that exit gradients will approach -1 by 2036 ( Figure 11) under business as usual. In contrast, we estimate that accreting wetlands similar to the East Pond on Twitchell Island near the levees will result in increasingly positive exit gradients (Figure 11).

DISCUSSION
We provided unique insights about temporally and spatially variable processes that affect vertical accretion rates and future impounded marsh development in the Delta. In contrast to previous modeling efforts, inorganic sediment and organic inputs were accounted for in more detail and substantial data were avail-able for model development, calibration and verification. There was good agreement between simulated and measured vertical accretion rates, bulk density and organic matter content.
Tidal marsh simulation (Franks Wetland) demonstrated that organically dominated deposition occurred during most of the last 4,700 years BP, and increased sediment deposition influenced accretion during and after the Gold Rush period. Also deep subsidence likely from gas withdrawal, influenced accretion during the last 80 years. This simulation provided insight for the Twitchell Island impounded marsh simulation. Using the Twitchell Island simulation results, we estimated time-frames to reach tidal elevations throughout the Delta. We also demonstrated the benefits of impounded marsh creation to levee stability.

Franks Wetland
Accretion during most of the last 4,700 years BP was in response to sea level rise of about 0.17 cm yr -1 . Accretion within the last several hundred years was affected by more mineral input than during the previous several thousand years (Figures 4 and 5). Peak sediment input occurred during the late 1800s to early 1900s. Within the past 80 years, deep regional subsidence affected accretion.
Geochemical and physical data (Figures 4 and 5) indicate that high sediment loads during the post-1850 gold mining period resulted in a reduction in the organic component of the peat and an increase in bulk density of the peat matrix. There was also likely movement of mineral material into the sediment profile below the 1850 horizon as evidenced by geochemical data presented in Alpers et al. (2008). The model simulated these changes using varying mineral input to the accreting sediment column, which resulted in changes with depth similar to those that were measured. The actual cause of the measured changes with depth was deposition of relatively highly dense mineral material moving in the Delta during the latter half of the 19th century. Deposition of Gold Rush sediments has been reported in marshes throughout the San Francisco Estuary (May 1999;Goman and Wells 2000;Reed 2002;Drexler et al. 2009aDrexler et al. , 2009bDrexler et al. , 2011. As indicated by the relatively high lead/titanium ratio at an elevation of about -65 cm (Alpers et al. 2008), there was likely downward movement of mineral material into the profile. We estimated that the sediments above about -59 cm MSL (within 82 cm of land surface) were deposited within the last 350 years.
There has been little previous attention to regional subsidence from gas withdrawal and its effects on vertical marsh accretion. The analysis presented in Rojstaczer et al. (1991) indicated that it has been a significant process resulting in 0.29 to 0.61 cm yr -1 of regional subsidence. Their analysis included 137 Cs data for Franks Wetland described above which indicated a subsidence rate of 0.4 cm yr -1 . Model results for Franks Wetland indicate marsh communities compensated for subsidence through more sediment and organic-matter accumulation. The literature indicates similar responses in other locations. For example, Atwater et al. (2001) and Cahoon and Lynch (1997) presented evidence for increased vertical accretion in response to sudden and gradual subsidence.

Twitchell Island
The Twitchell Island impounded marsh facilitated greater organic accumulation and faster accretion than tidal wetlands such as Franks Wetland.
Model results indicate that accretion rates will likely decrease during the next decade because of the consolidation of accumulating sediments. Four variables are pivotal to effective simulation of accretion in the varied Delta environments: inorganic and organic inputs, porosity and variables that control consolidation. Core data collected by Drexler et al. (2009a, and Miller and Fujii (2010) provided bounds for model inputs for porosity as well as inorganic and organic inputs and data from the Twitchell Island mesocosms provided some confidence in the values of k cons . Because accretion estimates are sensitive to porosity and variables that control consolidation, we estimated substantial variance in times to reach tidal elevations to account for uncertainty in input variables that affect consolidation. Consolidation of deposited organic and inorganic sediments will affect future accretion rates.
Decomposition of organic detritus results in structural rearrangement of organic particles and reduction of pore space which leads to consolidation (Hobbs 1986). The extent of consolidation of highly organic deposits is a function of the effective stress which depends on the buoyant support and weight of the deposits (Pizzuto and Schwendt 1997). Pizzuto and Schwendt (1997) estimated consolidation rates in peat deposits in Delaware of about 1 mm yr -1 during the last 6,000 years. Tornqvist et al. (2008) estimated consolidation rates of 5 mm yr -1 in highly organic Holocene deposits in the Mississippi Delta. Appendix A provides visual evidence for substantial physical change and consolidation during burial of sediments in the marsh (see Figures A-3, A-4).

Delta-Wide Vertical Accretion and Benefits to Levee Stability
The coupling of GIS analysis and the vertical accretion model provided a spatial distribution of times required for accretion to future sea level (Figure 9). Time-periods for accretion to sea level ranged from several decades to over 250 years. Shorter accretion times coincided with less subsided areas on the periphery of the Delta. Longer times coincided with more deeply subsided parts of the western and central Delta and central parts of islands.
Our results demonstrate the difficulty of developing solutions within the next several decades to reduce the risk of levee failure and to create habitat on deeply subsided islands (Figure 9). This is especially true in the western Delta where levee failure can result in a major disruption of the state's water supply (e.g., Suddeth et al. 2010). The Delta Risk Management Strategy (DRMS) Phase 1 February 2009Report (CDWR 2009 estimated the combined probability of levee failure and island flooding from earthquake, high-water flooding, and sunny-day levee failure for most of the deeply subsided central and western Delta islands ranges from 53% to over 84% during the next 25 years. The estimated economic costs of future multiple Delta levee failures were in billions of dollars. Levee failure or instability from seepage occurs when hydraulic gradients are large enough to erode or move levee internal and/or foundational materials. Continuing subsidence and sea level rise will increase seepage. Deverel et al. (2007) predicted that by 2050 seepage onto Twitchell Island will increase by 22% to 34% under a "business-as-usual" scenario. Future risk of levee failure will be reduced if entire islands or areas near levees are converted to impounded marshes for subsidence reversal from reduced cumulative hydraulic forces (CF) and hydraulic gradients. If wetlands cover all of Twitchell Island, Deverel et al. (2007) estimated that seepage will be reduced by 14% to 24% relative to present-day conditions. We used CF as an indicator of factors that affect levee stability. We estimated that impounded marsh accretion will cause CF to decrease to 87% of its current value in 50 years whereas subsidence would cause it to increase by 44% during the same time period (Figure 10). We also estimated that restoring wetlands near levees will result in increasing positive exit gradients during the next several decades whereas continued subsidence will result in exit gradients that approach critical values for levee foundational piping during the same time-period ( Figure 11). Though localized factors and hydraulic gradients also affect levee stability and require consideration, the risk of levee failure will drop substantially as a result of vertical accretion in impounded marshes because subsidence will stop, land-surface and island water levels will accrete, land-owners will no longer be required to deepen drainage ditches, and hydraulic gradients will decrease. Accretion in wetlands for subsidence reversal will likely provide a greenhouse gas benefit. We estimated long-term carbon sequestration rates for impounded marshes such as the Twitchell Island project of about 12 to 15 t C ha -1 yr -1 . Permanently flooded wetlands on organic soils will stop the loss of carbon dioxide (median = 4.6 t C ha -1 yr -1 ) from organic matter oxidation (Deverel and Leighton 2010).

SUMMARY AND CONCLUSIONS
As part of a study on peat accretion processes in the Sacramento-San Joaquin Delta, we simulated marsh vertical accretion at Franks Wetland and the Twitchell Island impounded marsh. We used a modified cohort-accounting model (Callaway et al. 1996), together with physical and chemical data collected during the study, and from the literature, to develop model inputs. For Franks Wetland, we simulated accretion since 4,700 calibrated years BP using 14 C data for model calibration. We also simulated recent accretion and compared estimates with 210 Pb data, lead, and titanium concentration data. For the Twitchell Island impounded marsh, we calibrated and verified the model with sedimentation-erosion table results and core data collected by the USGS from 1997 to 2008. For both islands, our model results compared favorably with measured accretion, bulk density, and organic matter content.
Accretion on Franks Wetland was predominantly organic with varying degrees of inorganic sediment input. From 4,700 to 350 years BP, accretion on Franks Wetland was about 0.12 cm yr -1 , consistent with accretion rates in tidal wetlands worldwide. Franks Island received relatively little inorganic sediment during this time period. Our simulations of recent accretion indicate heightened sediment influx during and after the onset of the California Gold Rush period.
At Twitchell Island, all organic and inorganic inputs stay in the confines of the marsh. This has resulted in greater organic inputs and average accretion rates as high as 9.2 cm yr -1 . Future simulations indicate that managed impounded marshes will accrete highly organic (over 90%), low bulk density material at long-term average rates of about 3 cm yr -1 . The extent of future consolidation is a key uncertainty in determining rates of future accretion.
Extension of the results to other Delta islands indicates that large areas of the periphery of the Delta could be restored to tidal elevations within 50 to 100 years by creating impounded marshes. Most of the central Delta would require 50 to 250 years to be restored. A large portion of the western Delta could be restored to sea level within 50 to 150 years (large areas on Sherman, Jersey, and Bethel islands, and small areas on Bradford, Twitchell, and Brannan islands and Webb Tract). Creation of impounded marshes on Delta islands can provide benefits to levee stability as demonstrated by cumulative force and hydraulic gradient calculations.