Understanding the Eco‐Geomorphologic Feedback of Coastal Marsh Under Sea Level Rise: Vegetation Dynamic Representations, Processes Interaction, and Parametric Sensitivity

A growing number of coastal eco‐geomorphologic modeling studies have been conducted to understand coastal marsh evolution under sea‐level rise (SLR). Although these models quantify marsh topographic change as a function of sedimentation and erosion, their representations of vegetation dynamics that control organic sedimentation differ. How vegetation dynamic schemes contribute to simulation outcomes is not well quantified. Additionally, the sensitivity of modeling outcomes to parameter selection in the available formulations has not been rigorously tested to date, especially under the influence of an accelerating SLR. In this study, we used a coastal eco‐geomorphologic model with different vegetation dynamic schemes to investigate the eco‐geomorphologic feedbacks of coastal marshes and parametric sensitivity under SLR scenarios. We found that marsh platform relief increased with SLR rate. The simulations with different vegetation schemes exhibited different spatial‐temporal variations in elevation and biomass. The nonlinear Spartina scheme presented the most resilient prediction with generally the highest marsh accretion and vegetation biomass, and the least elevation relief under SLR. But the linear Spartina scheme predicts the lowest unvegetated‐vegetated ratio. We also found that vegetation‐related parameters and sediment diffusivity, which were not well measured or discussed in previous studies, were identified as some of the most critical parameters. Additionally, the model sensitivity to vegetation‐related parameters increased with SLR rates. The identified most sensitive parameters may inform how to appropriately choose modeling representations of key processes and parameters for different coastal marsh landscapes under SLR and demonstrate the importance of future field measurements of these key parameters.

For their representation of vegetation-related processes, some modeling studies assumed static vegetation with a constant influence of vegetation on hydrodynamics and sedimentation (Allen, 1995;D'Alpaos et al., 2011;J. R. French, 1993;Mudd et al., 2009;Rogers et al., 2012;Schile et al., 2014;Stralberg et al., 2011;van Wijnen & Bakker, 2001). Other studies modeled more detailed vegetation-water-land interactions by considering the impact of vegetation density, height, and submergence on water flow and sediment transport (e.g., Da D'Alpaos et al., 2007;Duvall et al., 2019;Mudd et al., 2009Mudd et al., , 2004Temmerman et al., 2005). Morris et al. (2002) first proposed a clear relationship between marsh vegetation biomass and its depth below mean highest tide level based on the field observation on the coastal marsh in South Carolina, USA. Other studies extended this work to explicitly integrate quantitative representations for vegetation dynamics into coastal marsh evolution by assuming (1) a linear relationship between Spartina-dominant vegetation and its inundation condition (Belliard et al., 2015;D'Alpaos et al., 2007), (2) a nonlinear relationship between Spartina-dominant vegetation and its inundation condition (Kirwan & Murray, 2007;Mariotti & Fagherazzi, 2010), or (3) a linear relationship between multiple vegetation species and their inundation condition (Belliard et al., 2015;D'Alpaos et al., 2019D'Alpaos et al., , 2007Marani et al., 2013Marani et al., , 2004Silvestri et al., 2005). Many of these modeling studies evaluated the vulnerability of coastal marshes under SLR by using a lumped approach, where they treated coastal marshes as a single point or only focused on the marsh near the seaward boundary without an examination of the marsh spatial variation from the ocean to the upland in responding to SLR (D'Alpaos et al., 2011;J. French, 2006;Kirwan et al., 2010;Kirwan & Temmerman, 2009;Mudd et al., 2009;Temmerman et al., 2003;van Wijnen & Bakker, 2001). Other studies investigated the spatial and temporal variation of coastal marsh evolution under SLR Kirwan, Walters, et al., 2016;Marani et al., 2013;Ratliff et al., 2015). However, the response of coastal marsh evolution under SLR to varying representations of vegetation dynamic processes is still not well understood, especially the coevolution of coastal marsh elevation and vegetation. Furthermore, as the complexity and sophistication of these coastal models continues to increase, there is a critical knowledge gap in how sensitive model predictions are to model parameterizations under different SLR conditions. This knowledge is critical for developing effective model parameterizations and designing field studies to constrain those model parameters under different SLR scenarios. Currently, this knowledge gap limits our confidence in the application of these types of models to inform coastal wetland management and protection.
In this study, we used a coastal eco-geomorphic model with different vegetation dynamic representations to investigate the eco-geomorphologic feedbacks in coastal marshes under future SLR conditions to address the following two questions: 1. How will the selection of vegetation representations result in spatial and temporal differences in eco-geomorphologic outcomes of coastal marshes under SLR? 2. How will the different vegetation representations and different rates of SLR affect model parametric sensitivity?
To address these questions, we simulated the evolution of a one-dimensional coastal marsh transect using a well-established coastal eco-geomorphologic model from D' Alpaos et al. (2007). Specifically, under two commonly used future global mean SLR scenarios (SLR = 0.005 m/yr and SLR = 0.01 m/yr, corresponding to Representative Concentration Pathways (RCP) 4.5 and RCP 8.5 scenarios in Phase 5 of the Coupled Model Intercomparison Project (Spencer et al., 2016), we explored three different dependencies of vegetation biomass on elevation above mean sea level: linear and nonlinear formulations for the Spartina-dominant vegetation Mariotti & Fagherazzi, 2010;Morris et al., 2002) and the mixed-species linear function . After comparing the spatial and temporal variations of coastal marsh evolution under SLR with different vegetation equations, we used a global sensitivity approach to evaluate the sensitivity of eco-geomorphologic processes to model parameterizations spanning a wide range of the parameters.
The study begins by introducing process representation in Section 2, followed by model introduction, study site, experiment design, and model setting in Section 3. Then, we analyze the marsh evolution and model sensitivity under different rates of SLR, vegetation schemes, and maximum organic soil production rates in Section 4. Finally, we discuss the implications of this study for understanding the vulnerability of coastal marsh under SLR, guiding data-model integration, representativeness, and uncertainties.

Background: Process Representation in Eco-Geomorphologic Models
Eco-geomorphologic models represent topographic change of coastal marsh as the net balance of sediment erosion and deposition (Fagherazzi et al., 2012). Based on mass conservation, the spatially averaged dynamics of topographic elevation in a coastal landscape can be expressed as where z is the surface elevation relative to the mean sea level with the dimension of (L); t is time (T); p is the porosity of bed sediment; D and E represent local sediment deposition and erosion rates with the dimensions of (LT −1 ), respectively; R is the rate of SLR (LT −1 ).
where E shear is the erosion due to bed shear stress (LT −1 ). Erosion occurs when the bed shear stress (τ 0 ) exceeds the critical shear stress for erosion (τ e ), viz where α is the erosion coefficient. E break in Equation 2 is the erosion due to wave breaking (LT −1 ). According to Mariotti and Fagherazzi (2010), E break is a function of wave power dissipated by breaking: where β is the wave erosion coefficient; P is the wave power per unit area (WT −1 L −2 ); P cr is the threshold of wave power for wave erosion (WT −1 L −2 ); d is the spatial interval over which wave breaking occurs (L).
The sedimentation rate, D in Equation 1, is given by where D s is the inorganic sediment settling rate (LT −1 ), which is a function of settling velocity (w s ; LT −1 ; Cao et al., 2020), suspended sediment concentration (C), bed shear stress (τ 0 ) due to water flow (ML −3 T −2 ), and critical shear stress for sedimentation (τ d ; ML −3 T −2 ; Krone, 1962), namely, D t in Equation 5 is the inorganic sediment trapping rate due to the effect of vegetation canopy (LT −1 ), which can be represented by an empirical form where U is the water flow velocity (LT −1 ); ϵ is a capture efficiency of vegetation stems, h w is the water flow depth (L), and several vegetation characteristics, such as plant stem diameter (d s ), stem density (n s ), and vegetation height (h s ) (Mudd et al., 2004;Palmer et al., 2004). Additionally, D o in Equation 5 is the organic matter production rate (LT −1 ), which is a function of plant biomass, viz where K b is the maximum production rate of belowground organic material (LT −1 ), B is the aboveground plant dry biomass at the current time (ML −2 ); B max is the maximum vegetation biomass ( coastal marsh vegetation is controlled by several factors related to nutrient inputs (e.g., nitrogen and phosphorous) and soil environmental stress (e.g., oxygen availability, salinity, and sulfide concentration) . Morris et al. (2002) proposed a relation between vegetation biomass and the depth of the marsh surface below the mean highest tidal level based on observations at a coastal marsh in South Carolina, USA. Based on this relation, several empirical functions were derived to represent equilibrium vegetation biomass under different ponding conditions. The empirical function can be expressed as a linear  or a parabolic (Morris et al., 2002) function of salt marsh elevation relative to tide level. For the linear dependency, the lowland area with frequent flooding is more favorable for salt-tolerant and flood-tolerant species, such as Spartina alterniflora. The vegetation biomass is proportional to inundation depth. Quantitatively, the biomass equation can be written as (see the blue line in Figure 2) where B 1 is the time-averaged aboveground biomass density (ML −2 ); B max is the maximum biomass density (ML −2 ); MHTL represents the mean highest tide level (L); D biomax and D biomin are the highest and lowest depth below MHTL, respectively, which bounds the upper and lower limits of vegetation growth range . MHTL − D biomin and MHTL − D biomax represent the elevations of the upper and lower boundaries for vegetation growth (the dashed lines in Figure 2). Whereas, some mixed species on marshland prefer higher elevation region with less flooding and better aerated soil (see the orange line in Figure  2), namely where B 2 is the time-averaged aboveground biomass density for mixed species (ML −2 ; D' Alpaos et al., 2007). Besides these linear functions, a parabolic formulation describes that the plant biomass goes to zero when the marsh surface elevation reaches the upper (MHTL − D biomin ) or lower bound (MHTL − D biomax ), and the biomass reaches its peak at a certain elevation between MHTL − D biomin and MHTL − D biomax (see the yellow line in Figure 2): where B 3 is the time-averaged aboveground biomass density (ML −2 ; Morris, 2006); D is the ratio between MHTL − D biomin − z and D biomax − D biomin ; a, b, and c are fitting coefficients.
The representation of marsh hydrodynamics driven by tides and waves is also an essential part of eco-geomorphologic modeling because both erosion and sedimentation are fundamentally tied to surface water flow (Scheidegger, 1961). The shallow water equations, derived from the depth-integrated Navier-Stokes equations, have been widely used to compute hydrodynamics in coastal regions where the water horizontal length scale is much greater than the vertical length scale (Vreugdenhil, 2013). Specifically, the shallow water equations consist of two conservation equations: (1) conservation of mass and (2) conservation of momentum. Namely, in a one-dimensional (1-D) domain, where h is the water surface elevation = land surface elevation (z) + local water flow depth (η; L), thus h varies not only depending on the change in water depth, but also the simultaneous morphological change; u is the flow velocity (LT −1 ); g is the gravitational acceleration (LT −2 ); x is the spatial direction along the 1-D domain (L); C is the Chezy's friction coefficient.

Numerical Model
We ), except sediment erosion due to waves because the effect of waves in controlling the spatial and temporal variation of coastal marsh evolution was well studied by Duvall et al. (2019) and Mariotti and Fagherazzi (2010), and vegetation can significantly mitigate waves if the waves are not too strong, thus wave-induced erosion is not a focus in this study. We focused on conditions with regular semidiurnal tidal cycle and background SLR. For the representation of vegetation biomass, the original D-model included functions (e.g., Equations 9 and 10) that assume a linear relationship between annual averaged biomass and the elevation relative to mean sea level and considered different responses of Spartina and mixed vegetation species (see details in Section 2). To have a comprehensive understanding of the differences of the eco-geomorphologic feedbacks under different representations of vegetation dynamics, we incorporated the nonlinear function (e.g., Equation 11) into the D-model as well. For the computation of hydrodynamics, the D-model uses an approach similar to the kinematic-wave form that assumes a balance between water surface slope and friction in the momentum equation (Equation 13) Rinaldo et al., 1999). The detail of the hydrodynamic component is referred to the Text S1 and D' Alpaos et al. (2007), and the detail for the sediment transport component is referred to Section 2 above and D' Alpaos et al. (2007).

Numerical Experiment
We used a 1-D transect based on a marsh platform along the Delaware Bay, USA, as a prototype for our simulations (the black solid line in Figure 3c). Marsh surface elevation in the 1-D transect is at a level close to the MHTL (gray dashed line in Figure 3c), consistent with observations in Delaware Bay based on the CoNED coastal elevation database (Danielson et al., 2016) and National Oceanic and Atmospheric Administration (NOAA) tide observations (NOAA, 2001), which indicates that the landscape is at or close to an equilibrium state under the current sea level conditions .
This study simplifies the 1-D transect topography by using a linear interpolation of the observed topography (red line in Figure 3c) as the initial land surface elevation for the numerical experiments. The origin of the 1-D model domain is placed at the seaward boundary (x = 0), whereas the upland boundary is located at x = L. Water and sediment can only flow through the seaward boundary with zero flux flowing through the upland boundary. The current mean sea level (MSL) is at −0.13 m above North American Vertical Datum of 1988 (NAVD88), and the averaged tide amplitude is about 0.8 m based on the NOAA tide and current observation at station Cape May, NJ [8536110] (the red star in Figure 3b). We used a constant suspended sediment concentration (C 0 = 20 mg/L) at the ocean boundary (x = 0). The value of C 0 falls at the lower bound of the range of sediment concentration used in the previous coastal eco-geomorphologic modeling studies (e.g., Kirwan, Walters et al., 2016). Thus, this study makes a conservative prediction of coastal marsh change under SLR. However, a comparable numerical experiment with the same model settings but with a higher suspended sediment concentration (C 0 = 100 mg/L) was also conducted, and the results can be found in the supporting information (see Figures S2 and S3).
ZHANG ET AL. 6 of 23

10.1029/2020JF005729
In order to speed-up simulations to geomorphologically relevant time scales, the simulations adopted a morphological scaling factor (MSF, e.g., Lesser et al., 2004;Roelvink, 2006;Zhang et al, 2016), which effectively assumes that changes in the topographic profile over time scales smaller than the scaling factor do not appreciably affect the flow field and the eco-geomorphic dynamics. Hence, elevation change is computed offline by applying sediment fluxes determined in a tidal cycle, assumed to be constant for a period of time equal to the MSF. Thus, in this study, the simulations were run for 500 years (consistent with the simulation time in D'Alpaos et al., 2007 to make sure the landscape reaches an equilibrium state) with a spatial interval of 1 m and a time interval of 10 min for hydrodynamics in a single tidal cycle and MSF = 50 for the eco-geomorphologic change of 50 tidal cycles.
We designed several focused numerical experiments to characterize eco-geomorphologic feedbacks under different representations of vegetation dynamics and SLR scenarios for the future 500 years. We adopted two commonly used future global mean SLR scenarios from global climate model predictions, including (1) the relatively low SLR rate (0.005 m/yr) Ganju et al., 2020;Kirwan & Temmerman, 2009;Spencer et al., 2016) and (2) the relatively high rate of SLR (0.01 m/yr) (Ganju et al., 2020;Kirwan, Walters et al., 2016;Orson et al., 1985;Spencer et al., 2016). In addition, we considered three different representations of vegetation dynamic processes, such as the Spartina-dominant linear function, Spartina-dominant nonlinear function, and mixed species linear function. Also, in simulating vegetation organic soil production, we incorporated two different rates of maximum organic production rates: (1) K b = 0.003 m/yr, a commonly used maximum organic production rate under current climate (Langley et al., 2009;Morris et al., 2016) and (2) K b = 0.005 m/yr, a larger maximum organic production rate, reflects the increase of belowground biomass productivity under elevated atmospheric CO 2 in the future (Ratliff et al., 2015). Specifically, based on a comprehensive literature review, Ratliff et al. (2015) found that biomass productivity increased about 33% for a 400 ppm increase in atmospheric CO 2 . Here, we assumed a present-day maximum organic production rate of 3 mm/yr. Under the RCP 8.0 climate scenario (the most ambicis future CO 2 emission scenario), the CO 2 -equivalent levels was projected to exceed 1,200 ppm, which means an additional 800 ppm with respect to present (Hayhoe et al., 2017). This leads to a future maximum organic production rate equal to 0.00498 m/yr ≈ 0.005 m/yr. Scenario details are listed in Table 1. The parameters for these individual simulations are listed in the fourth column in Table 2.

Sensitivity Analysis
There are many sensitivity analysis approaches available to understand parametric sensitivity of model behavior (see Song et al., 2015 for a detailed review). In this study, we used a widely applied sensitivity analysis approach, the Fourier Amplitude Sensitivity Test (FAST) technique (Cukier et al., 1973;Xu & Gertner, 2008a. FAST is computationally efficient and can be used for both nonlinear and nonmonotonic relationships between parameters and model outputs (Xu & Gertner, 2011). FAST uses a periodic sampling strategy to assign a characteristic periodic signal for each parameter. Within FAST, a Fourier transformation is used to decompose the variance in model outputs into partial variance contributions by individual model parameters based on the assigned signals. The ratio of partial variance contributed by a specific parameter to the total variance of a model output is defined as the first-order sensitivity index to measure the importance of each model parameter. The FAST analysis has been incorporated into a software tool, the UASA High SLR rate (0.01 m/yr) and low K b (0.003 m/yr) Case 7 Case 8 Case 9 Low SLR rate (0.005 m/yr) and low K b (0.003 m/yr) Case 10 Case 11 Case 12 Abbreviation: SLR, sea-level rise.   Xu and Gertner (2008b) and provides a rigorous way of defining, executing, and analyzing experiments for model parametric sensitivity.
This study selected 11 common parameters that have been used in many coastal eco-geomorphologic models (see the list of the parameters in Table 2). Based on this selection, the UASA ToolBox was used to generate 1,100 groups of parameters for the model ensemble simulations to quantify the models' individual parametric sensitivities. The range of each parameter is estimated based on our literature survey and empirical knowledge. However, because there is not enough data to derive informative probability density distributions, we used a uniform distribution for our sensitivity analysis.
Model sensitivity is defined in terms of relevant quantitative metrics describing the final state of the system: (1) the difference between the MHTL and the elevation at the seaward boundary (MHTL minus elevation, hereinafter referred to as Depth_m), (2) the difference between minimum and maximum elevations (hereinafter referred to as elevation relief) from each ensemble simulation under different scenarios, (3) domain averaged sediment fluxes, (4) the vegetation biomass at the seaward boundary, and (5) the vegetation biomass at the upland boundary. Notably, the first metric, Depth_m, measures how the landscape elevation (at least the seaward boundary) responds to SLR. While the second metric (elevation relief) measures the difference of elevation at the seaward boundary and inland and possible inland depression on the 1-D marshland.

Topographic Change Across Individual Simulations
We first used the 12 individual simulations (cases 1 to 12 in Table 1) as examples to compare the elevation change under different vegetation equations for biomass estimation and SLR scenarios simulated by the D-model over 500 years (see Figure 4). The corresponding sediment fluxes at the end of the 500 years are illustrated in Figure 5. Domain-wide, the elevations in the cases with a higher maximum organic production rate (K b ) (the first column in Figure 4) were higher than the elevation in the lower K b cases (the second column in Figure 4). At the seaward boundary, the relative locations between elevation and MHTL near the seaward boundary in all vegetation-covered scenarios remained constant after 400 years' simulation (not shown at here) due to a balance between sediment fluxes ( Figure 5), which indicate that (1) (Figure 4b), where the simulated marsh surface elevation with the mixed-vegetation equation (hereinafter referred to as mixed-veg case) was equal to the MHTL because the organic production rate is equal to the SLR (see the sediment flux in Figure 5j).
Moving landward, the marsh elevations declined due to a decrease in sedimentation rate landward. Some of the marshland became totally submerged in water as the elevation was below the final mean sea level (final MSL indicated by the blue dashed lines in Figure 4). With a higher K b , a shorter portion of the marshland was below the final MSL because a higher K b resulted in a higher organic sedimentation rate, which dominantly contributed to the accretion rate at the upland area where inorganic sediment from the ocean was restricted to this region.
High SLR caused a larger elevation relief up to 5 m (Figures 4a and 4c), compared to the low SLR scenarios that only had a maximum elevation relief of 0.48 m in the higher K b case ( Figure 4b) and 2.85 m in the lower K b case (Figure 4d), respectively. Notably, in the lower SLR scenario with a higher K b , the elevation only slightly declined landward (solid lines in Figure 4b), which means that the accretion rate domain-wide can always keep pace with the rate of SLR as illustrated in Figures 5h-5j. Different vegetation schemes also highlight different influences on the topographic outcomes. Among all the vegetation cases, the Spartina-nonlinear case showed the highest final elevation and the least elevation relief due to the highest sedimentation rate throughout the domain, particularly due to organic soil production in the middle and upper portions of the transect. Elevation declined closer to the ocean boundary in the mixed-veg cases than the elevations with Spartina-dominant linear and nonlinear functions (hereinafter referred to as Spartina-linear case and Spartina-nonlinear case, respectively) (thick and thin black solid lines in Figure 4). Notably, in the mixed-veg case under low SLR and high K b (gray solid line in Figure 4b), the elevation reached a level similar to the MHTL. This was because the vegetation growth in the mixed-veg case is greater at lower inundation levels. Thus, vegetation continued growing even when the elevation was at the same level of the MHTL. In the Spartina-dominant cases, the Spartina-nonlinear cases showed the declines of elevation started closer to the seaward boundary than the elevation decline in the Spartina-linear cases (black solid lines in Figures 4a, 4c, and 4d).
The simulations with a higher suspended sediment concentration in the ocean (C o = 100 mg/L), a higher SLR rate (0.01 m/yr), and a higher K b (0.005 m/yr) showed similar profiles with the simulations under a lower suspended sediment concentration from the ocean: the marsh elevation near the ocean boundary was at the similar level with the future MHTL, and the marshland at the upland was drowning (see Figure S2) at the end of the 500 years' simulation. However, with a higher sediment supply from the ocean, we observed a longer portion of the marsh elevation (∼300 m from the ocean boundary) can keep pace with the increase of sea level, compared with the elevation profile in the lower sediment concentration cases. The high sediment concentration case also predicted less elevation relief than that in the low sediment concentration case due to a higher sediment supply from the ocean.
For the contributions of sediment fluxes to marshland accretion, in general, sediment settling rate contributed more than sediment trapping rate and organic production rate near the seaward boundary (light blue lines in Figure 5) in all vegetation-covered cases, except the mixed-veg cases (Figures 5j and 5n) where the organic production rate was higher than the other fluxes. This is because the mixed-veg case assumes that vegetation can grow better under lower inundation or no inundation conditions where vegetation organic production always plays a role in contributing to marsh accretion, but inorganic sediment settling contributes less due to limited delivery of sediment landward. Given that the elevation near the seaward boundary accreted faster than the inland area, the inundation depth near the seaward boundary was shallower than the inland, which provided a more favorable condition for mixed vegetation species to grow near the seaward boundary, resulting in a higher organic production rate there.
Moving landward, the inorganic sediment settling rate was still a dominant sediment flux contributing to the accretion rate, except the cases with a lower SLR rate and higher K b , where the organic production rate was dominant (purple lines in Figures 5h-5j).

Model Parametric Sensitivity From Ensemble Simulations
We explored the model parametric sensitivity represented by the ratio of individual parametric variance to the total variance from the ensemble simulations across different combination of parameters spanning wide ranges of their values (see Table 2).

Parametric Sensitivity for Topographic Change
For the sensitivity of modeled Depth_m (defined as MHTL minus elevation in Section 3.3) to parameterization (Figure 6a), vegetation-related parameters showed a larger influence on Depth_m under the higher SLR rate scenarios (e.g., the first three columns in Figure 6a). While, under the lower SLR rate scenarios, the sediment-related parameters, especially the "sediment concentration," were the dominant parameters (the last three columns in Figure 6a). For the different vegetation dynamic schemes, the mixed-veg cases were highly sensitive to the "maximum organic production rate" indicating that the Depth_m was highly dependent on the organic matter production rate regardless of the rates of SLR because some species in the mixed-veg cases can grow under more prolonged flooding condition, and the other species are adapted to less frequent and prolonged flooding condition, such that the vegetation processes can contribute to sedimentation in all conditions. While, in the Spartina-dominant cases, the vegetation can only grow under more prolonged flooding condition driven by SLR and tide. Thus, the parametric sensitivities in the Spartina-linear and -nonlinear cases (the first, second, fourth, and fifth columns in Figure 6a) were controlled by the inundation condition, the sediment settling, and vegetation processes and did not present a huge difference among parameters, compared to the mixed-veg cases. The relatively more sensitive parameters are "maximum organic production rate," "maximum biomass," "water depth for plant growth," "sed concentration," and "critical shear stress for deposition." Among these parameters, the "maximum organic production rate" was the most sensitive parameter in the high SLR scenario (the first and second columns in Figure 6a). In contrast, "sediment concentration" was the most sensitive parameter in the low SLR scenarios (the fourth and fifth columns in Figure 6a). This is because the high rate of SLR has a larger potential to cause a higher inundation condition by high tides, a favorable condition for Spartina to grow. Therefore, the vegetation effect had a larger contribution to sedimentation than the contribution from vegetation in the low SLR cases.
For the sensitivity of elevation relief, in the high SLR scenario (the first, second, and third columns in Figure  6b), the model simulations were more sensitive to "sediment diffusivity," an important parameter in the sediment diffusion equation that controls how much sediment could diffuse landward. "Tide amplitude" was also one of the most sensitive parameters in the Spartina-linear and -nonlinear cases. The vegetation-related parameters showed relatively low sensitivity under the higher SLR rate (the first and second columns in Figure 6b), which means that the elevation relief was more dependent on how much sediment ZHANG ET AL. 12 of 23 10.1029/2020JF005729 can transport landward and deposit under a high SLR rate. However, in the mixed-veg case, the "maximum organic production rate" along with "sediment diffusivity" were the most sensitive parameters, which reflects the tolerance of the growth of mixed vegetation species in different conditions. For the lower SLR scenario (the fourth to sixth columns), the vegetation-related parameters showed higher sensitivity, which means that the vegetation processes were more dominant to the change of elevation relief, especially for the Spartina-linear and -nonlinear cases under a low rate of SLR. The values of sensitivity for each parameter in each scenario can be found in Tables S1 and S2.

Parametric Sensitivity for Sediment Fluxes
The parametric sensitivities of sediment fluxes to model parameterization are similar to corresponding vegetation cases under both SLR scenarios. For example, the Spartina-linear cases (the first to third columns in Figure 7a and first to third columns in Figure 7b) under both the high and low SLR scenarios show a similar parametric sensitivity for each corresponding sedimentation processes.
Specifically, for the sediment settling process, all the cases (the first, fourth, and seventh columns in Figure  7a and the first, fourth, and seventh columns in Figure 7b) were most sensitive to the "maximum organic production rate," which may be because the organic production influences elevation changes that indirectly control sediment settling process. Besides the "maximum organic production rate," the model simulations were also sensitive to some sediment settling-related parameters, such as "sediment concentration," "settling velocity," and "critical shear stress for deposition," which are the key parameters directly control sediment settling process. For the organic soil production by vegetation, all the cases (the second, fifth, and eighth columns in Figure 7a and the second, fifth, and eighth columns in Figure 7b) were most sensitive to the "maximum organic production rate," the key parameter in organic soil production process. For the sediment trapping process, the sensitivity was almost evenly distributed for each parameter because sediment settling, sediment diffusion and advection, and vegetation all influence sediment trapping, however, the parameters of "sediment diffusivity" and "sediment concentration" that control the distribution of ZHANG ET AL.  Tables S3 and S4.

The Spatial and Temporal Variation of Vegetation Biomass From Individual Simulations
The different formulations for vegetation growth in response to inundation conditions (illustrated in Figure  2) lead to distinct patterns in biomass distributions and marsh response to tidal and SLR-induced flooding. Figure 8 showed the spatial variation of vegetation biomass at the end of 500 years in the simulations under the different vegetation dynamic schemes, rates of SLR, and K b s. In general, the spatial patterns of vegetation biomass corresponded to the marsh elevation profiles in Figure 4. For example, the locations of dramatic declines of vegetation biomass in the high SLR scenarios are well aligned with the topographical depression area in Figures 4a and 4c. In this low-lying region, the marsh elevations approach an unfavorable inundation condition for vegetation growth with high ponding water detrimental to vegetation growth. In contrast, the vegetation biomass even increased landward in the lower SLR and higher K b scenario because the entire domain kept pace with the SLR rate, and the inundation condition was still within the vegetation's growth range (see the elevation profiles in Figure 4b). Notably, the vegetation biomass of the mix-veg case reached its maximum biomass across the entire model domain (the gray dashed line in Figure 8). However, with the lower K b , the simulation shows an abrupt decrease when the marshland was submerged in water (Figure 8d), similar to the final biomass profile in the high SLR scenarios.
For the mixed-veg cases, despite the different locations of the abrupt decreases, they showed similar patterns under the lower K b (Figures 8c and 8d), but different responses under the higher K b (Figures 8a and  8b). In the low SLR condition with the higher K b (Figure 8b increased approximately linearly and then decreased to zero at further locations landward, compared with the mixed-veg cases. The Spartina-nonlinear cases showed a higher estimated vegetation biomass than the biomass in the Spartina-linear cases, but the biomass started to decrease to zero closer to the seaward boundary in the Spartina-nonlinear cases, which reflected the nature of the differences in the assumptions in the vegetation equations. In order to examine the temporal evolution of biomass across the marsh, we plotted the time series at three locations: the seaward boundary, and 100 m and 400 m landward of the boundary (Figure 9). Across the 12 simulation cases, the temporal evolution of biomass may be divided into three stages, though not all stages are presented at all locations or every scenario.
Rapid change characterizes the first stage. With the exception of mixed vegetation (Figures 9c, 9f, 9i, and 9l), all locations exhibited rapid increases in biomass for the first 100-200 years of the simulation. During the second stage, biomass continued to adjust but at significantly slower rates than the first stage. These adjustments are seen at the seaward boundary and 100 m locations in the Spartina-linear simulations under both SLR forcings and K b values (orange and cyan circles in Figure a, d, g, and j) and Spartina-nonlinear simulations (orange and cyan circles in Figures 9b, 9e, 9h, and 9k). A dramatic exception to the gradual adjustments in Stage 2 is the 400 m location in the Spartina-linear and -nonlinear rapid SLR scenarios (green circles in Figures 9a, 9b, 9d, and 9e) and the low K b scenarios (green circles in Figures 9d, 9e, 9j, and 9k) in which biomass rapidly drops back to a value of zero between 100 and 300 years.
The third stage is the period when a system enters a stable state or equilibrium state, which indicated that a new equilibrium or quasiequilibrium state was reached under the new rate of SLR. Examples of this stability include the limited changes in vegetation biomass near the seaward boundary and at the 100 m locations in all the cases. This is because the vertical accretion rate at these locations in all the cases always kept pace with the rates of SLR.

Parametric Sensitivity for Vegetation Dynamics
After analyzing the spatial and temporal variation of vegetation biomass change, we computed the sensitivity of biomass estimation at the seaward boundary and upland boundary to model parameterization based on the ensemble simulations. The biomass estimations were more sensitive to the vegetation-related parameters, especially the parameters of "maximum organic production rate" and "maximum biomass" (Figure 10). For the vegetation biomass at the upland (the second, fourth, and sixth columns in Figures 10a  and 10b), "maximum organic production rate" and "maximum biomass" were the two most dominant parameters that control the estimation of biomass. However, the vegetation biomass near the seaward boundary was also sensitive to sediment settling-related parameters. Specifically, in the higher SLR scenario, the biomass estimations near the seaward boundary (the first, third, and fifth columns in Figure 10a) were also sensitive to all the other parameters, except the parameters for erosion (e.g., "erosion coefficient"). In contrast, in the low SLR scenario, the most sensitive sediment settling-related parameters were only "sediment concentration" and "settling velocity" in the Spartina-dominant cases (the first and third columns in Figure  10b). The vegetation biomass estimation near the seaward boundary in the mixed-veg case was more sensitive to "maximum biomass" and "maximum organic production rate" than the rest parameters. The values of sensitivity for each parameter in each scenario can be found in Tables S5 and S6. ( Figure 4) and (2) the elevation landward declined and part of it drowned in water for the high SLR scenarios and low SLR with low K b . The elevation near the seaward boundary started to approach a new equilibrium state under the rising SL conditions around 100 years (e.g., the cyan circles in Figures 9a, 9b, 9d, 9e, 9g, 9h, 9j, and 9k), which was consistent with the findings in previous studies (D'Alpaos et al., 2011;Kirwan et al., 2008;Kirwan, Temmerman et al., 2016;Kirwan & Temmerman, 2009;Temmerman et al., 2003;van Wijnen & Bakker, 2001). This pattern of lower accretion rates in the interior of marshes has been previously documented in both modeling (D'Alpaos et al., 2019Kirwan, Walters et al., 2016;Langston et al., 2020;Marani et al., 2013;Mariotti, 2016;Ratliff et al., 2015;Thorne et al., 2018) and field studies (Friedrichs & Perry, 2001;Palinkas & Engelhardt, 2019;Schepers et al., 2017;Temmerman et al., 2003). Eventually, the interior marshland died-off and turned into water pools as shown in Figure 9 (especially the temporal change of biomass at the 400 m location). Schepers et al. (2017) also reported that the size of water pools would expand through time, influencing the connectivity between marshland and channels.

Coastal Marsh Vulnerability Under
Under climate change, if the maximum organic soil production rate (K b ) increases to a similar level as the rate (0.005 m/yr) used in this study due to the increase of temperature and CO 2 in the future, the spatial and temporal variations of vegetation biomass are relatively small and vary within the vegetation growth range (Figures 8b and 9g, 9h, and 9i) under the lower SLR rate (0.005 m/yr). Based on these results, a SLR of 0.005 m/yr does not appear to threaten the survival of coastal marsh systems characterized by these types of vegetation on a 500-year scale. However, for a K b rate commonly observed today (0.003 m/yr), all the SLR scenarios showed clear declines of surface elevation starting near the middle or upper of the domain (solid lines in Figure 4a) and continuing to the upland boundary illustrating that the accretion rate at the inland portion of the coastal marsh cannot keep pace with the future SLR rates. These inland areas turned into open water habitats with degradation and marsh vegetation mortality occurring after 200-300 years in these locations (Figures 8a and 9a and 9b), which may lead to the change of coastal marsh ecosystem functions and hydrological regime shift (Ganju et al., 2020).
The simulations above used a conservative sediment concentration rate from the ocean boundary (C 0 = 20 mg/L), which limited the delivery of sediment landward under the high SLR rate, resulting the drowning of upland marsh. However, in our simulations with a higher sediment concentration from the ocean (C 0 = 100 mg/L), more sediment entered the domain and improved the potential for survival of coastal marshland under a high rate of SLR. However, simulations with the higher sediment concentration delayed, but did not prevent upland submergence, which further demonstrated that coastal marsh is largely vulnerable under the high rate of SLR (0.01 m/yr) (see Figures S2 and S3). A microtidal regime (tidal range = 1.6 m) was used in the individual simulations. However, the elevation profile would be subject to change under various tidal conditions because different tidal conditions would cause different inundation conditions that would results in distinct marsh vertical accretion. The role of different tidal regime in controlling marsh evolution was also demonstrated by our ensemble sensitivity analysis, where we identified the "highest tide amplitude" as one of the most sensitive parameters for elevation relief (e.g., Figure 6b).

Marsh Vulnerability due to Vegetation Representation
The experimental cases with different vegetation schemes consistently predicted coastal marsh vulnerability under future SLR. Under a conservative sediment concentration from the ocean (C 0 = 20 mg/L), at the seaward boundary, marsh elevation accretion should keep pace with future SLR, regardless the rate of SLR and K b values. Landward, the inland part of the coastal marsh was resilient under the lower rate of SLR (0.005 m/yr) and simultaneously with the higher K b , but potentially vulnerable to collapse under high rate of SLR or with the lower K b .
Our simulations also highlighted marsh response to increased ponded water depth under future SLR. The mixed-veg scheme was the most resilient scenario under the lower SLR and with the higher K b (gray solid line in Figure 4b): the marsh accretion rate was equal to the SLR rate throughout the entire domain due to less inundation condition and high organic soil production rate. However, the mixed-veg scheme was the most vulnerable scenario under the higher SLR or with the lower K b (see the mixed-veg cases in Figures 4  and 8): the decline of marsh elevation started closer to the seaward boundary due to unfavorable high inundation conditions for vegetation growth. Except for the mixed-veg case under the lower SLR and higher K b , the Spartina-nonlinear scheme was the most resilient scheme in all cases-it predicted the largest elevation increases throughout the domain, the least elevation depression (Figure 4), and the highest vegetation biomass (Figure 8). These most resilient predictions from the Spartina-nonlinear treatment were attributed to the assumption of the nonlinear relationship between vegetation biomass and inundation condition. Vegetation biomass reaches its peak when the inundation depth is at the middle level or near the middle level of the vegetation growth range (defined by the D biomax and D biomin ) and does not have to be at the highest inundation level, compared with the Spartina-linear scheme. However, we also found that the elevation and vegetation biomass started to decrease closer to the seaward boundary in the Spartina-nonlinear case, compared with the Spartina-linear case, which implies that the Spartina-nonlinear case predicted a bit higher unvegetated-vegetated marsh ratio as defined in Ganju et al. (2017).
In addition, our simulation depicted the evolution of vegetation biomass with the evolution of marsh landscape (Figure 9), reflecting some of the plant life-history traits (Schwarz et al., 2018). The vegetation biomass of our studied marshland varied through different trajectories at the seaward boundary, mid-marshland, and the upland (Figure 9)

Implication to Data-Model Integration and Future Coastal Eco-Geomorphologic Modeling
Our sensitivity analysis captured the overall parametric sensitivity of the eco-geomorphologic processes in the model and highlighted how different representations of vegetation dynamics and SLR conditions affect the parametric sensitivity. We found that the "sediment concentration" and "tidal amplitude" are the most sensitive parameters for coastal marsh evolution, which are in agreement with the findings in prior studies Kirwan et al., 2010;Kirwan, Walters, et al., 2016;Temmerman et al., 2003). More importantly, this study also identified additional parameters that are highly sensitive for the spatial and temporal variations of key landscape characteristics, such as (1) the Depth_m (depth between MHTL and marsh elevation at the seaward boundary), (2) elevation relief, (3) averaged sediment fluxes, and (4) vegetation biomass near the seaward boundary and at the upland. These parameters include "sediment diffusivity," "maximum organic production rate," and "maximum biomass." Thus, this sensitivity analysis highlights the need for future modeling and field observations to better measure and parameterize these controls on marsh evolution.
In particular, our sensitivity analysis identified the parameter of "sediment diffusivity" as one of the most sensitive parameters for predicting marshland evolution, especially controlling elevation relief, which implies the importance of hydrodynamic process that brings water and sediment landward and back to ocean. Although the evaluation of coastal hydrodynamics is outside the scope of this study, a good representation of coastal hydrodynamics as a function of coastal boundary condition (e.g., tide and wave), topographic gradient, and vegetation effect (e.g., influencing surface roughness) is critical for predicting sediment budget accurately and is worth deeper investigation in future modeling studies (Best et al., 2018;Duvall et al., 2019).

Representativeness of the Model Simulations
In this study, we selected the parameter values and the rates of SLR that were widely used in previous modeling studies or were established in the literature from field measurements to ensure that the simulations were realistic and representative. Additionally, the formulations used to represent the dominant processes were selected from broadly used sedimentation, erosion, and vegetation dynamic equations. Thus, the individual simulations should reflect current model capabilities and formulations used to understand process interactions and marsh response to SLR. Based on the ensemble simulations, we generated a large number of parameter samples for the sensitivity analysis. Thus, the results of the sensitivity analyses reasonably reflected the overall sensitivity of the model processes over their physical parameter ranges.
To further demonstrate that the D-model appropriately captures the behavior of coastal evolution under SLR, we conducted some of the same simulations by using another well-established coastal eco-geomorphologic model developed by Mariotti and Fagherazzi (2010) Figures S4-S7). The simulations from the M-model also identified similar most sensitive parameters for different scenarios. For example, the sensitivity of Depth_m in the M-model to vegetation-related parameters also increased with the SLR rates. The most sensitive parameters for elevation relief under the higher SLR were also "sediment diffusivity," "sed concentration," and "Highest tide amplitude." The "maximum organic production rate" was a more dominant parameter for elevation relief under the lower SLR. Meanwhile, the most sensitive parameters for sediment fluxes and vegetation biomass were also the "maximum organic production rate" and "maximum biomass." Both the D-and M-models predicted the elevation relief of marshland under the higher SLR (0.01 m/yr) was up to 2.5-5 meters (Figures 4a and 4c and Figure S4). This is an accumulative effect of the different accretion rate between the marsh near the ocean boundary and the interior marsh over the 500-year scale, which reflected the model behaviors under different process representations, parameters, and external drivers. It also highlighted the transition zone between more resilient marshland near the river/ocean boundary and more vulnerable interior marsh. Future work is needed to validate the reality of this transition zone with a better field measurement of sedimentation, vegetation biomass, and other marsh accretion-related parameters in these areas.

Uncertainties and Future Work
We used two maximum organic production rates (0.003 and 0.005 m/yr) in the individual simulations in this study. The former one was adopted from the previous studies by Langley et al. (2009) and Morris et al. (2016), representing an averaged value of the maximum production rate under the present climate condition. The latter one was calculated based on Ratliff et al. (2015) and the adoption of the IPCC highest CO 2 emission scenario (Hayhoe et al., 2017), reflecting an increase of vegetation belowground organic production rate under a warming climate. The use of both rates illustrated model parametric sensitivity and the role of organic sedimentation in controlling future marsh elevation change. However, there are still some uncertainties that may affect the organic production rate under a warming climate. On the one hand, the models did not consider the increase of organic soil decomposition rate under the warming climate, which could be comparable or even higher than the organic production rate (Kirwan & Blum, 2011;Langley et al., 2009). On the other hand, it is possible that the soil organic decomposition rate will not increase much due to the constraint from soil aeration level, a factor controlling soil organic decomposition (Roner et al., 2016;. Therefore, the balance/imbalance between organic soil production and decomposition will be the key to a better understanding of organic accretion in coastal marsh evolution. Thus, future work should focus on better quantifying the organic accretion components (organic soil production and decomposition), as well as the drivers and limits that control these components.
Our sensitivity analysis showed the importance of "maximum biomass" and "organic production rate" for the prediction of marshland elevation changes. Within most of the current eco-geomorphologic models, they are fixed through time. However, future climate changes, higher temperature and CO 2 conditions might change the value of these parameters gradually in the evolution process. Therefore, to improve the prediction accuracy, it is critical to have process-based models that can incorporate the impact of a dynamic future climate on vegetation production and litter decomposition through time.
The simulations used a no-flow boundary condition at the upland boundary, which limits the water and sediment supply from uplands through upland surface and subsurface environments. An appropriate consideration of the hydrologic and geomorphologic connectivity with the upland region may improve the flexibility of our test model in realistically representing a wider variety of settings, in terms of the relevant hydrodynamic and sediment transport processes (Wohl et al., 2019;Zhang et al., 2018), especially for intertidal areas receiving water and sediment from both riverine and ocean sources (Gleichauf et al., 2014;Kirwan, Walters, et al., 2016;Wolfram et al., 2016;Yousefi Lalimi et al., 2020). Also, water and sediment fluxes from tidal channel to marshland were not considered in these 1-D simulations. Tidal channel can compensate for the spatial discrepancy in sediment accretion by routing water and sediment from upstream to the coastal area or from the ocean boundary to the upland (Belliard et al., 2016). At the seaward boundary, the models used constant sediment concentration in rivers/ocean, while variability in this concentration could contribute to the uncertainty in predictions of the accretion rate on coastal marshes. In addition, a more precise estimation of sediment concentration in the aquatic systems by using high resolution field measurements or a high-resolution, process-based coastal ocean model would improve the predictive capability of coastal marsh eco-geomorphologic models (Stumpf, 1983;Temmerman et al., 2003).

Conclusion
We used a coastal eco-geomorphologic model with different vegetation dynamic representations to investigate eco-geomorphologic feedbacks on the coastal marsh and changes in model parametric sensitivity under various future SLR conditions. We conducted model simulations by using a standard set of test cases with consistent model settings and parameters. This study explored coastal marsh evolution under SLR not only from the domain averaged features, but also from the spatial and temporal variations of key landscape characteristics, such as the elevation relief and biomass at the seaward boundary and upland. We found that evaluating the spatial and temporal coastal marsh evolution under different representations of vegetation dynamic process provides new insights to better understanding the uncertainty of predicting coastal marshes vulnerability facing future accelerating SLR from different process representations.
Qualitatively, the three vegetation dynamic schemes (Spartina-linear, Spartina-nonlinear, and mixed-vegetation linear equations) produce consistent evaluations of the vulnerability of the coastal marsh under high and low SLR rates. However, the Spartina-nonlinear scheme predicted the highest vegetation biomass and organic production rate, yielding the highest accretion rate and elevation, except for the mixed-veg case under the low SLR. The mixed-veg case represents the most resilient marsh type under low SLR with high K b , but is the most vulnerable case under high SLR. Except the mixed-veg case under the low SLR, all the Spartina-linear cases predicted the largest marsh extent and smallest open water area.
The sensitivity analysis study identified the parameters whose values most critically affect model outcomes under different SLR conditions. The parametric sensitivity of the eco-geomorphologic models (e.g., the Dand M-model used in this study) were not the same under the high and low SLR conditions. For example, the most sensitive parameter, such as the maximum organic production rate, in the simulation under the high SLR, was not the most sensitive parameter in the low SLR scenario. The differences in parametric sensitivity highlighted the importance of evaluating parametric sensitivity under different external drivers.
The identified most sensitive parameters can help inform how to appropriately model key processes in different coastal marsh landscapes under SLR and vegetation evolution. These identified key parameters under different climate change conditions can also serve to inform future field measurements studies.

Data Availability Statement
This study is a model-based study that uses two numerical models and synthetic experiments. The model source code are available in the model repositories from the Community Surface Dynamics Modeling System (CSDMS) https://csdms.colorado.edu/wiki/Model_download_portal. Model parameters are listed in the tables above. and Andrea D'Alpaos acknowledge support by Provveditorato for the Public Works of Veneto, Trentino Alto Adige, and Friuli Venezia Giulia, provided through the concessionary of State Consorzio Venezia Nuova and coordinated by CORILA in the framework of the Venezia 2021 Research Program.