Quantifying Slopes as a Driver of Forest to Marsh Conversion Using Geospatial Techniques: Application to Chesapeake Bay Coastal-Plain, United States

Coastal salt marshes, which provide valuable ecosystem services such as flood mitigation and carbon sequestration, are threatened by rising sea level. In response, these ecosystems migrate landward, converting available upland into salt marsh. In the coastal-plain surrounding Chesapeake Bay, United States, conversion of coastal forest to salt marsh is well-documented and may offset salt marsh loss due to sea level rise, sediment deficits, and wave erosion. Land slope at the marsh-forest boundary is an important factor determining migration likelihood, however, the standard method of using field measurements to assess slope across the marsh-forest boundary is impractical on the scale of an estuary. Therefore, we developed a general slope quantification method that uses high resolution elevation data and a repurposed shoreline analysis tool to determine slope along the marsh-forest boundary for the entire Chesapeake Bay coastal-plain and find that less than 3% of transects have a slope value less than 1%; these low slope environments offer more favorable conditions for forest to marsh conversion. Then, we combine the bay-wide slope and elevation data with inundation modeling from Hurricane Isabel to determine likelihood of coastal forest conversion to salt marsh. This method can be applied to local and estuary-scale research to support management decisions regarding which upland forested areas are more critical to preserve as available space for marsh migration.


INTRODUCTION
Salt marsh survival through the end of this century is threatened as sea level rise (SLR) continues to accelerate at rates unseen for two millennia (IPCC, 2013). In the context of this paper, salt marsh refers to "estuarine intertidal emergent wetlands" defined by Cowardin classification as areas with hydrophytes, tidally influenced, with at least sporadic access to ocean waters and dilution by freshwater runoff from land (Cowardin et al., 1979). Given the essential ecosystem services provided by marshes, such as habitat for key species, carbon sequestration, and protection from storms (Barbier et al., 2011;Möller et al., 2014;Narayan et al., 2017), it is crucial to understand how and where marshes will persist. As sea level rises, salt marshes respond by vertically accreting both organic and inorganic sediment (Morris et al., 2002;Kolker et al., 2009); where accretion rates are outpaced by sea level rise, marshes degrade and eventually drown (Reed, 2002). Marshes are additionally vulnerable in the lateral direction from wave erosion and sediment deficits that degrade the marsh at the seaward edge and can lead to marsh loss (Kirwan et al., 2016). Salt marshes can migrate inland where upland is available for conversion to marsh (Williams et al., 1999;Enwright et al., 2016), responding to sea level rise and counteracting marsh loss from those lateral processes. Marsh migration inland is well documented on the Atlantic and Gulf coasts of the United States (Williams et al., 1999;Smith, 2013;Raabe and Stumpf, 2016;Anisfeld et al., 2017;Schieder et al., 2018;Schieder and Kirwan, 2019), however, the mechanisms that control this process are less firmly established.
Marshes on the Atlantic seaboard commonly transgress into upland coastal forests and farmland (Kirwan and Gedan, 2019). The process of forest conversion to marsh occurs in two stages which take place on differing time scales. First, as sea level rises, upland ecosystems are more regularly inundated with saline water. The increase in soil salinity and saturation kills seedlings, saplings and understory vegetation making room for marsh vegetation to encroach (Brinson et al., 1995;Fagherazzi et al., 2019). In the second stage, extreme events such as hurricanes can kill the mature trees through uprooting and/or canopy loss. The loss of regenerative capacity from the first stage reduces the ability of coastal forests to be resilient to these extreme events (Kirwan et al., 2007;Fagherazzi et al., 2019), which are increasing in frequency due to climate change (Melillo et al., 2014). With the canopy removed and light penetrating to the forest floor, marsh vegetation moves into the available understory. Along the North American Atlantic and Gulf coasts, tree stumps surrounded by marsh vegetation, termed "ghost forests," identify locations where migration has occurred and is likely still occurring (Kirwan and Gedan, 2019).
Rates of marsh transgression into coastal forest range regionally from close to zero in New England (Field et al., 2016) to almost 7 m/yr in the Mid-Atlantic (Hussein, 2009). Slope, salinity, rate of sea level rise, frequency and duration of inundation, groundwater, and soil properties are some of the key factors that play a role in the likelihood and rate of migration (Brinson et al., 1995;Hussein, 2009;Field et al., 2016;Fagherazzi et al., 2019). When other factors are controlled for, variation in marsh migration rate on both local and regional scales has been attributed to the slope of land the marsh is transgressing over (Brinson et al., 1995;Hussein, 2009;Smith, 2013;Field et al., 2016). Low-lying coastal plains with a slope of ∼0.0004, such as those in Florida (Raabe and Stumpf, 2016), provide the greatest accommodation space for marshes to migrate as uplands are inundated due to sea level rise . By contrast, marshes along the glaciated New England coastline often abut steep bluffs which limit migration. As slope decreases, low-lying areas are inundated more frequently and for longer periods (Hussein and Rabenhorst, 2001). Additionally, there is a reduction of drainage area that accompanies decreasing slope which reduces freshwater input to the system (Hussein, 2009). The increase in salinity to the system, coupled with a reduction in freshwater inputs makes upland conditions more favorable for marsh migration in areas with low slopes. However, in certain field studies, the correlation between marsh migration rate and slope is weaker than expected (Schieder et al., 2018). Additional research is needed to better understand the interplay between slope and other factors influencing migration rates.
When considering marsh migration potential, it is most useful to know the slope of land perpendicular to the marshforest boundary as that is the direction of transgression. Current methods of determining the slope across the marsh-forest boundary rely on field measurements using real-time kinematic with global navigation satellite systems (RTK GNSS) or taking sediment cores to determine the slope of past transition zones (Hussein, 2009;Anisfeld et al., 2017;Schieder and Kirwan, 2019). These methods are time and resource intensive, which make regional assessments impractical. Using geospatially complete, high resolution elevation data, generated from remote sensing methods such as lidar, is a more appropriate option for regional analyses. However, this requires a new method to calculate slope across land use boundaries.
The new slope quantification methodology developed in this study repurposes the Digital Shoreline Analysis System (DSAS) (Himmelstoss et al., 2018) to cast transects perpendicular to marsh-forest boundaries and a digital elevation model to determine the average slope across each transect. We apply these geospatial techniques to identify and analyze marsh-forest boundaries in the coastal-plain surrounding the Chesapeake Bay, which has well-documented occurrences of marsh migration (Kirwan et al., 2007;Hussein, 2009;Schieder et al., 2018;Kearney et al., 2019;Schieder and Kirwan, 2019;Gedan et al., 2020;Burns et al., 2021). Seamless high-resolution elevation, land cover and land use datasets that improve the precision of this method are available for the region. Ongoing field studies of marsh migration in the region provide measurements of slope across the marsh-forest boundary in several locations to assess the accuracy of our results. We also discuss additional factors that contribute to marsh migration and demonstrate the use of a dynamic inundation model to estimate marsh migration likelihood. This method quantifies the geospatial variability of slope across the marsh-forest boundary that can be applied both locally and regionally to provide valuable information to land managers for prioritizing acquisition of coastal forests to allow for future marsh transgression.

MATERIALS AND EQUIPMENT Data
The method was developed using ArcDesktop version 10.6.1 and the Digital Shoreline Analysis Software (DSAS) version 5.0 (Himmelstoss et al., 2018). All data layers used in this study were converted to the same projected coordinate system, North American Datum (NAD) 1983 Universal Transverse Mercator (UTM) Zone 18N.

Salt Marsh
Salt marsh data were downloaded from the U. S. Fish and Wildlife Service (USFWS) National Wetlands Inventory (NWI) (U. S. Fish and Wildlife Service, 2018). NWI analysts identify and delineate wetland types with GIS technology using high altitude imagery with collateral data sources and field work (U. S. Fish and Wildlife Service, 2018). Imagery in the Chesapeake Bay region is primarily from the 2000s and 2010s (U. S. Fish and Wildlife Service, 2018). Marsh coverage in the coastal-plain surrounding the Chesapeake Bay was created by combining the Maryland and Virginia state marsh extents.
This study considers salt marshes to be estuarine intertidal emergent wetlands as classified by NWI following Cowardin classification (Federal Geographic Data Committee, 2013). Marine, lacustrine, and palustrine wetland categories are not considered as part of this study. Typical regularly flooded marsh vegetation in this region includes Spartina alterniflora and Schoenoplectus americanus (Perry et al., 2001;Schieder et al., 2018). Typical irregularly flooded marsh vegetation includes Spartina patens and Juncus romerianus (Perry et al., 2001). Phragmites australis was also noted in the dataset and has been identified as a transition species during forest to marsh conversion (Perry et al., 2001;Smith, 2013).

Coastal Forest
Forest coverage was obtained from the Chesapeake Conservancy's Chesapeake Bay High-Resolution Land Cover and Land Use Data Projects datasets (Chesapeake Conservancy, 2018a; Chesapeake Conservancy, 2018b). The land cover dataset is based on land cover conditions in the National Agriculture Imagery Program images from 2013/2014. The land use dataset was created from the land cover dataset that was then further modified using 13 ancillary datasets, such as zoning, parcel boundaries, and floodplains.
This study uses the "forest" class (value 8) from the land use dataset which consists of all standing trees and areas of tree harvest contiguous over one acre, and "tree canopy and shrubland" class (value 2) from the land cover dataset which consists of both deciduous and evergreen woody vegetation over 2 m tall, including stand-alone individuals and discrete clumps. Although there was no characterization made of tree species in these categories, loblolly pine (Pinus taeda) forests dominate uplands adjacent to marsh in this region (Kirwan et al., 2007).

Elevation
Elevation was defined using a topobathymetric digital elevation model (TBDEM) at 1-m horizontal resolution for the Chesapeake Bay region, which was obtained from the U. S. Geological Survey (USGS) Coastal National Elevation Database (CoNED) (Danielson and Tyler, 2016). The data layer was put together from over 369 different data sources such as topographic and bathymetric lidar point clouds, hydrographic surveys, side-scan surveys, and multi-beam surveys (Danielson and Tyler, 2016). Accuracy of the dataset is discussed in Results.

Study Area
We developed and tested our slope quantification method along the coast of Chesapeake Bay, United States ( Figure 1). The Chesapeake Bay is one of the largest estuaries in the world, draining 166,000 km 2 of watershed across 6 states and Washington, D.C. (Perry et al., 2001). The Bay is microtidal with a maximum mean tidal range of approximately 1 m near the mouth (Perry et al., 2001).
The mid-Atlantic coast of the United States, and in particular, Chesapeake Bay, are hotspots of relative sea level rise (sea level measured with respect to the surface of the Earth (IPCC, 2013)), due to a weakening of the Gulf Stream coupled with rapid regional subsidence from glacial rebound (Engelhart et al., 2009;Sallenger et al., 2012). As a result, historical rates of relative sea level rise in the Chesapeake Bay, measured by local tidal gauges, ranged from 3 to 6 mm/yr 1 compared to global mean sea level rise (sea level change with respect to ocean water volume) of 3.2 mm/yr from 1993 to 2010 (IPCC, 2013). However, estimates of modern sea level rise rates in Chesapeake Bay in 2011 were 4-10 mm/yr (Ezer and Corlett, 2012). High rates of sea level rise force ecosystems to adapt on a shorter time scale, which makes the Chesapeake Bay an ideal location for studying marsh-forest conversion on event-based and decadal timescales.
Due to the large geographic extent and density of marsh-forest boundaries in the study area, the analyses were performed in 10 tiles (2 columns x 5 rows) of the total extent (Molino et al., 2020).

Marsh-Forest Boundary
Marsh extent was based on the NWI categorization of "Estuarine Intertidal Emergent" wetlands. The original projection was NAD 1983 Albers and we converted it to UTM Zone 18N. The forest extent was determined by merging the areas classified as "forest" from the land use dataset and the areas classified as "tree canopy and shrubland" from the land cover dataset. Both raster datasets were projected from their original United States Contiguous Albers Equal Area Conic USGS version projection to UTM Zone 18N. The land cover and land use forest rasters were converted to shapefiles and merged to create a single forest layer. Any overlap between the wetlands layer and forest layer was assumed to be wetlands, as later verified by examining aerial images.
Forest polygons that were within 10 m of wetlands were considered potential migration zones for marshes and their boundaries were considered marsh-forest boundary. Interior holes in the forest polygons were removed to limit the selection to the exterior forest boundaries. These holes, areas without forest, were usually from a house or pond inside the forested area. Additionally, forest polygons with areas less than 900 m 2 were excluded from analysis. These polygons represented 1.5% of the total forest area considered in the study. While they were present across the study area, more were removed from the area to the east of Chesapeake Bay.
Once the forest polygons representing marsh-forest boundary were identified, multiple steps were taken to simplify the geometry created when the forest raster was converted to polygons in a previous step (Molino et al., 2020). Buffers were created inside the forest polygons to remove complex edge features, and then buffers were created back out to the original extent. Overlapping forest features created from the buffering were merged which removed dissecting roads. This was followed by a polynomial approximation with exponential kernel (PAEK) smoothing with a 30 m tolerance to reduce sharp angles along the boundaries. The sharp angles were often a remnant of the original forest raster file which approximates a forest edge using square pixels. The polygons retained their original size, with a simplified shape, which is necessary to reduce the file size and processing time in subsequent steps.

DSAS
We repurposed the USGS Digital Shoreline Analysis System (DSAS) version 5.0 to cast transects perpendicular across the marsh-forest boundary. Although predominantly used to compute changes in beach shorelines, DSAS has previously been used to cast transects at the marsh-forest boundary (Smith, 2013).
DSAS requires a baseline from which to cast perpendicular transects, and shorelines to determine the length of transects. As DSAS was created to study shoreline change of beaches, these terms refer to those features. However, in this study, we were interested in the slope values within 10 m of the marsh-forest boundary, therefore we created a "shoreline" and a "baseline" on opposite sides of the marsh-forest boundary. Specifically, we first generalized the marsh-forest boundary by smoothing with a PAEK algorithm with a 30 m tolerance and then created a 10 m offset outside and inside of the boundary to define a shoreline and a baseline, respectively. As these polyline features extend for the entire perimeter of the original forest polygon, they include areas of the forest which do not border marsh. Transects cast along these areas are removed in a later FIGURE 1 | Study area with the NWI estuarine intertidal emergent wetlands in teal and forest extent obtained from Chesapeake Conservancy Land Use and Land Cover datasets in tan. Inset has been rotated 90°clockwise from actual position and is the same location as Figure 2.
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 616319 step. We chose 20 m transects to increase the likelihood that this method captured the marsh-forest boundary in scenarios where the boundary based on the geospatial data might have been offset from the physical boundary due to inaccuracies in the original marsh and forest files, or due to the smoothing algorithm (Marsh-Forest Boundary).
Following the requirements in the DSAS User Guide (Himmelstoss et al., 2018), we set transect spacing to 30 m apart along the baseline with a smoothing distance of 0 m. DSAS requires that all resulting transect files be saved in an ArcGIS personal geodatabase. The resulting transect file is limited in size to approximately 65,000 transects per DSAS run. If the estimated number of transects (total forest polygon perimeter for a subregion divided by the transect interval) is larger than this number, the subregion must be divided into sub-subregions.
In general, the parameters used in this section of the study, such as transect length and distance between transects, were customized for the Chesapeake Bay coastal-plain and can be modified for other regions with different coastal morphology and marsh-forest characteristics.

Assignment of Slope to Transects
Once the transect file was created, we removed any transects cast further than 10 m from salt marsh using the estuarine intertidal emergent wetlands shapefile created in a previous step. We also removed any transects which were cast across artificial forest edges created when we split the polygons into subregions and sub-subregions as these did not reflect a true marsh-forest boundary.
Using the "Add Surface Information" tool in ArcMap, transects were then assigned an average slope in percent rise based on the TBDEM (Figure 2). A single slope for each transect was calculated from the weighted average of slopes between 1 m cells that the transect crosses. Transects with a slope value of exactly 0 were removed as this meant they spanned areas that had artificial elevation data. Transects that have very low slope are typically not exactly 0, therefore, a value of exact 0 reflects a transect located entirely over an area with a single fill value for elevation. The origin of the fill values in the TBDEM is discussed in Uncertainty and Completeness.
Finally, transects from each tile and subregion were merged into a single dataset that describes the geospatial distribution of slope across the Chesapeake Bay coastal-plain (Molino et al., 2020) (Figure 3).

Uncertainty and Completeness
Ultimately, the accuracy of the slope values is limited by the accuracy of the source data. The TBDEM is derived from 369 topographic and bathymetric datasets acquired between 1859 and 2014. The elevation data was obtained primarily through lidar, with a vertical accuracy of 15-20 cm in root mean square error (Danielson and Tyler, 2016). Additionally, the creators of the TBDEM used a fill value of -1 to denote areas where hydroflattening (i.e., assignment of a fill value in water bodies with no bathymetric data) occurred, which does not reflect the underlying topobathy (Danielson and Tyler, 2016). Areas of artificial elevation also occurred where the underlying source datasets for the TBDEM were extrapolated to a 1 m resolution from a lower resolution. We used a focal statistics tool in ArcMap to identify these areas and the impacted transects. Transects spanning areas with fill values represent 2.5% of the entire dataset. Transects that partially overlap areas with artificial elevation values from either of these processes have been left as part of the dataset for future updates as new TBDEM data FIGURE 2 | Location of marsh-forest boundary is the same as the inset in Figure 1. (A) Shoreline (yellow) and baseline (pink) input files for DSAS set 20 m apart, between which (B) transects, blue, are cast with a 30 m spacing (C) Removed extraneous transects cast between forest and non-marsh land types and assigned remaining transects a slope value from TBDEM (indicated by color gradient of green, low slope, to red, high slope).
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 616319 become available. However, these transects are not considered during our analysis of the slope dataset nor the assessment of marsh migration likelihood. Similarly, the accuracy of the marsh-forest boundary is based on the accuracy of the marsh and forest extents. The marsh dataset is based on aerial imagery primarily from the 2000s and 2010s. The forest datasets are based on aerial imagery from 2013/ 2014. There are undoubtedly land use and land cover changes that have occurred since then that are not captured in the datasets and thus resulted in transects that are not located exactly along the marsh-forest boundary. For many of these, the center point of the transect, which is the approximated marsh-forest boundary, is located over a creek or river. This results in a negative elevation value that is an artifact of the processing and does not reflect the actual elevation of the marsh-forest boundary along that transect. Additionally, there are likely inaccuracies in land classification; however, the land cover dataset showed approximately 90% overall accuracy with a 98% accuracy for the forest land cover category (Pallai and Wesson, 2017). We did not include transects with a negative elevation value at the midpoint in our analysis given that these commonly reflected artifacts from hydro flattening or processing steps. This impacted less than 10% of our transects.
Despite these challenges, this method of determining slope along the marsh-forest boundary provides the first geospatial product of slope values across the entire Chesapeake Bay coastal-plain. As more accurate information becomes available, the marsh-forest boundary and slope values themselves can be updated. Additionally, if extreme accuracy is required for localized analyses of specific forested areas, this method can easily be scaled down and used with site-specific elevation and land cover data.

Accuracy Assessments
We conducted several assessments to determine how the parameters set in DSAS influence the slope values of the transects. Given that the parameter values we chose were FIGURE 3 | Average slope (percent rise) along transects across marsh-forest boundaries throughout the Chesapeake Bay coastal-plain. Additionally, to assess how the vertical accuracy of the TBDEM influences the transect slope, we took elevation measurements in the field at three locations around Chesapeake Bay which represent low, medium, and high slope environments. At each site, a GNSS base station was set on a tripod over a temporary benchmark set up near the transect. The wooden stake used for the benchmark was driven into the marsh approximately 1 m, in a position with open skies to the south for good satellite coverage. The base station collected satellite readings continuously for approximately 4 h, while a GNSS rover and total station were connected to the base station in order to conduct an integrated survey. Using this setup, a topography survey was conducted perpendicular to the marshforest boundary that gives latitude, longitude, and elevation with an average precision of 0.0092 m in the horizontal, 0.0109 m in the vertical dimensions, and 0.0360 m in terms of elevation (averages based on fifteen surveys completed in the summer of 2019

RESULTS
This method successfully produced a dataset of slope values on an estuary-scale that can be used to determine geospatial patterns in slope essential to studies of marsh migration. This methodology produced a dataset that contains 217,200 transects which represent the slope values across over 6,500 km of marshforest boundaries in the entirety of the Chesapeake Bay coastal-plain. The median average slope of all the transects in the Chesapeake Bay region is 4.0% rise. Median is used to describe the slope dataset as the values follow a Rayleigh distribution. Over 90% of the slope values range from nearly 0 to 14% rise. However, less than 3% of transects have a slope value less than 1%. Slope values increase as the elevation of the marsh-forest boundary increases when the elevation values are sorted into 1-m bins (R 2 0.99, p < 0.001 from Pearson (r) score) ( Figure 4A). The elevation bin of 0-1 m has 81% of slope transects, the 1-2 m bin has 14% of transects, the 2-3 m bin has 3% of transects, the 3-4 m bin has 1.2% of transects, and the +4 m bin has 1.1% of transects ( Figure 4B). Transects with a marsh-forest boundary elevation of greater than 4 m are often associated with transects that cross man-made structures that divide these ecosystems. Additionally, there are also a small number of transects which have very high marsh-forest boundary elevations because the marsh abuts a steep slope or bluff on which there are trees or which acts as a divide between ecosystems. If the marsh and forest extents are off by a meter or so in these locations, and thus the marsh-forest boundary is off by a meter or so, the elevation of the marsh-forest boundary will be higher than expected. When the data are further examined within these bins, there are clusters of slope values at certain elevations. These typically represent a single geographic area within the Chesapeake Bay coastal-plain that has similar slope values at similar elevations. The slope values across the marsh-forest boundaries in the Aberdeen Proving Grounds, located in region 7 in Figure 5, are an example of a cluster of low slope values in the dataset.
The lowest slope values across the marsh-forest boundary are in region 7 which has median slope of 2.6% ( Figure 5). The highest slope values are located in region 6 which has a median slope of 6.3% ( Figure 5).
The results of the accuracy assessments showed minimal changes in slope values at the two test locations, represented by white triangles in Figure 5. In the first test, casting transects at resolutions of 20 and 30 m apart, changed the slope values by less than 0.5% ( Table 1). The second test, which shifted the placement of the transects by 15 m, altered the slope values by less than 2% in both test locations (Table 1). In the third test, forest polygons at each site were processed without the smoothing step: the difference in the resulting transect values were 0.3% on the western side and 2.7% on the Delmarva peninsula ( Table 1). The results of these tests lead to the conclusion that the parameters and processing steps do not strongly influence the slope value for transects produced by this method.
Additionally, in our assessment of how vertical accuracy of the TBDEM influences the transect slope, the average percent rise calculated for the transects generated in ArcGIS corresponded well with slope values of those taken at the corresponding field locations, represented by black circles (Figure 5). While there was a vertical offset in the data, which varied by site, the geospatial patterns held true, i.e., the transect with the lowest slope from the field data also had the lowest slope using the TBDEM. This vertical offset ranged from −0.2 m to 0.1 m at points along the transects, with the TBDEM on average 0.01 m higher than the field measurements. The difference in slope between the field data and TBDEM altered the percent rise of the transects between 0.7 and 1.1% rise.

Geospatial Patterns of Slope Across the Marsh-Forest Boundary
Using this slope quantification method, we have created an estuary-wide dataset of slope values across marsh-forest  Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 616319 8 boundaries in the Chesapeake Bay coastal-plain. Interestingly, less than 3% of transects cross boundaries with a slope value considered favorable for forest to marsh conversion. The majority of these transects are located on the Virginia portion of the Delmarva peninsula (region 10 on Figure 5), where salt marsh gently slopes upward for kilometers into wide expanses of coastal forest. Given the key role of slope in controlling marsh migration rates, this has significant implications for directing conservation efforts throughout the Chesapeake Bay coastal plain.
High slope values are common along the tributaries which segment the land to the west of Chesapeake Bay where marshes abut steep bluffs leading to coastal forest. In these areas, marsh migration is blocked by topographical highs, limiting migration inland. However, there are segments of marsh-forest boundary associated with wide expanses of gently sloping marsh where the slope increases rapidly and unexpectedly (>5.0% rise). Further examination of the aerial imagery (Google Earth Pro) and the elevation data often reveals a structure, such as a berm, dividing the marsh and forest ecosystems. In the Chesapeake Bay region, especially on agricultural lands, individual homeowners have historically built berms to protect their property from rising sea levels (Kirwan and Gedan, 2019). It is not uncommon in our dataset to see transects with unexpectedly high slopes across the marsh-forest boundary, where the forest abuts farmland on the opposite side. In these cases, berms were likely built to protect the farmland from inundation, and this artificially higher land provided favorable conditions for forest growth. As these features were put in by the individual landowners, there is no official record of them, and they are often difficult to identify from aerial imagery. A dataset of berms and other built barriers would be invaluable to an assessment of migration potential. Interestingly, this method offers the first step at identifying these structures for future work assessing barriers to marsh migration in the region.
Simple estimates of marsh migration rates throughout Chesapeake Bay can be calculated using the median slope (4.0%) from this dataset and modern sea level rise rates in Chesapeake Bay of 4-10 mm/yr (Ezer and Corlett, 2012). This predicts that salt marsh will transgress inland at rates of 0.1-0.25 m/yr. Given the wide variation of slope values discussed above, these rates will range significantly throughout the coastal plain. Migration rates have the potential to reach meters per year on the Delmarva peninsula, while no migration may occur where natural or man-made barriers exist inland. Our ultimate goal with the establishment of this method is to encourage analysis of marsh migration on an estuary-scale using slope data alongside other factors to assess marsh migration potential.

Application of Method to Estimate Marsh Migration Likelihood
Here we detail an example application of this method with the intent of supporting and spurring more comprehensive studies of marsh migration in this region. We analyzed storm inundation data for Hurricane Isabel, provided by the Advanced Circulation (ADCIRC) Prediction System developed by the University of North Carolina at Chapel Hill. Hurricane Isabel, which made landfall in September 2003, is considered to be one of the most influential tropical cyclones to impact central and eastern Virginia in 50 years (Beven and Cobb, 2004). The inundation from the storm extended throughout the Chesapeake Bay making it ideal to examine the influence on marsh-forest boundaries throughout the region. Storm surge in the southern portion of the Bay was 1.5-1.8 m above predicted astronomical tidal level. In the central portion of the Chesapeake Bay, storm surge was observed to be 0.9-1.5 m, with surge values reaching 1.8-2.4 m above predicted astronomical tidal levels in the upper reaches of the Bay (Beven and Cobb, 2004).
We extracted maximum storm inundation depth (meters) and duration of storm inundation (hours) at the midpoint of each slope transect from the ADCIRC datasets. Then these data were assessed in combination with the slope dataset to assign migration likelihood at each location. The values in all three datasets were divided into 5 categories, which were assigned a score of 0-4 corresponding to values that reflect an increase in migration likelihood (Table 2), with 0 being little to no influence on migration and 4 a large increase in migration likelihood.
The categories for the slope dataset were assigned using slopes associated with different geomorphologies. Slope values 0-0.1% represent gentle upland slopes typical of coastal plains, moderate upland slopes of 0.1-1% represent formerly glaciated coasts, and steep upland slopes of 1-20% are typical of active margin coasts (Kirwan et al., 2016). It should be noted that geomorphic processes beyond glaciation and tectonics, such as river incision, can cause steep slopes; we have maintained these categories to provide consistency with the previous assessment of slope on marsh expansion and so that this method will apply to geomorphic settings beyond Chesapeake Bay. Low slope environments provide greater accommodation space for migration than moderate and steep upland slopes of glaciated and active margin coasts . Additionally, low slope coastal environments are more frequently inundated and don't drain as well as high slope environments, both factors creating more favorable conditions for marsh migration inland over gentle upland slopes (Hussein and Rabenhorst, 2001;Hussein, 2009). This method of categorizing slope to assess likelihood of ecosystem changes has been previously applied to coastal environments (Pendleton et al., 2004;Torio and Chmura, 2013). However, we further subdivided the 1-20% category (steep upland slope) into two bins, 1-5% and 5-20% slope. A number of small-scale berms built by homeowners to protect their land from inundation were detected by our slope quantification method (Results); these structures range in slope from about 5 to 20%. While not all transects with a slope between 5 and 20% cross berms, this category acknowledges a reduction in migration potential as the slope increases past 5%, regardless of if a manmade barrier exists or not. Scores for all three variables at each point were added together and averaged for a final score of marsh migration potential that ranged from 0 to 4. The average final score of marsh migration potential is 1.1, indicated in red in Figure 6 as low potential for marsh migration. Marsh migration potential based on these factors does point to several areas of high migration potential (blue areas) throughout 2 | Categories for slope, maximum storm inundation, and storm inundation duration, and the assigned score for each. A score of 0 reflects the lowest likelihood of marsh migration while a score of 4 represents the highest likelihood based on these three factors. Chesapeake Bay coastal-plain ( Figure 6), which correlate well with the spatial patterns of migration rates measured throughout Chesapeake Bay coastal-plain in Schieder et al. (2018). In particular, the eastern coasts of regions 2 and 3 and the southern portion of region 9 ( Figure 5), have both higher measured migration rates by Schieder et al. (2018) and are areas of high migration potential in our study. Overall, Figure 6 demonstrates that there are relatively few locations in the Chesapeake Bay where marsh migration is very favorable; areas which are considered favorable to migration are clustered predominantly in region 9 with a few smaller areas in regions 2, 3, and 7. This example stands as one application of the method described to determine slopes across the marsh-forest boundary and how it might be used in ongoing research. We have identified several additional data layers worth examining for more in-depth future marsh transgression potential indexes: • Salinity: Modeling systems capable of nowcasting and forecasting salinity levels for individual estuaries would further inform inundation data. • Sea level rise: The National Oceanic and Atmospheric Administration (NOAA) has released global and regional sea level rise predictions for multiple probabilistic ranges which can be examined locally for estuary-specific studies (Sweet et al., 2017 -Caraballo et al., 2015). Maximum inundation, inundation duration and frequency can be obtained from both datasets to further inform migration potential.
Once acquired, further assessment of how these variables impact geospatial patterns of marsh migration and if there are certain threshold levels (for example, of salinity), which amplify migration rates, is needed.
Such work would also benefit from additional studies measuring marsh migration rates across the Chesapeake Bay coastal-plain to further corroborate our results.

Broader Utility of Method
The slope quantification method established in this study provides an opportunity for researchers to assess the slope across marsh-forest boundaries on a scale not previously possible with the time and resource limitations of field work. Researchers and organizations can potentially apply this method at smaller scales to improve understanding of marsh migration at a local level and support conservation decisions. Organizations can use local marsh and forest geospatial data, which are less time consuming to corroborate with field work and can alter the DSAS parameters to better fit the characteristics of the marsh-forest transition zones in that area. For example, Gedan et al., 2020 estimated probabilities of upland to marsh conversion based on elevation and distance to shoreline for Somerset Country, MD (region 9 in Figure 5). Our method would be well suited to quantify slope across land use type boundaries to add to these probability assessments. As new data become available, such as land cover/land use data from more recent imagery or an updated TBDEM, the methodology is easily replicable to update the marsh-forest boundary or the slope values. Additionally, land managers interested in buying forested areas for conservation efforts can incorporate this method to determine the most likely area for marsh migration inland. Those interested in marsh migration in the previous century could use historical photographs or topographic sheets to identify the marsh-forest boundary and assign slope data from older elevation datasets. Other major estuaries where marsh migration studies are underway (e.g., Narragansett Bay) can supplement their research efforts with this method. Further, the slope quantification method is not restricted to use at the marshforest boundary; this method can be applied to marsh boundaries with other upland environments such as agricultural lands.

CONCLUSION
To our knowledge, this slope quantification method and the published dataset are the first of its kind and provide critical information on marsh survival to a rapidly growing field of study. We recognize that the number of ground-truthed transects available for our accuracy assessment represents only three locations out of over 200,000 transects in our dataset. We acknowledge that this method is not a substitute for the accuracy of field measurements as it is limited by the precision and accuracy of the input datasets. However, this method does provide a reasonable and consistent estimate of slopes across an area too extensive for field work and establishes a starting point for more targeted studies of high interest locations. Further, the reproducible workflow presented in this study can be applied more locally, allowing for utilization of higher resolution datasets and application of numerous field measurements for verification purposes. Our goal is to have this method and its future applications better inform land managers and others interested in conservation to ensure coastal salt marshes continue to provide essential ecosystems services in the coming decades.

AUTHOR CONTRIBUTIONS
GM, NG, and JC designed the overall study; GM and ZD developed and applied the geospatial methods. AA developed the oceanographic model analysis. GM analyzed the results with aid from all authors. JC collected and provided field data for accuracy assessment of results. GM wrote the manuscript with feedback from all authors. All authors provided critical feedback and helped shape the research, analysis and manuscript.

FUNDING
Funding for this study was provided by the United States Geological Survey's Coastal/Marine Hazards and Resources Program and Ecosystems Mission Area.