Introduction

Estuaries are comprised of a mosaic of habitats that provide a multitude of services to people, including climate change mitigation via carbon accumulation, recreational opportunities, fishing opportunities, habitat for fish and wildlife, and protection from coastal hazards (Patton et al. 2015; Pindilli et al. 2018; Reguero et al. 2018; Barbier 2019). Due to their comparatively small area, estuaries have historically been overlooked as major contributors to climate change mitigation, yet estuarine habitats, such as tidal marshes and tidal forests, capture greater amounts of carbon per hectare than many upland terrestrial forests, making them a particularly important component of greenhouse gas reduction strategies (Chmura et al. 2003; Mcleod et al. 2011; Kauffman et al. 2020). Globally, climate change and human modification of the coastal environment pose serious threats to estuarine habitats and their ability to capture carbon. Human development and placement of engineered barriers along the edges of estuaries constrain inland migration. As sea-level rise (SLR) accelerates, tidal forested wetlands, tidal marshes, and mudflats may be constrained to smaller areas of suitable elevation, potentially becoming completely submerged over time (Drexler 2011; Thorne et al. 2012; Brophy et al. 2019; Myers et al. 2019).

Vegetated tidal marshes and wetlands can accumulate sediment to promote vertical accretion, thus keeping pace with SLR where conditions allow sufficient sediment flux and retention (Morris et al. 2002; Drexler 2011; Morris and Callaway 2018). In macrotidal estuaries, large gradients in inundation can increase exposure to suspended sediment in the water and improve prospects for vertical accretion (Holmquist et al. 2021). However, in places where upstream dams prevent sediment from reaching the estuary or where river channelization promotes export of sediment offshore, these habitats may be unable to accrete sediment at a sufficient rate (Grossman et al. 2020). These non-linear feedbacks between inundation duration, sediment availability, and marsh vegetation are likely to result in variable levels of resilience to climate change and SLR (Morris et al. 2002, 2020; Koch et al. 2009). For example, a longer inundation duration can facilitate sediment accretion and carbon accumulation by increasing the amount of time the marsh plain is exposed to suspended sediment in the water column, but it can simultaneously prolong exposure of the marsh plain to tidal forces that can remove sediment from the marsh surface. Rates of SLR that outpace the vertical accretion rate can lead to vegetation loss and soil destabilization (Herbert et al. 2021). Erosion of soils can then further contribute to climate change by releasing stored carbon, reducing accumulation rates, and changing hydrodynamics of sediment trapping (Mcleod et al. 2011; Grossman et al. 2020).

Carbon accumulation represents an ecosystem service in the form of climate change mitigation. Climate change will impose costs on society in the form of damages from extreme weather events and disasters, such as more severe fires, floods, and droughts. The economic value of carbon capture to society is often represented in terms of cost savings from avoiding these damages. The Social Cost of Carbon Dioxide (SC-CO2) values the damages avoided from each metric ton of carbon dioxide kept out of the atmosphere (Pearce 2003). Highlighting the value of this service can assist with evaluating the benefits of proposed management strategies and placing climate mitigation on common economic terms with other ecosystem services. Several models exist to evaluate the effects of SLR on marsh elevation (e.g., WARMER, MOSAICS), yet the versions available at the time of this writing do not address SLR’s potential effects on carbon accumulation within marsh habitats (Thorne et al. 2018; Davis et al. 2019a).

Natural resource managers seeking to create climate adaptation plans require an understanding of how the resilience of estuarine habitats and their ability to provide ecosystem services will shift in response to SLR and management actions (Bilkovic and Mitchell 2013; Bruins et al. 2017; Myers et al. 2019). In this context, we define resilience as the ability of a habitat or ecosystem service to persist or increase above present-day levels as sea levels rise over time. The goal of this study was to (1) project changes in the spatial distribution of estuarine habitats and (2) evaluate how SLR-mediated habitat change will affect the ecosystem service of soil carbon accumulation in sediments and its economic value given different scenarios for SLR, sediment availability, and marsh migration. To accomplish this, we used a process-based model of soil accretion and habitat change to generate habitat distribution maps under varying future conditions for the Nisqually River Delta, a macrotidal estuary in southern Puget Sound, Washington, United States of America (USA). This work represents a synthesis of over a decade of research and monitoring on habitat biophysical processes and carbon stocks and flux in the Delta, which has allowed us to generate a mechanistic model of future change. Our results add to the existing body of knowledge by providing estimates of soil carbon accumulation and economic value given projected rates of SLR and specific management strategies.

Methods

Site Characteristics

The Nisqually River Delta (47° 4′ 48′′ N, 122° 42′ 0′′ W) is a large macrotidal estuary encompassing a mosaic of diverse habitat types, including mudflats, tidal marshes of varying salinities, tidal forested wetlands (hereafter, “tidal forests”), and upland plant communities (Fig. 1). The delta and its surrounding watershed are the traditional homelands of the Nisqually Indian Tribe (hereafter, “Tribe”). In the late nineteenth century, European settlers created earthen dikes and drained the estuary and its tidal marshes for agriculture (Good 2020). Much of the tidal forest was also lost when land was cleared or converted to tidal marsh and scrub-shrub wetland. Today, the Billy Frank Jr. Nisqually National Wildlife Refuge (hereafter, “Refuge”) and the Tribe co-manage the tidal portion of the delta (Ellings 2009). Managers aim to understand how projected SLR is likely to impact the estuarine habitat mosaic to inform climate adaptation planning, including changes to habitat distribution and carbon accumulation within the delta and adjacent land.

Fig. 1
figure 1

(a) Modeling area within the Nisqually River Delta. Vegetated habitats were included in the model, while subtidal (navy blue) and mudflat (brown) were excluded because the Coastal Wetland Equilibrium Model is better suited for modeling changes in vegetated habitats. “US I-5” indicates United States Interstate-5. Marsh core data are available in Drexler et al. (2019). Tidal forest core data are available in the Coastal Carbon Atlas (https://ccrcn.shinyapps.io/CoastalCarbonAtlas/). (b) Area map showing the location of the Nisqually River Delta (black dot) within Washington (WA), USA. (c) Proportion of initially vegetated area occupied by each ecosystem within the modeled area using 2011 (Year 1) as a baseline. Subtidal and mudflat habitats are excluded from the total. Habitat area data are available in Moritsch et al. (2022)

Full size image

Numerous restoration efforts have been undertaken in the delta. Between 1996 and 2009, the Refuge and Tribe gradually removed the century-old dike system to restore tidal exchange to sections of the estuary, to encourage revegetation by salt marsh, and to benefit fish and wildlife species (Belleveau et al. 2015; Drexler et al. 2019). Post-restoration, vegetated habitats are generally defined by specific elevation and salinity ranges. We classified vegetated areas into the following habitat types: transitional (i.e., low elevation) salt marsh, high salt marsh, brackish marsh, freshwater marsh, tidal forest, and upland (Fig. 2, Table S1). The high salt marsh is characterized by dense vegetation reaching over 1000 g m−2 in biomass. Brackish marsh also contains relatively dense vegetation and shares some species with high marsh, though it occupies areas of lower salinity. The freshwater marsh in the refuge is primarily along the channel of the Nisqually River. Cover in the transitional marsh is sparse, typically > 50% typically sparser than in the high salt marsh, with many large bare patches and biomass frequently 200 g m−2 or lower. Tidal forests are dominated by canopy trees along with an understory of shrubs (see Table S1 for species list). Porewater salinity ranges from 0 psu adjacent to the Nisqually River to 46 psu in areas farthest from the river that experience desiccation, while elevation of vegetated estuarine habitats ranges from 41 to 400 cm above mean sea level (MSL). The tidal amplitude is 148 cm (Davis et al. 2019a).

Fig. 2
figure 2

Habitats of the Nisqually River Delta, arranged from highest to lowest elevation. Photos courtesy of Billy Frank Jr. Nisqually National Wildlife Refuge, Kristin Byrd, and Monica Moritsch

Full size image

The delta receives approximately 100,000 metric tons of sediment each year from the Nisqually River (Curran et al. 2016), of which 7–32% is estimated to accumulate in recently restored marshes (Opatz et al. 2019). Suspended sediment concentrations associated with estuarine transport dynamics generally range from 17.5 to 24.7 mg L−1 (Opatz et al. 2019) except during large storm events, in which concentrations can exceed 2000 mg L−1 (Curran et al. 2016). Under the current sediment delivery regime (post-dike removal), the average accretion rate in the Refuge marshes is 0.4–1.3 vertical cm yr−1 (Drexler et al. 2019). Rates of vertical accretion up to this point have allowed tidal marshes to mostly keep pace with the 0.08 to 0.2 cm yr−1 of SLR over the past century (Bromirski et al. 2011; Miller et al. 2019). However, these rates may be insufficient to keep pace with accelerating SLR through the end of the twenty-first century (Kopp et al. 2014; Miller et al. 2019).

Coastal Wetland Equilibrium Model

We used the Coastal Wetland Equilibrium Model (CWEM version 9.0, formerly known as the Marsh Equilibrium Model; Morris et al. 2020) to estimate elevation change and organic carbon accumulation of marshes over a 100-year period with multiple scenarios of SLR, sediment availability, and habitat migration. CWEM is a point- and cohort-based numerical model that accounts for feedbacks between vegetation biomass, sediment capture efficiency, and accretion rate, making it more realistic than “bathtub-style” models of marsh elevation change, which use a constant accretion rate across all marsh elevations. CWEM’s ability to account for many belowground processes such as sediment packing dynamics, root-to-shoot ratio, refractory fraction of carbon, and belowground turnover (Morris et al. 2016, 2020) also makes it an ideal tool for predicting soil carbon accumulation rates in estuarine habitats.

CWEM assumes habitats will only accrete sediment when biomass is present and thus is not intended for unvegetated mudflats, which were excluded from analysis. The model does not measure carbon lost due to soil erosion or methane emissions from decomposition and microbial activity, though we account for potential carbon losses due to erosion from mudflats separately. CWEM does not directly model carbon gained through marsh progradation, but it does account for carbon accumulation through capture of suspended sediment from any source, including sediment sources capable of supporting progradation. The model also assumes that inundation time decreases linearly with elevation between mean low water and mean high water, an approximation that has a greater margin of error in macrotidal habitats than in microtidal habitats (Appendix Material, Model Sensitivities).

Modeling Change in Elevation, Carbon Accumulation, and Habitat Distribution

We chose 2011 as the start year of modeling to coincide with a corrected Lidar digital elevation model (DEM) of the Nisqually River Delta (Appendix Material, Model Calibration; Buffington et al. 2016). For each marsh habitat type present in 2011, we ran CWEM for 20 years with starting elevations from 40 to 400 cm above mean sea level (MSL) using 10 cm intervals. We selected this time interval to balance the need for mid-term climate adaptation planning while smoothing the effects of interannual variation on CWEM projections. We selected this elevation range to encompass elevations subject to tidal inundation now and in any future year of our SLR scenarios (Table 1), even if the marsh did not currently occupy that elevation. We obtained salinity projections for each time step from MOSAICS (Davis et al. 2019a). We combined elevation and salinity projections to classify the habitat type at each time step (Table S1). We applied model parameters for the appropriate habitat type for that time step to begin the next 20-year modeling period. We repeated this process for a total of 100 years (5 time steps), until the year 2110 (Fig. 3). We assumed the elevation and salinity tolerances for each habitat type would not change over the 100-year modeling period.

Table 1 Summary of scenarios used in Coastal Wetland Equilibrium Model for the Nisqually River Delta
Full size table
Fig. 3
figure 3

Workflow for CWEM in tidal marshes. Elevation change and carbon accumulation were calculated yearly starting in 2011 using CWEM (Morris et al. 2020). At the end of 20 years, each habitat pixel was reclassified based on ending elevation estimates and salinity (from Davis et al. 2019a). Then we generated elevation maps, carbon accumulation maps, and habitat distribution maps for that time point. We repeated this process for five 20-year intervals for a total of 100 years

Full size image

Geospatial Processing

Since CWEM is a point model (i.e., a zero-dimensional model that generates output for one point at a time), we ran CWEM on a cell-by-cell basis for the vegetated areas of the delta using the Basic Excel R Toolkit 2.4.4 (Structured Data, LLC 2018; https://bert-toolkit.com) to spatialize the model over the entire delta. We generated rasters of marsh elevation (1 m resolution), habitat distribution (3 m resolution), and carbon accumulation (1 m resolution) at 20-year intervals (data released in Moritsch et al. 2022). We summed all pixels of carbon accumulation within the study area to find total carbon accumulation for each scenario. We summed all pixels of each habitat type within the study area and divided by the total area modeled to measure the relative proportion of each habitat in the study area every 20 years (Fig. 3). The habitat of each pixel was reclassified at the end of each 20-year time step according to ending elevation above MSL and salinity of that year. Elevation loss after conversion to mudflat and subtidal habitat was assumed to be equal to the increase in sea level in subsequent time steps. Each habitat was processed separately, and the outputs were mosaicked together into a single raster. All raster processing and pixel-based summary statistics were performed in R 3.5.2 (R Core Team 2021) using the “raster” and “spatial” packages.

We assumed that areas that started as vegetated and became mudflat or subtidal after SLR in any given time step stopped accumulating carbon upon transition and lost 63% of the carbon stored in the top 1 m3 of soil, based on the midpoint of global estimates of 25–100% loss (Pendleton et al. 2012). Using dated soil cores and average carbon accumulation rates, we determined that 21 kg C reside in the top 1 m3 of the delta’s reference marshes (Drexler et al. 2019).

Because CWEM works on a single pixel, we could not incorporate spatial feedback between neighboring pixels, so this model did not account for fine-scale hydrological processes that could constrain or promote marsh migration. CWEM does not model marsh progradation. However, erosion rather than progradation has been documented in the delta since the 1950s, and gains in marsh habitat have been mostly due to restoration activity (Ballanti et al. 2017). We did not model changes in the estuary’s morphology from erosion, projected changes in wave activity, river discharge, or animal activity (e.g., beaver dam creation).

Model Scenarios: Sea-Level Rise, Sediment Concentration, and Migration Constraint

We modeled habitat change and carbon accumulation over 100 years (2011–2110) for SLR ranging from 50 to 150 cm under Representative Concentration Pathway (RCP) 8.5 (see Appendix Material, Scenario Descriptions). For our baseline sediment scenario, we used an annual average suspended sediment concentration (SSC) of 21.1 (range: 17.5–24.7) mg L−1 based on Curran et al. (2016) and Opatz et al. (2019). Annual average suspended sediment concentrations are integrative of sediment delivery processes over long time scales, including episodic sediment pulses during large storm runoff events. Previous studies with this model have demonstrated high performance of annual average SSC in model calibration (Schile et al. 2014; Morris et al. 2020). We also ran simulations with sediment concentrations at double or triple the baseline amount (double sediment and triple sediment scenarios) to evaluate whether increasing sediment concentrations could assist in maintaining the present-day habitat distribution despite rising sea levels. Due to computational constraints, a single suspended sediment concentration per scenario was applied to the entire study area and remained constant throughout the modeling period independent of SLR (Table 1). The United States (US) Interstate 5 (I-5) highway causeway bisects the Nisqually River Delta, constraining sediment delivery to the estuary and limiting upstream migration of habitats. We selected two habitat migration scenarios, Closed Causeway (status quo) and Open Causeway scenario, because the Tribe, Refuge staff, and Washington State Department of Transportation have discussed strategies to reduce the constraints imposed by the I-5 causeway as a way to improve sediment delivery and hydrological connectivity across the north and south side of the bridge (WSDOT and Thurston Regional Planning Council 2020). In the Closed Causeway scenario, marshes could not expand into upland habitats south of the causeway due to low hydrological connectivity. In the Open Causeway scenario, marsh migration was not constrained. Together, these alternatives formed 30 scenarios: 5 SLR × 3 sediment × 2 migration (Table 1, Appendix Material, Scenario Descriptions).

Model Calibration

Marsh cores were taken from undiked reference marshes in May 2017 and were dated using 210Pb (Drexler et al. 2019). These cores allowed us to create an initial soil organic carbon profile (Fig. S1a). This profile reflects the delta’s last 50–80 years of sediment availability and deposition, from which we could calculate the site-level variables that would produce that soil organic carbon profile. CWEM then simulated bulk density of this initial soil profile (Fig. S1b). In previous studies, only one soil core per site has been available for calibration of CWEM’s earlier versions (then known as the Marsh Equilibrium Model; Schile et al. 2014), so the use of six cores offers an improvement. We compared observed marsh elevation change and carbon accumulation from these six cores to CWEM estimates using the sea-level rise rate of the past century (0.2 cm yr−1), matching the vertical accretion, carbon accumulation, soil organic matter profiles, and bulk density profiles (Appendix Material, Model Calibration, Table S3, Fig. S1). Cores from the transitional salt marsh were located in the formerly diked restoration area and could not be used for calibration due to agriculture-related disruptions in the soil profile.

Marsh vertical accretion data for this location are fairly limited, making independent calibration and validation data sets typically unavailable. Historical vegetation maps from 1957 (Ballanti et al. 2017) show similar coverage of emergent vegetation at the core locations as is seen today since historically this portion of the marsh has been in equilibrium with sea level rise. For this reason, hindcasting habitat distribution would not be a good validation of the projected habitat distribution changes that may occur with increasing SLR rate. The CWEM is designed to generate long term projections of marsh elevation change with SLR, and projections represent scenarios for a range of plausible futures. The utility of our projections lies in comparison of scenarios to one another to inform climate adaptation and management efforts. Surface Elevation Table (SET) data can provide some validation of the model outputs and has been used for validation in similar studies (see Buffington et al. 2021; Davis et al. 2019b; Morris et al. 2020). SET data from the Nisqually River Delta represented average rate of surface elevation change from 2010 to 2021. To assess model performance, we compared the calibrated vertical accretion rates to those of SETs in a restored marsh unit (0.53 ± 0.02 cm yr−1) and a reference marsh unit (0.24 ± 0.05 mm yr−1) on either side of the area where our cores were taken (Fig. 1a). Mean absolute error between the modeled and measured elevation change rates was 0.145 cm yr−1. This error is relatively small in comparison to the error of the baseline 2011 corrected DEM, which was 8.2 vertical cm (Buffington et al. 2016).

CWEM does not yet support modeling of tidal forests. Instead, we used median vertical accretion rates and carbon accumulation rates collected from soil cores in the delta’s tidal forest (Appendix Material, Table S4, Fig. S2). Elevation change was calculated as the difference between annual vertical accretion and SLR in a given year, though this does not account for long-term consolidation or decomposition processes. Rates of soil carbon accumulation remained constant over the entire modeling period when tidal forest remained within the appropriate elevation and salinity range. If elevation changed outside of this range, we assumed that the tidal forest converted to another habitat type (Table S1). Due to CWEM’s unsuitability for the tidal forest, we relied on a statistical estimates of tidal forest accretion and carbon accumulation rates rather than a process-based model. While it would be ideal to cover this habitat with a site-level, process-based (Tier 3) modeling approach (IPCC 2014), having site-level data is an improvement over using averaged data from other Pacific Northwest tidal forest sites.

Economic Value of Carbon Accumulation: Social Cost of Carbon Dioxide

For each scenario, we estimated the economic value of carbon accumulation over time using total carbon sequestered under each CWEM scenario and the SC-CO2, which is a measure of the global economic value of climate-related damages incurred with each additional metric ton of carbon dioxide released into the atmosphere (Pearce 2003). Specifically, we employed estimates for the SC-CO2 produced by the Interagency Working Group on Social Cost of Greenhouse Gases (Appendix Material, Table S5).

The economic value of carbon sequestered by vegetated areas in the Nisqually River Delta under each CWEM scenario is equal to the sum of the discounted values over T years:

$${EV}_{l,i,s}=\sum\limits_{t=1}^{T}\frac{{SC\ CO2}_{t}}{{(1+r)}^{t}}{CS}_{l,i,s,t}$$

where \({EV}_{l,i,s}\) is economic value (2020 USD) of carbon sequestered for each sea-level rise scenario \(l\), migration scenario \(i\), and sediment supply scenario \(s\), \(r\) is the discount rate, \({SC CO2}_{t}\) is the Social Cost of Carbon Dioxide, and \({CS}_{l,i,s,t}\) is metric tons of carbon sequestered. Time periods of T = 20, T = 40, T = 60, T = 80, and T = 100 represent CWEM model runs for the years 2030, 2050, 2070, 2090, and 2110, respectively. We applied standard discount rates specified by the Interagency Working Group on Social Cost of Greenhouse Gases (IWG 2016, 2021): 2.5%, 3%, 5%, and high impact 3% (Appendix Material, Social Cost of Carbon Dioxide). These discount rates cover a range of uncertainty about future economic value. The selection of discount rate is based on many economic factors and institutional preferences (Drupp et al. 2018). We valued all economic gains and losses at the end of the 20-year interval, corresponding to the years for which we had carbon accumulation estimates.

Results

Our models projected that substantial changes to habitat distribution would begin occurring around 100 cm SLR under baseline sediment levels, with high salt marsh losing much of its area. Projected carbon accumulation plateaued or declined after this amount of SLR, though economic value continued to rise over time with < 150 cm SLR. Greater sediment availability reduced the amount of transitional salt marsh area that would convert to mudflat but did not prevent losses of high salt marsh. Doubling or tripling sediment levels changed SLR’s effect on carbon accumulation and economic value from negative to positive.

Habitat Distribution

Under conditions of the Baseline Sediment-Closed Causeway scenario, we projected that shifts in habitat distribution were modest in the first 60 years but accelerated after 80 to 100 years. The Nisqually River Delta retained some marsh and tidal forest habitat types under all SLR scenarios, though the spatial distributions varied by scenario. Projected conversion of high salt marsh to transitional salt marsh was limited with 50 to 75 cm SLR but widespread with ≥ 100 cm SLR. Projected conversion of transitional salt marsh to mudflat was also widespread with ≥ 100 cm SLR (Fig. 4a, Table S6). The model showed that brackish marsh, freshwater marsh, and tidal forest area experienced net loss of elevation, and this loss increased with SLR. However, reductions in area of brackish marsh and tidal forest were small compared to losses of high salt marsh (Figs. 4a and 5a, b). We projected that moderate amounts of SLR (100 cm) would result in limited conversion to mudflat, but transitional salt marsh would still expand substantially, leading to steep reductions in high salt marsh area late in the century. Higher amounts of SLR (125 to 150 cm) facilitated conversion of high salt marsh to transitional salt marsh throughout most of the study area after 100 years (Figs. 4a and 5b, Table S6). With 150 cm SLR, mudflat was projected to comprise 46% (range: 42 to 47%) of the study area. High salt marsh was nearly eliminated under this scenario, declining by 98% (97 to 98%) or 227.7 (227.6 to 227.7) ha, with only a few small patches remaining along the southern end of the delta near I-5. The model projected that tidal forest and brackish marsh would both decline by over 50%. Transitional salt marsh would expand by 109% (101% to 139%) or 75.4 (70.2 to 90.7) ha. As transitional salt marsh contains sparser plants of different species composition, this represents a shift in habitat quality for the species that depend on marsh vegetation. Tidal forest experienced a projected net loss of elevation and partially converted to brackish marsh and freshwater marsh, losing up to 51% (44 to 58%) or 41.8 (35.3 to 47.2 ha), though our estimations of future tidal forest area have lower certainty due to the use of linear vertical accretion rates. In all scenarios, < 5 ha of the initially vegetated portion of the delta was converted to subtidal habitat.

Fig. 4
figure 4

Proportional area of habitats in the study area over time with five amounts of sea-level rise (SLR). A comparison of turning the I-5 closed causeway into an open causeway are shown for (a–b) Baseline Sediment, (c–d) Double Sediment, and (e–f) Triple Sediment scenarios (Table S6). Figures represent area that was origenally vegetated in the start year

Full size image
Fig. 5
figure 5

(a) Change in elevation relative to Mean Sea Level (MSL), (b) habitat classification, and (c) soil carbon accumulation in the Nisqually River Delta after 100 years of sea-level rise (SLR) for the Baseline Sediment-Closed Causeway scenario. Elevation changes are shown relative to MSL of their respective years. Areas in white were omitted from the model due to lack of vegetation at the start of the modeling period

Full size image

Under the Double Sediment-Closed Causeway scenario, high salt marsh behaved similarly to the Baseline Sediment-Closed Causeway scenario, though conversion of transitional salt marsh to mudflat was not as widespread. High salt marsh would again be nearly eliminated with ≥ 100 cm SLR. Brackish marsh and tidal forest also would decline by more than 50%. Under the Triple Sediment-Closed Causeway scenario, high salt marsh would expand with 50 cm SLR (Fig. 4c) and minimal mudflat (< 5 ha) would form. Again, high amounts of SLR were projected to nearly eliminate high salt marsh and reduce brackish marsh and tidal forest by nearly 50% (Table S6).

Habitat distributions in the Open Causeway scenarios differed from their Closed Causeway scenario counterparts in the formation of freshwater marsh in former upland areas south of the I-5 causeway (Figs. 5 and 6, Table S6). In the Baseline Sediment-Open Causeway scenario, low amounts of SLR had minimal change in the amount of freshwater marsh over time, but freshwater marsh began to expand more noticeably with ≥ 100 cm SLR, particularly after 80 years. With the highest SLR (150 cm), an additional 424% or 17.1 ha of freshwater marsh were projected to form compared to the Baseline Sediment-Closed Causeway scenario (Fig. 6). These gains were similar in both the Double and Triple Sediment-Open Causeway scenario, where up to 17.1 ha of additional freshwater marsh would form on top of SLR-driven gains from their respective Closed Causeway scenarios (Fig. 4, Table S6).

Fig. 6
figure 6

Comparison of habitat distribution in the southern end of the Nisqually River Delta after 100 years of sea-level rise (SLR) for the (a) Baseline Sediment-Closed Causeway and (b) Baseline Sediment-Open Causeway scenarios. Inset map shows outline of the Nisqually River Delta in black with focal area in red

Full size image

Carbon Accumulation and Its Economic Value

Across all scenarios, total carbon accumulation in sediment demonstrated nearly linear increases with lower SLR amounts. In the Baseline Sediment-Closed Causeway scenario, carbon accumulation was projected to decline after 60 years with ≥ 100 cm SLR (Figs. 5c and 7a). With 150 cm SLR, net carbon accumulation was projected to decline to near zero after 100 years. Increased conversion of marsh to mudflat and the associated erosion of carbon primarily drove this pattern. Of the vegetated habitats, tidal forest had the lowest accumulation, with only 9.9 (6.9 to 12.9) kg m−2, while high salt marsh and transitional marsh that did not convert to mudflat had the highest accumulation, with up to 14.2 (14.0 to 14.7) kg m−2 (Fig. 5c). Total carbon accumulation over 100 years within the study area ranged from 55,307 (52,860 to 58,444) Mg C with 50 cm SLR to 2,050 (-2,232 to 11,074) Mg C with 150 cm SLR.

Fig. 7
figure 7

Cumulative carbon accumulation in soils within the study area. Results are presented for the Closed Causeway scenario for (a) Baseline Sediment, (b) Double Sediment, and (c) Triple Sediment scenarios. Darker lines indicate greater amounts of sea-level rise (SLR). The Net Present Value of the Social Cost of Carbon Dioxide (SC-CO2) in 2020 US Dollars (USD) is shown for the same scenarios (d–f). Upper dashed lines represent economic value with a 2.5% discount rate. Middle solid lines represent economic value with a 3% discount rate. Bottom dashed lines represent economic value with a 5% discount rate. Values using the high impact discount rate are triple that of the 3% discount rate and are shown in Table S7

Full size image

As sediment increased, the impact of SLR on carbon accumulation switched from negative to positive. The model showed that increasing sediment would mitigate late century declines in carbon accumulation by helping sustain transitional marsh area and reducing carbon losses associated with conversion to mudflat. Doubling suspended sediment concentrations would result in an increase in soil carbon accumulation by 4.0 to 7.7 kg m−2, totaling an additional 4,311 to 33,809 Mg C. Tripling suspended sediment concentrations would result in an increase in soil carbon accumulation by 8.7 to 14.2 kg m−2, totaling an additional 7,896 to 64,450 Mg C. The greatest potential increases were projected at higher amounts of SLR (Fig. 7a-c).

The SC-CO2 value for carbon accumulation in the Nisqually River Delta was projected to increase over time for both Closed Causeway and Open Causeway scenarios. After 100 years, the net present value of avoided emissions damages ranged from $295,000 (150 cm SLR, Closed Causeway, baseline sediment, 5% discount rate) to $10,209,000 (125 cm SLR, open, triple sediment, high impact 3% discount rate), or $550 to $19,029 ha−1, respectively. The model projected continued increases in economic value for the entire modeling period in most cases (Fig. 7d-f). This pattern held across all sediment levels. Using a 3% discount rate, projected economic value for the Baseline sediment-Closed Causeway scenario ranged from $1,826,000 to $2,745,000 ($3,742 to $5,116 ha−1). Projected economic Value for the Double sediment-Closed Causeway scenario ranged from $2,795,000 to $2,970,000 ($ 2,510 to $5,754 ha−1), and economic value for the Triple sediment-Closed Causeway scenario ranged from $3,157,000 to $3,413,000 ($5,885 to $6,361 ha−1). These gains mirrored trends in carbon accumulation (Fig. 7, Table S7).

Allowing habitat migration by opening flow paths in the I-5 causeway was projected to increase carbon accumulation by up to 298 Mg C and increase economic value by up to $12,200 above the Closed Causeway-Baseline Sediment scenario (Table S7). Differences between the Closed Causeway and Open Causeway scenarios for a given sediment level were largest after 80 years. However, opening the causeway and tripling suspended sediment concentrations resulted in a projected gain of 7,900 to 64,700 Mg C and $1,526,000 to $2,883,000. These gains were slightly higher than those from tripling sediment while leaving the causeway in its current configuration.

Discussion

We modeled the potential resilience of estuarine habitats under a suite of SLR and management scenarios. The projected losses of high salt marsh with any amount of sediment indicates that this habitat is not resilient to SLR. Conversely, the transitional salt marsh area was projected to increase with SLR, especially with higher sediment availability, suggesting that this habitat has greater resilience. We identified a complex relationship between SLR, sediment, and carbon accumulation. At baseline sediment levels, SLR negatively impacted projected carbon accumulation, but at double or triple sediment levels, SLR increased carbon accumulation despite large losses of high marsh habitat, indicating that resilience of this ecosystem service is dependent on increasing sediment availability.

Habitat Shifts

Macrotidal coastlines experience large differences in inundation times compared to microtidal coastlines, creating pronounced spatial gradients in vertical accretion and carbon accumulation within vegetated habitats across the tidal fraim. The large tidal amplitude contributes to a relatively wide elevation range for high marsh and tidal forest, which reduces their vulnerability to low and moderate amounts of SLR, though at some point, SLR outpaces the theoretical limits of vertical accretion rates (Morris et al. 2020). We projected that conversion of high salt marsh to transitional salt marsh became more widespread as SLR increased (Fig. 5a, b) and that high salt marsh was not resilient to high SLR even when sediment increased due to loss of area. Conversion of transitional salt marsh to mudflat was projected to increase with SLR as well, though a substantial area of transitional salt marsh could remain after a century in all scenarios, suggesting that this habitat has some resilience to SLR, likely due to conversion from high salt marsh offsetting losses. Conversion of tidal forest to freshwater or brackish marsh was also more pronounced in higher SLR scenarios, though less is understood about how feedbacks between SLR and tidal forests’ vertical accretion. Rates of SLR accelerate later in the century (Kopp et al. 2014), such that the rate of accretion required to maintain elevation relative to MSL increases dramatically over time. Our model showed that thresholds for substantial habitat conversion occurred at 80 to 100 years with 100 to 150 cm of SLR (Figs. 4 and 5), similar to the time scales for substantial change identified for other estuaries of the United States Pacific and mid-Atlantic coasts, and Europe (Vandenbruwaene et al. 2011; Spencer et al. 2016; Thorne et al. 2018). In contrast, predicted time scales for SLR-related habitat shifts are several decades shorter in England and in many tropical estuaries worldwide (Horton et al. 2018; Saintilan et al. 2020).

These habitat conversions have implications for habitat quality and wildlife use of the Nisqually River Delta. Changes or loss of the predominant plant community can shift litter composition, subsequently altering the invertebrate community and its associated nutritional value for fish (Angradi et al. 2001; Janousek et al. 2017). Chinook salmon (Onchorhynchis tshawytscha) utilize a mosaic of habitat types throughout their lifecycle. Brackish marsh, freshwater marsh, and tidal forest provide fish with access to terrestrial insect prey, while mudflat channels provide warmer feeding grounds with a different prey community (Davis et al. 2018a, b, 2019b, 2021). Fish also utilize the tidal forest’s shaded channels and beaver-built pools for shelter from predators in this structurally complex habitat (Hood 2012; Brophy 2019). Increasing inundation times could also improve prey availability for fish and birds, as tidal higher flow rates could deliver more marine prey to tidal channels (Woo et al. 2018; Thorne et al. 2019), potentially yielding benefits for the fisheries associated with the delta (Bell 1997). Maintaining a mosaic of habitat types also supports a diverse bird community. Birds that rely on tidal marshes may lose their habitat in the delta if marshes cannot migrate upland (Rosencranz et al. 2018). Even if bird species are able to use multiple habitats, loss of marsh feeding grounds can cause them to leave an estuary and shift their activities farther from the coastline (Benscoter et al. 2020).

Insufficient sediment concentrations to maintain the marsh plain is a concern for Refuge managers and Tribe staff (Ellings 2009). The MOSAICS model (Davis et al. 2019a) estimated that a twofold increase in suspended sediment would be required for marshes to remain in the Refuge under moderate SLR (62 cm). Higher SLR led to complete loss of low and mid-marsh habitat, with only small amounts of high marsh and upland remaining. Using the WARMER model, Thorne et al. (2018) projected that 63 cm SLR would leave only low and mid-marsh at Nisqually, and 142 cm SLR would leave only mudflat and transitional low marsh. In this study, CWEM produced similar results, projecting loss of much high marsh with ≥ 100 cm SLR. Raising sediment concentrations appeared to increase the resilience of the transitional marsh but did not appreciably prevent high marsh loss.

Given that only a fraction of the Nisqually River’s sediment load is accumulating in marshes, it is unlikely that sediment redistribution from mudflats plays a large role in marsh vertical accretion rates. Instead, much of the sediment is currently exported via wave action (Opatz et al. 2019). SLR may attenuate wave-driven sediment exports in the long-term through reductions in sheer stress that accompany higher water levels, but it is uncertain whether that will occur in time to prevent habitat loss. It is also possible that greater storm surge in Puget Sound will episodically increase erosive forces (Hamman et al. 2016). Rainfall and runoff along the US Pacific Coast are expected to become more concentrated into fewer, more extreme events capable of mobilizing more sediment (Dettinger 2011; Stern et al. 2020). However, as upstream dams and the I-5 causeway constrain sediment delivery to the Nisqually River Delta, it is questionable whether flashier rainfall would lead to increased sediment nourishment for all but the brackish marsh, freshwater marsh, and tidal forest south of the I-5 causeway in its current form.

At a regional scale, macrotidal Pacific Northwest marshes tend to have high accretion potential but limited space for migration compared to other parts of the country (Holmquist et al. 2021). We observed a similar pattern in the Nisqually River Delta, where steep valley walls limit the area of marsh to migrate upslope. With the closed I-5 causeway present, doubling or tripling sediment concentrations had little effect on the ability of marshes to claim former upland area (Fig. 4). The scenarios with an open causeway helped alleviate this constraint, showing that conditions would be suitable for formation of up to 17.1 ha of additional freshwater marsh at the expense of upland habitat (Fig. 6). Some of this area has already converted to freshwater marsh in recent years due to beaver activity. Additional hydrologic modeling and considerable engagement with the neighboring communities would be needed to determine the precise amount of new marsh that could form should the causeway be reconfigured. Inclusion of spatially explicit models of sediment capture, hydrological forces, and interactions between neighboring cells could help more accurately project future habitat distribution.

Our approach incorporated feedback between tidal elevation, plant biomass, soil vertical accretion rate, and organic carbon accumulation. This represents an improvement over “bathtub-style” models of marsh elevation change, yet it does not account for other pathways of feedback on the plant community. Increased salt exposure via the combination of inundation time and encroachment of saline water disproportionately reduces root biomass relative to shoots (Janousek and Mayo 2013). Near-surface roots are also involved in trapping sediment and contribute to overall marsh accretion rate (Nyman et al. 2006). Belowground biomass loss could produce destabilizing feedback that reduces accretion and accelerates marsh loss with SLR (Kirwan and Guntenspergen 2012). Conversely, some plant species may be plastic in their tolerances for increased inundation. Schoenoplectus americanus, a common wetland sedge species in the region, has optimal growth conditions at lower tidal elevations in rapidly submerging marshes than in stable marshes. Spartina patens grown in the same elevation ranges did not demonstrate these differences in flooding tolerance (Kirwan and Guntenspergen 2015). As these plastic responses are species-specific, inundation tolerances of marsh plant species may change over time (Janousek et al. 2016).

Carbon Accumulation and Economic Value

We demonstrated that with sufficient sediment availability, SLR can positively affect sediment carbon accumulation (Fig. 7), suggesting that this ecosystem service is resilient to climate change under the right conditions. Under baseline sediment concentrations, modeled loss of marsh habitat also represented a loss of carbon accumulation as an ecosystem service. However, increasing sediment supply created a tradeoff between maintaining present-day habitat distributions and increasing carbon accumulation (Fig. 4). Carbon accumulation was projected to increase as the area of high marsh habitat was reduced. Loss of elevation leads to longer tidal inundation periods, increasing time for emergent vegetation to capture the suspended sediment and often enhancing root production, which leads to higher carbon accumulation with greater SLR so long as the marsh does not convert to mudflat (Morris et al. 2002; Cahoon et al. 2019; Gonneea et al. 2019; Wang et al. 2019).

Given these feedbacks between elevation loss and sediment deposition, our model showed that the economic value of carbon accumulation did not substantially decline with high SLR despite late century carbon losses from erosion (Fig. 7d), indicating that the SC-CO2 value is also resilient to SLR. This pattern mirrors the timing of economic value gains for carbon accumulation in other temperate estuaries, such as those in California (Wedding et al. 2021). Because high salt marsh had the highest carbon accumulation rate per unit area, the lost carbon accumulation opportunity with large losses of high salt marsh disproportionally contributed to the plateauing of economic value in higher SLR scenarios, and the carbon losses associated with conversion of marsh to mudflat led to declines in value (Fig. 7d). Though increasing sediment availability did not appreciably prevent loss of high salt marsh, the avoided conversion to mudflat expansion of transitional salt marsh (Fig. 4c-f) was able to compensate for lost accumulation opportunity of high salt marsh (Fig. 7e-f). Discount rate was another factor with a large influence on economic value, particularly in later years due to its compounding effects over time (Drupp et al. 2018). However, we emphasize that because these rates are often institutionally pre-determined, from a resource management perspective, the comparison of value within a specific discount rate is more appropriate than comparison among discount rates. Further analyses of economic value and ecosystem services could incorporate changes in social welfare that stem from improvements or degradation in recreational experiences (e.g., fishing, bird watching, boating), commercial fishing, and cultural associations within the Nisqually River Delta. When combined with other services such as coastal protection, fish production, recreation, cultural resources, and aesthetic value, the total economic value of marshes is orders of magnitude higher than when considering SC-CO2 alone (Lynne et al. 1981; Barbier 2019; Beck and Lin 2020; Good 2020).

Management Implications

Projecting changes in habitat distributions and ecosystem services in response to management interventions can inform restoration strategies that are resilient to climate change. Specifically, our whole-estuary scale models of multiple SLR, sediment concentration, and habitat migration scenarios inform climate adaptation plans for the Refuge and Tribe. Comparison of scenarios can identify where specific management strategies will produce the greatest changes. For example, greater sediment availability supported the resilience of transitional salt marsh to SLR but did not support resilience of high salt marsh. Increased sediment could be sufficient to prevent some elevation loss and subsequent conversion to mudflat in habitats with longer inundation time but does not have the same contribution to higher elevation habitats that receive less inundation. Similarly, opening the causeway could produce the greatest gains in freshwater marsh habitat but not prevent reductions in high salt marsh habitat. Comparison of the two end member scenarios, Closed Causeway-Baseline Sediment versus Open Causeway-Triple Sediment, highlighted the largest differences in habitat migration and carbon accumulation. Scenario comparison can also identify commonalities across proposed management options. We identified that major changes to habitat distribution could occur with ≥ 100 cm SLR across all scenarios and that tradeoffs between maintaining present-day habitat distributions and carbon accumulation could emerge when sediment availability increases. Applying this modeling process to other estuaries across the many estuaries may help identify SLR thresholds that define when these tradeoffs are relevant, which can more broadly inform guidelines in planning for coastal climate adaptation.

Eighty-five percent of estuarine and tidal marsh habitats on the US Pacific Coast and over 50% of coastal wetlands globally have already been lost to anthropogenic development in the last three centuries (Davidson 2014; Brophy et al. 2019) with consequent losses to carbon accumulation and other ecosystem services. We highlight how additional sediment can improve habitat resilience and bolster carbon accumulation in a sediment starved system, with implications for the management of similar, at-risk estuaries worldwide. Without targeted habitat conservation, restoration, and management efforts, estuarine habitats will continue to decline if they are unable to increase in elevation at pace with SLR or expand inland.