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.

1 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; Hayn et al., 2014). 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 the effects of hydrodynamics and light availability to model the growth of SAV. Everett et al. (2007) and Hipsey and Hamilton (2008) coupled the effects of chlorophyll and water to account for SAV variability, while Bissett et al. (1999a, b) used spectral underwater irradiance to model the light availability required for SAV growth. Carr et al. (2012a, b) developed a one-dimensional coupled hydrodynamics, sediment, and vegetation growth dynamics model. The model solved for vertical 1-D dynamics of SAV growth through a change in biomass that depended on water temperature, irradiance, and seagrass properties. Ganju et al. (2012) used a three-dimensional circulation model (Regional Ocean Modeling System; ROMS) coupled to a nutrient–phytoplankton–zooplankton–detritus (NPZD) eutrophication (water-column biogeochemistry model) developed by Fennel et al. (2006) and integrated the spectral light attenuation formulation (bio-optical model) provided by Gallegos et al. (2011). These models were linked to a benthic seagrass model to calculate seagrass distribution (Zimmerman, 2003) and applied on the temperate estuary of West Falmouth Harbor (del Barrio et al., 2014). While the model was able to capture the loss of SAV due to insufficient light, it did not include interactions with epiphytes or exchanges with water-column nutrients and gas pools. The hydrodynamic feedbacks (changes in currents and waves) and morphodynamic changes (sediment distribution) due to presence of SAV were also ignored. While these dynamic processes have significant implications for coastal ecosystem resilience, numerical models that allow for the two-way feedbacks between hydrodynamics, sediment transport, and SAV growth and nutrient cycling have generally been lacking.

Recently, Beudin et al. (2017) implemented the physical effects of SAV in a vertically varying water column through momentum extraction and vertical mixing as well as accounting for wave dissipation due to vegetation. These processes were implemented within the open-source Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) modeling system that couples the hydrodynamic model (ROMS), the wave model (Simulating WAves Nearshore; SWAN), and the Community Sediment-Transport Modeling System (CSTMS) (Warner et al., 2010). Through this effort, the COAWST framework accounted for the coupled seagrass–hydrodynamics interactions. The model reproduced the turbulent shear stress across the canopy interface and peaked at the top of the canopy similar to the observations of Ghisalberti and Nepf (2004, 2006). The presence of a seagrass patch led to a reduced shear-scale turbulence within the canopy and enhanced wake-scale generated turbulence. For more details on the impact of seagrass on hydrodynamics, readers are referred to Beudin et al. (2017). The inclusion of the physical effects of SAV on flow and sediment dynamics (Beudin et al., 2017) in COAWST allows us to develop a framework that results in dynamic growth of SAV using the temperature, nutrient loading, and light availability in the water column. Therefore, in this work, we implement a SAV growth model that 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.

2 Methods

2.1 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 (Warner et al., 2008). The COAWST framework utilizes ROMS for hydrodynamics with a wave model (SWAN) via the Model Coupling Toolkit (MCT), generating a single executable program (Warner et al., 2008). ROMS is a three-dimensional, free surface, finite-difference, terrain-following model that solves the Reynolds-averaged 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 dissipation. 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.

2.2 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.

Table 1 SAV growth model variable descriptions.

Table 2 SAV growth model input parameter descriptions and their values.

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.

2.2.1 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

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{\mathrm{s}}={\mathit{\mu }}_{\mathrm{s}}\min ({\mathit{\phi }}_{\mathrm{s}},{\mathit{\vartheta }}_{\mathrm{s}}).\end{array}$$

The maximum potential growth (μ_s) can be described as

$$\begin{array}{}\text{(2)}& {\mathit{\mu }}_{\mathrm{s}}={\mathit{\lambda }}_{\mathrm{s}}{\mathit{\vartheta }}_{\mathrm{s}}{g}_{\mathrm{s}}\exp \left[r_{\mathrm{s}}\left({\displaystyle \frac{\mathrm{1.0}}{T-T_{\mathrm{o}}}}\right)\right],\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\mathit{\lambda }}_{\mathrm{s}}=\mathrm{1}-{\left({\displaystyle \frac{\mathit{\alpha }}{{\mathit{\lambda }}_{\mathrm{s},\max }}}\right)}^{\mathrm{2}},\end{array}$$

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 bio-optical 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:

$$\begin{array}{}\text{(4)}& {\mathit{\phi }}_{\mathrm{s}}={\displaystyle \frac{\mathit{\theta }}{l_{\mathrm{s}}+\mathit{\theta }}},\end{array}$$

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

$$\begin{array}{}\text{(5)}& {\mathit{\vartheta }}_{\mathrm{s}}={\text{DIN}}_{\mathrm{wc}}+{\displaystyle \frac{n_{\mathrm{s},\mathrm{1}}{\text{DIN}}_{\mathrm{sed}}}{n_{\mathrm{s},\mathrm{2}}{\text{DIN}}_{\mathrm{sed}}+n_{\mathrm{s},\mathrm{1}}{\text{DIN}}_{\mathrm{sed}}}},\end{array}$$

where DIN_wc is the dissolved inorganic nitrogen concentration in the water column based on the sum of NH₄ (ammonium) and NO₃ (nitrate) in the water column, DIN_sed is the amount of dissolved inorganic nitrogen ($\text{DIN}={\mathrm{NH}}_{\mathrm{4}}+{\mathrm{NO}}_{\mathrm{3}}$) in the sediment bed layer, and n_s,1 is the half-saturation for nutrient limitation for SAV roots.

2.2.2 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

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}_{\mathrm{s}}={\mathit{\rho }}_{\mathrm{s}}{p}_{\mathrm{s}}\exp \left(r_{\mathrm{s}}T\right),\end{array}$$

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

$$\begin{array}{}\text{(7)}& {\mathit{\epsilon }}_{\mathrm{s}}=c_{\mathrm{s}}\exp \left(h_{\mathrm{s}}T\right),\end{array}$$

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.

2.2.3 Mortality

The mortality of SAV is computed separately for aboveground and belowground biomass, where aboveground mortality (ω_s,1) accounts for the sloughing of leaves and grazing in combination as

$$\begin{array}{}\text{(8)}& {\mathit{\omega }}_{\mathrm{s},\mathrm{1}}=(m_{\mathit{\alpha }}\mathit{\alpha }{)}^{\mathrm{2}},\end{array}$$

where m_α is the aboveground SAV mortality rate (sloughing).

Belowground mortality, ω_s,2, is a function of temperature and is given as

$$\begin{array}{}\text{(9)}& {\mathit{\omega }}_{\mathrm{s},\mathrm{2}}=\mathrm{0.01}\mathit{\beta }\exp \left(m_{\mathit{\beta }}T\right),\end{array}$$

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

$$\begin{array}{}\text{(10)}& {\mathit{\sigma }}_{\mathrm{s},\mathrm{1}}={\mathit{\rho }}_{\mathrm{s}}{k}_{\mathrm{1}},\end{array}$$

where k₁ 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

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_{\mathrm{s},\mathrm{2}}=\mathit{\beta }{k}_{\mathrm{2}},\end{array}$$

where k₂ 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).

2.2.4 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

$$\begin{array}{}\text{(12)}& {\mathit{\rho }}_{\mathrm{e}}={\mathit{\mu }}_{\mathrm{e}}\min ({\mathit{\phi }}_{\mathrm{e}},{\mathit{\vartheta }}_{\mathrm{e}}).\end{array}$$

The maximum potential growth of EPB (μ_e) can be described as

$$\begin{array}{}\text{(13)}& {\mathit{\mu }}_{\mathrm{e}}={\mathit{\lambda }}_{\mathrm{e}}{\mathit{\vartheta }}_{\mathrm{e}}{g}_{\mathrm{e}}\exp \left[r_{\mathrm{e}}\left({\displaystyle \frac{\mathrm{1.0}}{T-T_{\mathrm{e},\mathrm{o}}}}\right)\right],\end{array}$$

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

$$\begin{array}{}\text{(14)}& {\mathit{\lambda }}_{\mathrm{e}}=\mathrm{1}-{\left({\displaystyle \frac{\text{EPB}}{{\mathit{\lambda }}_{\mathrm{e},\max }}}\right)}^{\mathrm{2}},\end{array}$$

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

$$\begin{array}{}\text{(15)}& {\mathit{\phi }}_{\mathrm{e}}={\displaystyle \frac{\mathit{\theta }}{l_{\mathrm{e}}+\mathit{\theta }}},\end{array}$$

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

$$\begin{array}{}\text{(16)}& {\mathit{\vartheta }}_{\mathrm{e}}={\displaystyle \frac{n_{\mathrm{e}}{\text{DIN}}_{\mathrm{wc}}}{n_{\mathrm{e}}{\text{DIN}}_{\mathrm{wc}}+n_{\mathrm{e}}}},\end{array}$$

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.

2.2.5 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

$$\begin{array}{}\text{(17)}& {\mathit{\epsilon }}_{\mathrm{e},\mathrm{1}}={\mathit{\rho }}_{\mathrm{e}}{p}_{\mathrm{e}}\exp \left(r_{\mathrm{e}}T\right),\end{array}$$

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

$$\begin{array}{}\text{(18)}& {\mathit{\epsilon }}_{\mathrm{e},\mathrm{2}}=c_{\mathrm{e}}\exp \left(h_{\mathrm{e}}T\right).\end{array}$$

2.2.6 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:

$$\begin{array}{}\text{(19)}& {\mathit{\omega }}_{\mathrm{e}}=m_{\mathrm{e}}\mathit{\gamma },\end{array}$$

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

$$\begin{array}{}\text{(20)}& {\mathit{\tau }}_{\mathrm{e},\mathrm{1}}=z_{\mathrm{e},\max }[\mathrm{1.0}-\exp \left(-z_{\mathrm{e}}\right)],\end{array}$$

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

$$\begin{array}{}\text{(21)}& {\mathit{\tau }}_{\mathrm{e},\mathrm{2}}=\left({\displaystyle \frac{{\mathit{\omega }}_{\mathrm{s},\mathrm{1}}t}{\mathit{\alpha }}}\right),\end{array}$$

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⁻²) 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

$$\begin{array}{}\text{(22)}& l_{\mathrm{v}}=\mathrm{0.45}\left({\displaystyle \frac{\mathit{\alpha }}{\mathrm{120.0}+\mathit{\alpha }}}\right).\end{array}$$

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⁻²) is computed from a similar empirical formulation based on relationships in Krause-Jensen et al. (2000) and is computed as

$$\begin{array}{}\text{(23)}& n_{\mathrm{v}}=\mathrm{4.45}\mathit{\alpha }.\end{array}$$

2.3 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₃, NH₄), 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 three-dimensional, 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:

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial C}{\partial t}}+u{\displaystyle \frac{\partial C}{\partial x}}+v{\displaystyle \frac{\partial C}{\partial y}}+w_{\mathrm{d}}{\displaystyle \frac{\partial C}{\partial z}}={\displaystyle \frac{\partial }{\partial z}}\left(v{\displaystyle \frac{\partial C}{\partial z}}\right)+C_{\mathrm{source}},\end{array}$$

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 (Kalra et al., 2019) and in the current simulations, we utilized the multidimensional positive definite advection transport algorithm (MPDATA) scheme (Smolarkiewicz, 1984) that has been derived from the Lax–Wendroff (LW) family of schemes. The time-marching 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 (Shchepetkin and McWilliams, 2005, 2009). The predictor step calculates the tracer values that updates the momentum equations at an intermediate time step. At that point, the split-explicit algorithm is executed and the update of tracers is done using the corrector step after the new values of velocity are available. For more details of this algorithm, readers are referred to Shchepetkin and 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₃ (nitrate), NH₄ (ammonium), DO (dissolved oxygen), CO₂ (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:

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\partial {\text{DIN}}_{\mathrm{SAV}}}{\partial t}}=({\mathit{\delta }}_{\mathrm{s}}+{\mathit{\epsilon }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{s}})(\mathrm{1}-s_{\mathrm{f}})t+({\mathit{\epsilon }}_{\mathrm{e},\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{e},\mathrm{2}}-{\mathit{\rho }}_{\mathrm{e}})t,\end{array}$$

where $\frac{\partial {\text{DIN}}_{\mathrm{SAV}}}{\partial 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

$$\begin{array}{}\text{(26)}& s_{\mathrm{f}}=\mathrm{1}-\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp [-m_{\mathrm{f}}\left({\text{DIN}}_{\mathrm{wc}}-k_{\mathrm{f}}\right)]}}\right),\end{array}$$

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₄ (ammonium) and NO₃ (nitrate) in the water column. If net growth from SAV and epiphytes is negative, the net nitrogen regeneration is realized as NH₄ production in the water column ($\frac{\partial {\mathrm{NH}}_{\mathrm{4}}}{\partial t}=\frac{\partial {\text{DIN}}_{\mathrm{SAV}}}{\partial t})$. If there is net growth originating from SAV and epiphytes, the associated water-column uptake of DIN is apportioned between NO₃ and NH₄ relative to their availability in the water column via the following equations:

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial {\mathrm{NH}}_{\mathrm{4}}}{\partial t}}=\left({\displaystyle \frac{\partial {\text{DIN}}_{\mathrm{SAV}}}{\partial t}}\right)\left({\displaystyle \frac{{\mathrm{NH}}_{\mathrm{4}}}{{\text{DIN}}_{\mathrm{wc}}}}\right)\text{(28)}& {\displaystyle \frac{\partial {\mathrm{NO}}_{\mathrm{3}}}{\partial t}}=\left({\displaystyle \frac{\partial {\text{DIN}}_{\mathrm{SAV}}}{\partial t}}\right)\left({\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{{\text{DIN}}_{\mathrm{wc}}}}\right)\text{(29)}& {\displaystyle \frac{\partial \text{DO}}{\partial t}}=\left({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\delta }}_{\mathrm{s}}-{\mathit{\epsilon }}_{\mathrm{s}}+{\mathit{\rho }}_{\mathrm{e}}-{\mathit{\epsilon }}_{\mathrm{e},\mathrm{1}}-{\mathit{\epsilon }}_{\mathrm{e},\mathrm{2}}\right)t\text{(30)}& {\displaystyle \frac{\partial {\mathrm{CO}}_{\mathrm{2}}}{\partial t}}=\left({\mathit{\delta }}_{\mathrm{s}}+{\mathit{\epsilon }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{s}}+{\mathit{\epsilon }}_{\mathrm{e},\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{e},\mathrm{2}}-{\mathit{\rho }}_{\mathrm{e}}\right)t\text{(31)}& {\displaystyle \frac{\partial \text{LDeN}}{\partial t}}=\left({\mathit{\omega }}_{\mathrm{s},\mathrm{1}}+{\mathit{\omega }}_{\mathrm{e}}+{\mathit{\tau }}_{\mathrm{e},\mathrm{1}}\right)t\text{(32)}& {\displaystyle \frac{\partial \text{LDeC}}{\partial t}}=\left({\mathit{\omega }}_{\mathrm{s},\mathrm{1}}+{\mathit{\omega }}_{\mathrm{e}}+{\mathit{\tau }}_{\mathrm{e},\mathrm{1}}\right)t.\end{array}$$

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).

2.4 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.

Figure 1 Schematic showing the coupling of SAV growth module implementation in the COAWST model.

Figure 2 Empirical relationships between aboveground biomass and SAV shoot height for Zostera marina populations in polyhaline regions of Chesapeake Bay and Chincoteague Bay. DW refers to dry weight. Data are from Moore (2004) and Ganju et al. (2018).

Figure 3 Planform view of the idealized test domain simulation.

3 Model setup

3.1 Idealized test case

The implementation of the SAV growth model within the COAWST framework is first tested on an idealized domain. The test case consists of an idealized rectangular domain of 9.2 km width and 9.8 km length with a 1 m deep basin. The number of interior domain points is 90 in the x direction and 98 in the y direction with 10 vertical sigma layers. The resulting domain has a grid resolution of 100 m 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⁻², respectively. The vegetation density, height, diameter, and thickness are initialized to be 400 stems m⁻², 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⁻¹ 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⁻² with a 24 h period. The kinematic surface longwave radiation flux is set to zero (W m⁻²). The surface air temperature varies between 1.5 and 18.5 ^∘C over a yearly period. The surface solar 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₃ at the northern boundary are set to 35 psu and 20 mmol N m⁻³, respectively, and we impose a suspended sediment concentration of 0.5 g L⁻¹ 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.

3.2 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 air–sea 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, phytoplankton growth, and light attenuation are highest. The initial SAV properties include a plant height of 0.195 m, plant density of 110 stems m⁻², plant diameter of 0.001 m, and plant thickness of 0.0001 m. The vegetative drag coefficients C_D in the flow model and the wave model are set to 1 (typical value for a cylinder at a high Reynolds number). We utilize the SAV growth model parameters described in Table 1. The model setup for West Falmouth Harbor (Sect. 3.2) is simulated for 56 d beginning 2 July 2010 (Ganju et al., 2012).

Figure 4 Planform view of (a) depth-integrated suspended sediment concentration (SSC) and (b) light attenuation averaged over the last day of the simulation in the idealized domain.

Figure 5 Planform view of (a) aboveground biomass and (b) vegetation stem density averaged over the last day of the simulation in the idealized domain. The red dot and blue triangle represent two points that are located at 0.1 and 4.5 km into the SAV bed, respectively.

Figure 6 Time series of (a) light attenuation, (b) aboveground biomass, (c) net primary production of SAV (${\mathit{\rho }}_{\mathrm{s}}-{\mathit{\delta }}_{\mathrm{s}}-{\mathit{\epsilon }}_{\mathrm{s}}$), and (d) SSC in the bottom cell averaged every day from the two locations identified in Fig. 5a.

Figure 7 Magnitude of bottom stress (a) and depth-integrated SSC (b) at the end of the simulation plotted along the y axis of the idealized domain at two locations, including one outside (x=1.8 km; panel a) and one inside the SAV bed (x=4.8 km, panel b).

4 Results and discussion

4.1 SAV, sediment, and hydrodynamics in the idealized test case

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⁻¹ (and 0 g L⁻¹ 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⁻¹) 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⁻² (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⁻² 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⁻² 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.

Figure 8 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 red circle indicates the outer harbor (a, c) and the blue triangle indicates the inner harbor (b, d) points for time-series data in Fig. 9.

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 attenuation remains high (above 30 m⁻¹) 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⁻². (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 inside the SAV bed at y=4.5 km can be observed (Fig. 5c) by looking at the net primary production of SAV (${\mathit{\rho }}_{\mathrm{s}}-{\mathit{\delta }}_{\mathrm{s}}-{\mathit{\epsilon }}_{\mathrm{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⁻¹.

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⁻² is obtained at x=1.8 km (outside of the SAV bed) that corresponds to a depth-averaged SSC of 0.31 g L⁻¹. 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⁻² at the northernmost location and a corresponding SSC of 0.26 g L⁻¹), while in the region where the SAV growth has occurred, the bottom stresses are close to zero (i.e., from 2.4 km onwards).

Figure 9 Time series of (a) chlorophyll, (b) light attenuation, (c) near-bottom PAR, and (d) aboveground biomass from outer and inner harbor locations identified in Fig. 6.

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.

4.2 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⁻², 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.

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.

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 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.

Figure 11 Modeled AGB∕AGB_initial (aboveground biomass) distribution compared with field data showing seagrass coverage extent (solid red line). Values of $\text{AGB}/{\text{AGB}}_{\mathrm{initial}}>\mathrm{1}$ represent seagrass growth potential and those below 1 indicate potential seagrass decline on day 14 of the simulation.

4.3 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 $\text{AGB}/{\text{AGB}}_{\mathrm{initial}}>\mathrm{1}$, there is a potential for seagrass growth, and for $\text{AGB}/{\text{AGB}}_{\mathrm{initial}}<\mathrm{1}$, the conditions are unfavorable for seagrass growth. In Fig. 11, the model and field data show an 89 % agreement 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⁻² (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⁻²) and middle (14.9 to 37.4 g C m⁻²) harbors. The 2–19 July model range of 45.7 to 156.3 mmol N m⁻² across the middle and outer harbors is also consistent with annual mean Z. marina biomass (10–88 mmol N m⁻²) 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.

4.4 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.

5 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 three-dimensional 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 oxygen cycling. The implementation of the model is successfully tested in an idealized domain where the introduction of sediment in the water column (SSC) at one end of the domain provides suboptimal light conditions that cause SAV dieback in that region. The model was applied to the temperate estuary of West Falmouth Harbor, where simulations show the coupled effect of enhanced nitrate loading in the inner harbor leading to poor light conditions for the SAV to grow, thus modeling the physical effect of eutrophication leading to the loss of a SAV habitat. Among other applications, in the future, the model will be used assess the effects of sea-level-rise scenarios that limit light availability and potentially cause the loss of SAV habitats.