Scholarship, Research, and Creative Work at Bryn Mawr College Scholarship, Research, and Creative Work at Bryn Mawr College

: Elevation is a major driver of plant ecology and sediment dynamics in tidal wetlands, so accurate and precise spatial data are essential for assessing wetland vulnerability to sea-level rise and making forecasts. We performed survey-grade elevation and vegetation surveys of the Global Change Research Wetland, a brackish microtidal wetland in the Chesapeake Bay estuary, Maryland (USA), to both intercompare unbiased digital elevation model (DEM) creation techniques and to describe niche partitioning of several common tidal wetland plant species. We identified a tradeoff between scalability and performance in creating unbiased DEMs, with more data-intensive methods such as kriging performing better than 3 more scalable methods involving post-processing of light detection and ranging (LiDAR)-based DEMs. The LiDAR Elevation Correction with Normalized Difference Vegetation Index (LEAN) method provided a compromise between scalability and performance, although it underpredicted variability in elevation. In areas where native plants dominated, the sedge Schoenoplectus americanus occupied more frequently flooded areas (median: 0.22, 95% range: 0.09 to 0.31 m relative to North America Vertical Datum of 1988 [NAVD88]) and the grass Spartina patens , less frequently flooded (0.27, 0.1 to 0.35 m NAVD88). Non-native Phragmites australis dominated at lower elevations more than the native graminoids, but had a wide flooding tolerance, encompassing both their ranges (0.19, −0.05 to 0.36 m NAVD88). The native shrub Iva frutescens also dominated at lower elevations (0.20, 0.04 to 0.30 m NAVD88), despite being previously described as a high marsh species. These analyses not only provide valuable context for the temporally rich but spatially restricted data collected at a single well-studied site, but also provide broad insight into mapping techniques and species zonation.


INTRODUCTION
Saline and brackish marshes are vegetated ecosystems with specialized species that are adapted to live within the intertidal zone.These marshes provide myriad ecosystem services, such as fisheries habitat (Barbier et al. 2011), protection of coasts from storms (Doughty et al. 2017), and carbon storage (Chmura et al. 2003, Ouyang & Lee 2013, Howard et al. 2017).The wetlands are vulnerable to direct degradation (Brophy et al. 2019), and their resilience or vulnerability to sea-level rise (SLR) is uncertain, although they have some capacity to accrete in equilibrium with SLR (Kirwan & Megonigal 2013).Vulnerability or resilience depends on a system's capacity to capture suspended sediment (Morris et al. 2002), the ability of its plants to form soil mass by belowground root addition (Nyman et al. 2006), and the space available for upslope transgression (Kirwan et al. 2016).
Models used for forecasting marsh elevation's dyna mic response to SLR all require accurate elevation as an initial condition.These models include the Marsh Equilibrium Model (MEM; Morris & Bowden 1986, Morris et al. 2002, 2016, Morris 2006), MEM variations that include more detailed spatial analyses (Schile et al. 2014), Hydro-MEM (Alizad et al. 2016, 2018, Bacopoulos et al. 2019), Wetland Accretion Rate Model for Ecosystem Resilience (WARMER; Swanson et al. 2014, Thorne et al. 2018), models by Kirwan and Mudd (Kirwan & Mudd 2012), and the SLR Affecting Marshes Model (SLAMM; Park et al. 1989).Commonly used pointbased models such as the MEM, WARMER, and Kirwan-Mudd models also require information on plant tolerance to flooding, for which relative position in the tidal range is used as a proxy (Janousek et al. 2019).Small relative elevation changes can affect phenomena such as vegetation zonation along gradients as much as a few centimeters (Castillo et al. 2000).Information to parameterize these models comes from plants grown under controlled conditions at varying tidal elevations, so-called marsh organ experiments (Kirwan & Guntenspergen 2015, Langley et al. 2013, Mozdzer et al. 2016), as well as high resolution GPS surveys of elevation and biomass (Schile et al. 2014) or plant presence and absence (Thorne et al. 2018).
Monitoring multiple effects of global change on tidal wetland processes, including tidal elevation, has been a major focus of the Global Change Research Wetland (GCReW), located in Kirkpatrick Marsh in the Chesapeake Bay estuary, Maryland (MD), USA.Various wetland grasses (C 4 ), sedges (C 3 ), and the invasive Phragmites australis (Cav.)Trin.ex Steud have been the subject of studies investigating the physiological tolerances and niches of species using experimental, marsh organstyle, floo ding treatments at GCReW (Langley et al. 2013, Mozdzer et al. 2016).GCReW has not yet had an extensive elevation survey (Nelson et al. 2017), which limits the ability to verify that the modern elevations represented in the marsh organ experiments are representative of the landscape.
The most scalable techniques for producing digital elevation models (DEMs) use light detection and ranging (LiDAR); however, LiDAR data are prone to bias in tidal wetlands because laser pulses often do not fully penetrate low, dense vegetation or standing dead biomass and because wet soil can absorb some light, leading to obfuscated returns after post-processing (Hladik et al. 2013, Medeiros et al. 2015, Buffington et al. 2016).These biases yield modeled elevations that are higher than those measured using ground-based techniques.Elevation errors are particularly problematic in coastal wetland systems since subtle vari ations in elevation result in different plant communities and functions (Alizad et al. 2020).
There are methods for generating unbiased wetland DEMs which range in complexity from simple adjustments based on height to data-intensive statistical correction models.For a national scale assessment, Holmquist et al. (2018) corrected the elevation bias by subtracting a single average offset estimated from a literature review, for a contiguous USA-wide probabilistic area estimate of coastal lands.This single offset approach did not account for spatial variation in the bias term, so national-scale products could be enhanced with more localized analysis.Some studies correct elevation zones using maps of vegetation community distribution (Hladik et al. 2013, Schile et al. 2014).Buffington et al. (2016) applied a multivariate linear regression to correct bias in tidal wetlands using the LiDAR elevation estimate and the Normalized Difference of Vegetation Index (NDVI), a measure of vegetation greenness mapped from National Agriculture Imagery Program (NAIP) imagery.This method is called the LiDAR Elevation Adjustment using NDVI (LEAN) algorithm.Other studies have avoided LiDAR altogether, instead spatially interpolating a 3-dimensional surface among point-based surveys (Thorne et al. 2013(Thorne et al. , 2018)).While such surveys are essential for sea-level resiliency assessments, more extensive field-based collection of elevation data is costly and damages sensitive sites, so there is a need for techniques that can be applied across larger areas with a minimal number of calibration points.
Furthermore, despite the intensive, long-term experiments at GCReW, no previous study has documented the natural history of important vegetation− elevation relationships at the site.Previous studies have focused on 3 major communities that dominated the site when experiments were first started in 1987.Schoenoplectus americanus (Pers.)Volkart ex Schinz & R. Keller (formerly Scirpus olneyi) occupies more frequently inundated portions of the marsh and a C 4 grass mix of Distichlis spicata (L.) Greene, and Spartina patens (Aiton) Muhl.dominates less frequently inundated portions of the marsh.Marsh organ experiments at GCReW have been able to constrain a minimum and maximum elevation tolerance for S. americanus, but have not constrained the upper growing limit for C 4 grasses (Langley et al. 2013).This lack of a documented upper limit for high-marsh grasses is not unique to GCReW; it has also been observed for marshes on the opposite shore of Chesapeake Bay along Blackwater Bay, MD (Kirwan & Guntenspergen 2015).An invasive lineage of the common reed P. australis dominates major areas of the site.This species has been present in the mid-Atlantic since before 1910 and has expanded greatly since the 1960s (Rice et al. 2000, Saltonstall 2002), expanding its area at GCReW 22-fold since 1972 (McCormick et al. 2010).Documenting the relationships between elevation and P. australis dominance could provide insight into its status as an invasive species, as well as provide information on how it may affect marshes in the future.Finally, there are far more than 3 species present at the marsh.For example, the native shrub Iva frutescens L. has not been integrated into any global change experiments despite being a dominant component of common vegetation communities at GCReW (Byrd et al. 2018).
Given that there are multiple methods for creating unbiased DEMs with differing levels of data requirements, and that the extensive experimental data at GCReW need site-wide elevation to add context to and provide comparison for marsh organ studies, our goals were 3-fold.First, to survey elevation and relative plant cover of the entire marsh complex at 20 × 20 m resolution and to document the distribution of elevation relative to tidal datums from a single year.Second, to use these elevation data to generate 4 different unbiased DEMs utilizing different techniques, then assess their accuracies and precisions and quantify the tradeoffs between performance and scalability.We hypothesized that there would be a tradeoff between scalability and performance, with more labor-intensive methods providing better performance.Third, to quantify wetland plant presence and relative dominance in terms of observed elevation/inundation niches for different plant species.We hypothesized that dominant plant communities partition elevation niches, with S. americanus being located in more frequently flooded areas and C 4 grasses located in less frequently flooded areas.

Site description
GCReW is located at the Smithsonian Environmental Research Center in Edgewater, MD, USA.The site is one of the most studied tidal wetlands in the world and is the source of many impactful insights on relationships among plant traits, microbial activity, and soil biogeochemical responses to global change factor interactions.GCReW occupies about half of a larger 22 ha marsh complex known locally as Kirkpatrick Marsh (Fig. 1).The site has been commonly referred to as a high marsh (Langley et al. 2013, Lu et al. 2019), but this status as a high marsh has not been explicitly quantified relative to tidal flooding across the whole site.The site is brackish throughout the year but has seasonal and interannual variance in salinity due to freshwater input from precipitation.The system is microtidal, with a greater diurnal tidal range of 0.44 m, as measured at a nearby tide gauge in Anna polis, MD (NOAA 2020a) over the last datum period (1983− 2001).In the Chesapeake Bay, relative SLR was 3.67 ± 0.20 mm yr −1 from 1928−2019 (NOAA 2020b).A portion of this is eustatic SLR, and a portion is attributable to subsidence, 1.3−1.5 mm yr −1 in the region according to GPS (Karegar et al. 2016).The climate zone is classified as temperate with no dry season and hot summers.

Tidal datum transformations
We summarized elevation data relative to tidal datums that were calculated from the nearby Annapolis tidal gauge located 13 km from GCReW.We calculated a suite of tidal datums for 2016 from complete, verified, 6 min tide gauge data.We used the 'ftide' function in the R package 'tideharmonics' (Stephenson 2016) to calculate a set of 60 harmonic constituents.We then predicted tides at the same 6 min resolution as the data using nearest neighbor analysis to isolate high and low tides, then differentiated higher high tides from lower high tides, and lower low tides from higher low tides.We used a 52 tide moving window and reclassified the highest high tide as spring higher high water and the lowest predicted tide of the lunar period as the lower low spring tide.
We calculated mean sea level (MSL), highest observed tide (HOT), and lowest observed tide (LOT) by calculating the mean, minimum, and maximum of the observed data set.We calculated mean high water (MHW), mean higher high water (MHHW), mean higher high water spring tides (MHHWS), mean low water (MLW), mean lower low water (MLLW), and mean lower low water spring tides (MLLWS) by taking the average water levels of the classified predicted water levels.We quantified the highest and lowest astronomical tide (HAT and LAT) from the minimum and maximum of the predicted water levels.

Improvements to geodetic control network
We strengthened the geodetic control network property-wide by performing 4 long occupations of fixed monuments with JAVAD company total station GPS units, one of which was used as the base station for the current study.Two of these occupations were done on geodetic markers on the north side of the Rhode River on monuments installed by the National Geodetic Survey and the Smithsonian Institution.For the 2 survey-grade monuments on the north side of the Rhode River, GPS units were mounted on a 2 m tripod on a 0.25 m adapter with a 0.025 m attachment.
In wetlands, soil elevation is prone to dynamic changes because of sediment consolidation, expansion of the rooting zone by root growth, and contraction of the rooting zone by soil organic matter decomposition.Therefore, geodetic monuments need to be attached to rods driven several meters into the ground, hitting bedrock (Cahoon et al. 2002).At GCReW, we occupied 2 such deep-rod monuments, one near a tidal creek near a flume and marsh organ experiments, and the other closer to the upland interface near the sediment elevation table (SET) used to monitor net-elevation change in an experimental plot (Langley et al. 2009).For the flume location, the total station GPS unit was mounted on a 1.5 m tripod with a 0.25 m adapter and a 0.025 m attachment.The equipment was centered on the monument on a pointed tip and was leveled using a bubble level built into the tripod.At one SET (referred to as SET #17 in the supplemental data release; Holmquist et al. 2021), the unit was mounted on a 2 m surveygrade pole and a 0.025 m attachment.Instead of being centered on a monument point, it was screwed into a SET receiver using a customized 0.339 m adapter.The unit was supported by 2 bipods and was levelled using a bubble level attached to the top of the adapter with adhesive.All elevation measurements were made relative to GEOID12b and are reported in m relative to the North American Datum of 1988 (NAVD88).All occupations occurred over approximately 48 h.Total station GPS data were post-processed using OPUS Projects software (Mader et al. 2012).

Real-time kinematic GPS survey
Two Trimble R8 real-time kinematic (RTK) GPS receivers were used for elevation surveys of the marsh.The first was deployed as a base station attached to a SET.The base station was attached to the SET receiver by a customized mount of known height (0.339 m), and a 2 m long survey-grade pole extension rod.The unit was leveled with a bubblelevel attached with adhesive.The survey rod was stabilized using 2 bipods with their feet stuck into the adjacent soil.The base station used the known elevation from the total station occupation (mean ± SE: 0.362 ± 0.016 m NAVD88) to correct and constrain elevations measured using the roving receiver.
The RTK GPS roving receiver was mounted on a 2 m survey-grade pole and outfitted with a flat rounded bottom about 6 cm in diameter that prevented the unit from sinking into soft or muddy soil surfaces.We generally considered the surface of the marsh to be the point where the unit came to rest under its own weight.In cases where the rover bottom sank below the soil surface under its own weight, we held the rod at the soil surface at the point of initial resistance and noted the adjustment in the elevation plot notes.The nominal accuracy of RTK GPS is 0.04 m.Regularly, at the end and/or beginning of a workday, we measured the same spot on a boardwalk adjacent to the base station and found the standard deviation of the repeated measurements to be 0.006 m, less than the nominal accuracy of the instrument.
We completed 2 separate elevation surveys, one of the vegetation plots measured as part of the Smithsonian's Tennenbaum Marine Observatory Network (TMON), and the other, a 20 × 20 m grid over most of Kirkpatrick Marsh.The 37 TMON plots at Kirkpatrick Marsh and GCReW are arranged in 3 transects spanning vegetation and elevation gradients (Fig. 1).We surveyed each corner of the 1 × 1 m plots on 22 June 2016.
For the 20 × 20 m grid survey, we created an outline of the contiguous section of the wetland using the most recent (at the time) image basemap in ArcMap v.10.6 (ESRI, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, SSGS, AeroGRID, IGN and the GIS User Community; accessed 2016) and polygons delineated by the National Wetlands Inventory (NWI; US Fish and Wildlife Service 2014).We digitized the outlines of major features on the marsh such as boardwalks and experimental chambers.At the time of this survey, some experiments and new boardwalks had not been captured in the high-resolution Google Earth Ima gery, so we surveyed the boundaries of these experiments and boardwalks.In ArcGIS, we created a 20 × 20 m grid using the gridlines function, converted the grid to center points, then subsetted those points that intersected the marsh and omitted those within a 5 m buffer of site infrastructure or the TMON plots (Fig. 1).
We surveyed the entire 525 point grid between 26 July and 15 August 2016.If an elevation sample point did not intersect the marsh as mapped because of trees, creeks, or marsh edges, the point was moved to nearest access, and we noted this change.In some places, we opportunistically surveyed the edges of the marsh if they were either near a sampling point or if the edge of the marsh was not visible in the imagery used to plan the survey.In one case, a large section of P. australis-dominated marsh was growing under trees and was not apparent in the imagery used to plan the survey (Fig. 1).RTK GPS data were post-processed using Trimble Business Center v.3.7.

Producing unbiased elevation maps
We implemented 4 methods for creating an unbiased DEM for GCReW.These included (1) applying a single average offset from a literature review of coastal wetland studies that performed bias corrections to LiDAR-based DEMs (Holmquist et al. 2018), (2) using the LEAN algorithm to correct the LiDARbased DEM (Buffington et al. 2016), (3) applying vegetation community-specific corrections to correct the LiDAR-based DEM, and (4) using spatial interpolation among points instead of the LiDAR-based DEM.For the 3 LiDAR-based techniques, we used an Anne Arundel County LiDAR-based DEM created from LiDAR flown in March 2011.The resulting 1 m resolution DEM was produced by the Sanborn Map Company for county-wide FEMA flood map moderniza-tion.The LiDAR flight and our RTK GPS survey were 4 yr apart, but we assumed that any accretion or subsidence over that time would be within the nominal error of the RTK GPS.
We used the LEAN technique to correct the vertical bias in LiDAR-based DEMs due to dense vegetation coverage (Buffington et al. 2016).LEAN uses ground measurements of elevation, typically from RTK GPS, to calculate the vertical bias in LiDAR data sets.Initial LiDAR elevation and the NDVI ( = [NIR − Red] / [NIR + Red]; red wavelength = 608−662 nm, near infrared [NIR] wavelengths = 833−887 nm) from NAIP were used in a multivariate linear regression model to estimate the bias.NDVI was assumed to be a metric of vegetation height and density.We used all RTK GPS marsh elevation points from the 20 m grid survey to calibrate LEAN.The 4 band NAIP imagery used to calculate NDVI was acquired 24 July 2015.To estimate model skill, we used a 10fold cross validation procedure in which the model was iteratively fit with 90% of the elevation data and tested against the 10% hold out.
We used a digitized map of GCReW vegetation zones to test applying cover-based elevation bias corrections.The map we used consists of the 6 major vegetation cover types: an Iva frutescensdominated community, I. frutescens and Phragmites australis mixed community, a Schoenoplectus americanus and Spartina patens mixed community, a P. australis-dominated community, a S. americanus-dominated community, and a S. patensdominated community.Community boundaries were delineated by trained personnel using highresolution 2010 summer image ry in Google Earth.Polygons were assigned one of 6 vegetation types using the interactive supervised classification tool in ArcGIS which uses the maximum likelihood method, trained using 26 total ground-based reference points that were collected in 2014.This map, created for a separate unpublished project, and the associated data release (Lu & Williams 2021) compares the vegetation communities in 2010 to those delineated from an image taken in 1973 (Jenkens 1973).
We used ANOVA to test if there were significant differences between measured and uncorrected mapped LiDAR-based elevation from the 20 m grid survey, using the mapped vegetation community in the 2010 map as an independent variable and LiDAR offset as the dependent variable.If the relationship was significant, then we could create a new DEM by subtracting those vegetation offsets from the mapped surface.
We performed empirical Bayesian kriging using ArcGIS Pro v.2.4 applied to all marsh surface data from the 20 m grid survey.We mapped both the median prediction and standard error of the kriged surface using a power semivariogram and a 100 maximum point standard circular neighborhood with a radius of 272 m.

Intercomparing unbiased elevation mapping strategies
We compared the 4 techniques by validating each of the 4 maps with the TMON elevation plot measurements, which were independent of all vegetation height correction strategies.We calculated the bias as mean error and precision as unbiased root mean square error (RMSE').We normalized both of these metrics (shown by *) to the standard deviation of the reference data set (σ r ).We also calculated the total RMSE, the sum of squares of mean error and RMSE'.Normalized total RMSE (RMSE*) provides a convenient performance threshold.When a strategy's RMSE* is less than 1, it performs better than simply applying the average from the reference data set (Jolliff et al. 2009).We then ranked the 4 strategies based on our interpretation of their scalability, meaning the feasibility of applying the strategy to different sites with a larger extent.

Vegetation survey
In conjunction with the grid-based elevation survey, we also recorded plant community presence, absence, and relative abundance.At each elevation plot, we estimated ordinal Braun-Blanquet scores for major vegetation community types (Braun-Blanquet 1932, Snedden & Steyer 2013).The area of the estimation was approximately a 1 m circumference circle around the center of the elevation point estimated as the length of the extended arm of the technician recording elevation and vegetation.Braun-Blanquet scores ranged from 0−5, with 0 indicating trace amounts of cover (<1% of the aerial view of the plot); 1 = 1−5%; 2 = 5−25%; 3 = 25−50%; 4 = 50−75%; and 5 = 75−100%.The major communities we recorded included S. americanus, the C 4 mix of Distichlis spicata and S. patens, the shrub I. frutescens, and the invasive reed P. australis.We recorded minor plant species as well.To enhance consistency, the same technician recorded all elevation and vegetation measurements.

Tidal datums for 2016
The full suite of tidal datums calculated for 2016 from the Annapolis tide gauge are listed in Table 1.The MSL was 0.092 m.Other tidal datums were as follows: MHW = 0.194 m NAVD88, MHHW = 0.278 m, and MHHWS = 0.340 m.Observed water levels ranged from −0.847 to 0.861 m.

Elevation distribution of the GCReW
Summary statistics of the elevation distribution measured in the 20 m grid survey are presented in tion centered around its median (0.221 m NAVD88), which skewed higher than its mean (Fig. 2).Overall, the distribution had a positive skew with a long tail on the lower end of the elevation gradient.
We classified the elevation survey points according to regularity of tidal flooding at different levels.We found that 37.1% of the marsh was below the MHW line (Fig. 2), receiving semidiurnal tides.We classified the majority of the marsh concentrated between the MHHW and MHW water lines (40.5%) as being inundated once per day.About one-fifth (19.7%) of the marsh was between the MHHWS and MHHW lines, being inundated once daily to twice a month.The remaining 2.7% was classified as above the MHHWS line.

Digital elevation mapping
We used our ground-based elevation data to test the performance and scalability of the 4 methods for producing an unbiased DEM.The mean error of the 2011 Anne Arundel County DEM at GCReW was 0.195 ± 0.114 (n = 567).For comparison, the weighted mean and SD of errors compiled from 19 sites across 3 studies (Hladik et al. 2013, Medeiros et al. 2015, Buffington et al. 2016) is 0.172 ± 0.206 (n = 19 349).
Because we detected bias (0.195 overestimation of elevation; Fig. 3A), we compared 4 techniques for producing an unbiased DEM product.The first method assumed a constant bias of 0.172 m estimated from previously collected data, and subtracted it from the LiDAR-based DEMs elevation (Fig. 3B).This simple calculation required no new data collection, and it provided a base for evaluating more complex correction methods.
Elevations corrected with the LEAN algorithm (Fig. 3C) had an RMSE of 0.081 m (estimated by 10-fold cross validation), a 28.8% improvement over the original LiDAR measured root mean square error (RMSE).Mean error from cross validation was −1.69 × 10 −5 ± 0.014 m.
A third correction method estimated separate correction factors for each of 6 vegetation communities.This independent and previously unpublished vegetation map included C 3 sedges, C 4 grasses, mixes between C 3 and C 4 communities, Phragmites australisdominated communities, Iva frutescens-dominated communities, and an Iva−Phragmites mixed community.The observed elevation bias was significantly different among the 6 communities (ANOVA, p < 0.0001).The largest LiDAR elevation biases were for the C 3 sedges (0.251 m; Fig. 4) and P. australis (0.237).The lowest bias was for the C 4 grass communities (0.120 m).We used the community-specific observed biases to produce a bias-corrected map (Fig. 3D).
Finally, we tested the more intensive technique, interpolating elevation among surveyed points using empirical Bayesian kriging (Fig. 3E).The standard error of the interpolated surface ranged from 0.0196 to 0.175 m.

Elevation mapping method comparison
The biases in the resulting post-processed DEMs were markedly lower than the original bias and ranged from 0.009 to −0.004 m NAVD88 for the simple single offset and the kriged maps, respectively.The RMSE' and RMSE of the 4 models did show general improvement at the cost of scalability, with the RMSE' of the kriged map being lowest (0.065 m) followed by the LEAN-corrected DEM (0.071 m).Despite being more scalable, the LEAN model performed better than the site-specific vegetation correction.The constant bias correction produced the poorest results in terms of RMSE' and RMSE.Model performance statistics were normalized to σ r (represented by *; Fig. 5).LEAN, site-specific vegetation correction, and empirical Bayesian kriging all produced maps that performed better than applying the average of the reference data set, whereas the simple, single offset did not.The value of σ r was greater than that of the modeled data set (σ m ) for LEAN and the kriged map, which means that the models produced maps that are flatter than the original surface.
We decided based on our experience that the order from most scalable to least was the single average offset, the application of LEAN, the site-specific vegetation community mapping, and kriging of elevation points.We plotted these ranks versus RMSE* to evaluate scalability versus performance tradeoffs.

Species presence and dominance by elevation zone
We summarized the Braun-Blanquet cover data in 2 ways: by vegetation presence and by relative vegetation dominance.Overall, we documented 11 plant  2).When analyzing the plant coverage for relative dominance, 8 species were observed to be dominant or codominant in plots somewhere on the marsh.The C 4 grasses were most often dominant (30.9% of observations), with P. australis second, I. frutescens third, and S. americanus fourth (Table 2).Spartina cynosuroides, T. latifolia, K. virginica, and S. alterniflora were less prevalent but had instances of high relative dominance in isolated subhabitats on the marsh.
Along the tidal elevation gradient, the native graminoid species cluster around different elevations, with the C 4 grasses dominating less frequently flooded areas between the MHHW and MHHWS lines where it floods a few times a month (Fig. 6).Schoenoplectus americanus dominated the more frequently floo ded zone between the MHHW and MHW lines where there are daily tides.Iva frutescens and P. australis were both much more dominant in the low marsh than the high marsh and had wider elevation distribution in general.Iva frutescens had a median growing elevation similar to that of S. americanus, but it had a much wider elevation distribution and dominated at elevations lower in the tidal range.Phragmites australis also occupied a relatively wide range of tidal elevations, occurring at both higher and lower elevations than the native marsh vege tation and being more dominant in frequently flooded areas.The latter observation is based on interquartile ranges of the elevation distributions, with I. fru tes cens and P. australis having wider ranges in and dominance (0.100 and 0.104 m, respectively) compared to the other 2 major marsh building species (C 4 = 0.100 and 0.070 and S. americanus = 0.100 and 0.076 m).
Of the minor species dominant in less than 20 instances, there were a few notable vertical and hor- ) on the x-axis, and total RMSE as the sum of squares (circles).Each metric has been normalized to the standard deviation of the reference data set so that a RMSE* value of 1 represents a performance threshold that models need to be lower than to be considered efficacious.RMSE' is artificially signed by whether the SD of the modeled (m) or reference (r) data sets are larger.(B) RMSE* relative to ordinal values ranking strategies from most to least scalable izontal spatial trends.Spartina cynosuroides occurred along a wide flooding gradient and was clustered in a few locations, including on raised berms on the marsh edge and on a gradual slope on the marsh's landward edge.Typha latifolia was clustered in continuously floo d ed depressions adjacent to sloped sections of the marsh's landward edge.Spartina alterniflora was rare but was occasionally dominant lower in the tidal range than any other species we documented.

DISCUSSION
We performed an extensive survey of elevation and vegetation at the GCReW.Our comparison of multiple mapping strategies motivated new recommendations on issues of scaling in tidal wetland elevation mapping.The distribution of elevation across the entire site, the distribution of plant species with respect to elevation, and our field observations suggest several useful future studies linking the physiology and ecology of species to mechanisms of marsh collapse and expansion.

Elevation mapping techniques: tradeoffs between performance and scalability
Maps of elevation relative to tidal range are essential tools for applying process-based models of marsh resilience to SLR and identifying trends in vulnerability and resilience.However, interference of live vegetation with LiDAR introduces elevation bias in an ecosystem where plant communities are separated by centimeters.Previous studies have developed vegetation correction methods with different spatial scales and different foci.In this study, we compared more scalable techniques, such as the single average offset, to techniques that require in crea s ing amounts of sitespecific ecologic and geographic data.
One of the goals of this study was to compare multiple mapping strategies for tidal marsh DEM bias correction.All the techniques, except the constant bias correction, outperformed simply using the mean bias of the reference data, meaning that the methods were efficacious.While the choice of correction technique depends on the research questions being asked, this study highlights limitations that resear chers should consider with regard to using a constant bias correction.
For the 3 efficacious bias correction techniques, there was a general tradeoff between performance and scalability.Spatial interpolation using kriging had the best performance overall, followed by the LEAN-corrected DEM, and then by estimating separate corrections for different vegetation communities.The LEAN-corrected DEM provided balanced tradeoffs between performance and scalability.For instance, Buffington et al. (2016) recommend a minimum of ~120 points site −1 to calibrate the statistical model, provided the survey captures variation in NDVI and elevation.Freely available high-resolution imagery from the NAIP program can be used for developing the NDVI input.Because it is based on imagery, the LEAN method has the advantages of being less labor intensive, being less invasive to the site, and not needing detailed site-specific vegetation maps.
Although the LEAN algorithm and kriging had the best overall performances, comparing the variances of the modeled and reference data suggests that extra consideration and guidance are needed when utilizing these techniques.Both LEAN-and krigingproduced maps had less variance than the reference data, meaning the maps are flatter than the measured marsh surface.This type of error could have implications for certain types of process modeling in which microtopography is a factor (Wu et al. 2011, Frei et al. 2012).

Plant elevation distributions and ramifications for vulnerability and resilience to SLR
This study is the first at GCReW to quantify elevation niche spaces occupied by marsh-building species.The native marsh species are typically presented as 2 functional groups: C 4 grasses and C 3 sedges.The elevation distributions of these groups overlap by 92% but still partition across the tidal elevation gradient, with sedges dominating between the MHW and MHHW line and grasses dominating above the MHHWS line.Other species did not have the same level of dominance in specific elevation zones, but still exhibited spatial clustering.For example, Typha latifolia was dominant in pools adjacent to sloping marsh edges, and Spartina cynosuroides was clustered at the open water edge and along landward edges of the marsh.Such spatial patterns are not explained by the concept that elevation is the dominant determinant of species distribution in tidal marshes, but instead may have more to do with salinity, anoxia, and ecological interactions.
Although not dominant at our site, we did find isolated patches of S. alterniflora, common in other studies of Chesapeake Bay marshes (e.g.Monie Bay Component of Chesapeake Bay Maryland National Estuarine Research Reserves; NERRS 2019), relatively low in the tidal frame (Kearney et al. 1994).Future studies could investigate whether future SLR and salinity changes could affect the dominant cover types leading S. alterniflora to be more dominant, as it is at other sites on the US Atlantic coast (Crosby et al. 2017).
At GCReW, Phragmites australis and Iva frutescens were dominant outside of their expected elevational ranges, which is an interesting observation that has implications for marsh resilience to SLR and merits further study.Previous research has reported that P. australis invasion is typically limited to the high marsh and mid-marsh zones (Bertness et al. 2002, Chambers et al. 2003), but we found that P. australis is dominant in both the lower and higher elevation landward edge.At GCReW, P. australis increased its distribution 22-fold between 197822-fold between and 200722-fold between (McCormick et al. 2010)), and it continues to expand unrestrained by its position in the tidal frame.In tidal wetlands, P. australis expansion into lower elevation zones is thought to be constrained by salinity, sulfide, and inundation time (Chambers 1997, Chambers et al. 2003).We found that P. australis at GCReW dominates tidal wetlands beyond the physiological thresholds thought to limit its establishment and distribution.This suggests that P. australis may have a broader ecological range than previously reported and greater overlap with elevation niches of native plants at GCReW (Kettenring et al. 2012) Regardless of the mechanism, the dominance and expansion of P. australis into low elevations may have unanticipated benefits.Phragmites australis is an ecosystem engineer that fixes 2−3 times more carbon than native plant communities at GCReW under current or near-future conditions (Caplan et al. 2015).Such increases in productivity can raise ecosystem blue carbon storage (Davidson et al. 2018).Once established, P. australis attenuates storm surge waves as effectively as native communities (Leonard et al. 2002, Theuerkauf et al. 2017), and builds soils and raises surface elevation more than native plant communities (Windham & Lathrop 1999, Rooth et al. 2003).This unique combination of higher productivity, high carbon storage, similar levels of wave attenuation, and greater surface elevation gain helps to maintain or increase the resilience of tidal wetlands that may otherwise succumb to accelerating SLR.Historical rates of relative SLR are 3.67 mm yr −1 measured at the Annapolis tide gauge (NOAA 2020b).
The dominance of P. australis in high elevation zones at the marsh−forest ecotone may also influence the development of future tidal wetlands undergoing transgression.We observed P. australis growing underneath the forest canopy at the forest−marsh ecotone (Fig. 1).This is consistent with a finding in Delaware Bay estuary that 60% of forests subject to SLR are being replaced by P. australis (Smith 2013), and with reports of it invading newly forming tidal wetlands under ghost forests (Kirwan & Gedan 2019).Future research should consider the ramifications of P. australis encroachment on adjacent uplands for efforts to conserve native plants and animal habitat.
We found I. frutescens lower in the marsh system than previously reported.Iva frutescens has been characterized as a high marsh plant dominating the landward edge above the S. patens zone on the US Atlantic Coast, with reduced photosynthesis even in the absence of competitors when transplanted to the low marsh (Bertness et al. 1992).In Louisiana, I. frutescens has been observed occupying more frequently flooded regions than in Rhode Island, having overrun S. alterniflora in created marshes (Owens et al. 2007).GCReW also shows a peak of I. frutescens dominance above MHW, but a substantial portion of I. frutescens distribution is below the MHW line (Fig. 6).These locations tend to be adjacent to tidal creeks, where I. frutescens either co-occurs with P. australis or Schoenoplectus americanus or has little to no understory.Future studies should investigate why I. frutescens grows lower in GCReW's tidal range and determine if this observation is unique to GCReW.We have 3 hypotheses.First, I. frutescens may have more phenotypic plasticity or genotypic variation than previously recognized.Second, I. frutescens could truly be a high marsh plant, but sections of the marsh could be collapsing due to a combination of SLR and muskrat herbivory destroying the understory plant community.The shape of the I. fru tescens elevation distribution (Fig. 6) has a peak in the high marsh between the MHHW and MHHWS lines and between the peaks of the 2 native graminoids, but it also has a fat tail skewing towards the low marsh and subtidal elevation zones.We would expect to see this type of a skewed distribution if understory loss in some sections led to subsidence from loss of aerenchymous tissue and increased organic matter decay (Chambers et al. 2019).Third, there is growing evidence on the effects of facilitation and positive interactions in tidal marshes (Zhang & Shao 2013, He et al. 2013).Given that both Phragmites and I. frutescens are found outside their predicted niches, it is possible that positive interactions between these 2 species are expanding their ranges.Future studies could investigate how I. frutescens and its understory are responding to SLR and herbivory, how soil integrity interacts with changes in plant community composition (Wigand et al. 2018), and how soil integrity relates to susceptibility to loss.To forecast marsh collapse events and the extent of upland wetland transgressions (Kirwan & Gedan 2019, Gedan & Fernández-Pascual 2019), we need to develop a mechanistic understanding of how plant and plant−animal interactions change the distributions of plant communities and the feedbacks among biogeochemical processes that govern elevation change.Gathering this data will require new monitoring and experiments across the full elevation range of the marsh.

CONCLUSIONS
This study is the first extensive vegetation and elevation survey of the Kirkpatrick Marsh, which contains the GCReW, one of the most intensively studied wetlands in the world.In the literature, the site is typically discussed as a high marsh.Our study supports this designation because the surveyed elevation of the majority of marsh area fell between the MHW and MHHWS tidal datum.However, there is also a substantial portion of the marsh that is at low elevation and a skewed distribution towards zones that flood twice daily or are constantly inundated.We used the vegetation and elevation data we collected to compare 4 DEM strategies.We found 3 methods with promising results for performance and scalability, with the LEAN method providing a balance between scalability and performance.However, the better-performing strategies produced maps that were flatter than the survey data.We observed a diverse array of plant communities.The distributions of the most frequent communities were strongly partitioned along elevation gradients.The C 3 sedge community more frequently occupied parts of the marsh that flood daily while the C 4 grass community occupied parts that flood a few times per month.We documented a wide inundation range for the reed Phragmites australis, which may promote its success as an invasive species.Future research should further document and model interactions between flooding and plant distribution at the margins of the marsh and near adjacent forests on the landward side of the marsh.We particularly recommend more extensive monitoring of P. australis along the marsh−forest ecotone and more study of elevation and soil integrity where Iva frutescens occurs in frequently flooded areas.To conclude, this analysis adds value to past studies at the GCReW, will help plan future experiments for this site or other similar sites, and provides data on plant community distribution across elevation gradients to support cross-site synthesis and integration of monitoring data and/or model results.

Fig. 1 .
Fig. 1.Marsh extent and boundaries of existing experiments including C 3 and C 4 warming experiments, CO 2 enrichment crossed with community treatments, CO 2 crossed with nitrogen enrichment, and CO 2 crossed with nitrogen enrichment for Phragmites australis (Phrag.)at the Global Change Research Wetland, Maryland, USA, as well as the 20 × 20 m (n = 525) elevation and vegetation survey grid and the Tennenbaum Marine Observatory Network (TMON) plots (n = 37).Inset maps: site in relation to the state of Maryland and the USA.Basemap source: ESRI

Fig. 2 .
Fig. 2. Probability density distribution of elevation of the Global Change Research Wetland in 2016 (n = 525).Elevation is relative to the North American Vertical Datum of 1988 (NAVD88) measured during the 20 × 20 m elevation survey with the relationship to tidal datums from 2016: MLLW: mean lower low water; MLW: mean low water; MSL: mean sea level; MHW: mean high water; MHHW: mean higher high water; MHHWS: mean higher high water spring; HAT: highest astronomical tide

Fig. 3 .
Fig. 3. Elevation (in m) relative to North American Vertical Datum of 1988 (NAVD88) for the (A) uncorrected digital elevation model (DEM) and maps showing an unbiased DEM for various techniques including (B) using constant bias adjustment, (C) light detection and ranging (LiDAR) Elevation Correction with Normalized Difference Vegetation Index (LEAN) results, (D) vegetation zone corrected DEM, and (E) kriged DEM, built independent from LiDAR.The elevation color scheme is maximized to show differences among modeled tidal flooding regimes.In addition to average site elevation, the maps reveal variability in elevation and the presence or absence of a raised marsh edge in the north of the image

Fig. 5 .
Fig. 5. Performance metrics for unbiased digital elevation mapping strategies for the Global Change Research Wetland in 2016 (n = 148) including the uncorrected digital elevation model (Simple Offset), the light detection and ranging (LiDAR) Elevation Correction with Normalized Difference Vegetation Index (LEAN), the vegetation zone corrected DEM (Veg.Class), and the kriged DEM, built independent from LiDAR.(A) Target diagram showing bias on the y-axis, unbiased root mean square error (RMSE') on the x-axis, and total RMSE as the sum of squares (circles).Each metric has been normalized to the standard deviation of the reference data set so that a RMSE* value of 1 represents a performance threshold that models need to be lower than to be considered efficacious.RMSE' is artificially signed by whether the SD of the modeled (m) or reference (r) data sets are larger.(B) RMSE* relative to ordinal values ranking strategies from most to least scalable

Fig. 6 .
Fig. 6.Probability density distribution of elevation relative to the North American Vertical Datum of 1988 (NAVD88) at the Global Change Research Wetland in 2016, measured across the 4 most dominant vegetation types: C 4 grasses (C4), Iva frutescens (IVFR), Phragmites australis (PHAU), and Schoenoplectus americanus (SCAM), with the relationship to tidal datums in 2016: MLLW: mean lower low water; MLW: mean low water; MSL: mean sea level; MHW: mean high water; MHHW: mean higher high water; MHHWS: mean higher high water spring; HAT: highest astronomical tide

Table 2 .
The elevation gradient surveyed at GCReW ranged from −0.247 to 0.444 m.Mean (± SD) elevation was 0.209 ± 0.092 m.The distribution of eleva-

Table 1 .
Tidal datums at the 2016 Annapolis tide gauge (NAVD88, North American Vertical Datum of 1988)