Modeling Tidal Freshwater Marsh Sustainability in the Sacramento–San Joaquin Delta Under a Broad Suite of Potential Future Scenarios

In this paper, we report on the adaptation and application of a one-dimensional marsh surface elevation model, the Wetland Accretion Rate Model of Ecosystem Resilience (WARMER), to explore the conditions that lead to sustainable tidal freshwater marshes in the Sacramento–San Joaquin Delta. We defined marsh accretion parameters to encapsulate the range of observed values over historic and modern time-scales based on measurements from four marshes in high and low energy fluvial environments as well as possible future trends in sediment supply and mean sea level . A sensitivity analysis of 450 simulations was conducted encompassing a range of porosity values, initial elevations, organic and inorganic matter accumulation rates, and sea-level rise rates. For the range of inputs considered, the magnitude of SLR over the next century was the primary driver of marsh surface elevation change. Sediment supply was the secondary control. More than 84 % of the scenarios resulted in sustainable marshes with 88 cm of SLR by 2100, but only 32 % and 11 % of the scenarios resulted in surviving marshes when SLR was increased to 133 cm and 179 cm, respectively. Marshes situated in high-energy zones were marginally more resilient than those in low-energy zones because of their higher inorganic sediment supply. Overall, the results from this modeling exercise suggest that marshes at the upstream reaches of the Delta—where SLR may be attenuated—and high energy marshes along major channels with high inorganic sediment accumulation rates will be more resilient to global SLR in excess of 88 cm over the next century than their downstream and low-energy counterparts. However, considerable uncertainties exist in the projected rates of sea-level rise and sediment avail-ability. In addition, more research is needed to constrain future rates of aboveground and belowground plant productivity under increased CO 2 concentrations and flooding.


INTRODUCTION
Tidal freshwater marshes rely on accumulation of both inorganic sediment and organic matter on the marsh plain, which results in vertical accretion (cm per year), to maintain their position in the tidal frame (Reed 2000;Neubauer 2008;Drexler et al. 2009a;Drexler 2011). In contrast to salt marshes, because of their location at the landward end of estuaries, most tidal freshwater marshes are strongly influenced by fluvial processes and have muted tidal ranges (Odum 1988;Day et al. 1995). In addition, there is less robust feedback between the elevation of the marsh surface and marsh accretion in tidal freshwater marshes than in salt marshes (Kirwan and Guntenspergen 2012), making them particularly susceptible to drowning from relative sea level rise (Reed 1995).
Although recently there has been considerable work on peat formation processes in tidal freshwater marshes (e.g., Neubauer et al. 2002;Neubauer 2008;Drexler et al. 2009a;Drexler 2011), there has been limited application of marsh surface elevation models to project the future sustainability of such marshes. Craft et al. (2009) applied the Sea Level Affecting Marshes Model (SLAMM) to marshes in Gulf coast estuaries to evaluate the changes in areal extent of marshes along a salinity gradient. Within the SLAMM model, relative sea-level change is computed as the sum of the historic eustatic trend, the site-specific rate of elevation change from subsidence and isostatic adjustment, and the accelerated sea-level rise depending on the future scenario chosen. Vertical accretion and all the processes therein are not treated in a mechanistic fashion (Clough et al. 2010), so SLAMM cannot be used to infer processes at a localized scale. In the Sacramento-San Joaquin Delta of California (hereafter, the Delta), which is the most landward part of the San Francisco Estuary, Deverel et al. (2008Deverel et al. ( , 2014 hindcast vertical accretion in Franks Wetland, using an adaptation of the Callaway et al. (1996) model. The model applied constant or temporally variable mass accumulation rates to recreate peat accretion histories, but lacked the dynamic feedback between inundation and the magnitude of organic matter accumulation that approximates the natural self-regulating process responsible for sustaining marshes. Swanson et al. (2013) further adapted the Callaway et al. (1996) model in the Wetland Accretion Rate Model of Ecosystem Resilience (WARMER) to incorporate the appropriate dynamic feedback among vertical accretion processes, marsh surface elevation, and temporally variable sea-level rise for San Francisco Bay tidal salt marshes, but did not adapt or apply the model to tidal freshwater marshes of the Delta. Schile et al. (2014) examined the future resiliency of four marshes across the San Francisco Estuary using the Marsh Equilibrium Model coupled with digital elevation models; however no tidal freshwater sites in the Delta were included.
The Delta was once a 1,400 km 2 tidal marsh region, but, during the past 150 years, over 90% of the marshes have been drained for agriculture and the remaining natural ecosystems have been subject to a wide range of anthropogenic impacts (Atwater 1982;Cloern et al. 2011). The relict marshes provide critical habitat for flora and fauna (Foley et al. 2005). Future changes in global climate and anticipated operational and engineered changes to water diversions may further stress Delta ecosystems. In particular, recent work suggests that decreasing runoff and suspended sediment concentrations combined with rising sea level, water temperature, and salinity will present challenges to both natural and managed systems within the San Francisco Estuary, including the Delta (Cloern et al. 2011). Absent from these analyses are the implications for the long-term sustainability of tidal freshwater marshes in the Delta.
In this paper, we report on the application of WARMER to tidal freshwater marshes in the Delta in order to explore the future sustainability of these marshes under a broad range of possible conditions. We carried out a sensitivity analysis that includes 450 different future scenarios encompassing multiple values for each of the five input parameters: porosity, initial elevation, and rates of organic accumulation, inorganic accumulation, and sea-level rise. We derived inputs from field measurements in tidal freshwater marshes in the Delta, analysis of peat cores, and anticipated trends in sediment supply and sea-level rise in the region. Of the 450 scenarios, 120 scenarios correspond to marshes that receive high amounts of sediment because of their position in high energy, hydrogeomorphic zones, 90 scenarios correspond to marshes in low-energy hydrogeomorphic zones, and 240 scenarios that capture all of the other possible combinations of model parameters. In carrying out all these model runs, our main goal was to quantitatively evaluate the relative sensitivity of the model to the five input parameters.

Field Sites
The Delta region (Figure 1) has semi-diurnal, micro-tidal tides with a normal range of approximately one meter (Shlemon and Begg 1975;Atwater 1980). The climate in the Delta is characterized as Mediterranean, with cool winters and hot, dry summers (Thompson 1957). Mean annual precipitation is approximately 36 cm, but actual yearly precipitation varies from half to almost four times this amount. Over 80% of precipitation occurs from November through March (Thompson 1957).
The seasonality of precipitation and the timing of major storms, in particular, are the main drivers of suspended sediment in the two main rivers, the Sacramento and San Joaquin, which feed the Delta (Ingram et al. 1996;Wright and Schoellhamer 2004). Unlike salt marshes in the lower estuary, which receive sediment largely from tidal fluxes, sediment in the Delta is received predominantly from upstream riverine sources. Approximately two-thirds of the sediment that enters the Delta is deposited in its marshes, sloughs, and mudflats and, of this sediment load, more than 80% derives from the Sacramento River versus the San Joaquin River (Wright and Schoellhamer 2005).
Beginning in the mid-1800s, the Delta was largely drained for agriculture (Thompson 1957). This massive reclamation effort resulted in its current configuration of over 100 islands and tracts surrounded by 2,250 kilometers of man-made levees and 1,130 kilometers of waterways (Prokopovich 1985). Within the channels, small relatively undisturbed, remnant marshes remain, which not only provide critical habitat, but can also be used to study natural processes.
We chose tidal freshwater marsh sites to incorporate the various hydrogeomorphic settings and salinity regimes of the Delta (Figure 1, Table 1). We selected sites from high-energy, hydrogeomorphic environments such as the confluence of the Sacramento and San Joaquin rivers to more quiescent, low-energy environments, such as distributaries of the San Joaquin River. The four sites-Browns Island, Franks Wetland, the Tip of Mandeville Tip, and Bacon Channel Island-are all relatively undisturbed marshes, which, in contrast to the majority of the Delta, were not drained and converted to agricultural production during the late 1800s to early 1900s.
Browns Island is an oligohaline marsh (mean salinity ranges from 0.5 to 5 ppt) located at the confluence of the Sacramento and San Joaquin rivers . Because it is located at the ecotonal western border of the Delta and is one of the only remaining large expanses of marsh in the Delta, we have included it in this study and will from here on categorize it along with the following freshwater marshes for the sake of simplicity. The Tip of Mandeville Tip is a freshwater marsh situated along a main channel of the San Joaquin River. Franks Wetland is a freshwater marsh located along a distributary of the San Joaquin River behind a natural breakwater and adjacent to a permanently flooded farmed island. Bacon Channel Island is a fresh water marsh situated along a distributary of the San Joaquin River.
For the purposes of this paper, the four study sites can be grouped into two main categories based on Channel Island. The high-energy sites are located in either main channels or at the confluence of rivers and are the recipients of much greater amounts of sediment and have more variable sediment inputs than the low-energy sites, which are located in quiescent tributaries (Drexler 2011). Rates of peat formation-including vertical accretion, organic matter accumulation, and inorganic sediment accumulationare higher in high-energy sites than in low-energy sites. Low-energy sites rely more on organic matter accumulation than on inorganic sediment accumulation to maintain elevation of the marsh plain in the tidal frame (Drexler 2011).

Peat Coring
We used peat cores from each of the marsh study sites to estimate net mass accumulation (both organic and inorganic) rates for recent (the past 50+ years), as well as the entire lifetimes of the marshes in order to include the range of feasible mass accumulation rates. For Browns Island, we estimated recent and long-term (millennial) mass accumulation rates from two peat cores (Browns Island High cores) that Callaway et al. (2012) collected in the high marsh. These cores were approximately 0.3 km from the peat cores collected for the study by Drexler et al. (2009). Methods for the collection and analysis of the Browns Island High peat cores are provided in Callaway et al. (2012), and the data for the top ~30 to 40 cm of these cores are included in Table A1 (see Appendix A). Details concerning the collection and analysis of peat cores (one core per site encompassing the entire peat profile) from Browns Island, Franks Wetland, the Tip of Mandeville Tip, and Bacon Channel Island are provided in Drexler (2011), with the exception of methods for 210 Pb and 137 Cs dating of the top meter of each core, which follow.

Laboratory Analyses
We measured total 210 Pb, 226 Ra, and 137 Cs simultaneously by gamma spectrometry as described in Baskaran and Naidu (1995), Fuller et al. (1999), and Van Metre et al. (2004). We counted subsamples of dried peat samples using a high-resolution, intrinsic germanium well detector, gamma spectrometer. Samples were sealed in 7-mL polyethylene scintillation vials and placed in the detector bore-hole or well, which provides near-4π counting geometry. The supported 210 Pb activity, defined by the 226 Ra activity, was determined on each sample from the 352 KeV and 609 KeV gamma emission lines of the short-lived daughters 214 Pb and 214 Bi daughters of 226 Ra, respectively. We accounted for self-absorption of the 210 Pb 46 KeV gamma emission line using an attenuation factor calculated from an empirical relationship between self-absorption and bulk density that was developed based on the method of Cutshall et al. (1983). Self-absorption of the 214 Pb, 214 Bi, and the 661.5 KeV 137 Cs gamma emission lines was negligible for the well detector. We determined detector efficiency from National Institute of Standards and Technology (NIST) traceable standards. NIST and International Atomic Energy Agency reference materials were counted monthly to check detector calibration. We calculated the uncertainty in the measured activity to be typically within ±10% at the one standard deviation level based on the random counting error of sample and background spectra. The measured activities of replicate analysis of material from the same sample agreed to within ±15%. We decaycorrected measured activities of 137 Cs for the period between sample collection and analysis.
We found the overall magnitude of unsupported 210 Pb in the sediment profile to be very low, resulting from the very low atmospheric delivery of unsupported 210 Pb in the region (Fuller and Hammond 1983). Because low unsupported activities at depth often were within the uncertainty of the measurements, the total inventory of unsupported 210 Pb and its maximum depth were poorly defined, resulting in large associated errors in determining mass accumulation using the Constant Rate of Supply method (Robbins 1978). Thus, we applied the Constant Flux: Constant Sedimentation rate (CF:CS) method to estimate the mean (net) inorganic sediment accumulation rate at each of the study sites (Oldfield and Appleby 1984). We did not use the 137 Cs peak for dating, per se, but only to verify the general position of the 1963 6 horizon determined with 210 Pb dating. The 137 Cs peak was deeper than the 1963 horizon calculated using the CF:CS ( 210 Pb) method for two of the five Delta peat cores, and shallower for one of the two Browns Island cores.

WARMER
The Wetland Accretion Rate Model of Ecosystem Resilience (WARMER) is a 1-D cohort model of wetland vertical accretion based on the Callaway et al. (1996) model, hereafter the Callaway model. The Callaway model calculates the change in the elevation of the marsh plain relative to mean sea-level (MSL) at a point that represents the marsh surface as a function of accretion of inorganic and organic matter accumulation minus compaction, decay, and increases in relative sea level within each annual cohort. The Callaway model was previously applied in the Delta to project future rates of vertical accretion and carbon sequestration in impounded marshes on subsided Delta islands (Deverel et al. 2008(Deverel et al. , 2014. Swanson et al. (2013) adapted the Callaway model for San Francisco Estuary salt marshes by including temporally variable SLR, elevation dependent organic matter accumulation, and site-specific inorganic sediment accumulation functions. We set parameters in WARMER with site-specific field surveys, sediment core analysis, and inundation patterns. The elevation of the marsh surface, E, at time t (yr) relative to local MSL is calculated as where E(0) is the initial elevation relative to MSL, SLR(t) is the sea level at time t relative to the initial sea level and V i (t) is the volume per unit area, or height, at time t, of the cohort formed during year i. The total volume of an individual cohort is the sum of the mass of water, calculated from the porosity of the cohort, sediment, and organic matter content divided by the cohort bulk density.
In this study, we adapted WARMER for the slightly brackish and freshwater tidal marshes of the Delta by modifying the organic matter accumulation function to the elevation range and productivity of the dominant vegetation types within the marshes and allowing for temporally constant or variable inorganic sediment accumulation. Model parameters have been defined from previous (Drexler 2011) and new peat core analyses and from previous marsh surface elevation modeling (Deverel et al. 2008). Each of the model inputs is discussed in detail below. The model scenarios have been designed to explore a range of possible future conditions to assess the range of conditions that are favorable for marsh sustainability, not to forecast the evolution of marsh surface elevation at any single site. Mass accumulation rates and porosities for the low-and high-energy hydrogeomorphic settings are identified and the results of scenarios for each of these settings are discussed as well.

Initial Surface Elevation and Sea-Level Rise
We chose initial marsh surface elevations of 20, 30, and 40 cm relative to MSL to span current elevations of the marsh plain commonly found in the Delta (Table 1). We did not include the elevation of Browns Island (~50 cm) because it is higher than other marshes because of its position at the western boundary of the Delta, and, therefore, does not represent the region as a whole.
For the purposes of this paper, we considered SLR estimates for the Delta ranging from 43 to 179 cm from 2000 to 2100. This range is based on the SLR projections from two studies. The first is the semiempirical method of Vermeer and Rahmstorf (2009) from the International Panel on Climate Change's Fourth Assessment Report (IPCC 2007), which produced global SLR scenarios increasing from 81 to 179 cm from 2000 to 2100, with noticeable accelerations in the rate of change over the latter half of the century. The second is a downscaled range of projected SLR along the Pacific coast also based on the semi-empirical method of Vermeer and Rahmstorf (2009), which ranges from 43 to 160 cm . One high, two moderate, and one low sea-level rise scenario for 2000 to 2100 were included in the modeled scenarios based on the 43 to 179-cm range of projected SLR. Additionally, because the upper tidal reaches of the Delta may have attenuated SLR, we included a scenario that is 50% of the lowest sce-nario as a lower bound of potential SLR for the Delta region.
The nonlinear sea-level rise curves applied here were developed based on the limits of SLR described above, and the guidance proposed by the U.S. Army Corps of Engineers (USACE 2009)  (2) The five SLR curves for 21.5, 43, 88.3, 133.6, and 179 cm included in model scenarios are shown in Figure 2. and Schoellhamer 2005). At Mallard Island, which is situated just west of Browns Island, landward dispersive suspended sediment flux was only 20% of the total seaward flux from 1995 to 2003(McKee et al. 2006. Suspended sediment concentration in the rivers is highest during the winter storm season when large storms raise sediment loads and temporarily increase the depth and duration of inundation in Delta marshes. Suspended sediment concentration in the Sacramento River, the dominant source of sediment to the Delta, has decreased by 50% since the late 1950s as the remainder of high-sediment loads from historic gold mining operations has worked its way through the river system and the supply of sediment below dams has decreased (Wright and Schoellhamer 2004;Schoellhamer et al. 2013).
For WARMER, we used the two future scenarios for sediment supply to the Delta developed by Cloern et al. (2011): (1) a constant supply at present rates and (2) a continued decrease at the recent trend of -1.6% per year. We determined the constant and initial rates of net sediment accumulation from peat core analyses as described above and from Drexler et al. (2009b). This direct approach precludes calculation of net sediment accumulation from inundation, suspended-sediment concentration, settling velocity, and erosion, which vary in time and space and are, therefore, uncertain. In the cases of decreasing sediment accumulation, the mass of sediment accumulated per unit area per year at time t, M s (t), is equal to where M s (0), is the initial mass accumulation rate per unit area [M T -1 L -2 ].
The model was run for scenarios of constant high, moderate and low inorganic sediment accumulation described by the maximum, median, and minimum of the measured rates in Table 2. In addition to these three scenarios, the model was run for two scenarios with sediment supply decreasing at a rate of -1.6% per year. The lowest inorganic sediment accumulation scenario was not included in the temporally variable scenarios because the initial value was already extremely small.

Inorganic Sediment Accumulation
Unlike salt marshes in the lower estuary, which receive sediment largely from tidal fluxes, sediment in the Delta is received predominantly from upstream via the Sacramento and San Joaquin rivers (Wright

Organic Matter Accumulation
The magnitude of organic matter accumulation in salt marshes in the San Francisco Estuary is small relative to inorganic matter accumulation ) and is not a sensitive parameter in marsh surface elevation modeling in saline and brackish pickleweed (Sarcocornia pacifica) dominated marshes (Stralberg et al. 2011;Swanson et al. 2013). In contrast to salt marsh vegetation, plant communities in the tidal freshwater marshes of the Delta occupy the full ~1 m tidal range, are highly productive, and are major contributors to the building of peat (Atwater and Hedel 1976;Orr et al. 2003;Drexler et al. 2009bDrexler et al. , 2009aDrexler 2011;Watson et al. 2011).
WARMER was adapted by Swanson et al. (2013) to include dynamic feedback between marsh surface elevation-a proxy for tidal inundation regime-and organic matter accumulation rate, which is subsequently divided into aboveground and belowground organic accumulation. Such an adaptation was challenging in this current adaptation of WARMER, because little is known about the marsh surface elevation-organic matter accumulation feedback function in Delta marshes. Furthermore, no data were available on aboveground and belowground primary productivity in natural, fully tidal, freshwater marshes in the Delta. Morris et al. (2002) quantified a plant productivity function for Spartina alterniflora, which is the dominant salt marsh species along the east and gulf coasts of the United States. Recent research has shown that this function determined by Morris et al. (2002) is also the general shape of plant productivity for Schoenoplectus americanus and S. patens found in Maryland (Kirwan and Guntenspergen 2012). Because S. americanus is found in the Delta and is similar to the more broadly distributed S. acutus, we adapted the Morris et al. (2002) function to the elevation range of Delta vegetation; we then used it to estimate organic matter accumulation in the selected sites as an approximation for the marsh surface elevationorganic matter accumulation feedback function. As in Swanson et al. (2013), the parabolic function has been fit to the elevation range of the marsh vegetation in the Delta, which is from mean lower low water (MLLW) to mean higher high water (MHHW). This elevation range (also known as the great diurnal tide range) is consistent with the range of elevation for Delta marsh vegetation described at Bethel Island in the Delta by Atwater and Hedel (1976) and Watson and Byrne (2009).
The parabolic equation describing M o , the annual mass of organic matter accumulated per unit area [M L -2 T -1 ], at a given elevation, z, is: where a and b are constants with units of [M L -4 T -1 ] for aboveground and belowground production, respectively, fit to the measured organic matter accumulation rates in the surface layer (0 to 2 cm) of each peat core at the elevation of the marsh surface where the core was collected. This approach was used in Swanson et al. (2013) and a similar procedure was carried out by Deverel et al. (2008Deverel et al. ( , 2014 to simulate vertical accretion at Franks Wetland in the Delta. We are aware that the use of annual aboveground and belowground primary productivity data would have been preferable, however such data were not available for the dominant marsh species in the Delta and collection of such data was beyond the scope of this project. We determined the organic matter accumulation rate at the elevation of the core location from the inorganic sediment accumulation rate, M s , and the ratio of organic material to inorganic sediment in the surface layer of each sediment core, We determined organic matter accumulation rates for the entire lifetimes of each of the marsh sites using the radiocarbon age-dating of the peat cores (Drexler et al. 2009a). The peak organic matter accumulation rate for each curve, M o (MSL), is reported in Table 2. We used the minimum, maximum, and median peak organic matter production values to define the three modeled organic matter accumulation curves shown in Figure 3. We used a symmetrical great diurnal tide range of 1 m for all model runs as representative of the Delta environment where these marshes are located.
The total organic mass accumulated at each timestep was divided between aboveground and belowground productivity using a root:shoot ratio of 1:1 after Deverel et al. (2008). Aboveground organic mass was assumed to accumulate in the surface cohort. Belowground organic mass accumulation was distributed exponentially through the modeled soil column.

Compaction
Compaction and decomposition functions of WARMER follow the Callaway model and we set parameters according to the calibration by Deverel et al. (2008). We determined compaction of marsh sediment by a rate of decrease in porosity from the highly porous surface sediments from each peat core and a lower limit of porosity. Porosity is variable through the peat column and the peat is not highly compactible (Drexler et al. 2009b). For this reason we approximate the porosity of new surface material as the mean porosity plus the standard deviation. Because low porosities are often associated with very high sediment loads and the sediment supply to these marshes is low and possibly decreasing, the mean measured value represents a more realistic value for the limit of compaction. These values from Drexler et al. (2009b) are summarized in Table 3. The highenergy sites have lower porosity than the low-energy sites and there is reasonable agreement within these hydrogeomorphic classifications to use one representative porosity profile (surface and lower limit) for each hydrogeomorphic setting to define model inputs.
The rate of decrease, r, in the porosity of a given cohort is a function of the density of all of the material above that cohort such that: where ρ b is the density of the material above a cohort and k 1 is a calibration constant defined by Deverel et al. (2008).
matter accumulation (five rates), and SLR (five rates) described above for a total of 450 runs (Table 4). A copy of the source code text (see Appendix B), which was compiled using Compaq Visual Fortran, and sample input files for the 88.3 cm SLR suite of scenarios (see Appendix C) are included as separate files associated with this paper. Model inputs as defined for each hydrogeomorphic setting are shown in Table 5. Of the 450 scenarios, 120 scenarios correspond to high-energy sites, 90 scenarios correspond to low-energy sites, and the remaining 240 capture all other possible combinations of model parameters. The discrepancy between the numbers of high-and low-energy scenarios arises from the inclusion of temporally decreasing sediment accumulation for the high and medium accumulation rates, but not the smallest inorganic accumulation rate, which corresponds to the low-energy sites.
We then evaluated the model results to examine the timing and magnitude of marsh surface elevation changes for each input parameter and each hydrogeomorphic zone. A major objective of the study was to quantitatively evaluate the sensitivity of the model to the five input parameters. This evaluation was complicated by the fact that there were between two to five values for each input parameter, and for some input parameters there was no response in marsh surface elevation regardless of the parameter value, whereas for other parameters there was a significant response in elevation as the input parameter increased in value. Because of these complexities, we used the following approach to quantify the relative control of the various input parameters over the resulting final marsh surface elevations. For each input parameter p, we assigned percentile bins of each input value based on the number of values for that input. Marsh surface elevation results from the modeling exercise were grouped into the same bins. We calculated the fraction of marsh surface elevation results in the expected percentile bin, f p , was then calculated. If the results were independent of the parameter p (e.g., if the marsh surface elevation distributions for different values of an input parameter were identical), then f p = 1/n p , where n p is the number of input values for parameter p. If the results were completely sensitive to and varied monotonically

Decomposition
In WARMER, decomposition is assumed to follow a three-stage process where the youngest organic material less than one year old decomposes at the fastest rate, the organic matter one to two years old decays at a moderate rate, and the organic matter greater than two years old decays at the slowest rate. The percentage of refractory organic material is determined by dividing the mean percent organic matter at the top of a core by the mean percent organic matter in the bottom 4 cm of the core. The fraction of nonrefractory organic material, F M o , in each age class remaining after each time step is calculated as where d is the depth of the cohort and coefficients C (i) and k decomp (i) are coefficients for each age class, i, and are determined empirically by measuring the remaining mass of organic material in decomposition bags at varying times after initial placement. For this study, we used values for Typha spp., Schoenoplectus spp., and Phragmites australis under similar climates compiled by Deverel et al. (2008).

Model Application
WARMER was adapted as described above and applied for all combinations of porosity (two values), initial marsh surface elevation (three values), organic matter accumulation (three rates), inorganic  with the input parameter, f p = 1. For example, if the results were proportional to p, the percentile bin for a particular result would be identical to the percentile bin of that value of p. In order to compare the fraction of scenarios that fell within the expected bin for complete sensitivity between parameters with differing numbers of values, we adjusted f p by 1/n p to calculate a control parameter, C p : which is a value between zero, indicating no control, and one, indicating complete control, which can be used to compare the relative control across parameters with different numbers of input values (see Figure A1 in Appendix A for a detailed example). Once we identified most sensitive input parameter, we repeated the calculation for isolated suites of scenarios for each value of the most sensitive parameter to identify the next most sensitive parameter.

RESULTS
The WARMER model results show that marsh surface elevations relative to MSL in the Delta will decrease for most of the 450 scenarios. As an example, the 150 runs with an initial elevation of 30 cm are plotted in Figure 4 to demonstrate the variability in elevation for the scenarios. Eight scenarios with 21.5 or 43 cm SLR-the largest sediment input-and high porosity have accretion beyond the upper limit of modeled organic matter accumulation and up to the upper limit of inorganic sediment accumulation. At 1 m above MSL, further inputs to the marsh surface are eliminated, leading to the imposed maximum elevation that is observed for several scenarios (Figure 4). Subsequent decreases in elevation from SLR, compaction, and decay may lead to additional sediment input, but the modeled marsh surface elevation is not allowed to exceed 1 m. The decrease of relative elevation of the marsh surface in most scenarios accelerates in the latter half of the century when SLR accelerates appreciably (Figure 4).
Though most scenarios lead to some loss of marsh surface elevation, modeled marsh surfaces maintain a position within the range of marsh plant growth (MLLW-MHHW) for a large number of the 450 scenarios ( Figure 5). Of the 450 runs, 65.5% finish in 12 2100 with a marsh surface elevation above MLLW ( Figure 5A). The effects of most inputs appear additive and SLR is the primary control on final marsh surface elevation ( Figure 5A, Table 6).
The total number of scenarios in which marshes change from high marsh to low or drowned marsh increases as SLR increases. The suite of scenarios with a total of 88 cm of SLR terminate with a median elevation in the low marsh range, MSL-MLLW, and, of these, 84% of terminate above MLLW (Figures 6  and 7). Just 32% of the scenarios with 133 cm of SLR terminate with modeled marsh surface elevation above MLLW, indicating that 88 cm of SLR may represent a critical threshold of SLR for marsh survival (Figures 6 and 7). Five percent of scenarios transition to low marsh from initial high marsh elevations in 50 years with the median SLR (88 cm in 2100). Drowning of marshes (marsh surface elevation below MLLW) occurs only in the last 15 years of the modeled century for 88 cm of SLR but begins as early as 60 years in the modeled century for 169 cm of SLR and only 11% of marshes do not drown with this SLR forcing ( Figure 6).
For scenarios with the smallest final SLR (21 cm), which is quite similar to SLR over the previous century (17 cm, Equation 2), all 90 scenarios remain higher than MSL ( Figure 5A). Thus, given nearly historical forcing, the model indicates that existing marshes are sustainable. These results show that WARMER is reasonably reproducing current accretion processes in the study area.
As demonstrated by the sensitivity parameter, C p , SLR was the primary control on final marsh surface elevation for the suite of 450 scenarios considered (C SLR = 0.36; Table 6). When we repeated the calculation for each suite of 90 model runs for each SLR scenario independently, sediment supply was the secondary control on final marsh surface elevation (mean C sediment = 0.41; Table 6). The small differences in porosity and initial elevation did not result in large differences in final marsh surface elevation ( Figure 5A). When we repeated the sensitivity calculation a third time, isolating input values of each SLR and sediment supply (25 suites of 18 model scenarios each), we could not clearly define a single tertiary control on final marsh surface elevation. Instead, we found organic matter to be a sensitive parameter for only the smaller sediment inputs. Porosity, on the other hand, was a sensitive parameter for the higher sediment inputs scenarios (Table 6). We occasionally found that initial marsh surface elevation was a sensitive input parameter at low and moderate sediment inputs.
The high-and low-energy hydrogeomorphic zones follow the general trend of decreasing final marsh surface elevation with increasing SLR ( Figure 5B). The distribution of final marsh surface elevations for low-energy sites for a given SLR is lower and narrower than for high-energy sites. The low-energy sites have lower accumulation rates of both organic and inorganic material, and the higher porosity values associated with the peat cores from these quiescent environments does not offset the slower accretion when compared to the high-energy scenarios. . Several scenarios reach an elevation of 100 cm around year 2080, above which no further inputs of sediment or organic matter were allowed in the model. About 55% of all scenarios finish within the elevation zone of plant growth (50 cm below to 50 cm above MSL, which corresponds to MLLW to MHHW, respectively). The distribution of final elevation at year 2100 for (A) each input parameter where the inputs are ordered by rank and for (B) high and low energy scenarios for each modeled SLR. The number of model runs that define each distribution is plotted on the right axis. Each input is arranged by rank, but the magnitude of inputs is not reflected by the linear spacing along the horizontal axis. The elevation zone of plant growth is highlighted in green.

DISCUSSION
In this paper, we used 450 different scenarios to explore the future sustainability of tidal freshwater marshes in the Sacramento-San Joaquin Delta. For most of the scenarios, marsh surface elevations decreased over the modeled period (2000 to 2100). Despite the anticipated elevation loss, modeled marsh surfaces maintained a position within the elevation range of plant growth (MLLW-MHHW) for 65.5% of the runs (Figure 4). Scenarios for which marshes changed from high marsh to low or drowned marsh (marsh surface elevations < MLLW) increased as SLR increased. The two highest increases in SLR-133 cm and 169 cm-68% and 89%, respectively, terminated with marsh drowning by the end of the modeling period ( Figure 6). For these high-SLR scenarios, only those with the greatest inorganic sediment accumulation rate remained sustainable (Figure 7). In addition, model scenarios for high energy sites indicated that the high vertical accretion in these dynamic environments makes them marginally more resilient than their low-energy counterparts ( Figure 5).
More than 84% of the scenarios resulted in sustainable marshes with 88 cm of SLR by 2100, but only 32% and 11% of the scenarios resulted in sustainable marshes when SLR was increased beyond 88 cm. This strongly suggests that ~88 cm of SLR is a critical threshold for marsh sustainability for the given range of inorganic and organic accumulation parameters used in this modeling exercise. Overall, these results suggest that marshes at the upstream reaches of the Delta where SLR may be attenuated and high-energy marshes in major channels with large inorganic sediment accumulation rates will likely be more sustainable than their downstream and low-energy counterparts, particularly at SLR values of greater than 88 cm over the next century.
Currently, lack of ample marsh vertical accretion data precludes independent calibration and validation of this and other marsh sustainability models under use in the San Francisco Estuary (e.g., Deverel et al. 2008Deverel et al. , 2014Stralberg et al. 2011). Nevertheless, for the smallest final SLR projection (21 cm), which is nearly equal to SLR over the previous century Proportion and timing of marsh surface elevation changes for each SLR scenario. Each plot is for a suite of 90 scenarios for a given SLR. Drowned marsh has a surface elevation below MLLW, low marsh is from MLLW to MSL, and high marsh has a surface elevation above MSL.
than mean sea level, indicating that existing marshes are sustainable ( Figure 5A). This demonstrates that WARMER is reasonably reproducing current vertical accretion processes in the study area. Swanson et al. (2013) previously demonstrated good model performance of WARMER by comparing (1) soil profiles of percent organic matter and bulk density generated from 200 years of 2.1 mm yr -1 SLR to data from marshes in the San Francisco Estuary and (2) computed equilibrium marsh surface elevations and known elevations of marsh plant growth to existing marsh plant elevation ranges.
The sensitivity analysis of WARMER and calculation of control parameters showed that SLR is the primary driving force in marsh sustainability and also the parameter with the largest uncertainty ( Figure 5, Table 6). Like the Krone (1987) model, marsh surface elevation projections from WARMER reflect the shape of the SLR curve ( Figure 4). In addition, the results demonstrated that not only is the magnitude of SLR over the next century important, but the acceleration of SLR in the latter half of the century is also critical to marsh sustainability ( Figure 6). Simulations of marsh surface elevations resulted in a similar response to SLR over the next century as projected for salt marsh surface elevation changes in the San Francisco Estuary (Stralberg et al. 2011;Swanson et al. 2013). SLR and inorganic sediment accumulation appear to be dominant factors in determining tidal marsh sustainability in both freshwater ( Figure 5) and saltwater tidal marshes in the San Francisco Bay-Delta system (Reed 2002;Stralberg et al. 2011;Swanson et al. 2013;Schile et al. 2014). It is important to note, however, that uncertainty in the input parameters could be large and the corresponding uncertainty in marsh surface elevation increases with time in these scenarios. Sediment deposition was shown to be the second most important control parameter leading to resilient marshes under scenarios of future SLR (Table 6; Figure 7). Inorganic sediment is clearly an important part of vertical accretion in the Delta, both on millennial (Drexler et al. 2009a) and recent time scales (Reed 2002). Currently, the Delta traps approximately two-thirds of the sediment it receives from upstream (Wright and Schoellhamer 2005); however, sediment transport in the watershed appears to have waned in recent years because of the depletion of the pulse of gold mining sediments as well as the effects of dam-building (Wright and Schoellhamer 2004;Cloern et al. 2011;Schoellhamer et al. 2013). The magnitude and timing of changes in future sediment supply as well as the rate of deposition on the marsh plain are key uncertainties in projecting marsh sustainability. Because marsh sustainability is sensitive to an uncertain future sediment supply, our strategy was to bracket the likely future supply by simulating likely upper (no further decrease) and lower (constant rate of decrease) bounds.
Organic matter accumulation was expected to be a primary control on marsh surface elevation and sustainability in Delta marshes because of the history of organic soil (peat) formation in the Delta (Drexler et al. 2009a;Drexler 2011) and the dynamic feedback between organic matter productivity and elevation loss for high tidal marsh (Morris et al. 2002;Kirwan and Guntenspergen 2012). In addition, recent work by Schile et al. (2014)

Figure 7
Marsh habitat defined by the median elevation of the final marsh surface for 18 scenarios for each paired SLR and sediment scenario. High marsh is defined as marsh elevation above MSL and low marsh is from MLLW to MSL. Low marsh is more vulnerable to drowning because of a positive feedback between elevation loss and reduced organic matter accumulation. Drowned marsh scenarios terminate with a relative elevation below MLLW.
contain highly productive vegetation were projected to be more resilient than more saline sites with lower productivity. This finding suggests that organic matter accumulation is crucial for the long-term sustainability of low salinity marshes. However, the range of organic matter accumulation rates calculated for the Delta marshes did not exceed inorganic sediment deposition and SLR in controlling marsh surface elevation in the suite of model scenarios. It is possible that we may have underestimated organic matter accumulation by using measured organic matter accumulation rates in the surface layer (0 to 2 cm) of peat cores to approximate annual primary production. However, a potentially greater problem in applying parameters for this input is a major lack of data regarding the full potential productivity of Delta marshes under conditions of climate change and increased flooding. Marshes have been shown to have increased belowground plant productivity under higher atmospheric CO 2 concentrations (Langley et al. 2009). Furthermore, belowground plant productivity may increase at a higher rate than aboveground productivity in response to increased inundation (Kirwan and Guntenspergen 2012;Nyman et al. 2006). In addition, changes in salinity regime, increases in temperature, and nutrient loading may also affect plant productivity, yet these factors are currently not included in most marsh sustainability models (Schile et al. 2014).
Alternatively, a possible explanation for why organic matter accumulation did not rank as high as SLR and inorganic sediment accumulation in the sensitivity analysis of input parameters may be the parameterization of decomposition in WARMER. The decomposition function may oversimplify the belowground biogeochemical processes that favor organic matter accumulation during periods of low sediment supply and/or may simply overestimate the decomposition of refractory organic material. Working in restored, impounded Delta marshes, Miller and Fujii (2010) recorded carbon sequestration rates between 1,200 to 3000 g C m -2 yr -1 . This suggests that the upper limit for organic matter accumulation (which can be estimated as roughly double the rate of carbon sequestration or 2,400 to 6,000 g C m -2 yr -1 , Mitsch and Gosselink [2000, p. 157]) is likely much higher than what we used in the application of WARMER (maxi-mum 200 to 1700 g C m -2 yr -1 , Table 2) and may only be known when marshes are subject to increased inundation. Clearly, to better constrain this and other marsh sustainability models, more data are needed on the feedback between flooding and annual primary productivity (both above ground and below ground) as well as on long-term decomposition rates for a range of marsh species.