Coastal wetlands are likely to lose productivity under increasing rates of sea‐level rise (SLR). This study assessed a fluvial estuarine salt marsh system using the Hydro‐MEM model under four SLR scenarios. The Hydro‐MEM model was developed to apply the dynamics of SLR as well as capture the effects associated with the rate of SLR in the simulation. Additionally, the model uses constants derived from a 2‐year bioassay in the Apalachicola marsh system. In order to increase accuracy, the lidar‐based marsh platform topography was adjusted using Real Time Kinematic survey data. A river inflow boundary condition was also imposed to simulate freshwater flows from the watershed. The biomass density results produced by the Hydro‐MEM model were validated with satellite imagery. The results of the Hydro‐MEM simulations showed greater variation of water levels in the low (20 cm) and intermediate‐low (50 cm) SLR scenarios and lower variation with an extended bay under higher SLR scenarios. The low SLR scenario increased biomass density in some regions and created a more uniform marsh platform in others. Under intermediate‐low SLR scenario, more flooded area and lower marsh productivity were projected. Higher SLR scenarios resulted in complete inundation of marsh areas with fringe migration of wetlands to higher land. This study demonstrated the capability of Hydro‐MEM model to simulate coupled physical/biological processes across a large estuarine system with the ability to project marsh migration regions and produce results that can aid in coastal resource management, monitoring, and restoration efforts.

1 Introduction The fate of coastal wetlands may be in danger due to climate change and sea‐level rise (SLR) in particular. Identifying and investigating factors that influence the productivity of coastal wetlands may provide insight into the potential future salt marsh landscapes and to identify tipping points [Nicholls et al., 1999; Nicholls, 2004]. Salt marsh systems play an important role in coastal protection by attenuating waves and providing shelter and habitat for various species [Daiber, 1977; Halpin, 2000; Moller et al., 2014]. A better understanding of salt marsh evolution under SLR supports more effective coastal restoration, planning, and management [Bakker et al., 1993]. Plausible projections of global SLR, including its rate, are critical to effectively analyze coastal vulnerability [Parris et al., 2012]. In addition to increasing local mean tidal elevations, SLR alters hydrologic circulation patterns and sediment transport, which can affect the ecosystem as a whole and wetlands in particular [Nichols, 1989]. Studies have shown that adopting a dynamic modeling approach is preferred over static modeling when conducting coastal vulnerability assessments under SLR scenarios. Dynamic modeling includes nonlinearities that are unaccounted for by simply increasing water levels. The static or “bathtub” approach simply elevates the present‐day water surface by the amount of SLR and projects new inundation using a digital elevation model (DEM). On the other hand, a dynamic approach incorporates the various nonlinear feedbacks in the system and considers: interactions between topography and inundation that can lead to an increase or decrease in tidal amplitudes and peak storm surge; changes to tidal phases and timing of maximum storm surge; and modifications of depth‐averaged velocities (magnitude and direction) [Hagen and Bacopoulos, 2012; Atkinson et al., 2013; Bilskie et al., 2014; Holleman and Stacey, 2014; Cheng et al., 2015; Passeri et al., 2015; Bilskie et al., 2016b; Passeri et al., 2016]. Previous research has investigated the effects of sea‐level change on hydrodynamics and coastal wetland changes using a variety of modeling tools [Wolanski and Chappell, 1996; Liu, 1997; Hearn and Atkinson, 2001; French, 2008; Leorri et al., 2011; Hagen and Bacopoulos, 2012; Hagen et al., 2013; Valentim et al., 2013]. Several integrated biological models have been developed to assess the effect of SLR on coastal wetland changes. They employed hydrodynamic models in conjunction with marsh models to capture the changes in marsh productivity as a result of variations in hydrodynamics [D' Alpaos et al., 2007; Kirwan and Murray, 2007; Temmerman et al., 2007; Hagen et al., 2013]. However, these models were designed for small marsh systems or overly simplified complex processes. Alizad et al. [2016] coupled a hydrodynamic model with the Marsh Equilibrium Model (MEM) to include the biological feedback in a time stepping framework that incorporates the dynamic interconnection between hydrodynamics and a marsh system by updating elevation and bottom roughness in each time step. Lidar‐derived DEMs are a generally accepted means to generate accurate topographic surfaces over large areas [Medeiros et al., 2011; Bilskie and Hagen, 2013]. While the ability to cover vast geographic regions at relatively low cost makes lidar an attractive proposition, there are well‐documented topographic errors, especially in coastal marshes, when compared to Real Time Kinematic (RTK) topographic survey data [Hsing‐Chung et al., 2004; Hladik and Alber, 2012; Medeiros et al., 2015]. The accuracy of lidar DEMs in salt marsh systems is limited primarily because of the inability of the laser to penetrate through the dense vegetation to the true marsh surface [Hladik and Alber, 2012]. Filters, such as slope‐based or photogrammetric techniques applied in post processing, are able to reduce but not eliminate these errors [Kraus and Pfeifer, 1998; Lane, 2000; Vosselman, 2000; Carter et al., 2001; Hicks et al., 2002; Sithole and Vosselman, 2004; James et al., 2006]. The amount of adjustment required to correct the lidar‐derived marsh DEM varies on the marsh system, its location, the season of data collection, and instrumentation [Montane and Torres, 2006; Yang et al., 2008; Wang et al., 2009; Chassereau et al., 2011; Hladik and Alber, 2012]. In this study, the lidar‐derived marsh table elevation is adjusted based on RTK topographic measurements and remotely sensed data [Medeiros et al., 2015]. Sediment deposition is a critical component in sustaining marsh habitats [Morris et al., 2002]. The most important factor that promotes sedimentation is the existence of vegetation that increases residence time within the marsh allowing sediment to settle out [Fagherazzi et al., 2012]. Research has shown that salt marshes maintain their state (i.e., equilibrium) under SLR by accreting sediments and organic materials [Reed, 1995; Turner et al., 2000; Morris et al., 2002; Baustian et al., 2012]. Salt marshes need these two sources of accretion (sediment and organic materials) to survive in place [Nyman et al., 2006; Baustian et al., 2012] or to migrate to higher land [Warren and Niering, 1993; Elsey‐Quirk et al., 2011]. The Apalachicola River (Figure 1), situated in the Florida Panhandle, is formed by the confluence of the Chattahoochee and Flint Rivers and has the largest volumetric discharge of any river in Florida, and drains the second largest watershed in Florida [Isphording, 1985]. Its wide and shallow microtidal estuary is centered on Apalachicola Bay in the northeastern Gulf of Mexico (GOM). Apalachicola Bay is bordered by East Bay to the northeast, St. Vincent Sound to the west, and St. George Sound to the east. The river feeds into an array of salt marsh systems, tidal flats, oyster bars, and submerged aquatic vegetation as it empties into Apalachicola Bay, which has an area of 260 km2 and a mean depth of 2.2 m [Mortazavi et al., 2000]. Offshore, the barrier islands (St. Vincent, Little St. George, St. George and Dog Islands) shelter the bay from the GOM. Approximately 17% of the estuary is comprised of marsh systems that provide habitats for many species, including birds, crabs, and fish [Halpin, 2000; Pennings and Bertness, 2001]. Approximately 90% of Florida's annual oyster catch, which amounts to 10% of the Nation's, comes from Apalachicola Bay. In Franklin County, the home to the bay, and 65% of workers are or have been employed in the commercial fishing industry [FDEP, 2013]. Therefore, accurate assessments regarding this ecosystem can provide insights to environmental and economic management decisions. Figure 1 Open in figure viewer PowerPoint Study area and marsh organ locations. (a) Location of the Apalachicola estuarine system in Florida. (b) The Apalachicola River, Bay, and other locations including the transects for assessing velocity variations in the estuarine system. (c) The marsh organ experimental sets in the estuarine system. It is projected that SLR may shift tidal boundaries upstream in the river, which may alter inundation patterns and remove coastal vegetation or change its spatial distribution [Florida Oceans and Coastal Council, 2009; Passeri et al., 2016]. Apalachicola salt marshes may lose productivity with increasing SLR due to the microtidal nature of the estuary [Livingston, 1984]. The primary objective of this study is to assess the response of the Apalachicola salt marsh for various SLR scenarios using a high‐resolution Hydro‐MEM model for the region. The Hydro‐MEM model was developed and previously tested for a mesotidal estuary in the Timucuan salt marsh system (northeast Florida) [Alizad et al., 2016]. The Timucuan marsh is a relatively small system (∼60 km2) and is not affected by major river flow influence. Migration of the Timucuan salt marsh was found to be minimal under SLR due to the hard infrastructure along the river banks and some of the outer extents of the system. This research builds upon the work of Alizad et al. [2016] and examines how the Hydro‐MEM model performs in a microtidal marsh system, such as the Apalachicola marsh, and how the marsh may migrate upland. Due to the nature of the microtidal environment and the riverine influence from the Apalachicola River, it is critical to ensure the marsh platform elevations used in Hydro‐MEM are accurate. Therefore, the marsh platform elevations were ground‐truthed and corrected based on RTK‐GPS (Global Positioning System) surveys. Furthermore, this model incorporates river inflow (in addition to astronomic tides) to better capture the hydrodynamics of the system which leads to a more accurate simulation of biomass density. The results obtained include inundation area, a new comparison for the Hydro‐MEM model using remotely sensed data, and the demonstration of potential marsh migration by the year 2100 for four SLR scenarios. To our knowledge, no other salt marsh model has presented a capability to dynamically simulate marsh migration. Finally, we present biomass productivity results for the Apalachicola region for the years 2020, 2050, 2080, and 2100 under the four SLR scenarios to illuminate how and when the marsh loses productivity and converts to open water. In summary, this paper applies the Hydro‐MEM model in a microtidal estuary by adjusting the marsh platform and incorporates river inflow to address (1) inundated area extent, (2) marsh productivity change, and (3) marsh migration under high SLR scenarios.

2 Methods 2.1 Hydro‐MEM The integrated Hydro‐MEM model was used to assess the response of the Apalachicola salt marsh under SLR scenarios. The hydrodynamic and biologic components are coupled and exchange information at discrete coupling intervals (time steps). Incorporating multiple feedback points along the simulation timeline via a coupling time step permits a nonlinear rate of SLR to be included in the model. This contrasts with other techniques that apply the entire SLR amount in one single step. This allows for a more realistic integration of the physical and biological components associated with the hydrodynamics of the salt marsh systems over time. The coupling time step, considering the desired level of accuracy and computational expense, depends on the rate of SLR and governs the frequency of information exchange between the hydrodynamic and biological models [Alizad et al., 2016]. The model was initialized from the present‐day marsh surface elevations and sea level. First, the hydrodynamic model computes water levels and depth‐averaged currents and yields astronomic tidal constituent amplitudes and phases. From these results, mean high water (MHW) was computed and passed to the parametric marsh model, and along with field/laboratory analyzed biomass curve parameters, the spatial distributions of biomass density and accretion were calculated. Accretion was applied to the marsh elevations and the bottom friction was updated based on the biomass distribution and passed back to the hydrodynamic model to initialize the next time step. Additional details in the Hydro‐MEM model can be found in Alizad et al. [2016]. The next two sections provide specific details to each portion of the Hydro‐MEM model framework, the hydrodynamic model and the marsh model. 2.1.1 Hydrodynamic Model Hydrodynamic simulations were performed using the two‐dimensional, depth‐integrated, ADvanced CIRCulation (ADCIRC) finite element‐based shallow water equations model to solve for water levels and depth‐averaged currents. An unstructured mesh for the Apalachicola estuary was developed with focus on simulating water level variations within the river, tidal creeks, and daily wetting and drying across the marsh surface. The mesh was constructed from manual digitization using recent aerial imagery of the river, distributaries, tidal creeks, estuarine impoundments, and intertidal zones. To facilitate numerical stability of the model, with respect to wetting and drying, the number of triangular elements within the creeks was restricted to three or more to represent a trapezoidal cross‐section [Medeiros and Hagen, 2013]. Therefore, the width of the smallest captured creek was approximately 40 m. For smaller creeks, the computational nodes located inside the digitized creek banks were assigned a lower bottom friction (that is equal to creeks value) to allow the water to flow more readily, thus keeping lower resistance of the creek (compared to marsh grass vegetation). The mesh extends to the 60° west meridian (the location of the tidal boundary forcing) in the western North Atlantic and encompasses the Caribbean Sea and the GOM [Hagen et al., 2006]. The triangular elements range from a minimum of 15 m within the creeks and marsh platform to 300 m in Apalachicola Bay to 136 km in the open ocean. The source data for topography and bathymetry consisted of the lidar‐derived DEM provided by the Northwest Florida Water Management District (http://www.nwfwmdlidar.com/) and Apalachicola River surveyed bathymetric data from the U.S. Army Corps of Engineers, Mobile District. Medeiros et al. [2015] showed that the lidar topographic data for this salt marsh contained high bias and required correction to increase the accuracy of the marsh surface elevation and facilitate wetting of the marsh platform during normal tidal cycles. The DEM was adjusted using biomass density estimated by remote sensing; however, a further adjustment was also implemented. The biomass adjusted DEM (Figure 2b) reduced the bias of the lidar‐derived marsh platform elevation (Figure 2a) from 0.65 to 0.40 m when compared against a sample of 229 RTK spot elevations acquired by the authors in this marsh. This remaining 0.40 m of bias was removed by further lowering the DEM at the southeastern salt marsh shoreline where the vegetation is densest and linearly decreasing the adjustment value to zero moving upriver (i.e., to the northwest). Figures 2a and 2b show the uncorrected and corrected topography of the marsh system. Figure 2a illustrates an example of why this method was necessary in order to properly capture the cyclical tidal flooding of the marsh platform; without removing this remaining 0.40 m bias, the majority of the platform was incorrectly specified to be above MHW and was unable to become wet in the model during high tide. This would have resulted in an unrealistic wetting and drying regime for the marsh, with consequences in sediment supply and biomass density prediction as well. Figure 2 Open in figure viewer PowerPoint Topographic model input of the Apalachicola estuarine system and the elevation change along the transect “T” shown in (a). Color bar elevations are referenced to NAVD88 in meters. (a) Lidar data elevation without any correction; (b) adjusted marsh platform elevation. Bottom friction was included in the model as spatially varying Manning's n and were assigned based on the National Land‐Cover Database 2001 [Homer et al., 2004; Bilskie et al., 2015]. Three values for low, medium, and high biomass densities in each land cover class were defined as 0.035, 0.05, and 0.07, respectively [Arcement and Schneider, 1989; Medeiros, 2012; Medeiros et al., 2012]. At each coupling time step, using the computed biomass density and hydrodynamics, Manning's n values for specific computational nodes were converted to open water (permanently submerged due to SLR) or if their biomass density changed. In this study, the four global SLR scenarios for the year 2100 presented by Parris et al. [2012] were applied; low (0.2 m), intermediate‐low (0.5 m), intermediate‐high (1.2 m) and high (2.0 m). The coupling time step was 10 years for the low and intermediate‐low SLR scenarios and 5 years for the intermediate‐high and high SLR scenarios. This protocol effectively discretizes the SLR curves into linear segments that adequately capture the projected SLR acceleration for the purposes of this study [Alizad et al., 2016]. Thus, the total number of hydrodynamic simulations performed in this study are 60: 10 (100 years divided by 10 coupling steps) for the low and intermediate‐low scenarios and 20 (100 years divided by 20 coupling steps) for the intermediate‐high and high scenarios. The hydrodynamic model is forced with astronomic tides at the 60° meridian (open ocean boundary of the Western North Atlantic Tidal model domain) and river inflow for the Apalachicola River. The tidal forcing is comprised of time varying water surface elevations of the seven principal tidal constituents (M 2 , S 2 , N 2 , K 1 , O 1 , K 2 , and Q 1 ) [Egbert et al., 1994; Egbert and Erofeeva, 2002]. The river inflow boundary forcing was obtained from the United States Geological Survey (USGS) gage near Sumatra, FL in the Apalachicola River (USGS 02359170). The mean discharge value from 38 years of record for this gage (http://waterdata.usgs.gov/nwis/uv?02359170) was 670.17 m3/s as of February 2016. The discharge at the upper (northern) extent of the Apalachicola River is controlled at the Jim Woodruff Dam near the Florida—Georgia border and was considered to be fixed for all simulations. Two separate hyperbolic tangent ramping functions were applied at the model start, one for tidal forcing and the other for the river inflow boundary condition. The main outputs from the hydrodynamic model used in this study were the amplitude and phase harmonic tidal constituents, which were then used as input to the Hydro‐MEM. This hydrodynamic model has been extensively validated for astronomic tides in this region [Bilskie et al., 2016a; Passeri et al., 2016]. 2.1.2 Marsh Model B) (g m−2 yr−1) using a parabolic curve [Morris et al., 2002 (1) The parametric marsh model employed was the Marsh Equilibrium Model (MEM), which quantifies biomass density () (g myr) using a parabolic curve []: The biomass density curve includes composited data for three different marsh grass species of Spartina cynosuroides, Juncus romerianus, and Spartina alterniflora. The curve was divided into left and right branches that meet at the maximum biomass density point, resulting in an asymmetric biomass response curve. The left and right curve coefficients used were a l = 197.5 g m−3 yr−1, b l = −9870 g m−4 yr−1, c l = 1999 g m−2 yr−1, and a r = 326.5 g m−3 yr−1, b r = −1633 g m−4 yr−1, c r = 1998 g m−2 yr−1, respectively. They were derived using field bioassay experiments, commonly referred to as “marsh organs,” which consist of planted marsh species in an array of PVC (Polyvinyl Chloride) pipes cut to different elevations in the tidal frame [Morris, 2007]. The pipes each contained approximately the same amount of biomass at the start of the experiment and were allowed to grow based on their natural response to the hydroperiod [Morris et al., 2013]. They were installed at three sites in the Apalachicola estuary (Figure 1c) in 2011 and inspected regularly for 2 years by qualified staff from the Apalachicola National Estuarine Research Reserve (ANERR). Morris et al., 2002 Y) (cm/yr) during the study time (dt) was calculated by incorporating the amount of inorganic accumulation from sediment load (q) and the vegetation effect on organic and inorganic accretion (k): (2) The accretion rate in the coupled model included the accumulation of both organic and inorganic materials and was based on the MEM accretion rate formula denoted by []. The total accretion (d) (cm/yr) during the study time (d) was calculated by incorporating the amount of inorganic accumulation from sediment load () and the vegetation effect on organic and inorganic accretion (): where the parameters q (0.0018 yr−1) and k (1.5 × 10−5 g−1 m2) [Morris et al., 2002] were used to calculate the total accretion at each computational point and update the DEM at each coupling time step. The accretion was calculated for the points with positive relative depth, where average mean high tide is above the marsh platform [Morris, 2007].

3 Results 3.1 Hydrodynamic Results The MHW results for future SLR scenarios predictably showed higher water levels and altered inundation patterns (Figure 3). The warmer colors in Figure 3 indicate higher water levels; however, please note that the range of each scale bar was adjusted in order to capture the variation in water level for each scenario. The water level change in the creeks was spatially variable under the low and intermediate‐low SLR scenario, which shows the nonlinear response to SLR and the water level increase is different than the amount of SLR (Figures 3a–3c). However, the water level in the higher scenarios was more spatially uniform (Figures 3d and 3e). In the low SLR scenario, the water level changed from 10 to 40 cm in the river and creeks, but the wetted area remained similar to the current condition (Figures 3a and 3b). The water level and wetted area increased under the intermediate‐low SLR scenario and some of the forested area became inundated (Figure 3c). Under the intermediate‐high and high SLR scenarios, all of the wetlands were inundated, water level increased by more than a meter, and the bay conditions extended to the uplands (Figures 3d and 3e). Figure 3 Open in figure viewer PowerPoint Mean high water results for projected SLR scenarios for 2100 for (a) current conditions; (b) low sea‐level rise (SLR) scenario (0.2 m); (c) intermediate‐low SLR scenario (0.5 m); (d) intermediate‐high SLR scenario (1.2 m); and (e) high SLR scenario (2 m). Note that the scale bar colors are not consistent. The table in Figure 3 shows the inundated area for each SLR scenario in the Apalachicola region including the bay, rivers, and creeks. Under the intermediate‐low SLR, the wetted area increased by 12 km2, an increase by 2%. However, the flooded area for the intermediate‐high and high SLR scenarios drastically increased to 255 and 387 km2, which is a 45% and 68% increase in the wetted area, respectively. These results dynamically affect biomass density and salt marsh migration. 3.2 Biomass Density Biomass density is a function of both topography and MHW and is herein categorized as low, medium, and high. Figure 4b shows the biomass density map for the current condition in the lower part of the East Bay islands. The black boxes in the map (Figures 4a and 4b) draw attention to areas of interest for comparison with satellite imagery [Medeiros et al., 2015]. The three boxes qualitatively illustrate that the model reasonably captured the spatial patterns of low, medium, and high productivity in these areas. These areas were selected because their extents are relatively consistent in both the satellite‐derived and model predicted biomass density images. A close inspection of Figure 4 reveals that some land areas shown in the satellite‐derived biomass density image (Figure 4a) is covered with water (blue color) in Figure 4b. This is an artifact of the land cover‐based masking procedure used in Medeiros et al. [2015]. The simulated categorized (low, medium, and high) biomass density (Figure 4b) was compared with the biomass density derived from satellite imagery over both the entire satellite imagery coverage area (Figure 4a, n = 54,686 pixels) as well as over the highlighted areas (n = 6411 pixels). The confusion matrices for these two sets of predictions are shown in Table 1. Figure 4 Open in figure viewer PowerPoint Medeiros et al., 2015 Biomass density maps focused on the wetland area in the islands between the Apalachicola River and the East Bay classified into regions of low, medium, and high biomass represented by red, yellow, and green, respectively and blue shows the wet regions. (a) The remotely sensed biomass densities [] with black boxes highlights three selected regions for comparison. (b) Biomass density predicted by Hydro‐MEM under current conditions and the black boxes for comparison. Table 1. Confusion Matrices for Current Biomass Density Levels (Number of Pixels) Based on Remotely Sensed and the Hydro‐Marsh Equilibrium Model Predictions Predicted True Entire Coverage Area Highlighted Areas Low Medium High Total Low Medium High Total Low 4359 7112 4789 16,260 970 779 576 2325 Medium 1764 10,331 7942 20,037 308 809 647 1764 High 2129 8994 7266 18,389 396 728 1198 2322 Total 8252 26,437 19,997 54,686 1674 2316 2421 6411 These confusion matrices indicate a true positive rate over the entire coverage area of 26.8%, 51.5%, and 39.5% for low, medium, and high, respectively, and an overall weighted true positive rate of 40.12%. For comparison, a neutral model [Gardner and Urban, 2007] where the biomass category was assigned randomly (stratified to match the distribution of the true data derived from satellite imagery) produced true positive rates of 30.5%, 36.8%, and 32.8% for low, medium, and high, respectively, and an overall weighted true positive rate of 33.6%. For the highlighted areas, the true positive rates were 41.7%, 45.8%, and 51.6% for low, medium, and high, respectively, and an overall weighted true positive rate of 46.4%. To show the riverine discharge effect on biomass productivity, biomass density derived from a tides‐only scenario was compared to simulated biomass derived from a simulation that included tides and river discharge. Under the tides‐only scenario, much of the forested land (the islands between the Apalachicola River and East Bay and extending upstream to the area near Brothers River (Figure 1b)) around the river and creeks was incorrectly shown to be covered by salt marsh (Figure 5b). By including river inflow, the marsh area was correctly represented (Figure 5a) as well as the forested region (shown as no marsh area), which was not correct when not including river discharge. Figure 5 Open in figure viewer PowerPoint River inflow effect on biomass productivity for current conditions. (a) Simulated biomass density with tidal and river discharge forcings and (b) simulated biomass density with only tidal forcing. Temporal changes in the biomass density are demonstrated in the four columns (a–d) of Figure 6 for the years 2020, 2050, 2080, and 2100. The rows from top to bottom show simulated biomass density for the low, intermediate‐low, intermediate‐high, and high SLR scenarios, respectively. Under the low SLR scenario (top row), both the spatial distribution of vegetated marsh and its productivity show almost no change from the present time to year 2020 (Figure 6a) and an expansion of vegetated habitat marsh into upland by 2080, but from 2080 the trend reversed. However, for the higher SLR scenarios (second to fourth row from top), the medium and low biomass density became more dominant and the distribution of open water generally increased over time at the expense of vegetated marshland (Figures 6b–6d). Figure 6 Open in figure viewer PowerPoint Temporal changes in biomass density under future sea‐level rise (SLR) scenarios. Biomass density is categorized into low, medium, and high productive regions represented by red, yellow, and green, respectively, blue shows open water, and unclassified upland area is tan. The pie charts are based on the number of computational nodes. The columns (a–d) show biomass density for the years 2020, 2050, 2080, and 2100, respectively, and the rows from top to bottom display the results for the low (0.2 m), intermediate‐low (0.5 m), intermediate‐high (1.2 m), and high (2 m) SLR. The sequence of change in the marsh landscape is complex. For example, in the intermediate‐low SLR scenarios (second row from top), there was first a trend of marsh loss by the year 2050 (Figures 6a and 6b), and then by the year 2080 the medium and low productivity marsh regions was produced, but from 2080 to 2100 salt marsh areas lost productivity. The loss of marsh habitat by year 2080 was most significant in the mid‐high and high SLR scenarios (third and fourth rows from top). In both cases, vegetated marsh area was largely replaced by open water by 2100 (Figure 6d). Inland marsh migration was forecasted to compensate for drowned marsh at times in some scenarios. Marshland lost productivity in some areas while other, formerly upland areas gained productivity. Productivity increased in regions in the lower part of the islands between the Apalachicola River and East Bay, while most of the other marsh areas lost productivity or converted to open water near Lake Wimico. Under the intermediate‐low SLR scenario (second row of Figure 6d), the upper part of the islands became flooded and most of the salt marshes lost productivity while some migrated to higher lands. Salt marsh migration was more evident in the intermediate‐high and high SLR scenarios (third and fourth row of Figure 6d). The productive band around the extended bay under higher scenarios implied the possibility for productive wetlands in those regions. It is shown in the intermediate‐high and high‐SLR scenarios (third and fourth row) from 2020 to 2100 (column [a] to [d]) that the drowned marsh initiated around regions of open water and extended into low productivity areas. Under higher SLR, marsh migration extended to its topographic limits.

4 Discussion The salt marsh response to SLR is dynamic and strongly depends on both MHW and geomorphology of the marsh system. The hydrodynamic results for the SLR scenarios are a function of the rate of SLR and its magnitude, the subsequent topographic changes resulting from marsh platform accretion, and the change in flow resistance induced by variations in biomass density. This was demonstrated in the low and intermediate‐low SLR scenarios, where the water level varied in the creeks, on the marsh platform and where the flooding started from no marsh productivity regions (Figures 3b–3c). In the intermediate‐high and high SLR scenarios, the overwhelming extent of inundation dampened the impact of topography and flow resistance and the new hydrodynamic patterns were generally dependent on SLR magnitude. In both extended and highlighted areas of landscape (Figure 4), as shown in Table 1, differentiating between medium and high biomass classifications was problematic. This may be attributed to saturation in satellite imagery based coverage where at high levels of “greenness” and canopy height, the satellite imagery can no longer differentiate subtle differences in biomass density, especially at fine spatial resolutions. As such, the incorrectly classified cells are primarily located throughout the high‐resolution part of the model domain where the spacing is 15 m on average. This error may be mitigated by coarsening the resolution of the biomass density predictions using a spatial filter; predictions at this fine of a resolution are admittedly ambitious. This would reduce the “speckled” appearance of the predictions and likely lead to better results. Ultimately, both methods of biomass prediction rely on imperfect models. Nevertheless, it is encouraging that with the adjusted marsh platform and river inflow forcing in the hydrodynamic simulations, the model results demonstrated good agreement with the remotely sensed biomass data (Figures 4a and 4b). If we focus on the three black boxes highlighted in the Figures 4a and 4b, in the first box from left, the model predicted the topographically lower lands as having low productivity whereas the higher lands were predicted to have medium productivity. Higher lands near the bank of the river (middle box in the Figure 4a) were correctly predicted as a low productivity (Figure 4b). The right box in Figure 4a also depicts areas of both high and mixed productivity patterns. The model also predicted the vast expanse of area with no productivity. In the highlighted areas, the agreement between biomass classification by the Hydro‐MEM and by remote imagery was 46.4% of the cells being correctly identified. This indicates the promise of the method to capture the biomass density distribution in key areas. However, as discussed above, both methods of classifying biomass rely on imperfect models, and possibly the best metric of agreement is the similarity in total vegetated area. Other sources of error that decrease the prediction accuracy likely originate from the remotely sensed, experimental, and elevation data. The remotely sensed biomass density is primarily based on Interferometric Synthetic Aperture Radar (IfSAR) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data which have inherent uncertainty [Andersen et al., 2005]. In addition, the remotely sensed biomass density estimation was constrained to areas classified as Emergent Herbaceous Wetland (95) by the National Land Cover Dataset 2011. This explains the difference in coverage area and may have excluded areas where the accuracy of the proposed biomass density prediction model was better. Other errors are likely that originate from variability in planting, i.e., survivorship and growth unrelated to the elevation that is minimized, but unavoidable, in any field experiment of this nature. These data were used to train the satellite imagery based biomass density prediction method and errors in their computation, while minimized by sample replication and outlier removal, may have propagated into the coverage used as the ground truth for comparison (Figure 4a). Another source of error occurs in the topography adjustment process in the marsh system. The true marsh topography is highly nonlinear and any adjustment algorithm cannot include all of the nonlinearities and microtopographic changes in the system. Considering all of these factors and the capability of the model to capture the low, medium, and high productive regions qualitatively as well as the difficulty to model the complicated processes within salt marsh systems, the Hydro‐MEM was successful in computing marsh biomass density. The tides‐only scenario demonstrated drastically lower water levels in the river as well as its tributaries and distributaries due to the lack of inflow from upstream. The water level in the lower creeks branched from East Bay remained the same as the base condition because tidal forcing dominates water level in this area. As a result, the spatial distribution of the biomass density on the islands and near the lower creeks followed a pattern similar to the real case scenario (Figures 5a and 5b). However, the results of tides‐only and the real case scenarios differed in the upper part of the islands and areas near the river due to differences in computed MHW. It should be noted that our simulation used a yearly average for inflow boundary condition assuming a constant flow rate controlled by Jim Woodruff Dam, which is located upstream. In order to improve the predictions and correctly identify additional regions with limited or no salt marsh productivity, future versions of the Hydro‐MEM will include inter‐annual variability in flow induced by the effects of climate change on extreme events [Wang et al., 2013; Cheng et al., 2015; Hovenga et al., 2016] and their complex impact on geomorphology, salinity, and sediment transport on wetlands productivity [Cahoon, 2006; Greening et al., 2006; Turner et al., 2006]. Lastly, investigating the temporal changes in biomass density helps to understand the interaction between the flow physics and biology used in the coupled Hydro‐MEM. In the year 2020 (Figure 6a), under the low (linear projected) SLR case, accretion aids the marshes to maintain their position in the tidal frame, thus enabling them to remain productive. However, under higher SLR scenarios, the magnitude and rate of inundation prevented this adaptation and salt marshes began to lose their productivity. The sediment accretion rate in the SLR scenarios also affects the biomass density. The accretion and SLR rate associated with the low SLR scenario in the years 2050 and 2080 (first row from top in Figures 6b and 6c) increased productivity. The higher water levels associated with the intermediate‐low SLR scenario lowered marsh productivity and generated new impoundments in the upper islands between the Apalachicola River and East Bay (second row of Figures 6b and 6c). Under the intermediate‐high and high SLR, the trend continued and salt marshes lost productivity, drowned or disappeared altogether. It also caused the disappearance or productivity loss and inundation of salt marshes in other parts near Lake Wimico. The inundation direction under the intermediate‐low SLR scenario began from regions of no productivity in upland, extended over low productivity areas, and stopped at higher productivity regions due to an increased organic and inorganic accretion rates and larger bottom friction coefficients. Under intermediate‐high and high SLR scenarios, the migration to higher lands was apparent in areas that have a topographic profile that accommodates regular flooding when sea level rises. In the year 2100 (Figure 6d), the accretion and SLR rate associated with the low SLR scenario increased productivity near East Bay and produced a more uniform salt marsh. The productivity loss near Lake Wimico is likely due to the increase in flood depth and duration. For the intermediate‐high and high SLR scenarios, large swaths of salt marsh were converted to open water and some of the salt marshes in the upper elevation range migrated to higher lands. The area available for this migration was restricted to a thin band around the extended bay that had a topographic profile within the new tidal frame. If resource managers in the area were intent on providing additional area for future salt marsh migration, targeted regrading of upland area projected to be located near the future shoreline is a possible measure for achieving this. The results from this research are critical to properly define open water areas that are created from marshes drowning under SLR and then to parameterize bottom friction for storm surge modeling [Bilskie et al., 2016b]. Local storm surge is nonlinearly increased (relative to the SLR, particularly under intermediate‐high and high SLR) when the open water regions are expanded and bottom roughness is modified so that marsh areas that are likely to be submerged have a reduced roughness.

5 Conclusions The Hydro‐MEM was used to simulate the wetland response to four SLR scenarios in Apalachicola, Florida. The model coupled a two‐dimensional depth‐integrated hydrodynamic model and a parametric marsh model to capture the dynamic effect of SLR on salt marsh productivity. The parametric marsh model used empirical constants derived from experimental bioassays installed in the Apalachicola marsh system over a 2‐year period. The marsh platform topography used in the simulation was adjusted to remove the high elevation bias in the lidar‐derived DEM. The average annual Apalachicola River flow rate was imposed as a boundary condition in the hydrodynamic model. The results for biomass density in the current condition were validated using remotely sensed‐derived biomass density. The water levels and biomass density distributions for the four SLR scenarios demonstrated a range of responses with respect to both SLR magnitude and rate. The low and intermediate‐low scenarios resulted in generally higher water levels with more extreme gradients in the rivers and creeks. In the intermediate‐high and high SLR scenarios, the water level gradients were less pronounced due to the large extent of inundation (45% and 68% increase in inundated area). The biomass density in the low SLR scenario was relatively uniform and showed a productivity increase in some regions and a decrease in the others. In contrast, the higher SLR cases resulted in massive salt marsh loss (conversion to open water), decreased productivity and migration to areas newly within the optimal tidal frame. The inundation path generally began in areas of no productivity in upland, proceeded through low productivity areas, and stopped when the local topography prevented further progress.

6 Future Considerations One of the main factors in sediment transport and marsh geomorphology is the velocity variation within tidal creeks. These velocities are dependent on many factors including water level, creek bathymetry, bank elevation, and marsh platform topography [Wang et al., 1999]. The maximum tidal velocities for four transects, two in the distributaries flowing into the Bay in the islands between Apalachicola River and East Bay (T1 and T2), one in the main river (T3), and the fourth one in the tributary of the main river far from main marsh platform (T4) (Figure 1) were calculated. The magnitude of the maximum tidal velocity generally increased with increasing sea level (Figure 7). The rate of increase in the creeks closer to the bay was higher due to reductions in bottom roughness caused by loss of biomass density. The magnitude of maximum velocity also increased as the creeks expanded. However, these trends leveled off under the intermediate‐high and high SLR scenarios when the boundaries between the creeks and inundated marsh platform were less distinguishable. This trend progressed further under the high SLR scenario, where there was also an apparent decrease due to the bay extension effect. Figure 7 Open in figure viewer PowerPoint Maximum velocity variation with time under the low, intermediate‐low, intermediate‐high, and high SLR in four transects shown in Figure 1 Under SLR scenarios, the maximum flow velocity generally, but not uniformly, increased within the creeks. Transect 1 was chosen to show the marsh productivity variation effects in the velocity change. Here, the velocity increased with increasing sea level. However, under the intermediate‐low and intermediate‐high SLR scenarios, there were some periods of velocity decrease due to creation of new creeks and flooded areas. The new creeks were created using the complex physical processes considering the changes in marsh platform elevation, Manning's n, tidal flow velocity, etc. within the ADCIRC element of the Hyro‐MEM model. This effect was also seen in transect 2 for the intermediate‐high and high SLR scenarios. Under the low SLR scenario in transect 2, a velocity reduction period was generated in 2050 as shown in Figure 7. This period of velocity reduction is explained by the increase in roughness coefficient related to the higher biomass density in that region at that time. Transect 3, because of its location in the main river channel, was less affected by SLR‐induced variations on the marsh platform but still contained some variability due to the creation of new distributaries under the intermediate‐high and high SLR scenarios. All of the curves for transect 3 show a slow decrease at the beginning of the time series indicating that tidal flow dominated in that section of the main river at that time. This is also evident in transect 4 located in a tributary of the Apalachicola River. The variations under intermediate‐high and high SLR scenarios in transect 4 are mainly due to the new distributaries and flooded area. One exception to this is the last period of velocity reduction occurring as a result of the bay extension. The creek velocities generally increased in the low and intermediate‐low scenarios due to increased water level gradients, however, there were periods of velocity decrease due to higher bottom friction values corresponding with increased biomass productivity in some regions. Under the intermediate‐high and high scenarios, velocities generally decreased due to the majority of marshes being converted to open water and the massive increase in flow cross‐sectional area associated with that. These changes will be considered in the future work of the model related with geomorphologic changes in the marsh system.

Acknowledgments This research is partially funded under Award No. NA10NOS4780146 from the National Oceanic and Atmospheric Administration (NOAA) Center for Sponsored Coastal Ocean Research (CSCOR), and the Louisiana Sea Grant Laborde Chair. Partial support was provided by NSF DEB 1052636. The computations for the hydrodynamic model simulations were performed using the STOKES Advanced Research Computing Center (ARCC) at University of Central Florida, the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation Grant Number ACI‐1053575 [Towns et al., 2014], and High Performance Computing at Louisiana State University (LSU) and the Louisiana Optical Network Initiative (LONI). The authors would like to extend our thanks appreciation to ANERR staffs, especially Mrs. Jenna Harper for their continuous help and support. The statements and conclusions do not necessarily reflect the views of NOAA‐CSCOR, Louisiana Sea Grant, STOKES ARCC, XSEDE, LSU, LONI, ANERR or their affiliates. All data for this paper are properly cited and referred to in the reference list.