Development of a Submerged Aquatic Vegetation Growth Model in a Coupled Wave-Current-Sediment- Transport Modeling System (COAWST v3.4)

0 Abstract The coupled biophysical interactions between submerged aquatic vegetation (SAV), hydrodynamics (currents and 10 waves), sediment dynamics, and nutrient cycling have long been of interest in estuarine environments. Recent observational studies have addressed feedbacks between SAV meadows, current velocity, sedimentation, and nutrient cycling and suggest SAV are ecosystem engineers whose growth can be self-reinforcing. 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 model COAWST. In this study, we extend the 15 COAWST modelling framework to account for dynamic changes of SAV and associated epiphyte biomass. Modelled SAV biomass is represented as a function of temperature, light, and nutrient availability and 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 water column and sediment bed properties, hydrodynamics, and sediment transport (i.e., a two-way feedback). We demonstrate the 20 behavior of these modelled processes through application to an idealized domain, then apply the model to a eutrophic harbour where SAV dieback is a result of anthropogenic nitrate loading and eutrophication. These cases demonstrate an advance in the deterministic modelling of coupled bio-physical processes and will further our understanding of future ecosystem change.

1 Introduction 25 Submerged aquatic vegetation (SAV), or seagrasses, are rooted vascular plants that inhabit sediments of estuaries and coastal waters, with a wide global distribution.SAV are important primary producers in shallow environments, provide habitat for a number of aquatic organisms, can slow water velocities and dampen wave energy to trap particulate material (Carr et al., 2004), and can alter biogeochemical cycles through oxygenation of sediments (Larkum et al., 2006).The positive impact of ecosystem services provided by SAV presence has been 30 well-studied (Hemminga and Duarte, 2000, Nixon et al., 2001, Terrados and Borum, 2004. and 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 owing 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 andNeckles, 1999).
The complex interactions between light availability, nutrient loading, SAV dynamics, hydrodynamics, and 5 sediment transport can be investigated using numerical modelling tools.Few modelling efforts have attempted to couple the effects of hydrodynamics and light availability to model the growth of SAV.Everett et al., 2007;Hipsey and Hamilton, 2008 coupled the effects of chlorophyll and water to account for SAV variability while Bissett et al., 1999aBissett et al., , 1999b used spectral underwater irradiance to model the light availability required for SAV growth.Ganju et al., 2012 used a three-dimensional circulation model (ROMS) coupled to a Nutrient Phytoplankton Zooplankton 10 Detritus (NPZD) eutrophication (water column bio-geochemistry 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 et al., 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 15 nutrient and gas pools.The hydrodynamic feedbacks (change 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 20 column through momentum extraction, vertical mixing as well as accounting for wave dissipation due to vegetation.These processes were implemented within the open source COAWST (Coupled-Ocean-Atmospheric-Wave-Sediment Transport) modelling system that couples the hydrodynamic model (ROMS), the wave model (SWAN) and the Community Sediment Transport Modelling System (CSTMS) (Warner et al., 2010).The SAV modelling was based on modifications to the flow field resulting from three-dimensional drag, in-canopy wave-induced 25 streaming, and production of turbulent kinetic energy in the hydrodynamics model (ROMS), along with energy dissipation and resultant hydrodynamic feedback from the wave model (SWAN).The presence of SAV in the COAWST model modifies the wave characteristics primarily by wave energy dissipation resulting from the work done by drag force on SAV stems and secondarily by influencing the water level and current fields.The change in current fields alter the bed friction, thus affecting the sediment dynamics.The inclusion of the physical effects of 30 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 related to biomass that includes the above ground (stems and shoots), below ground (roots and rhizomes) biomass and epiphyte biomass.The biomass equations are based upon previous 35 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 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 section 2; in section 3, the model setup for the idealized domain and a realistic simulation of West Falmouth Harbor, MA are described; in section 4, we present the results from the two model configurations along with a discussion of limitations of the current modelling work and in section 5, we 10 summarize our work and outline areas of future research.

Inclusion of SAV effect on flow and sediment dynamics in the numerical model
The integration of the model to account for the effects of SAV on flow within the COAWST numerical 15 modelling system (Warner et al., 2008) have been described in Beudin et al. (2017) and summarized in Kalra et al. (2017).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, the Community Sediment Transport Modelling System (CSTMS) (Warner et al., 2010) is used to track the transport of suspended-sediment and bed load transport under the action of 20 current and wave-current forcing.The model can represent an unlimited number of user defined sediment classes.

Implementation of SAV growth model
With the inclusion of SAV affecting the flow and sediment dynamics in the COAWST model, in this work we extend the modelling framework and implement coupling with an SAV growth model.The SAV growth model is 25 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 above ground biomass (AGB) that consists of stems or shoots, and the below ground 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 30 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 are modelled as a function of primary production, respiration, and mortality.The remaining section describes the model equations used to simulate AGB, BGB and EPB and how they were implemented within a three-dimensional framework.The default input parameters 35 required by the following model equations are described in Table 1. (1) The maximum potential growth () can be described as: 5 where   is a self-shading parameter that accounts for crowding and self-shading within the SAV canopy,  accounts for SAV's maximum growth fraction,  is the active SAV respiration coefficient,  is the temperature in water column,   is the user defined optimum temperature that allows for species-specific sensitivities to temperature.The self-shading parameter,   used in Eq. 3 is calculated by setting a maximum aerial biomass of 10 SAV (Madden and Kemp 1996), thereby making growth rates density-dependent and is defined as: where  is the above ground SAV biomass and  − accounts for the maximal SAV biomass.
The availability of photosynthetically active radiation () for SAV leaves in the bottom cell is simulated using a 15 bio-optical model (Gallegos et al. 2009, del Barrio et al. 2014).While the bio-optical model generates predictions of light available across the spectrum within PAR, the light availability (  )used to compute primary production (Eq. 1) is obtained through traditional photosynthesis-irradiance curves based on total PAR used to represent SAV growth responses to light: where  is the half-saturation for light limitation for SAV and  refers to photosynthetically available radiation that is obtained from the bio-optical model.
The nutrient limitation (  ) 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: 25 where   is the dissolved inorganic nitrogen concentration in the water column,   is the amount of dissolved inorganic nitrogen in the sediment bed layer, and   is the half-saturation for nutrient limitation for SAV roots.The active respiration term is defined as: where   is the primary production term (Eq.1),  is the maximum fraction of photosynthesis available for respiration,  is the SAV's active respiration coefficient,  is the temperature in the water column.
The basal respiration term is defined as: 5 where  is the maximum fraction of SAV below ground biomass (BGB) that is respired,  is the SAV basal respiration coefficient for both AGB and BGB,  is the temperature in the water column.

Mortality:
The mortality of SAV is computed separately for above-ground and below-ground biomass, where 10 AGB mortality accounts for the sloughing of leaves and grazing in combination as: where   is the above ground SAV mortality rate (sloughing).
Below ground mortality,   , is a function of temperature and is given as: where   is the below-ground 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 20 species differences or regional differences in the physiology of particular species.
The epiphyte biomass (EPB) is computed similarly to SAV biomass by simulating EPB as a function of primary production, respiration, and mortality (e.g., grazing).

Primary production (𝑝𝑝 𝐸𝑃𝐵 ):
The primary production of EPB depends on the maximum potential growth rate 25 (  ) and a limitation between light (  ) and nutrient (  ) availability, as: The maximum potential growth of EPB (  ) can be described as: where   is the self-shading parameter that accounts for spatial limits on the epiphyte population,   accounts for epiphyte's maximum growth fraction,   is the  is the temperature in water column,  − is the user defined optimum temperature that allows for species-specific sensitivities to temperature.  is calculated by setting a maximum aerial biomass of EPB, thereby making growth rates density-dependent similar to the SAV growth rate, as: 35 where  is the epiphyte biomass and  − is the maximum epiphyte biomass.
The light availability (  ) used to compute primary production (Eq.10) is obtained through traditional photosynthesis-irradiance curves used to represent epiphyte growth response to light, as: 5 where   is the half-saturation for light limitation for epiphytes and  is the photosynthetically available radiation obtained from the bio-optical model.
The nutrient limitation (  ) 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: 10 where   is the amount of dissolved inorganic nitrogen in the water column,   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 15 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   is the primary production term (Eq.1),   is the maximum fraction of photosynthesis for 20 epiphytes,   is the epiphyte's active respiration coefficient,  is the temperature in the water column.
The basal respiration term is defined as: where  is the maximum fraction of epiphyte biomass that is respired,  is the epiphyte basal respiration and  is the temperature in the water column.25 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: where   is the epiphyte mortality rate.
The loss of epiphyte biomass due to grazing (grz  ) modelled using an Ivlev function can be described as: where  − is the maximum grazing rate on epiphytes and   is the grazing coefficient on epiphytes.
The reduction of epiphyte biomass due to the SAV sloughing loss is computed as: where   is the above ground mortality term described in Eq. 8, is the time step size in per day units and AGB is the above ground biomass.
The above ground biomass (AGB) computed in the SAV growth model is utilized to obtain SAV shoot height 5 (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 (  ) is related to AGB 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-10 dimensional models have used similar formulations (e.g., Cerco and Moore, 2001 for Chesapeake Bay).and CDOM was previously integrated (Gallegos et al. 2009, del Barrio et al. 2014) into the BGCM model.Along with the light attenuation model, the effects of algal respiration, seagrass kinetics and diel oxygen dynamics were also added to BGCM model.The SAV growth model described in Section 2.2 interacts dynamically with BGCM model to simulate SAV growth.25

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 modelled SAV community exchanges nutrients.detritus, dissolved oxygen, and dissolved inorganic carbon with the water-column BGCM.Changes in SAV biomass, and canopy characteristics also impacts hydrodynamics, wave dynamics and sedimentary dynamics (resuspension-transport) (as described in Section 30 2.1), that 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 between hydrodynamic-biological coupling.Figure 1 describes the coupling process between different modules schematically.a width of 1.8 km, centered in the domain (Figure 3).The ROMS barotropic and baroclinic time steps are 0.05 s and 1 s respectively.The bed roughness is set to   =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 (CD) is set to be 1 (typical value for a 10 cylinder at high Reynolds number).The imposed surface wind speed is 3 m/s from the north to induce a wave field.

Idealized test case
The surface air pressure is initialized as 101.3 kPa.The kinematic surface solar shortwave radiation is set to an amplitude of 500.0W/m 2 with a 24-hour period.The kinematic surface longwave radiation flux is set to zero (W/m 2 ).
The surface air temperature varies between 1.5 ℃ to 18.5 ℃ over an 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 15 be zero.The bulk flux parameterizations in COAWST for surface wind stress and surface heat flux are based on the COARE code (Fairall et al. (1996a(Fairall et al. ( , 1996b) ) and 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 hours.Northern boundary conditions include a water temperature variation between 1.5 ℃ to 18.5℃ over an yearly period.Salinity and NO3 at the northern boundary are set to 35 psu and 20 mmol N/m 3 respectively, 20 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 6h 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 days and the model output is averaged over each day.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 30 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 35 (Ganju et al., 2012 anddel 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 The vegetative drag coefficients   in the flow model and the wave model are set to 1 (typical value for a cylinder at high Reynolds number).We utilize the SAV growth model parameters described in Table 1.The model setup for West Falmouth Harbor (Section 3.2) is simulated for 56 days, beginning 2 July 2010 (Ganju et al., 2012).10

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, SSC was imposed at the northern open boundary at concentrations of 0.5 g/L (and zero g/L within the bed), resulting in a decline in SSC as one moves 15 towards the southern boundary (Fig. 4a).This distribution of SSC input results in an increase in light attenuation (Kdpar=30.0m -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 modelled SAV region between the open boundary in the north and about 2.4 km into the SAV bed.
Consequently, these sub-optimal light conditions in the northern 2.4 km of the SAV bed cause AGB to decrease from 20 its initial value of 90.0 millimoles N/m 2 to 0.0 millimoles N/m 2 (Fig. 5a).Boundary effects associated with SSC inputs are substantially muted in the region between 2.4 km and 8.0 km within the SAV bed (Figs.4&5), where in-bed SSC concentrations are much lower than those outside the bed at the same distance from the boundary.As a consequence, where AGB biomass increases from its initial value of 90.0 millimoles N/m 2 to 320.0 millimoles 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 25 height, where SAV density increased from its initial value of 400 stems/m 2 to of 1400 stems/m 2 owing to favourable 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 30 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) reveal that close to the northern boundary (y=0.1 km), the daily averaged light attenuation remains high (above 30 m -1 ) over the 60-day 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 millimoles N/m 2 to 0.0 millimoles N/m 2 .(Fig. 5b).This decay in AGB over the 60-day 35 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 (  −   −   ).At this location (y=4.5 km), the SAV growth rate increases over the 60-day period while it remains close to zero 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. 5 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 favourable 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 10 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 60 th 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 0.49 N/m 2 of peak bottom stress is obtained at x=1.8 km (outside of the SAV bed) that corresponds to a depth-averaged SSC of 0.4 g/L.On the other hand, the transect within the SAV bed (Fig. 7b) shows that the region where SAV dieback 15 has occurred (between 0.0 km to 2.4 km) corresponds to increased bottom stresses (0.45 N/m 2 at the north most location) while the region where the SAV growth has occurred, the bottom stresses are close to zero (i.e. from 2.4 km and onwards).
The simulation of the idealized domain demonstrates the capability of the modelling framework to perform two-way feedbacks between hydrodynamics, sediment and biological dynamics.The SSC input in the northern 20 boundary affects the light attenuation in the domain and causes SAV dieback close to the northern boundary.The SAV grows in the region where favourable 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 act as bottom sediment stabilizers by reducing SSC. 25

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 30 near-bottom PAR (Fig. 8c).Peak AGB exceeds 150 millimoles N m -2 , while seagrass is nearly eliminated in the inner harbor, and completely eliminated 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 days demonstrate the diurnal variability, as well as the rapid loss of AGB in the inner harbor due to the local nitrate loading, 35 phytoplankton proliferation, and degraded light climate.The sizeable diurnal variability in AGB (Fig. 9d) appears to be an artifact of production/respiration formulations that are based on seasonal responses to environmental forcing, The modelling 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 section 3.2 with no nitrate loading, to simulate a hypothetical scenario where the groundwater input has no influence from the 5 wastewater treatment plant (unimpacted 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, as well as smaller increases in the outer harbor (Fig. 10d).This implementation represents an incremental improvement to the prior modelling exercise (Ganju et al., 2012 anddel 10 Barrio et al., 2014), because the interaction 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.

Limitations of SAV growth model and Future Work
While this modelling approach represents an advance in modelling coupled biophysical processes in estuaries, there are limitations that must be addressed in future work: 1.The modelling of SAV dieback/growth scenarios may require long-term simulations on decadal timescales.However, the short model time step limits the duration of such simulations.The time step size is of the order of 20 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 time scales at this juncture.
2. The biomass equations described in Section 2.3 are formulated for seasonal time scales and are being used in the model implementation at every ocean model time step.This leads to large daily variations in above and below 25 ground biomass that do not likely 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 30 applicable to a wider variety of locations.

Conclusions
The present study adds to the open source COAWST modelling framework by implementing a SAV growth model.Based on the change in SAV biomass (above ground, below ground) and epiphyte biomass, SAV density and 35 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 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 sub-optimal light conditions that causes 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 10 harbour leading to poor light conditions for the SAV to grow; thus modelling the physical effect of eutrophication leading to the loss of SAV habitat.Among other applications, in 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 habitat.

Code availability 15
The implementation of the SAV growth model has been implemented in the Coupled Ocean Atmosphere Waves .COAWST is an open-source community modeling system maintained by John C. Warner (jcwarner@usgs.gov)and distributed through the USGS code archival repository.It is available for download on https://code.usgs.gov/coawstmodel/COAWST.The COAWST distribution files contain source code derived from ROMS, SWAN, WRF, MCT and SCRIP, along with the Matlab code, examples and a User's 20 Manual.

Data availability
The model solutions can be obtained here: http://geoport.whoi.edu/thredds/catalog/clay/usgs/users/tkalra/journal_veg_growth/catalog.html 25 All the files can be accessed through the HTTPserver on this link.The input files are available on the above link for the idealized domain simulation in the folder: inputfiles_idealized_domain.zip and for the West Falmouth Harbor simulation in the folder: inputfiles_wfh_domain.zip.

30 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.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-271Manuscript under review for journal Geosci.Model Dev. Discussion started: 29 March 2019 c Author(s) 2019.CC BY 4.0 License. −ℎ = ( SAV stem density   , (in stems/m 2 ) is computed from a similar empirical formulation based on relationships inKrause-Jensen et al., 2000 and is computed as:  = 4.45   (21)152.3Interactions of SAV with Water-Column Biogeochemistry Model (BGCM model)The SAV growth model is built to interact dynamically with the water-column biogeochemistry model within the COAWST modelling framework.We utilize one of the existing biogeochemical models (BGCM model) developed byFennel et al., 2006 that accounts for nutrients (NO3, NH4), phytoplankton and zooplankton biomass, and detritus.The spectral irradiance model that provides the light attenuation in response to chlorophyll, sediment, 20 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-271Manuscript under review for journal Geosci.Model Dev. Discussion started: 29 March 2019 c Author(s) 2019.CC BY 4.0 License.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 are 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 the vertical (this varies with water level).A rectangular vegetation bed extends from the north boundary of the domain southward 8 km, with 5

25 3. 2
Realistic test case: West Falmouth Harbor, Massachusetts, USA del Barrio et al. ( Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-271Manuscript under review for journal Geosci.Model Dev. Discussion started: 29 March 2019 c Author(s) 2019.CC BY 4.0 License.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 5 plant height of 0.195 m, plant density of 110 stems/m 2 , plant diameter of 0.001 m, and plant thickness of 0.0001 m.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-271Manuscript under review for journal Geosci.Model Dev. Discussion started: 29 March 2019 c Author(s) 2019.CC BY 4.0 License.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.

15
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-271Manuscript under review for journal Geosci.Model Dev. Discussion started: 29 March 2019 c Author(s) 2019.CC BY 4.0 License.light predictions obtained from coupled hydrodynamics (temperature), bio-geochemistry (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 modelled sediment transport and fate.The resulting modelling framework provides a two-way coupled SAV-biogeochemistryhydrodynamic and morphodynamic model.This allows for the simulation of the dynamic growth and mortality of 5

Figure 2 :
Figure 2: Empirical relationships between above ground biomass and SAV shoot height for Z. marina populations in polyhaline regions of Chesapeake Bay and Chincoteague Bay.Data from Moore et al. 2004 and Ganju et al. 2018.

Figure 3 :
Figure 3: Planform view of the idealized test domain simulation.

Figure
Figure (a)

Figure 5 :
Figure 5: Planform view of (a) above ground biomass and (b) vegetation stem density averaged over the last 5 day of the simulation in the idealized domain.

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

Figure 10 .
Figure 10.Change in outcomes between impacted and non-impacted scenario (nitrate loading scenariono loading scenario).Difference in mean over 22 days of (a) depth-averaged chlorophyll, (b) light attenuation, (c) near-bottom PAR, and (d) peak above ground biomass at day 14 of the simulation.