Localized Scenarios and Latitudinal Patterns of Vertical and Lateral Resilience of Tidal Marshes to Sea‐Level Rise in the Contiguous United States

Coastal wetlands have two dimensions of vulnerability to sea‐level rise (SLR), a vertical one, in cases where SLR outpaces their capacity to vertically accrete, and a lateral one, in cases where they are restricted from migrating inland by topography and land use. We conducted a meta‐analysis of accretion rates, standardized our analysis by using only 137Cs based estimates, and used model intercomparison to generate a vertical resilience index, a function of local SLR, tidal range, and tidal elevation category for the tidal wetlands of the contiguous US. We paired the vertical resilience index with a lateral resilience index made up of elevation, water level, and land cover maps, then projected them both into the future using localized SLR scenarios. At the regional scale, the vertical resilience index predicts changes from marsh aggradation to submergence for the coastal US Mid‐Atlantic, Southeast, and portions of the Northeast by 2100. At the sub‐regional scale, there is a geographic tradeoff between vertical and lateral resilience with more northerly wetlands vulnerable to the lack of suitable proportional area to migrate into, and more southerly wetlands vulnerable to accretion deficit. We estimate between 43% and 48% of the existing contiguous US wetland area, almost entirely located in watersheds along the Gulf of Mexico and Mid‐Atlantic coasts, is subject to both vertical and lateral limitations. These vertical and lateral resilience indices could help direct future research, planning, and mitigation efforts at a national scale, as well as supplement more processed informed approaches by local planners and practitioners.

geological conditions interact with eustatic SLR to produce differing rates of relative sea-level rise (RSLR)

Materials and Methods
In this study, we built then analyzed indices of tidal wetland vertical and lateral resilience to RSLR for CONUS. Then we examined the implication of geographic patterns in environmental drivers by calculating those indices using localized RSLR scenarios. First, we designated an area of interest for the study. Second, we designed the vertical resilience index by using a linear modeling framework with vertical accretion rate as the dependent variable. We built the calibration dataset for vertical accretion by performing a meta-analysis of dated soil cores, focusing in the end on 137 Cs-dated cores. We build the environmental covariate dataset by querying databases and parsing the original sources for additional information. Third, for the lateral resilience index, we performed a geographic information system (GIS) analysis summarizing the relative area and types of land available for tidal wetland transgression by combining publically available elevation, water-level, and land cover data. Fourth, we compare both indices run into 2100 using RSLR scenarios by Kopp et al. (2014).

Defining an Area of Interest
To define our area of interest we mapped all CONUS coastal watersheds containing tidal wetlands based on the National Wetlands Inventory (NWI) (U.S. Fish and Wildlife Service, 2014). We downloaded NWI data for each oceanic coastal state including Washington, Oregon, California, Texas, Louisiana, Mississippi, Alabama, Florida, Georgia, South Carolina, North Carolina, Virginia, Maryland, Delaware, New Jersey, Pennsylvania, New York, Rhode Island, Connecticut, Massachusetts, New Hampshire, and Maine. We included any NWI polygons with tidal modifier "E" or "M," indicating estuarine or marine systems, or with both a "P," indicating palustrine system, and "R," "S," "T," or "V," which indicate palustrine tidal systems. For watershed units, we used intermediate watershed boundaries (Hydrologic Unit Code Level 8 [HUC8]; United States Department of Agriculture-Natural Resources Conservation Service (USDA-NRCS), the United States Geological Survey (USGS), and the Environmental Protection Agency (EPA), 2015). We selected HUC8 units that overlapped NWI tidal wetlands and manually removed watershed units that were selected because of obvious outliers caused by mistakes in NWI coding.

Vertical Resilience Index Formulation
We calculated vertical resilience (VR) as modeled net accretion (Equation 1), RSLR subtracted from modeled accretion. Values greater than 0 indicate aggradation, while values less than 0 indicate submergence.

 VR Modeled Accretion -RSLR
(1) We used a series of linear models to determine which environmental variables to include in the vertical resilience index. We also investigated the effect of time itself on accretion since longer periods of time could potentially decrease observed accretion because of sediment compaction and belowground carbon loss. RSLR, greater diurnal tidal range (GT), watershed level flow weighted average suspended sediment concentration (FWA-SSC), inundation class, and time were all treated as independent variables and accretion rate as the dependent variable. RSLR, GT, FWA-SSC, and time were all treated as an additive, with inundation status interacting with them. We used linear modeling functionality in base R (R Development Core Team, 2014). In a series of iterative steps we eliminated the least significant parameter until all parameters were significant (p < 0.05). We used Akaike's Information Criterion for small sample sizes (AICc) to intercompare models, quantifying tradeoffs between performance and parsimony (Barton, 2015). Since frequently inundated marshes respond to SLR more dynamically than infrequently inundated marshes , we used the linear model parameters for frequently inundated marshes to create a vertical resilience index.
Because there were significant differences between methods, and 1963 to core-collection-date time-scale overlaps with wide-scale water-level gauging, we screened the sediment accretion estimates to only include those based on 137 Cs dated sediment cores. In addition, such 137 Cs accretion spans a long enough timeframe to provide average rates for decadal periods rather than short records which may be highly influenced by an unusual annual event. We excluded data from mangroves and other forested wetlands, submerged aquatic vegetation beds, and bay bottoms. If a source included freshwater wetland data as well as estuarine data, we included freshwater data points, provided that they were identified as "coastal" or "tidal" according to the source or NWI maps (Section 2.1).
We investigated the potential effect of differing regional suspended sediment concentrations (SSC), hypothesizing that regional SSC should correlate positively with accretion. We replicated the FWA-SSC calculations outlined by Weston (2014) and extended the analysis to all USGS gauges within watersheds that HOLMQUIST ET AL.
10.1029/2020EF001804 4 of 24 contained tidal wetlands. In order to calculate watershed average FWA-SSC (P80154), we downloaded SSC from the USGS web portal (USGS, 2016), specifying gauges with more than 50 discrete records in all states that have coastal wetlands. For each gauge, we also downloaded average daily discharge (P00060) data from the same tide gauge over the entire gauge period. We used the "R" Package "hydrostats" (Bond, 2015) to calculate base flow and base flow index for the entirety of the daily discharge record and calculated FWA-SSC as the average of all SSC taken on days when baseflow was greater than 20, meaning baseflow made up greater than 20% of the discharge (Weston, 2014) (Table S2; https://doi.org/10.25573/serc.12921746.v1). Despite slight differences in methodology between Weston's (2014) analysis and our own, the two analyses correlate nearly 1:1 (R 2 = 0.97) for those records which they have in common. Because of this correlation, we integrated FWA-SSC values for gauges included in Weston's 2014 analysis but not in ours. We calculated the average FWA-SSC for all watersheds that had one or more gauges falling within their borders using Spatial Join in ArcGIS (Esri Inc, 2017).
For each 137 Cs dated core in the accretion dataset, we first determined the nearest NOAA tide gauge (Permanent Service for Mean Sea Level, 2016) with an overlapping time frame with the core, with a 5-year tolerance, using MapTools in R (Bivand & Lewin-Koh, 2016). We then downloaded monthly mean water level data from NOAA servers (https://tidesandcurrents.noaa.gov/) using the package downloader (Chang, 2015). We calculated RSLR as the slope of a linear model with time as the independent variable and tidal elevation as the dependent variable using the linear model in R (R Development Core Team, 2014). For all cores, we also assigned the GT range from the nearest NOAA tide gauge (NOAA, 2017) with a reported GT range and complete data over the most recent datum epoch using a spatial join in ArcGIS Pro (Esri Inc, 2017; Table S3; https://doi.org/10.25573/serc.12921746.v1).
As  did, we attempted to control for differential responses of accretion to RSLR between more and less frequently inundated marshes. However, we encountered complexity and lack of consistency in terminology in our expanded literature survey. Therefore, we wrote a decision tree attempting to classify as much of the meta-analysis dataset as possible into two categories, higher energy, more mineral dominated systems, as frequently inundated marshes, and lower energy, less mineral dominated systems, as infrequently inundated marshes. According to the decision tree we classified coring locations based on the following criteria, in order: explicit descriptions as relatively frequently or infrequently inundated; if not that then further parsing of the relative tidal elevation and plant community descriptions high, mid, and low marsh; if not that then gleaning additional context from coring location descriptors; if not that then the presence of indicator plant species described as present and dominant at the coring locations.
First, if the study described sites as "frequently" or "infrequently inundated" we utilized those definitions. Second, an initial ANOVA showed that high and mid marshes had significantly lower accretion rates than low marshes, but not from each other (p < 0.001). So, for studies that compared "high," "mid" and "low" elevation marshes we categorized both high and mid marsh elevations as "infrequently inundated," and low elevations as "frequently inundated." Third, for other studies that compared marsh position or inundation regime in other ways we classified marshes that were described as "shoreward" (Lagomasino et al. 2013), "levee" (Craft, 2007;Hatton et al., 1983), or "streamside" (Craft et al., 1993;DeLaune et al., 1978) as "frequently inundated" and opposed to marshes described as "inland" (DeLaune et al., 1978), "back marsh" (Craft et al., 1993;Lagomasino et al. 2013;Sturdevant et al., 2002), or "plains" (Craft, 2007) as "infrequently inundated." For studies that compared natural to impounded sites, we always classified impounded marshes as "infrequently inundated" (Bryant & Chabreck, 1998;Sturdevant et al., 2002). We classified "natural" sites according to additional vegetation descriptors.
Inundation classifications and the criteria we used to determine them are recorded in Table S3 (https://doi. org/10.25573/serc.12921746.v1). We recognize that there is a degree of subjectivity that is unavoidable in this schema, and further discuss implications in our discussion section under Section 4.2. Limitations and Caveats. Because the assignment of inundation regimes was based partially on context gleaned from original papers which had a subjective element, we tested the effect of eliminating the inundation term.

Lateral Resilience Index Formulation
We calculated lateral resilience (LR) index as the ratio of the land cover amenable to marsh migration in the SLR zone, to modern-day coastal wetland cover below modern-day mean higher high water spring (MHHWS).

 LR Amenable Migration Area / Tidal Wetland Area
(2) Values one or higher indicate amenable areas are at or above a one-to-one replacement level given lateral wetland migration. Values between zero and one indicate amenable area are less than replacement level for given lateral wetland migration.
We define the SLR zone as the land below MHHWS at 2100 given an RSLR scenario. We assumed that amenable landcover types in the SLR zone included currently nontidal wetlands, uplands, and bare land. We excluded developed land, cultivated land, and open water. For tidal wetland area, we included all brackish and saline wetlands (referred to throughout as estuarine wetlands) and any freshwater (referred to throughout as palustrine wetlands) below the modern day MHHWS elevation.

Lateral Resilience Index Input Datasets
To calculate the LR index, we used three types of datasets, one for land cover, one for land elevation, and one for MHHWS elevation. For landcover, we used the 1996-2010 Coastal Change Analysis Program (C-CAP) (NOAA, 2013) to map major coastal landcover types to parse freshwater (a.k.a. palustrine) wetlands. We downloaded and mosaiced 1996 to 2010 C-CAP data for all states with coastal wetlands. We simplified the 2010 C-CAP's 22 land classes into tidal wetlands, nontidal wetlands, water, developed, cultivated, upland, and bare land (Table 1).
HOLMQUIST ET AL.

Table 1 How Coastal Change Analysis Program (C-CAP) Land Cover Classes Were Simplified and Treated as Candidates for Lateral Marsh Migration in our Lateral Resilience Analysis
To compare broad regional trends in future modeled accretion to the potential future inland marsh migration we utilized publicly available digital elevation models (DEMs), spatially extrapolated high spring tide levels from a previous study (Holmquist, Windham-Myers, Bernal, et al., 2018), and localized SLR scenarios from Kopp et al. (2014). For elevation, we downloaded DEMs associated with the NOAA SLR viewer as well as Louisiana, Maryland State, and San Joaquin Delta data sources. The DEMs are all the same as those used to map coastal lands by Holmquist, Windham-Myers, Bliss, et al. (2018) and are listed in supplemental information therein. All DEMs were converted to meters relative to the North American Vertical Datum of 1988 (NAVD88), if not already in those units, and had variable pixel size ranging from 1 to 10 m. DEMs were resampled to 30 × 30 m resolution using a "maximum area" rasterization in ArcGIS in order to spatially match the resolution of C-CAP.
DEMs went through two transformations, from meters relative to current day NAVD88 to meters relative to MHHWS. We used a NAVD88 to MHHWS conversion layer used to determine the extent of "coastal lands" for a U.S.-wide coastal wetland carbon inventory (Holmquist, Windham-Myers, Bernal, et al., 2018), and assumed an average of 17.3 cm of positive elevation bias for wetland surfaces, introduced by dense vegetation interfering with Light Detection and Ranging penetration (Holmquist, Windham-Myers, Bernal, et al., 2018). The MHHWS line was also used to classify Palustrine Wetlands into finer tidal and nontidal classifications using the elevation of MHHWS as a cutoff point, Palustrine wetlands below are classified as tidal, above are classified as nontidal. We shifted the datum a second time, converting meters relative to modern MHHWS to meters relative to MHHWS 100 years in the future. (3)

Sea-Level Rise Scenarios
We applied both the vertical and lateral resilience metrics by extrapolating tide gauge data, or predictions using GIS summarized by both state boundaries and intermediate watershed units. We calculated the vertical resilience index at the location of tide gauges (Table S4; https://doi.org/10.25573/serc.12921746.v1), then extrapolated them spatially using GIS. Similar to our comparison of accretion to RSLR for the core dataset, we queried the NOAA Permanent Service for Mean Sea Level dataset and calculated RSLR over a deliberate and consistent time frame. The linear modeling exercise resulted in RSLR, GT, and inundation regime as being significantly correlated with accretion. We calculated RSLR from 1996 to 2010, the same time frame as the national-scale coastal land cover change data used in the LR index. We selected gauges within a 1-year tolerance of that time frame and had at least 66% data completeness. This resulted in 75 different tide gauges to model net accretion in recent decades. We also modeled the accretion rate using likely 2000-2100 RSLR data for the same gauges (Kopp et al., 2014). We used median estimates from probabilistic local SLR models based on low, medium, and high impact climate change scenarios known as representative concentration pathways (RCPs). We used RCPs 2.6, RCP 4.5, and RCP 8.5, which assumed 2.6, 4.5, and 8.5 watts m −2 of radiative forcing by 2100 (Meinshausen et al., 2011). We assumed an average linear increase in sea level.
In addition to using localized RSLR scenarios to estimate vertical resilience ratios, we also calculated the LR index by extrapolating mean total RSLR predicted under RCP's 2.6, 4.5, and 8.5, raising the current MHH-WS boundary by those amounts, and summarizing the land cover types in the SLR zone (Equation 3). We define the SLR zone as the land between the current MHHWS boundary and the future MHHWS boundary given an RSLR scenario.
We extrapolated prediction surfaces from gauges using Empirical Bayesian Kriging in ArcGIS Pro (Esri Inc, 2017). Empirical Bayesian Kriging extrapolates point data to a surface using by modeling spatial autocorrelation between points and accounts for uncertainty in the semivariogram used to model spatial autocorrelation (Gribov & Krivoruchko, 2020). Settings that are needed to replicate the calculation are unmodified from generic inputs and include an output cell size of 300 m, no data transformation, a power semivariogram model, 100 maximum local points in each model, a local model area overall factor of 1, and 100 simulated semivariograms. For the search neighborhood, we used a standard circular neighborhood, 15 max neighbors and 10 min, 1 sector, an angle of 0, and a radius of 15.
We report VR, LR, and wetland area for intermediate-scale watershed units (Table S5;  We calculated regional and national summaries from watershed level statistics for both vertical and lateral resilience in two ways, by wetland area, and by wetland percentage. For each watershed, we classified vertical resilience as a binary; wetlands and watersheds were classified as either vertically vulnerable if VR was greater than or equal to 0. For LR we summarized wetland area and watersheds using a weighted sum. If the LR was less than 1, then the total wetland area was multiplied by LR to estimate the area of resilient wetlands, and (1-LR) to estimate the area of vulnerable wetlands. If LR index was ≥1, then all the wetland area in the watershed was classified as laterally resilient.

Linear Model Results
We fit five linear models into two subsets of data. First, we tested models 1-3 (Table 2) using a dataset that included all data which could be categorized by inundation frequency and had watershed-level FWA-SSC data (n = 94). Neither FWA-SSC (Table 2; model 1) nor time (Table 2; model 2) was significantly correlated with accretion rate. In model 3, which had only inundation class, RSLR, GT, and interactive effects between inundation class, both RSLR and GT were each significant (p < 0.01). Model 3 had a lower AICc score than Models one or 2 (Table 2) and had an overall p-value of <0.001 and an R 2 of 0.54 (Table 2). For this model RSLR (in mm yr −1 ) and GT (in meters) were both positively correlated with accretion for frequently inundated marshes, and RSLR was positively correlated with accretion in infrequently inundated marshes. Accretion was significantly higher for infrequently inundated marshes.
Because model 1 indicated that FWA-SSC did not have a significant impact on accretion, we created models 4 and 5 using a larger subset of data in which cores could be categorized as frequently or infrequently inundated (n = 269). Once again time was a nonsignificant input variable (model 4) and removing it reduced the AICc score (model 5; Table 2). Overall model 5 had an R 2 of 0.33 and a p-value <0.001 ( Figure 2; Table 2) and indicated that inundation class, RSLR, GT, and interactive effects between inundation class contributed to variability in vertical accretion rates (p < 0.001).
For the subset of data without FWA-SSC and with inundation regimes (n = 269), eliminating the inundation term still resulted in an overall significant model (p < 2.2e−16), although the model explained less of the variance in the data (R 2 value decreased from 0.33 to 0.23), and resulted in a less parsimonious model (AICc increased from 1,204 to 1,222; Table 2). When the inundation term was removed, GT was no longer a significant parameter (p = 0.92). For the larger total dataset (n = 478), the same model form was significant (p < 2.2e−16), and R 2 increased from 0.23 to 0.30; GT was still not a significant parameter without the inundation term given the larger dataset.

Vertical Resilience Index Results
We predicted net accretion for frequently inundated marshes as a metric with which to broadly compare resilience to RSLR across major geographic gradients. VR is described by Equation 4, in which RSLR is in mm yr −1 and GT is in m (Figure 2).
At 75 tide gauges, predicted net accretion varied regionally between 1996 and 2010 based on GT and RSLR (Figure 3). The most vulnerable marshes would be those associated with the Grand Isle, Louisiana gauges where RSLR was 7.4 mm yr −1 and GT was 0.323 m (  When median local 2100 projections of RSLR based on RCP 2.6 (Kopp et al., 2014) were applied in place of 1996-2010 values, and GT was assumed to be static, many locations switched from aggrading to submerging ( Figure 3). Nine out of 10 gauges from the mid-Atlantic, and six of the more southern Northeast sites switch from aggrading to submerging under RCP 2.6. Under the more extreme RCP 8.5 values, portions of the Southwest and the Southeast, as well as more of the Northeast, switch from aggrading to submerging (Figures 4a, 5a and 6a). In the northernmost portion of the Atlantic Coast, as well as much of the Northwest, macrotidal conditions and projected RSLR is low, and Equation 3 predicts net accretion will remain positive.
In our analysis, we see that threshold beyond which marshes fail to accrete relative to RSLR is regionally variable, depending on GT, and in many cases is already crossed, or will likely be crossed by 2100. If we set VR to zero, and rearrange Equation 4, we get Equation 5, in which we solve for RSLR for a specific GT in meters. This is an RSLR threshold for when net accretion would from positive to negative.

Lateral Resilience Index Results
In addition to regional vulnerability due to RSLR and GT gradients, there is also substantial regional variation in the area available for marsh migration. Within the 198 watersheds for which there was adequate DEM data, "palustrine" wetlands were the dominant low-elevation land cover class, dominating 100 watersheds in total (Figure 3). "Palustrine" wetlands are a prominent feature along the CONUS South Central coast, as well as the Southeast and the mid-Atlantic (Figures 4b, 4c, 5b, 5c, 6b and 6c). Marsh migration potential is highest along the coast of North Carolina (Figure 5b) because low elevation "palustrine" areas are extensive in proportion to current tidal wetland areas.
Watersheds with relatively low migration potential values can broadly be classified into three different conditions. First, they can represent areas where the SLR zone is dominated by "cultivated" or "developed" lands. Southern California and much of the North Eastern Seaboard between New Jersey and New Hampshire are examples of the "developed" dominated classes (Figures 4b, 4c, 5b, 5c, 6b and 6c). Second, they can represent watersheds that do have adjacent undeveloped inland migration space, but that space is limited Scatterplot showing the variables used in the vertical resilience index, relative sea-level rise, inundation category, and tide range. The frequently inundated datapoints were used to create the index. The dashed line shows a one to one ratio between relative sea-level rise and accretion rate, and which data used to calibrate the index are aggrading, above the line, and submerging, below the line.
by the relief of adjacent slopes, such as the North-Eastern Coast of Maine ( Figure 6b). Third, they can represent watersheds in which inland migration space exists but in relatively low proportion to the area of tidal wetlands. This occurs in Louisiana where the tidal wetland area within some watersheds is very high (Figure 5b).
Regions varied generally from north to south in the proportion of land that was available to accommodate lateral wetland migration. The northwest (Oregon and Washington) had the least, 21% under RCP 4.5. This was followed by the Southwest (CA) with 21% and Northeast (from Maine in the North to New Jersey in the South) with 23% ( Table 3). The South Central (Louisiana, Texas, and the rest of the Gulf Coast), Southeast (Atlantic Florida south-north to North Carolina), and the mid-Atlantic all had much higher proportions of land available for wetland migration. The mid-Atlantic zone had 43%, the Southeast had 44%, and the South Central zone had the most with 50%, under RCP 4.5.

Paired Vertical and Lateral Resilience Metrics at Watershed Scale
At a watershed scale, the trends in vertical and lateral resilience indices were negatively correlated (Figure 7), meaning that watersheds that were more likely to be vulnerable to losing elevation to SLR, were also more likely to have adjacent areas available to accommodate marsh migration. Likewise, marshes that are more likely to have inland migration space available are less likely to be able to vertically accrete at pace with SLR.
For individual watersheds, both indices follow latitudinal and regional trends. The Pacific Coast is generally below a lateral replacement level of one to one, and above a mean net accretion rate of zero. In the Northwest, wetlands are the most resilient vertically and have the least available migration space. The coast of California is similarly restricted in terms of lateral migration, but also loses VR from north to south (Figure 7). The California coast is near a mean net accretion rate of 0, meaning vertically vulnerable in RCP 8.5, with mean accretion rates in some watersheds lower than the rate of RSLR. The Northeast similarly is below replacement level laterally, and vertically resilient in the north, which declines to the south. For the Southeast, mid-Atlantic, and south-central zones, vertical resilience declines southward, but LR increases. Regionally the highest proportions of doubly vulnerable wetlands are in the Mid-Atlantic under RCP 2.6 (59%), and in the Northeast (58%) under RCP 4.5 (Table 3).  vulnerable wetland in the Southeast and Northeast increased substantially from RCP 2.6 to 8.5; lateral and vertical vulnerability in the mid-Atlantic drops slightly from RCP 4.5 to 8.5.

Discussion
The goal of this paper was to synthesize existing data to extract spatial trends of vertical and lateral tidal marsh resilience to SLR in the contiguous U.S. and in closely adjacent Canada. We use the most extensive synthesis of accretion rates compiled thus far allowing us to test the inter-comparability of sediment accretion using a standardized method ( 137 Cs dated sediment cores) for developing the VR index. This meta-analysis adds new insight into the relationship between observed accretion, RSLR, and tidal range along the U.S. coast and reveals some striking regional and latitudinal differences. Paired with the GIS analysis of lateral resilience, our results suggest that at a broad spatial scale both current and future tidal wetland vulnerability is, and will be, spatially variable, with accretion response falling off southward, because RSLR, tidal range, and the space available for potential marsh migration are spatially correlated with latitude.
The fact that accretion was not sensitive to FWA-SSC was contrary to our hypothesis and we think necessitates some additional examination of the modeling literature, as well as the dataset itself. First, modeling work by Vinent et al. (2019) posits that for micro-tidal estuaries, accretion in interior marshes is insensitive to FWA-SSC and that organic matter accretion is the dominant driver of total accretion. In our core dataset, the higher range FWA-SSC values we record come from the Mississippi River delta region of HOLMQUIST ET AL.

10.1029/2020EF001804
12 of 24 the Gulf of Mexico, where micro-tidal conditions could make accretion relatively insensitive to the potential input. Second, modeling work by Morris et al. (2016) argues that low SSC, high productivity, and resulting organic-dominated soil accretion is typical of Atlantic Coast conditions.  showed that the majority of sampled marsh soil volume in the U.S. is dominated by organic matter.
On closer examination, two cores in our review dataset had inconclusive accretion rates, in which the actual accretion rate was higher than could be captured by the core depth. One core came from the Mississippi delta, in Louisiana (DeLaune et al., 1992), a region that has high SSC and a prograding delta (NOAA, 2013). The other came from a mineral-dominated low marsh setting near San Diego, California (Weis et al., 2001). We hypothesize that we did not see a significant relationship between FWA-SSC and accretion because of a combination of geomorphic, spatial, and sampling phenomenon: that there are negative spatial correlations between FWA-SSC and GT, that wetlands dominated by mineral accretion are less frequent than those dominated by organic accretion, and that it is harder to measure accretion using 137 Cs dating in faster accumulating mineral dominated systems because the 1963 peak is more likely to be below the maximum depth of the barrel style corer, typically used for this type of sampling.

Regional Trends in Resilience and Vulnerability
The two metrics we use to compare marsh resilience or vulnerability across CONUS-scale gradients, modeled net accretion and marsh migration potential, indicate a nuanced and potentially complex future for marsh vulnerability. DeLaune and White (2012) posited that the high rates of RSLR experienced by the Louisiana Gulf Coast, the highest historically observed in the CONUS, may be the closest analogs we have to 2100 projected SLR. For the Gulf coast, our model (Equation 1) predicts continued vertical deficits under both RCP 2.6 and 8.5. RSLR scenarios. But, as Kopp et al. (2014) point out, these RCP scenarios are based on trends in past anomalies that may be too pessimistic because fluid extraction, and associated submergence may decline in the future (Intergovernmental Panel On Climate Change, 2014). Although we find that the South Central zone is most vulnerable to elevation loss, the area has lowlands primarily dominated by land cover class types favorable to migration ("palustrine," "upland") and few watersheds dominated by "cultivated" or "developed" land. Nevertheless, because some ocean-facing watersheds in Louisiana contain extensive existing wetlands, inland migration potential is low-that is the limited area of nonwetland land available to migration is not proportional to the existing coastal wetland area. A few of these ocean-facing watersheds with low inland migration potential could be underestimating migration potential if migration can occur beyond the watershed boundaries that make up the fundamental unit of this analysis.
The reason for the southward decline in vertical accretion capacity remains unknown, but may have to do with several intercorrelated geographic patterns including higher rates of historic and projected RSLR from north to south, and decreasing tidal range from north to south. Future RSLR is projected to approach and exceed the rates (7.4 at Grand Isle, Louisiana;1996   averages because the loss of mass in ice sheets results in an overall weakening of their gravitational fingerprint (Engelhart & Horton, 2012;Kopp et al., 2014;Tamisiea & Mitrovica, 2011). The continued weakening of the AMOC will also contribute to RSLR (Sallenger et al., 2012). Many of these future potential vulnerable areas have low elevation land cover that is dominated by "palustrine" or "upland," and two mapped watersheds being dominated by "developed" land. Migration potential is relatively low, likely because "cultivated" land is often secondarily dominant in these watersheds, or because of topographic constraints. Compared to the Mid-Atlantic the Southeast is both less vulnerable to RSLR and to coastal squeeze. Tidal HOLMQUIST ET AL.

10.1029/2020EF001804
14 of 24 Thousands of ha (% of regional wetland area)  ranges are wider and projected RSLR is lower south of Cape Hatteras. There are also extensive low lying "palustrine" dominated areas in high proportion to the current marsh area, exceeding 100% in some cases.
Southwest and Northeast projections differ from those of the mid-Atlantic and Southeast in that they are less vulnerable to losing their capacity to accrete, but more vulnerable to coastal squeeze (Figure 7). These zones have both lower projected RSLR and wider tidal frames than other regions. This was observed by the recently published MARS index (Raposa et al., 2016). The Southwest and Northeast each have extensive low-lying areas predominantly dominated by development (Figures 4c and 6c), which would presumably halt any potential marsh inland migration (Linhoss et al., 2015). Some areas, such as Southern California, not only have small tidal wetland areas but are also the most potentially limited by lateral migration barriers (Figures 4b and 7). The Southwest contains a high proportion of low-lying agricultural land, especially around the Sacramento Delta, and in Northern California (Figure 4c). If RSLR is higher than the median likely scenarios predicted by RCP 8.5, then many Southwestern marshes assumed resilient by other metrics, could be particularly vulnerable to coastal squeeze, as seen in the recently published WARMER model which predicts extreme habitat losses for the Pacific Coast under high SLR scenarios . The CONUS Northwest is the least vulnerable to losing accretion capacity of all regions; migration potential is also high for the Northwest compared to the CONUS Southwest and Northeast.
The median GT of our calibration dataset (1.6 m) results in an RSLR threshold of 6.7 mm/yr (Equation 5).
For the sake of comparison, this threshold is only slightly higher than and estimates by Morris et al. (2016). They used a first-principles soil formation modeling approach and estimate that total steady-state accretion HOLMQUIST ET AL. should not be greater than 5 mm yr −1 given typical Gulf Coast and Atlantic conditions. Different regions likely have variable RSLR thresholds because of latitudinal and another spatial variability in tidal properties.
We hypothesize that the trends that we see in both vertical and lateral vulnerability are ultimately driven by latitudinal geographic patterns in postglacial geology, tidal ranges, and human settlement patterns. In more southerly regions, there is a more ongoing isostatic rebound from the last glacial period, meaning that there are higher projected rates of RSLR. Southerly latitudes also have smaller tidal ranges, thus more vertically vulnerable. More northerly regions have a less isostatic rebound, and some regions are isostatically uplifting, meaning slower projected RSLR. More northerly regions also have larger tidal ranges, leading to more VR.
Similarly, lateral vulnerability also follows latitudinal patterns. We hypothesize that in the further South Central region, lower topographic relief leads to more of a gradation between tidal and nontidal wetlands meaning that there is more amenable land to facilitate inland marsh migration. In the more Northern and Western regions, we hypothesize that higher relief land leads to less of a gradation. Also in more Northerly Latitudes, we hypothesize much more low elevation land is devoted to development and agriculture.

Limitations and Caveats
There are limitations to this analysis to be addressed in future research using process-based modeling and vetted site-appropriate data Swanson et al., 2014;Thorne et al., 2015). First, the net accretion model presented here was calibrated using modern data and extrapolated out into the future using projected SLR. In other words, future conditions may be out of the boundaries of our calibration dataset. For example, the coast of Maine has a combination of high projected SLR and wide tidal frames; most of our historical high SLR data is from areas with small tidal ranges. Although the linear modeling here and assumptions we have made may be viewed as simplistic, it allows the use of current rates of accretion for a first-cut resilience and vulnerability assessment that serves as a comparative metric, or an a priori hypothesis for future work. This could be approached next by using a targeted study applying process-based modeling over local to continental-scale gradients.
Similar to our analysis of the vulnerability of wetland accretion, our analysis of marsh migration is based on some key assumptions, and is meant to present a comparative metric based on current conditions. The analysis does not take into account artificially built tidal restrictions present now, or in the future. It also assumes a direct one-to-one replacement of noncultivated, nondeveloped land cover types with tidal wetlands as sea level rises. However, salinity intrusion into freshwater systems and wetland plant competition with upland species are both biogeochemically and ecologically complex; more research is needed to determine whether or not marsh migration into adjacent uplands and palustrine wetlands can keep pace with edge erosion and pool expansion .
Our meta-analysis found that 137 Cs dating is by far the most common dating technique available in the literature, with 480 out of 1,004 observations (Table S1; https://doi.org/10.25573/serc.12921746.v1). Drexler et al. (2018) point out that the utility of 137 Cs dating will have a limited lifetime, that quality of dating is likely geographically variable because of differing rates of fallout in different regions, and that peaks are prone to error because of erosion, redeposition, and sediment mobility. The 1963 137 Cs peak quality can be highly uncertain because of the lack of published raw profiles and lab errors and will become less reliable in the future as the peak decays exponentially with a half-life of 30.17 years. Our study reanalyzes data collected from 1978 to 2011 (Table S2), so the fact that the quality of 137 Cs is predicted to decline in the future is not as much an issue as it would be if we were producing new data. As far as depositional issues and lab errors, one of the advantages of synthesis studies is that combining data from multiple sources and depositional environments allows us to extract trends despite lab or site-specific effects. The regression modeling accretion as a function of RSLR, GT, and relative flooding frequency had an R 2 of 0.33, which leaves room for future studies exploring other factors affecting accretion. We recommend future studies apply hierarchical models to disentangle the effects of environmental drivers from data quality and site-level variability (LeBauer et al., 2013). Our study did not synthesize full profiles of 137 Cs but could have benefited from them if they were widely available. Future studies should use best practices in archiving 137 Cs information including full radioactivity profiles and lab errors (Drexler et al., 2018). Finally, we propose that we researchers need a way of propagating uncertainty in 137 Cs dating, similar to what we have 14 C (Blaauw & Christen, 2011) and 210 Pb dating (Aquino-López et al., 2018). Doing so could help us disentangle environmental signals from data quality issues as 137 Cs dating using the 1963 peak inevitably becomes less reliable into the future.
Despite the caveats of using the present and recent past to infer the future, such scenario development has value. Our analysis, based upon the largest such dataset of linked marsh accretion, RSLR, tidal range, and FWA-SSC data, provides evidence that tidal wetlands of the CONUS, and particularly the South Central zone are vulnerable to current rates of SLR. Many more in the CONUS Southeast, and mid-Atlantic, such as the Chesapeake Bay, will likely be vulnerable to even conservative estimates of future RSLR. Although tidal wetlands in the CONUS Southwest and Northeast are comparatively less vulnerable to losing their capacity to accrete relative to SLR, they are more vulnerable to coastal squeeze because of topographic constraints, and the dominance of developed and/or cultivated land in many watersheds.

Future Research Directions
Tidal wetlands are vital ecosystems that serve as nurseries for fishery species (Boesch & Turner, 1984), carbon sinks (Mcleod et al. 2011), and are home to endangered and endemic species (Rosencranz et al., 2019). Previous efforts analyzing their resiliency to RSLR (Enwright et al., 2016;Thorne et al., 2018), as well as current monitoring programs (Steyer et al., 2003), have largely been regionally focused. However, national-scale action such as funding, guidance, and goal setting, will likely be needed to ensure wetland persistence and function into the future as sea levels rise.
In our linear modeling of vertical vulnerability, we learned that the inundation regime is important for two reasons. Without it the model explains 10% less of the variance in the data and GT is not a significant parameter. The fact that GT's role in accretion is conditional on separate consideration of inundation regime Points to a potentially important physical link between the two that a more advanced process model could potentially elucidate. Inundation regime based on a somewhat subjective list of criteria. The problem we had was that there is not a consistent definition of relative inundation across the broad array of studies analyzed. Future studies could have a more nuanced integration by using a separate data model for inundation class, having multiple researchers code inundation status, and integrating more uncertainty for less straightforward classifications. Futures studies could also attempt at a more unified concept of tidal inundation statuses such as dimensionless tidal elevation (Janousek et al., 2016) or flooding time (Hickey, 2019). A continuous driver as opposed to a binary definition would not be as affected by regional definitions or judgment calls made by data interpreters.
What we cannot learn from this dataset and linear models is what marsh collapses resulting from accretion deficit will occur, whether they will be gradual or sudden, instigated by disturbance events such as hurricanes, and what the early stages look like. We hypothesize that monitoring with sediment elevation tables (SETs; Cahoon et al., 2002) could catch early stages of collapse by observing persistent accretion deficit in more frequently flooded, more mineral-dominated sections of the marsh. However, conclusive evidence supporting this would take an expanded database, additional monitoring, and possibly a more process-rich modeling approach. We should also anticipate places and times where thresholds are passed or as it could mean marsh collapse events could spatially and temporally cluster such as following El Niño events (Goodman et al. 2018), or correlated to the lunar nodal cycle (Peng et al., 2019).
Throughout the text, we have referred to vertical and lateral resilience indices because we think the metrics presented here are far from mature enough to be thought of as a forecast, used for actionable decision making and land use planning. For future research, we recommend more extensive data-model integration and the development toward a forecasting system (Dietze et al., 2018) by integrating process-and data-models, using the statistical modeling framework here for model intercomparison or model averaging, and further developing the indices presented here to inform management goals.
First, we recommend a forecasting system that integrates environmental drivers, process models, and data models. Drivers from this study that could be integrated include considerations of regional topography, low elevation land cover class, and localized RSLR scenarios. Process models could include members of the Marsh Equilibrium Model (MEM; Morris et al., 2002) family and WARMER . A more sophisticated data model could account for differences in data quality across time and between studies (Drexler et al. 2018). Bayesian hierarchical models can propagate uncertainties in environmental drivers, parameters, and models themselves (Dietze et al., 2018;Hobbs & Hooten, 2015) Second, near-term forecasting (Dietze et al., 2018), which is increasingly being used as a tool for promoting actionable decision support and collaboration between scientists and decision-makers, benefits from Bayesian models and from model selection (Conn et al., 2018) and model averaging techniques. The vertical vulnerability index, which is a simple linear model could be used to performance test more complex models against.
Third and finally the indices here could be used to inform goals for restoration planning. VR index, which is simply a modeled accretion deficit, could be used to set goals and both sediment and financial budgets for sediment diversion (Snedden et al., 2007) or sediment augmentation projects (Ulibarri et al., 2020). The two components of LR index-topography and low elevation land cover--could help planners determine acreage goals for securing marsh migration space through negotiation of conservation easements, restrictive covenants, and outright purchases (Field et al., 2017).  observed that marsh vulnerability can be overstated if dynamic feedback between marsh elevation and inundation is not taken into account, and potential inland migration of marshes is ignored. Our analyses support the existence of dynamic feedback because marsh accretion keeps pace with RSLR up to a point. Our analysis supports estimates by Morris et al. (2016) that estimates 5 mm yr −1 as a typical upper limit for accretion in most Atlantic and Gulf Coast marshes are already near there. Using a median GT results our VR equation outputs an RSLR threshold of 6.7 mm yr −1 . RSLR threshold may be spatially variable, lower in places with smaller tidal ranges and higher in places with larger tidal ranges. We see much of the Gulf Coast already passed their threshold, and much of the Atlantic Coast passing it by 2100 under even modest RSLR scenarios. However, our analyses also take into account the magnitude and geography of projected RSLR, comparative trends in low elevation land cover, as well as the fact that microtidal marshes are additionally vulnerable to SLR. Given our results, we believe that marsh sustainability into the future should not be generally assumed. In fact, if the amassed data and indices we present are representative, we expect future marsh vulnerability to be spatially variable with many currently stable marshes switching from aggrading to submerging, and some marshes failing to migrate inland because of topography and a lack of unoccupied low-elevation land cover.

Acknowledgments
This study was supported by the Department of Interior establishing grant #Y561461:03 for the Southwest Climate Science Center and USGS research grants #G12AC20505, #G13AC00083, #G14AP00178. James R Holmquist was supported by the Smithsonian Institution, as well as grants from NASA Carbon Monitoring System (#NN-H14AY67I) and the National Science Foundation (#1655622) while writing and revising this manuscript.