Development of a submerged aquatic vegetation growth model in the Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST v3.4) model

The coupled biophysical interactions between submerged aquatic vegetation (SAV), hydrodynamics (currents and waves), sediment dynamics, and nutrient cycling have long been of interest in estuarine environments. Recent observational studies have addressed feedbacks between SAV meadows and their role in modifying current velocity, sedimentation, and nutrient cycling. To represent these dynamic processes in a numerical model, the presence of SAV and its effect on hydrodynamics (currents and waves) and sediment dynamics was incorporated into the open-source Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) model. In this study, we extend the COAWST modeling framework to account for dynamic changes of SAV and associated epiphyte biomass. Modeled SAV biomass is represented as a function of temperature, light, and nutrient availability. The modeled SAV community exchanges nutrients, detritus, dissolved inorganic carbon, and dissolved oxygen with the water-column biogeochemistry model. The dynamic simulation of SAV biomass allows the plants to both respond to and cause changes in the water column and sediment bed properties, hydrodynamics, and sediment transport (i.e., a two-way feedback). We demonstrate the behavior of these modeled processes through application to an idealized domain and then apply the model to a eutrophic harbor where SAV dieback is a result of anthropogenic nitrate loading and eutrophication. These cases demonstrate an advance in the deterministic modeling of coupled biophysical processes and will further our understanding of future ecosystem change.


Introduction
Submerged aquatic vegetation (SAV), or seagrass, includes rooted vascular plants that inhabit sediments of estuaries and coastal waters, with a wide global distribution. SAV involves important primary producers in shallow environments and provides a habitat for a number of aquatic organisms; it can slow water velocities and dampen wave energy to trap particulate material (Carr et al., 2010), as well as alter biogeochemical cycles through oxygenation of sediments (Larkum et al., 2006). The positive impact of ecosystem services provided by SAV presence has been well studied (Hemminga and Duarte, 2000;Nixon et al., 2001;Terrados and Borum, 2004;McGlathery et al., 2007;. The growth of SAV is dependent upon light availability at the leaf surface, which is a function of light attenuation in the water column and the biomass of epiphytic algae growing on SAV stems. During the last several decades, the loss of SAV has accelerated due to anthropogenic pressures (Kennish et al., 2016) or natural causes such as storms (Hamberg et al., 2017). One of the dominant factors of SAV loss is eutrophication through nutrient loading, exemplified by increased phytoplankton growth and epiphytic growth on vegetation. This results in a reduction of light availability (Burkholder et al., 2007), causing a loss of SAV habitat (Cabello-Pasini et al., 2003;Short and Neckles, 1999).
The complex interactions between light availability, nutrient loading, SAV dynamics, hydrodynamics, and sediment transport can be investigated using numerical modeling tools. Few modeling efforts have attempted to couple dynamically changes the SAV properties (stem density and height). The growth of SAV is modeled as biomass which includes the aboveground biomass (stems and shoots), belowground (roots and rhizomes) biomass, and epiphyte biomass. Individual biomass equations described in the implementation of SAV growth model (Sect. 2.2) are based on previous SAV biomass models (primarily Madden and Kemp, 1996), some of which have been previously implemented to simulate growth conditions for SAV in three-dimensional numerical model simulations (e.g., Cerco and Moore, 2001). The change in biomass leads to a spatial and temporal variation of SAV density and height. With the inclusion of the SAV growth model, SAV can experience growth or dieback while contributing and sequestering nutrients from the water column (modifying the biological environment) and subsequently affect the hydrodynamics and sediment transport (modifying the physical environment). Conversely, a change in the physical environment, for instance, the amount of sediment in the water column, can decrease light availability and cause SAV dieback, leading to reduced wave attenuation, increased sediment resuspension, and a further decrease of light availability.
We demonstrate the two-way biophysical coupling framework as follows: the SAV growth model and integration into COAWST are discussed in Sect. 2; in Sect. 3, the model setup for the idealized domain and a realistic simulation of West Falmouth Harbor, Massachusetts, are described; in Sect. 4, we present the results from the two model configurations along with a discussion of limitations of the current modeling work; and in Sect. 5, we summarize our work and outline areas of future research.

Inclusion of SAV effect on flow and sediment
dynamics in the numerical model Beudin et al. (2017) implemented the parameterizations that accounted for the presence of SAV within a coupled hydrodynamic and wave model within the open-source COAWST numerical modeling system . The COAWST framework utilizes ROMS for hydrodynamics with a wave model (SWAN) via the Model Coupling Toolkit (MCT), generating a single executable program . ROMS is a three-dimensional, free surface, finitedifference, terrain-following model that solves the Reynoldsaveraged Navier-Stokes (RANS) equations using the hydrostatic and Boussinesq assumptions (Haidvogel et al., 2008). The transport of turbulent kinetic energy and generic length scale is computed with a generic (GLS) two-equation turbulence model. SWAN is a third-generation spectral wave model based on the action balance equation (Booij et al., 1999). In ROMS, the presence of SAV extracts momentum, adds wave-induced streaming, and generates turbulence dis-sipation. Similarly, the wave dissipation due to vegetation modifies the source term of the action balance equation in SWAN. All of these subgrid-scale parameterizations account for changes due to vegetation in the water column extending from the bottom layer to the height of the vegetation. SWAN only accounts for wave dissipation due to vegetation at the bottom layer. The coupling between the two models occurs with an exchange of water level and depth-averaged velocities from ROMS to SWAN and wave fields from SWAN to ROMS after a desired number of time steps. The vegetation properties are separately input in the two models at the beginning of the simulations. Through these changes, the SAV can affect the bottom stress calculations that determine the resuspension and transport of sediment, providing a feedback loop between SAV-sediment dynamics-hydrodynamics and wave dynamics. To account for sediment dynamics, CSTMS (Warner et al., 2010) is used to track the transport of suspended-sediment and bed load transport under the action of current and wave-current forcing. The model can represent an unlimited number of user-defined sediment classes.

SAV growth model
The SAV growth model is primarily based upon a previous growth model developed and implemented in Chesapeake Bay by Madden and Kemp (1996). The model simulates the temporal dynamics of aboveground biomass (AGB) that consists of stems or shoots, and the belowground biomass (BGB) that consists of roots or rhizomes. In addition to AGB and BGB, epiphytic algal biomass (EPB) is simulated to account for reductions in light availability to plant leaves due to shading of SAV leaves by epiphytes under high nutrient loading conditions. AGB, BGB, and EPB are simulated as total biomass per unit area, with nitrogen as the currency for biomass. Changes in AGB and BGB pools are simulated as a function of primary production and respiration, mortality (e.g., grazing), and nitrogen exchange through the seasonal translocation of nitrogen between roots and shoots. EPB is modeled as a function of primary production, respiration, and mortality. The remaining section describes the source terms that calculate the evolution of AGB, BGB, and EPB, denoted by α, β, and γ , respectively. Tables 1 and 2 describe model variables and input parameters along with their equivalent symbols used in the source code.

Primary production (ρ s )
The primary production of α depends on the maximum potential growth rate (µ s ) and downward deviations from this maximal rate resulting from light (ϕ s ) and nutrient (ϑ s ) availability as ρ s = µ s min(ϕ s , ϑ s ). (1) The maximum potential growth (µ s ) can be described as where λ s is a self-shading parameter that accounts for crowding and self-shading within the SAV canopy, g s accounts for SAV's maximum growth fraction, r s is the active SAV respiration coefficient, T is the temperature in the water column, and T o is the user-defined optimum temperature that allows for species-specific sensitivities to temperature. The self-shading parameter, λ s , used in Eq.
(3) is calculated by setting a maximum aerial biomass of SAV (Madden and Kemp, 1996), thereby making growth rates density dependent, and is defined as where α is the aboveground SAV biomass and λ s,max accounts for the maximal SAV biomass. The availability of photosynthetically active radiation (PAR) represented by mathematical symbol θ for SAV leaves in the bottom cell is simulated using a bio-optical model (Gallegos et al., 2011;del Barrio et al., 2014). While the biooptical model generates predictions of light available across the spectrum within PAR, the light availability (ϕ s ) used to compute primary production (Eq. 1) is obtained through traditional photosynthesis-irradiance (PI) curves based on total PAR used to represent SAV growth responses to light: where l s is the half-saturation for light limitation for SAV, and θ refers to photosynthetically available radiation that is obtained from the bio-optical model. This simplified PI formulation, which has been applied in previous SAV models (Madden and Kemp, 1996;Zaldívar et al., 2009;Jarvis et al., 2014), is applied so that a general and flexible SAV growth response is available for users in a wide variety of environments with different species. More complex approaches can easily be applied (e.g., Zharova et al., 2001;Carr et al., 2012b). The nutrient limitation (ϑ s ) required in Eq. (1) to compute primary production represents the fact that rooted plants can obtain nutrients from both sediments (as in Madden and Kemp, 1996) and the water column and is defined as ϑ s = DIN wc + n s,1 DIN sed n s,2 DIN sed + n s,1 DIN sed , where DIN wc is the dissolved inorganic nitrogen concentration in the water column based on the sum of NH 4 (ammonium) and NO 3 (nitrate) in the water column, DIN sed is the amount of dissolved inorganic nitrogen (DIN = NH 4 +NO 3 ) in the sediment bed layer, and n s,1 is the half-saturation for nutrient limitation for SAV roots.

Respiration
SAV respiration terms are partitioned into active and basal respiration, where the active respiration term represents respiration that is dependent on the photosynthesis rate, and the basal rate represents maintenance respiration rate. The active respiration term is defined as where ρ s is the primary production term (Eq. 1), p s is the maximum fraction of photosynthesis available for respiration, r s is SAV's active respiration coefficient, and T is the temperature in the water column. The aboveground basal respiration term is defined as where c s is the maximum fraction of SAV BGB that is respired, h s is the SAV basal respiration coefficient for both AGB and BGB, and T is the temperature in the water column.

Mortality
The mortality of SAV is computed separately for aboveground and belowground biomass, where aboveground mor-tality (ω s,1 ) accounts for the sloughing of leaves and grazing in combination as where m α is the aboveground SAV mortality rate (sloughing). Belowground mortality, ω s,2 , is a function of temperature and is given as where m β is the belowground SAV mortality rate. Additional terms include that modify the AGB and BGB include the seasonal exchange (translocation) of root material (nitrogen) quantified as a fraction of primary production and the translocation of BGB to AGB which represents the seasonal translocation of nitrogen from roots to stems as the plants initially emerge in spring. Each of these terms is initiated on a specified day of the year (Madden and Kemp, 1996) and can be altered to account for species differences or regional differences in the physiology of particular species. Translocation of nitrogen from leaves to roots/rhizomes (storage) is modeled as a continuous response to SAV primary production (ρ s ) and is given by defining σ s,1 (translocation of aboveground biomass to belowground biomass) as where k 1 is a downward translocation coefficient. The translocation from roots/rhizomes to leaves (upward translocation) is modeled as a simple linear function of belowground biomass (denoted by β) that begins after a user-defined threshold temperature is crossed and is given by defining σ s,2 (translocation of belowground biomass to aboveground biomass) as where k 2 is a upward translocation coefficient. The EPB is computed similarly to SAV biomass by simulating EPB as a function of primary production, respiration, and mortality (e.g., grazing).

Primary production (ρ e )
The primary production of EPB depends on the maximum potential growth rate (µ e ) and a limitation between light (ϕ e ) and nutrient (ϑ e ) availability as ρ e = µ e min(ϕ e , ϑ e ).
The maximum potential growth of EPB (µ e ) can be described as µ e = λ e ϑ e g e exp r e 1.0 where λ e is the self-shading parameter that accounts for spatial limits on the epiphyte population, g e accounts for epiphyte's maximum growth fraction, r e is the active epiphyte respiration coefficient, T is the temperature in the water column, and T e,o is the user-defined optimum temperature that allows for species-specific sensitivities to temperature. λ e is calculated by setting a maximum aerial biomass of EPB, thereby making growth rates density dependent, similar to the SAV growth rate, as where EPB is the epiphyte biomass and λ e,max is the maximum epiphyte biomass.
The light availability (ϕ e ) used to compute primary production (Eq. 10) is obtained through traditional PI curves used to represent epiphyte growth response to light as where l e is the half-saturation for light limitation for epiphytes, and θ is the photosynthetically available radiation obtained from the bio-optical model. The nutrient limitation (ϑ e ) required in Eq. (1) to compute primary production for epiphytes depends only on the nutrients in the water column and is a traditional algal form (e.g., Monod model) given as where DIN wc is the amount of dissolved inorganic nitrogen in the water column, and n e is the half-saturation for nutrient limitation for epiphytes.

Respiration
Epiphyte respiration terms are partitioned into active and basal respiration, where the active respiration term represents respiration that is dependent on the photosynthesis rate, the basal rate represents the maintenance respiration rate. The active respiration term is defined as where ρ e is the primary production term (Eq. 1), p e is the maximum fraction of photosynthesis for epiphytes, r e is the epiphyte's active respiration coefficient, and T is the temperature in the water column. The basal respiration term is defined as

Mortality
The mortality of epiphytes depends on mortality and grazing of algal cells, as well as losses associated with SAV sloughing (which effectively removes epiphytes from a cell). The mortality term is given as a simple linear form: where m e is the epiphyte mortality rate. The loss of epiphyte biomass due to grazing (τ e,1 ) modeled using an Ivlev function can be described as where z e,max is the maximum grazing rate on epiphytes and z e is the grazing coefficient on epiphytes. The reduction of epiphyte biomass due to the SAV sloughing loss is computed as where ω s,1 is the aboveground mortality term described in Eq. (8), t is the time step size in per day units, and α refers to the aboveground biomass.
The AGB computed in the SAV growth model is utilized to obtain SAV shoot height (meters) and stem density (stems m −2 ) to allow for the biomass model (AGB) to be translated into variables input into the SAV-hydrodynamic coupling. The shoot height (l v ) is related to AGB (denoted by α) as The relationship is based on measurements of Zostera marina in Chincoteague Bay and Chesapeake Bay (Fig. 2) but is consistent with relationships for Z. marina determined elsewhere (Krause-Jensen et al., 2000). Other three-dimensional models have used similar formulations (e.g., Cerco and Moore, 2001, for Chesapeake Bay). SAV stem density n v (in stems m −2 ) is computed from a similar empirical formulation based on relationships in Krause-Jensen et al. (2000) and is computed as

Integration of SAV growth model with the water-column biogeochemistry model (BGCM)
The SAV growth model is built to interact dynamically with the water-column biogeochemistry model (BGCM) within the COAWST modeling framework. We utilize one of the existing BGCM models developed by Fennel et al. (2006) that accounts for nutrients (NO 3 , NH 4 ), phytoplankton and zooplankton biomass, and detritus. The BGCM model in the current simulations solved for 12 state variables. The spectral irradiance model that provides the light attenuation in response to chlorophyll, sediment, and colored dissolved organic matter (CDOM) was previously integrated (Gallegos et al., 2011;del Barrio et al., 2014) into the BGCM model. The BGCM model was implemented within the hydrodynamic component of COAWST model, ROMS. ROMS is a threedimensional, free surface, terrain-following numerical model that solves finite-difference approximations of the RANS equations using the hydrostatic and Boussinesq assumptions (Chassignet et al., 2000;Haidvogel et al., 2000). ROMS is discretized in horizontal dimensions with a curvilinear orthogonal Arakawa C grid (Arakawa and Lamb, 1977). Each state variable is calculated based on the tracer transport equation with tracer concentrations calculated at the grid cell centers as follows: where C is the tracer quantity, t is time, x and y are the horizontal coordinates, and z is the vertical coordinates. u and v are the horizontal components of current velocity, with w d being the sinking velocity for tracers such as detritus. v is the turbulent diffusivity coefficient and C source is the tracer source/sink term, which represents the net effects of all sources and sinks in this representation. There are several choices of advection schemes for tracer advection available in COAWST  and in the current simulations, we utilized the multidimensional positive definite advection transport algorithm (MP-DATA) scheme (Smolarkiewicz, 1984) that has been derived from the Lax-Wendroff (LW) family of schemes. The timemarching scheme for tracers involves a predictor-corrector step using the leapfrog-trapezoidal methods. The 3-D tracer equations are solved at a different and shorter time step than the depth-integrated 2-D barotropic equations. The integration between the baroclinic mode and barotropic mode is performed using a split-explicit time step approach McWilliams, 2005, 2009 McWilliams (2005, 2009). The vertical tracer diffusion terms are solved using a fourth-order centered scheme (Shchepetkin and McWilliams, 2005). The vertical advective fluxes are computed using the piecewise parabolic method (Colella and Woodward, 1984). The vertical terms utilize a backwards Euler method for time marching.
The changes in water-column variables (dissolved and particulate nitrogen, dissolved oxygen, dissolved inorganic carbon) due to the SAV growth model occur locally at the bottom cell through the source terms (C source ) that affect six state variables in the BGCM model: NO 3 (nitrate), NH 4 (ammonium), DO (dissolved oxygen), CO 2 (carbon dioxide), LDeN (labile detrital nitrogen), and LDeC (labile detrital carbon). The change in these state variables based on the SAV growth model is as follows: where ∂DIN SAV ∂t is the net impact of SAV and epiphyte growth on water-column nitrogen concentrations, and s f decides the portioning of nutrient uptake between sediment and the water column using a logistic function and is defined as where m f and k f are constants and equal to 0.2 and 15.0, respectively, and DIN wc (dissolved inorganic nitrogen) is calculated as a sum of state variables NH 4 (ammonium) and NO 3 (nitrate) in the water column. If net growth from SAV and epiphytes is negative, the net nitrogen regeneration is realized as NH 4 production in the water column ( ∂NH 4 ∂t = ∂DIN SAV ∂t ). If there is net growth originating from SAV and epiphytes, the associated water-column uptake of DIN is apportioned between NO 3 and NH 4 relative to their availability in the water column via the following equations: ∂LDeN ∂t = ω s,1 + ω e + τ e,1 t (31) ∂LDeC ∂t = ω s,1 + ω e + τ e,1 t.
All the source terms in Eqs. (25) and (27)-(32) are solved using the SAV growth model described in Sect. 2.2. In Eqs. (30) and (32), these terms are converted to moles of carbon from moles of nitrogen assuming a fixed (and user-defined based on local data) C : N ratio in SAV tissue (we assumed a C : N of 30).

Two-way feedback from SAV to hydrodynamics, waves, sediment dynamics, and biogeochemistry
The addition of the SAV growth model leads to the biological evolution of SAV properties based on temperature, light, and nutrient availability. The modeled SAV community exchanges nutrients, detritus, dissolved oxygen, and dissolved inorganic carbon with the water-column BGCM. Changes in SAV biomass and canopy characteristics also impact hydrodynamics, wave dynamics, and sedimentary dynamics (resuspension-transport). By lowering the current speed and attenuation of wave flow, the reduction in bed shear stresses in the vegetation canopy reduces sediment resuspension, thereby altering sediment transport in the model (as described in Sect. 2.1), which uses the feedback to control light availability and, in turn, potential seagrass biomass production. This methodology of including the SAV growth model enables the COAWST framework to have a two-way feedback in hydrodynamic-biological coupling. Figure 1 describes the coupling process between different modules schematically.  by 100 m in horizontal and 0.1 m in vertical (this varies with water level) directions. A rectangular vegetation bed extends 8 km southward from the north boundary of the domain, with a width of 1.8 km, centered in the domain (Fig. 3). The ROMS barotropic and baroclinic time steps are 0.05 and 1 s, respectively. The bed roughness is set to z o = 1.5 mm. The k − ε turbulence model is implemented following the GLS method (Warner et al., 2005). The initial AGB, BGB, and EPB in the vegetation bed are set to be 90, 15, and 0.01 mmol N m −2 , respectively. The vegetation density, height, diameter, and thickness are initialized to be 400 stems m −2 , 0.19 m, 1.0 mm, and 0.1 mm, respectively. The vegetative drag coefficient (C D ) is set to be 1 (typical value for a cylinder at a high Reynolds number). The imposed surface wind speed is 3 m s −1 from the north to induce a wave field. The surface air pressure is initialized as 101.3 kPa. The kinematic surface solar shortwave radiation is set to an amplitude of 500.0 W m −2 with a 24 h period. The kinematic surface longwave radiation flux is set to zero (W m −2 ). The surface air temperature varies between 1.5 and 18.5 • C over a yearly period. The surface so-lar downwelling spectral irradiance just beneath the sea surface is set following Gregg and Carder (1990). The cloud fraction is set to be zero. The bulk flux parameterizations in COAWST for surface wind stress and surface heat flux are based on the Coupled Ocean-Atmosphere Response Experiment (COARE) code (Fairall et al., 1996a, b;Liu et al., 1979). The model is forced by oscillating the water level on the northern boundary with a tidal amplitude of 0.25 m and a period of 12 h. Northern boundary conditions include a water temperature variation between 1.5 and 18.5 • C over a yearly period. Salinity and NO 3 at the northern boundary are set to 35 psu and 20 mmol N m −3 , respectively, and we impose a suspended sediment concentration of 0.5 g L −1 as well. The northern boundary condition for tracers is a radiation condition with nudging on a 6 h timescale. For both flow and tracer fields (physical and biological), the western and eastern boundaries have a gradient condition and the southern boundary is closed. The model setup for the idealized domain is simulated for 60 d and the model output is averaged over each day.

Realistic test case: West Falmouth Harbor, Massachusetts, USA
In their study, del Barrio et al. (2014) used an offline coupling of the COAWST model with a bio-optical seagrass model (Zimmerman et al., 2003) to study the influence of nitrate loading and sea-level rise on seagrass presence/absence in West Falmouth Harbor, Massachusetts, USA. Nitrate concentrations in groundwater exceeded 200 µM due to a wastewater treatment plant in the watershed; however, recent mitigation is expected to eliminate the nitrate load in the future. The model of del Barrio et al. (2014) used the biogeochemical results to generate spectral irradiance fields which were then passed to the bio-optical model. While useful for investigating the interaction between phytoplankton dynamics, light climate, and potential seagrass coverage, that model did not account for the interaction of seagrass with water column and sediment nitrogen pools or hydrodynamics. Therefore, we tested the fully coupled hydrodynamic, biogeochemical, and vegetation model using the same hydrodynamic and biogeochemical model setup (Ganju et al., 2012;del Barrio et al., 2014) but with the full vegetative interaction implemented. Briefly, the model is forced with tides at the western boundary, groundwater and nitrate loading at the eastern boundary, and solar irradiance at the airsea boundary. Further details on the model setup are given by Ganju et al. (2012) and del Barrio et al. (2014). The hydrodynamic and biogeochemical (e.g., chlorophyll concentrations, light attenuation) results were assessed in those studies. In this work, we test the ability of the coupled model to reproduce the present-day spatial pattern of seagrass presence, with growth and persistence expected in the outer harbor and dieback in the inner harbor, where nitrate loading, phyto- Simulations of the coupled hydrodynamic-biogeochemical-SAV model revealed the integrated nature of estuarine dynamics in response to submerged macrophytes. In these simulations, suspended sediment concentration (SSC) was imposed at the northern open boundary at concentrations of 0.5 g L −1 (and 0 g L −1 within the bed), resulting in a decline in SSC as one moves towards the southern boundary (Fig. 4a). This distribution of SSC input results in an increase in light attenuation (K dpar = 30.0 m −1 ) in the region close to the northern boundary (0.0 km), while background conditions prevail in the southern reaches (Fig. 4b). In Fig. 4b, SSC input from the northern boundary causes a decrease in light availability within the modeled SAV region between the open boundary in the north and about 2.4 km into the SAV bed. Consequently, these suboptimal light conditions in the northern 2.4 km of the SAV bed cause AGB to decrease from its initial value of 90.0 to 30.0 mmol N m −2 (Fig. 5a). Boundary effects associated with SSC inputs are substantially muted in the region between 2.4 and 8.0 km within the SAV bed (Figs. 4 and 5), where in-bed SSCs are much lower than those outside the bed at the same distance from the boundary. As a consequence, AGB biomass increases from its initial value of 90.0 to 150.0 mmol N m −2 over the course of the simulation. Increases in SAV biomass within the bed during the simulation led to increases in SAV density and height, where SAV density increased from its initial value of 400 to 810 stems m −2 due to favorable light conditions from y = 2.4 km to y = 8.0 km. Thus, the model captured the role of SAV in resisting SSC transport into the bed, allowing for greater light availability and an increase in growth rates and biomass accumulation. The temporal evolution of SAV biomass in response to the SSC input at the northern boundary further emphasizes the self-stimulating role of SAV in the idealized simulations. A comparison of model simulations at two locations within the initially described SAV bed of the idealized domain (indicated in Fig. 5a and corresponding to y = 0.1 km and y = 4.5 km from the northern boundary) reveals that close to the northern boundary (y = 0.1 km), the daily averaged light at-  tenuation remains high (above 30 m −1 ) over the 60 d period (Fig. 5a). At y = 0.1 km, the increased light attenuation in the northern location corresponds to the lack of light availability and this causes a decay of AGB from its initial value of 90.0 to 30.0 mmol N m −2 . (Fig. 5b). This decay in AGB over the 60 d period at y = 0.1 km (SAV dieback) contrasts sharply with the AGB increases inside the SAV bed at the southern location (y = 4.5 km), where light attenuation is lower because sediments have not penetrated the SAV bed, allowing for higher SAV growth rates. The higher SAV growth rate in-side the SAV bed at y = 4.5 km can be observed (Fig. 5c) by looking at the net primary production of SAV (ρ s − δ s − ε s ). At this location (y = 4.5 km), the SAV growth rate increases over the 60 d period while it keeps decreasing in the northern location (y = 0.1 km). Due to the higher SAV growth inside the SAV bed (y = 4.5 km), the SSC in the bottom cell remains low (Fig. 5d) and at y = 0.1 km due to the SAV dieback; the sediment concentration in the water column stays high and above 0.25 g L −1 .   As mentioned above, the SSC input on the northern boundary of the idealized domain causes a region of suboptimal light conditions that lead to the SAV dieback, while the SAV growth occurs in the remaining bed where favorable light conditions exist. The effect of change in SAV density and height on the hydrodynamics and morphodynamics at the end of the simulation can be demonstrated by using the same idealized domain. To this end, two transects are chosen that are along the length of the SAV bed and extend from the northern boundary towards the southern boundary. The transects are chosen inside at x = 1.8 km (outside of the SAV bed) and at x = 4.8 km (inside of the SAV bed). The depth-integrated SSC and bottom stresses averaged on the 60th day in the transect (Fig. 7a) outside of the SAV bed show that the profile of bottom stress follows the distribution of SSC along the transect. In Fig. 7a, a peak bottom stress of about 0.2 N m −2 is obtained at x = 1.8 km (outside of the SAV bed) that corresponds to a depth-averaged SSC of 0.31 g L −1 . On the other hand, the transect within the SAV bed (Fig. 7b) shows that the region where SAV dieback has occurred (between 0.0 and 2.4 km) corresponds to increased bottom stresses (0.13 N m −2 at the northernmost location and a corresponding SSC of 0.26 g L −1 ), while in the region where the SAV growth has occurred, the bottom stresses are close to zero (i.e., from 2.4 km onwards).
The simulation of the idealized domain demonstrates the capability of the modeling framework to perform two-way feedbacks between hydrodynamics, sediment, and biological dynamics. The SSC input in the northern boundary affects the light attenuation in the domain and causes SAV dieback close to the northern boundary. The SAV grows in the region where favorable light conditions exist. The SAV dieback leads to increased bottom stresses, while the growth of SAV leads to a decrease in bottom stresses, illustrating the fact that the SAV acts as bottom sediment stabilizer by reducing SSC.

SAV growth in West Falmouth Harbor
The present-day simulation of seagrass dynamics reproduces the patterns of chlorophyll (via phytoplankton), light attenuation, and near-bottom PAR simulated by del Barrio et al. (2014). Nitrate loading from shoreline point sources led to increased phytoplankton growth indicated by increased chlorophyll and light attenuation in the landward, northeast portion of the harbor (Fig. 8a, b), while bathymetric controls in the deeper central basin led to decreased near-bottom PAR (Fig. 8c). Peak AGB exceeds 100 mmol N m −2 , while seagrass presence begins towards decline in the inner harbor and in the central basin as expected. Intertidal areas around the periphery of the harbor are devoid of AGB due to the enforced masking of areas with intermittent wetting and drying.
Time series of these parameters (Fig. 9) from selected outer and inner harbor locations over the first 22 d demonstrate the diurnal variability, as well as the rapid loss of AGB in the inner harbor due to the local nitrate loading, phytoplankton proliferation, and degraded light climate. The sizable diurnal variability in AGB (Fig. 9d) appears to be an artifact of production/respiration formulations that are based on seasonal responses to environmental forcing, rather than diurnal responses to solar irradiance. Future modifications could attenuate this variability by utilizing daily averaged environmental forcing or modifying the frequency of biomass updating.
The modeling framework developed in this work can be used to create hypothetical scenarios to estimate future environmental responses. For example, we ran the model setup of West Falmouth Harbor described in Sect. 3.2 with no nitrate loading to simulate a hypothetical scenario where the groundwater input has no influence from the wastewater treatment plant (non-impacted past or future scenario). The elimination of nitrate loading results in negligible changes in Figure 10. Change in outcomes between impacted and non-impacted scenarios (nitrate loading scenario -no loading scenario). Difference in mean over 22 d of (a) depth-averaged chlorophyll, (b) light attenuation, (c) near-bottom PAR, and (d) peak aboveground biomass on day 14 of the simulation. the outer harbor but greatly reduces chlorophyll and light attenuation in the inner harbor (Fig. 10a, b), while increasing the near-bottom PAR (Fig. 10c). Peak AGB responds to the decreased chlorophyll and increased light attenuation with an increase in the inner harbor (Fig. 10d). This implementation represents an incremental improvement to the prior modeling exercise (Ganju et al., 2012;del Barrio et al., 2014), because the interactions between SAV and the nitrogen pools are explicitly accounted for. For example, this model can now be used to test how changes in seagrass coverage influence nitrogen retention within the estuary or export to the coastal ocean. Further, the introduction of seagrass kinetics will allow for investigation of water-column oxygen budgets with and without seagrass under present and future scenarios.

Model evaluation in West Falmouth Harbor
In order to qualitatively evaluate the seagrass growth model, we have compared the modeled results with observations by del Barrio et al. (2014) that measured the extent of seagrass coverage in West Falmouth Harbor (red outline in Fig. 11). The field data are only available for the northern region of West Falmouth Harbor where the model-data comparisons are performed. The model results are compared by extracting the peak AGB on 14th day of the simulation and normalized with the initial aboveground biomass. The ratio of AGB/AGB initial is considered as a representative of seagrass growth. We assume that for AGB/AGB initial > 1, there is a potential for seagrass growth, and for AGB/AGB initial < 1, the conditions are unfavorable for seagrass growth. In Fig. 11, the model and field data show an 89 % agreement Figure 11. Modeled AGB/AGB initial (aboveground biomass) distribution compared with field data showing seagrass coverage extent (solid red line). Values of AGB/AGB initial > 1 represent seagrass growth potential and those below 1 indicate potential seagrass decline on day 14 of the simulation.
to determine the seagrass growth or dieback. The western region of outer harbor shows seagrass growth potential and agrees with the extent that the seagrass coverage is observed. In the eastern region, the field data show no seagrass coverage and the model also predicts potential seagrass dieback. The model predicts seagrass dieback because of nitrate loading from shoreline point sources that leads to increased chlorophyll and light attenuation (Fig. 8a, b). The model and observations do not compare well in the central basin of outer harbor where the model shows seagrass dieback potential, while the field data show the presence of seagrass. In the central basin, the field data show the presence of seagrass while its density remains low in this region. On the other hand, the modeled seagrass suffers dieback due to the bathymetric controls in the deeper central basin (decreased near-bottom PAR; Fig. 8c).
Direct estimates of aboveground SAV biomass have also been recently made in West Falmouth Harbor (Hayn et al., unpublished data). Although these measurements were not made during the same year as our simulations (measurements in 2006, 2007, 2013; model 2010), the mean aboveground biomass measured in the outer harbor of 49.5 (21 June-6 July 2006), 45.3 (6-19 June 2007), and 41.5 g C m −2 (15-19 July 2013) is consistent with the range of model simulations during a comparable period (2-19 July) in the outer (28.1 to 51.1 g C m −2 ) and middle (14.9 to 37.4 g C m −2 ) harbors. The 2-19 July model range of 45.7 to 156.3 mmol N m −2 across the middle and outer harbors is also consistent with annual mean Z. marina biomass (10-88 mmol N m −2 ) reported in nearby shallow systems on Cape Cod (Hauxwell et al., 2003) assuming a literature-based average that aboveground SAV biomass is 1.5 % N. The range in the model is computed based on the minimum and maximum values of AGB during the 18 d simulation period.

Limitations of SAV growth model and future work
While this modeling approach represents an advance in modeling coupled biophysical processes in estuaries, there are limitations that must be addressed in future work: 1. The modeling of SAV dieback/growth scenarios may require long-term simulations on decadal timescales (Carr et al., 2018). However, the short model time step limits the duration of such simulations. The time step size is of the order of seconds (typical of 3-D ocean models) and this, combined with the fact that the presence of SAV in the hydrodynamic model further limits time step size (due to hydrodynamic stability constraints), overall limits the applicability of the model to be utilized from monthly to annual timescales at this juncture.
2. The biomass equations described in Sect. 2.3 are formulated for seasonal timescales and are being used in the model implementation at every ocean model time step. This leads to large daily variations in above-and belowground biomass that likely do not occur in the environment, although diel variations on SAV growth have been measured in situ (Kemp et al., 1987). Hence, with the current formulations, the output from the biomass model needs to be analyzed as a daily averaged quantity.
3. The current implementation of the SAV growth model is limited to only one SAV species. However, it should be extended to include multiple SAV species to investigate competition under variable salinity and to make the model applicable to a wider variety of locations.

Conclusions
The present study adds to the open-source COAWST modeling framework by implementing a SAV growth model. Based on the change in SAV (aboveground, belowground) biomass and epiphyte biomass, SAV density and height evolve in time and space and directly couple to threedimensional water-column biogeochemical, hydrodynamic, and sediment transport models. SAV biomass is computed from temperature, nutrient loading, and light predictions obtained from coupled hydrodynamics (temperature), biogeochemistry (nutrients), and bio-optical (light) models. In exchange, the growth of SAV sequesters or contributes nutrients from the water column and sediment layers. The presence of SAV modulates current and wave attenuation and consequently affects modeled sediment transport and fate. The resulting modeling framework provides a two-way coupled SAV-biogeochemistry-hydrodynamic and morphodynamic model. This allows for the simulation of the dynamic growth and mortality of SAV in coastal environments in response to changes in light and nutrient availability, including SAV impacts on sediment transport and nutrient, carbon, and