Multi-Scale Inﬂuence of Flexible Submerged Aquatic Vegetation (SAV) on Estuarine Hydrodynamics

: Bottom friction is an important process in coastal and estuarine environments because it can reduce wave heights and moderate tidal currents. When modeling large systems with spatially varying hydraulic properties, bottom friction values are commonly derived from land use classiﬁcation products. However, estimation of bottom friction for vegetated areas can be more challenging due to the complicated and time-varying geometry of the roughness elements. This is particularly true of ﬂexible, buoyant submerged aquatic vegetation (SAV) species, such as seagrasses and kelps, that deform under waves and currents. In this study we incorporate a dynamic friction model that includes the temporal variation in SAV drag forces into a depth-integrated coupled circulation-wave model. In vegetated areas, the bottom friction is continuously updated based on plant geometry, water depth, and combined wave-current velocities. Taking a multi-scale approach, we use the model to investigate the impact of SAV dynamics on both the localized and the integrated bay-wide hydrodynamics of a riverine and tidally inﬂuenced estuary. First, we investigate SAV modiﬁcation of velocity ﬁelds and its implications for sediment transport and circulation pathways. Then, we show how SAV can modify tidal behavior throughout the estuary.


Introduction
The confluence of tides, wind, waves, and streamflow make estuaries one of the most hydrodynamically unique ecosystems [1]. The combination of terrestrial and oceanic forcings produces a wide array of commercially valuable services, including navigation, sustenance, and recreation [2,3]. Nearly all of these services depend on the balance between oceanic and terrestrial hydrodynamics, which dictates fluid and sediment fluxes. The ecological landscape of an estuary has also been shown to influence estuarine dynamics, such as reduction in tidal currents [4][5][6], tidal prism modification [7], as well as increased sediment retention [8]. Thus, quantifying changes to estuarine services under current and future climatologies requires accurate prediction of flow modification by vegetation canopies.
Restoration of natural and nature-based features (NNBFs) for ecological benefits and/or engineering services further necessitates understanding of the broader influence of vegetation on estuarine-scale hydro-and sediment dynamics. Due to the spatial scale variation between the study environment (O ∼ km) and the variability in vegetation features of interest (O ∼ cm-m), quantifying the sensitivity of estuarine hydrodynamics to vegetation conditions remains a challenge, especially for the case of buoyant, flexible submerged aquatic vegetation (SAV) species. We seek to fill a methodological gap to efficiently and accurately quantify drag forces from flexible SAV in the flow field for regional hydrodynamic computational models.

Estuarine Hydrodynamics
At the oceanic boundary, astronomical tides play a leading order role in controlling estuarine processes. As they propagate inland, the tides interact with increasingly shallower bathymetry, resulting in a combination of amplification (shoaling), reflection off the side walls, and dissipation (amplitude decay) by bottom friction and turbulent mixing [9][10][11]. Bottom friction and nonlinear interactions also produce higher-order harmonic overtides that result in time-averaged tidal asymmetry [12,13]. In shallow and/or long estuaries dominated by bottom friction, flood tide energy dissipates with distance into the estuary and reduces the intensity of the ebbing tide [10,14,15]. Such systems are referred to as flood-dominant and result in landward mean sediment transport [15,16]. In contrast, deeper channels (via reduced drag) and shorter channels (via resonance) tend to be ebbdominant and lead to mean sediment export as tidal energy is amplified [14,15,17]. Flood and ebb dominance can vary throughout a tidal basin depending on local bathymetry and proximity to river discharge [15,18]. Strong riverine discharge at the landward boundary can also de-correlate tide duration asymmetry and peak current asymmetry because riverine discharge opposes flood tide propagation while simultaneously enhancing ebb tide current velocities [16,18,19].
Enhanced bottom friction due to vegetation has been shown to shift tidal asymmetry [20][21][22]. Although increased bottom friction can lead to decreased tidal amplification, there is some dependence on estuarine shape that determines whether a system with enhanced friction remains hypersynchronous (tidal amplification with up-channel propagation) or becomes hyposynchronous (tidal dampening with up-channel propagation) [19]. Reduced tidal current velocities due to bottom roughness may also shift flow pathways and estuarine circulation. Additionally, reduction in tidal and wave-driven velocities by SAV habitats can have impacts at the local scale as well, including a decreased tendency for sediment transport (e.g., reduction in shoreline erosion) and increased shelter for fishery juveniles [23,24].

Vegetation-Induced Bottom Friction
The presence of vegetation in coastal environments reduces the momentum of the ambient flow, predominantly by overcoming a drag force acting on the plant's stems and leaves [25]. Approximating a single plant as a rigid cylinder, the drag across the plant can be formulated using a quadratic drag law 1 2 C d ρbl v U 2 , where C d is the drag coefficient, ρ is the fluid density, U is the flow velocity, and b and l v are the plant's width (or diameter, depending on plant morphology) and height, respectively [26]. Analytical solutions of vegetation-induced drag have been incorporated into coastal circulation models [6,27] and spectral wave models [28][29][30], which has enabled the study of their control on hydrodynamics across regional scales. Formulations of vegetation drag that are a function of plant geometry minimize the need to empirically parameterize drag forces. However, variation in natural plant geometry from the rigid cylinder approximation complicates the determination of C d [26]. Various studies have shown that C d values for vegetation can be predicted from the flow Reynolds number (Re) [31][32][33], the Keulegan-Carpenter number (KC) [34,35], and vegetation canopy geometry [34,36], but multi-factor models of C d for vegetation in combined wave-current flows have not yet converged to a singular approach [37,38]. Even more so than the drag coefficient, accurate calculation of drag forces has been shown to depend strongly on the height and vertical distribution of the vegetation of interest [4,39,40].
In contrast to calculating the vegetation 'sink' term by explicitly modeling drag over individual stems, an alternative approach can be derived using hydraulics, whereby vegetation is interpreted as a form of surface roughness in aggregate. Using this more implicit form of bottom friction to estimate vegetation-induced drag forces is particularly advantageous for depth-averaged circulation models because it can be adapted to all land cover types and it enables large scale parameterization efficiently [8,[41][42][43][44][45]. Mapping of Manning's n bottom roughness coefficients from land use and land cover (LULC) datasets has become a common methodology for determining Manning's n across a regional study environment [46][47][48]. However, the flexibility afforded by landscape-to-roughness mapping comes at a cost of precision and accuracy as bottom roughness can vary considerably at scales smaller than available LULC geospatial products, which also generally do not capture subaquatic features.
While there are significant advantages and disadvantages to both the explicit and implicit methods of modeling momentum loss due to flow-vegetation interactions, neither approach accounts for the variability of flow resistance by SAV due to its flexibility and deflection under waves and currents [49][50][51]. In order to fully capture the dynamic feedback between fluid drag forces and vegetation motion, three-dimensional (3D) or quasi-3D modeling techniques are required. To date, only a handful of community coastal models are available to study 3D vegetation-induced drag, turbulence, and wave streaming effects for deformable explicit vegetation stems [27,52], but they are substantially more computationally intensive than 2D depth-averaged (2DH) models. Furthermore, it is critical for more open-source models to rely on physics-based parameterizations to improve estimation of ecosystem services for large regions and to improve inter-model comparisons [53].
We explore the importance of sub-grid vegetation hydrodynamics to estuarine-and regional-scale circulation modeling. We modify commonly used, open-source circulation and wave models to allow flexible SAV to be parameterized as time-varying bottom friction. We then apply the model to an estuary with large SAV-inhabited tidal flats to assess how sensitive hydrodynamic quantities are to vegetation parameterization schemes. Rather than relying on empirical tuning, our dynamic parameterization scheme pulls directly from leading developments in the field of SAV behavior under waves and currents. We then apply the new method to our study site in Coos Bay, Oregon to demonstrate canopy-scale control on velocity fields, circulation patterns, and tidal characteristics.

Methods
First, we present our estuary of study. Then, we describe in detail the numerical modeling methods and forcing conditions used, paying greatest attention to the novel dynamic friction routine we developed. Finally, we briefly describe the numerous analytical methods used to quantify hydrodynamic statistics of interest from modeled output, as they are well documented in the literature elsewhere. We calculate both local and integrated quantities to illustrate the breadth of influence of the newly developed dynamic friction routine.

Study Site-Coos Bay Estuary
We center our numerical sensitivity study on the Coos Bay Estuary (CBE), located on the Oregon (USA) coast ( Figure 1). The CBE is one of the largest estuaries on the US Pacific coast [54,55]. The bathymetry of the estuary has evolved considerably since the Euro-American displacement of native tribes for the extraction of local resources [56]; work by Eidam et al. [55] and Eidam et al. [57] illustrates how anthropogenic modifications to the CBE's bathymetry have changed tidal and sediment dynamics over time. Presently, the estuary is 44 km 2 in size, approximately half of which is intertidal [54,[58][59][60]. Wide tidal flats support approximately 6.9 km 2 of SAV habitat [24,61], making the CBE one of the most SAV-rich Oregon estuaries and an ideal site for our numerical experiment. Selection of the CBE is also particularly apt due to the fact that it has been the focus of several recent studies [55,57,60,62].
In the CBE, both terrestrial-and oceanic-derived forcings exhibit a strong seasonal signal. The bay is fed by several freshwater sources with seasonally varying flowrates, peaking in winter months and lowest in summer. The predominant source of freshwater is the Coos River with an average winter (December, January, February; DJF) flowrate of 86.6 m 3 /s and average summer (June, July, August; JJA) flowrate of 9.87 m 3 /s. Additional flows from Palouse Slough to the north and Isthmus Slough to the southwest of the Coos River supply between 1.64 m 3 /s and 0.11 m 3 /s (DJF and JJA seasonal averages) and between 3.44 m 3 /s and 0.10 m 3 /s, respectively [63,64]. Winter peaks in streamflow also coincide with stronger offshore wave intensities, with significant wave heights H s of 3.01 m and peak wave periods T p of 12.6 s oriented nearly shore-normal; in the summer, waves arrive northwesterly (NW) with H s of 1.55 m and T p of 9.39 s (NDBC Station 46229).

Numerical Model
In this study, we embed a dynamic bottom friction routine into the open-source finite element coastal circulation model ADCIRC coupled with the spectral wave model SWAN [30,66,67]. The depth-averaged mode of ADCIRC has been used for numerous regional studies of estuarine hydrodynamics as well as for modeling impacts of compound storm events [60,68]. When coupled with SWAN, ADCIRC is highly computationally efficient and parallelizable, making ADCIRC+SWAN an excellent model choice for compoundforced estuary modeling [60,67]. Our study on the sensitivity of modeled velocities and tidal constituents to dynamic bottom friction builds on existing literature on the bottom friction sensitivity of the ADCIRC and coupled ADCIRC+SWAN models [47,48,[69][70][71].
For our numerical sensitivity experiment, we rely on a previous calibration of the computational mesh performed by Parker et al. [60] because it was found to be a faithful representation of the dominant physical drivers of hydrodynamics within the CBE, including the compound interaction of tides, waves, wind, and streamflow (R 2 > 0.97 for a month-long validation period). Topo-bathymetric elevations are sourced from a blend of DEMs from National Oceanic and Atmospheric Administration and LiDAR surveys taken by the Oregon Department of Geology and Mineral Industries. The mesh is made up of about 30,000 computational nodes, and triangular grid cell sizes range from 10 m (small channels) to 150 m throughout the bay itself. Outside of the bay, computational cell size increases to over 3 km at the tidal forcing boundary (located 133 km offshore).

Boundary Conditions
Water surface elevations at the offshore boundary are prescribed by the eight dominant tidal harmonics for the region; harmonic amplitudes, frequencies, and equilibrium arguments are extracted from the TPXO tidal database [72] using the OceanMesh2D MATLAB toolkit [73]. Offshore wave, wind, and streamflow conditions are modeled on seasonally averaged winter (DJF) conditions (Table 1). In addition to tides, the open ocean boundary nodes are forced with a temporally and spatially uniform JONSWAP spectrum parameter-ized by historical wave observations collected at NDBC Station 46229. Spatially uniform wind fields are applied over the domain to maintain wave heights as they propagate from the open boundary towards the bay.
Wind magnitudes and directions were determined from multiple data sources, including observations at the North Bend Airport [74,75] and long-term averaged wind reanalysis products by NCEP NCAR [76]. Wind velocities are assumed to be spatially uniform over the domain but do oscillate diurnally due to their dependence on the solar heat flux. Finally, constant streamflow rates are defined as the DJF average of the median monthly flow rates at each of the riverine boundaries; flow-duration statistics are derived from linear multiple regression analysis of streamflow records between 1906 and 2005 developed by USGS StreamStats [63,64].

Dynamic Bottom Friction
Using time-varying computed current velocities, wave statistics, and model SAV properties, we compute a dynamically varying Manning's n value in SAV habitats, or n SAV , at each time step for every node classified as SAV. A schematic of the ADCIRC+SWAN information exchange and the dynamic friction routine is shown in Figure 2. The arrows between the solid blue boxes in Figure 2 represent the tight coupling between ADCIRC, which calculates depth averaged velocities (U, V) and water surface elevation (η), and SWAN, which passes to ADCIRC wave radiation stresses (S xx , S xy , S yx , S yy ). To calculate the combined wave-current drag on our computational SAV, the time-varying root mean square (rms) wave velocity U w and mean direction φ w are also passed from SWAN to ADCIRC.
The dynamic friction routine inserted in both ADCIRC and SWAN ( Figure 2) uses a physically representative vegetation-induced bottom roughness formulation proposed and validated by Baptist et al. [77]. Developed using data-driven genetic programming, the roughness formulation is trained on results of a k − ε turbulence model of submerged canopy flow [78]. Assuming the bottom roughness of the bare bed is negligible relative to the roughness caused by the vegetation, the Chezy coefficient C is defined as where m is the density of vegetation within the canopy (blades/m 2 ), h is the water depth, g is the acceleration due to gravity, and κ is the von Karman coefficient. Because the formulation for C by Baptist et al. [77] depends on SAV height, we account for nonstationary plant posture by calculating the effective height l e of the deflected plant, which can be predicted from two dimensionless parameters, the Cauchy number (Ca) and the buoyancy parameter (B) [79]. Ca is used to quantify the posture-dependent ratio of drag forces to plant stiffness. The relationship between the effective height and Ca is modulated by B. In summary, these two dimensionless numbers quantify the response of a plant stem or blade to the flow as a function of its total resistance to motion (buoyancy and stiffness). The expression for l e developed and validated by Luhar and Nepf [79] is and the dimensionless parameters Ca and B are defined as where E is the modulus of elasticity, I is the moment of inertia about the rotational axis (equal to bt 3 /12), ρ v is the density of plant material, and t is the stem or blade thickness. For our study, the numerical dynamic vegetation is modeled after the SAV species Zostera marina (Z. marina). Because previous field surveys of Z. marina in the CBE have been limited [24], we defined the geometric and mechanical properties of our SAV mimic using the mean values of properties published in surveys of Z. marina in PNW estuaries [59,[80][81][82][83][84][85][86]. Values of the necessary Z. marina parameters are an additional input to the dynamic friction routine (dashed yellow box in Figure 2); the specific values used for our numerical experiment are shown in Nomenclature.
Although ADCIRC allows the bottom boundary to be parameterized with C, we use the Manning's n formulation, which can be found from the Chezy coefficient via n = h 1/6 /C. By substituting the equation of effective height l e (Equation (2)) for l v in Equation (1), we can translate sub-grid flow-vegetation interactions into a framework easily interpreted by friction-based numerical models. We will rely on the Luhar and Nepf [79] model of plant deflection that is a function of Ca for the combined wave-current flow magnitude, U wc , where U 0 is the near-bed wave orbital velocity equal to √ 2U w and φ c and φ w represent the approach angle of the current and wave velocities, respectively (φ c = atan2(V, U)) [87]. Thus, our approach retains dependence on plant buoyancy as well as the wave environment in the event that U w >> U. Based on Luhar and Nepf [79], we assume that C d for a deflected plant stem is adequately described by the drag over a rigid cylinder with equivalent height of the deflected stem; essentially, once the effective plant height is known, calculations of drag based on a rigid cylinder assumption are valid. In this work, we assume a constant C d of 1.0 [88]. Finally, we impose a maximum n SAV value of 0.10, which prevents n SAV from becoming unrealistically high [89].
Because bottom friction in SWAN is parameterized with the wave bottom friction factor, f w , ADCIRC first converts Manning's n values to an equivalent roughness height, z 0 [67,90]. SWAN converts z 0 to f w and then an equivalent friction coefficient, C b , which modulates the bottom friction dissipation term in the wave action balance [91]. Although SWAN also has an explicit vegetation module that adds an additional sink term based on formulations of drag over rigid cylinders [28][29][30], we are restricting the current work to not include additional vegetation modules at this time to maintain consistency across ADCIRC and SWAN, though its inclusion could be an important element of future research [92]. Our approach differs from other recent advances in vegetation drag for large-scale hydrodynamic modules because we account for both vegetation blade mechanics and deflection under waves and currents. For example, Moki et al. [93] calculate a spatially varying Manning's n that is only dependent on SAV leaf morphology but does not account for plant deflection because their focus is on spatial variability of bottom roughness rather than temporal variability. Zhao and Chen [94] do incorporate a dynamic Manning's n into a wave model, but their formulation is based on an iterative solution for n and the deflected vegetation height l e [95].

Experimental Design
We examine the effect of a dynamic n SAV on the hydrodynamics of the CBE by performing the same 40-day simulation with a number of bottom friction treatments in SAV habitats ( Table 2), each of which is compared with a control simulation where no SAV is included. In the control, SAV habitats are assigned a Manning's n equivalent to that of a bare sediment bed (0.02); outside of SAV habitats for the control and all other simulations, n is assigned based on LULC classification [60,96]. The first SAV parameterization scheme used, referred to as the status quo treatment (SQ), represents the common approach of assigning a single, empirically derived value of n SAV based on LULC classification (equal to 0.035 for vegetated intertidal flats) [41,70,97,98]. The simplicity of the SQ treatment contrasts in complexity with the dynamic parameterization scheme, which applies the dynamic friction routine ( Figure 2) to all nodes within the SAV habitats. This dynamic treatment is used for two different coverages of SAV, one with historical coverage [24] and one with potential coverage [55]. The former is referred to as the dynamic historic (DH) treatment, and the latter, the dynamic potential (DP) treatment, represents dynamic motion of SAV, should SAV be present in all locations designated as potential habitat based on light penetration thresholds (z SAV = −1 m < MLLW < 0.5 m [55]), which nearly doubles modeled SAV habitats from the historical extents. Finally, to isolate the effect of temporal and spatial variability in n SAV on hydrodynamic observations, we re-perform the simulation with time-and space-averaged values of n SAV calculated during the DH simulation. For the static varying (SV) parameterization scheme, the value of n SAV at any given SAV node is equivalent to the time-averaged value of n SAV at the same node during the DH simulation. While the SV treatment of n SAV is spatially varying and temporally static, the static uniform (SU) parameterization is both temporally static and spatially uniform. All SAV nodes are given the identical value of n SAV for the SU treatment, equal to the time-and space-averaged value of n SAV across all SAV nodes in the DH simulation.

Depth-Averaged Hydrodynamics
Depth-averaged horizontal current velocities U, V are reported by ADCIRC at 15-min sampling intervals. Spectral wave statistics calculated by SWAN are also reported at 15-min intervals, including H s , U w , T p , and φ w . Using the simulated wave and current velocity fields, we first calculate the effect of SAV flow attenuation on the likelihood of localized sediment transport. Then, we use the simulated results to calculate the effect of SAV coverage on large-scale estuarine hydrodynamics, specifically tidal distortion and spatial patterns of flow circulation.

Sediment Fluxes
Non-cohesive sediments are predicted to be mobilized when the dimensionless transport number (Shields number, θ) exceeds the critical Shields number θ c , defined as where D * is the dimensionless grain size, d 50 is the median grain size, s is the specific gravity of the material (equal to the ratio of the sediment to fluid density), and ν is the kinematic viscosity of water [99,100]. Spatially varying values of d 50 are interpolated to computational element centroids within the bay using samples collected by Eidam et al. [55]. Areas within the CBE made up of fine bed material (d 50 < 0.06 mm) are excluded from the analysis [99]. The tendency for transport was calculated for all sandy areas of the bay from the simulated current and wave velocities following [99,100], In Equation (8), C b represents the bed friction drag coefficient that varies with bed material type and relative depth [99].

Circulation
Passive tracer methods are commonly used to quantify an estuary's residence time, also known as its turn-over time or flushing time [101,102]. Numerical implementations of particle tracking have also been used to investigate circulation and dispersal as an alternative to in situ drifter deployments, particularly for large estuaries with complex geometries. Advected pathways (latitude, longitude) of neutrally buoyant numerical "drifters" have been shown to be well-correlated with physical observations and field deployments [103][104][105]. We use a similar technique to examine patterns of flow circulation throughout the CBE by tracking the pathlines of randomly dispersed numerical drifters. Simulated current velocities were interpolated spatially (from computational nodes to the instantaneous drifter locations) and temporally (from the 15-min reporting interval to a 1.5-min particle tracking interval) to track the movement of particles throughout the bay using a fourth-order Runge-Kutta method [106]. A total number of 600 particles are used to qualitatively investigate how increased SAV coverage can shift pathways of fluid transport (and thus, potentially sediment, nutrient, or contaminant transport) throughout the CBE.

Tidal Velocity and Duration Asymmetry
Although we are not modeling sediment transport directly, we can qualitatively predict the mean direction of sediment fluxes in and out of the estuary by characterizing the CBE's tidal dynamics. Quantitatively, the ebb or flood dominance of an estuary can be determined from peak velocity asymmetry or from tide duration asymmetry. Peak velocity asymmetry can be calculated as the cubed ratio of the average peak ebb velocity U ebb to the average peak flood velocity U f lood , which can vary throughout the estuary. Tide duration asymmetry is often quantified using the phase difference between the dominant constituent (M2 for the CBE) and its nonlinear overtide (M4) to determine whether the estuary is flood-dominant (2ω M2 − ω M4 = 0 • − 180 • , shorter rising tide) or ebb-dominant (2ω M2 − ω M4 = 180 • − 360 • , shorter falling tide) [10,12,16]. The ratio of dominant to nonlinear overtide amplitudes (e.g., A M4 /A M2 ) indicates the magnitude of tidal asymmetry, where A is the tidal constituent amplitude; A M4 /A M2 increases from a value of 1.0 with increasing tidal asymmetry [10,12,16].

Characterization of Dynamic n SAV
The inclusion of the dynamic Manning's n routine results in temporal and spatial variation of the Manning's n coefficient at the SAV nodes. As shown by the normalized probability distribution functions in Figure 3a, n SAV varied significantly across nodes. For each vegetated node, the value of n SAV varied over the model simulation between roughly 0.04 (mean minimum value) and 0.09 (mean maximum value). For comparison, blue dashed and solid lines at n SAV = 0.02 and 0.035, respectively, are included in Figure 3a to indicate commonly used values for open water (as used for the Control) and the SQ treatment, respectively (Section 2.5) [41,70,97,98]. A tidally phase-averaged value of n SAV was determined for each vegetated node by splitting the time series of the computed dynamic Manning's n by tidal phase. When we then take the mean (across SAV nodes) phase-averaged n SAV and compare it with the phase-averaged water surface elevation η, we see that the two signals are well-correlated ( Figure 3b). Values of n SAV are greatest when the tide stage is at extremes (low tide and high tide) and lowest during mid-tide when tidal current velocities are expected to be highest; the asymmetry in n SAV peak width during high and low tide closely matches the asymmetry in the phase-averaged tidal current magnitude (not shown). The strong dependence of n SAV on local flow velocities is further illustrated in Figure 3c, where the time-averaged mean and standard deviation σ of n SAV for each vegetated node is shown to decrease with increasing rms combined wave-current velocity U wc,rms [87]. Where U wc,rms is greater than approximately 0.2 m/s, the relative contribution of the waves to the total local velocity increases, which indicates that exclusion of waves from the simulation would underestimate the deflection of SAV and thus overestimate its dissipative capacity in areas where wave velocities are nontrivial (Figure 3d). In Figure 3c,d, the relative strength of the waves is represented as the wave kinetic energy fraction KE w /KE wc , which represents the ratio of kinetic energy (proportional to U 2 ) by waves alone KE w to the flow's total kinetic energy by combined waves and currents KE wc . Waves are said to dominate flow energy when values of KE w /KE wc exceed approximately 0.75 [107].

SAV Modification of Wave and Current Velocities
Spatial distributions of the change in combined wave-current velocities ∆U wc during ebb tide due to the inclusion of SAV relative to the control are shown in Figure 4, where blues indicate a decrease in U wc,ebb when vegetation is included and reds indicate an increase. For all treatments, the characteristic U wc,ebb , which is defined as the value of U wc averaged over instances of local minima in dη/dt across all tidal cycles, is shown to decrease inside the areas designated as SAV habitats (outlined in grey) as compared to the control. As the value of n SAV used to parameterize SAV increases, reduction in U wc,ebb becomes larger. When large static values of n SAV are used (SU and SV, Figure 4c,d), distinct areas of increased U wc,ebb are seen in between SAV canopies. While this flow acceleration around canopies is shown for the SQ and DH cases (as well as DP, not shown), the funneling effect is less intense because values of n SAV for those treatments are much lower during peak ebb currents. When tidal velocities are at a maximum, n SAV is as low as 0.04 across much of the estuary for the dynamic treatments (DH and DP), which is very close to the static value of n SAV used for SQ treatment (0.035).

Sediment Fluxes
The change in total area of the sediment bed estimated to be mobilized (Equations (6)-(8)) due to the inclusion of SAV is shown in Figure 5. Inside SAV habitats, the greatest reduction in total area subject to sediment transport results from the SU treatment of SAV, and the smallest reduction results from the SQ treatment (Figure 5a). The DH, DP, and SV treatments reduced likelihood of sediment mobility inside SAV canopies relatively equally. Flow acceleration between SAV canopies during the high-Manning static treatments (SU and SV) is shown to increase the tendency for sediment motion in unvegetated regions of the CBE, as compared to the no-vegetation control, as well as compared to the DH and DP simulations (Figure 5b). In contrast, when SAV is allowed to deflect during the DH and DP treatments, n SAV is reduced during peak flows. The reduction in n SAV prevents excessive flow resistance that accelerates flow between patches and increases θ to above θ c during the SU and SV treatments.

Circulation
In addition to quantifying SAV's effect on localized hydrodynamics, we use a numerical drifter analysis to illustrate the cumulative impact of SAV on flow and circulation patterns in the CBE. Across the broad tidal flats flanking the main channel, numerical drifter pathlines are more narrowly concentrated in between SAV canopies when present (for all SAV treatment types), as opposed to a wider distribution of paths when SAV was not included in the model. Figure 6 shows the discretized pathlines of the drifters overlying elevation contours and SAV node locations in the CBE's South Slough. The regions of the South Slough identified with arrows in Figure 6 demonstrate where SAV canopies can enhance or reroute dominant flow paths. As in the main channel, the presence of SAV along the deepest channel of South Slough narrows particle flow paths (denoted by solid grey arrow). For the DP simulation, SAV habitat in the southern-most tip of South Slough blocks the accumulation of drifters at the same location when SAV is absent during the Control and DH simulations (solid white arrow). The wider SAV coverage in the DP treatment also appears to trap particles outside of the main channel by creating low-velocity pockets (dashed arrows).

SAV Modification of Tidal Dynamics
The presence of SAV was shown to decrease the amplitude of the dominant tidal constituent A M2 and its overtide A M4 (Figure 7), for all SAV treatments excluding the SQ treatment. As shown in Figure 7a, the amplitude of the M2 constituent A M2 was reduced by approximately the same amount for all along-channel distances dx; the DH and DP treatments reduced A M2 by approximately 10 cm, while the reduction by the SV and SU treatments was slightly less at approximately 6 cm. Dissipation of the overtide M4 due to the presence of SAV increased with increasing dx (Figure 7b). As with A M2 , the DH and DP treatments led to a larger reduction in A M4 from the control. The SV and SU treatments also reduced A M4 along the main channel. In contrast, the SQ case led to an increase in A M4 for dx greater than 20 km. Outside the main channel, reduction in A M4 is further enhanced within SAV habitats.
Contours quantifying the relationship between the M2 and the M4 overtide are shown in Figure 8. The spatial variation in A M4 dissipation shown in Figure 7b leads to spatial variability in the nonlinearity factor A M4 /A M2 (Figure 8a-c). The presence of SAV appears to reduce A M4 /A M2 relative to the control, especially within the first 5 km of the mouth of the CBE and in some of the shallowest portions of SAV habitats. The presence of SAV also led to a shift in the M2 constituent phase ω M2 and the M4 constituent phase ω M4 . As compared to ω M2 for the control, the M2 tide arrives 87 • earlier for the SQ treatment, roughly 45 • later for the DH and DP treatments, and 10 • -20 • later for the SU and SV treatments. Shifts in ω M2 were nearly uniform over the CBE. In contrast, shifts in ω M4 vary spatially over the domain; greater lags in ω M4 due to the presence of SAV are observed outside of SAV patches. The spatial variability in ω M4 due to SAV leads to spatial variability in 2ω M2 − ω M4 (Figure 8d-f). As with the A M4 /A M2 , the phase lag 2ω M2 − ω M4 is shown to decrease slightly within the first 5 km, indicating a slight potential for flood dominance within the first third of the main channel, as well as for the SAV-inhabited tidal flats on the eastern side of the main channel.  Following Eidam et al. [55], tidal ebb-flood dominance can also be determined by the cubed ratio of the average peak ebb velocity U ebb to the average peak flood velocity U f lood , where values greater than 1 indicate ebb dominance and values less than 1 indicate flood dominance. Values of U 3 ebb /U 3 f lood in Figure 9c are shown for nodes within the main channel with increasing distance dx from the inlet. Although the presence of vegetation is shown to increase U ebb and U f lood approximately 5-20% as compared to the control between dx = 20-25 km (a and b, respectively), no significant change in U 3 ebb /U 3 f lood is observed in this section of the main channel. Excluding the region of dx = 20-25 km, the presence of SAV is shown to slightly decrease U ebb and U f lood for all treatments but the SQ treatment. Again, these slight modifications do not lead to an appreciable change in U 3 ebb /U 3 f lood .

Discussion
The following subsections are organized thematically to provide a holistic overview of the impact of SAV and the dynamic friction routine on estuarine processes. Flexible SAV is shown to influence hydrodynamics across a range of scales in the CBE. We also discuss the applicability and limitations of our findings, as well as ecological parameter uncertainty.

Spatial and Temporal Dependence of SAV Attenuation
Velocity reduction by SAV varies considerably throughout the CBE. Discontinuous SAV canopies significantly reduce combined wave-current velocities within the canopies themselves, but are shown to create a potential for flow acceleration and sediment mobility in bare areas neighboring and between canopies (Figures 4 and 5). The variation in velocity attenuation is not only due to the discontinuity of SAV habitat, but also the spatial variability in hydrodynamics of the CBE. The horseshoe shape of the CBE, combined with strong riverine, wave, and wind forcing, results in a non-uniform distribution of U wc,rms . The spatial and temporal variation in U wc,rms has important ramifications for accurate estimation of bottom friction by SAV because drag across flexible blades is limited by the intensity of local hydrodynamics. In areas where combined wave-current velocities are high, SAV deflection is also high, resulting in a reduction in n SAV (Figure 3c). Unlike static treatments, dynamic treatments of SAV bottom friction do not require a priori determination of n SAV for habitats across a spatially varying flow field. Furthermore, without including the influence of oceanic and wind waves on the total flow velocity, SAV deflection would be underestimated in unsheltered and near-mouth regions of the CBE (Figure 3d). Hydrodynamic patterns contribute to the spatial variability in velocity attenuation by SAV, and our results support the need for modeling both mean flow and wave statistics to accurately parameterize SAV-induced bottom friction in numerical models.
Because estuarine hydrodynamics additionally evolve on diurnal time scales, bottom friction by SAV also evolves temporally (Figure 3b). Time-varying SAV stem motion by naturally variable flow conditions necessitates the inclusion of dynamic deflection of SAV in numerical bottom friction routines. As shown in Figure 4, excluding SAV dynamics and assuming static behavior can overestimate flow dissipation within SAV canopies, increasing flow velocities adjacent to the canopies by as much as 5 cm/s. This overprediction of funneling around SAV habitats can overestimate the likelihood of sediment transport (Figure 5b). In compound estuaries like the CBE, variable streamflow rates, wave conditions, and wind fields additionally contribute to temporal variability in wave-current velocity fields. It is worth noting that there may be some contexts in which assuming temporal uniformity may be acceptable as long as the spatial variation remains. For example, flow attenuation and likelihood of sediment transport inside SAV canopies was shown to be nearly identical for the DH, DP, and SV treatments (Figure 5a).

Integrated Influence of SAV Attenuation
Particle tracking and tidal constituent analyses give an integrated view of flow attenuation by SAV across discontinuous habitats. As indicated in Figure 6, SAV meadows can lead to stalled flows outside the reach of the main channel. While not a precise metric of residence time, the apparent pooling of numerical drifters in sections of the CBE far from the main channel suggests that SAV can reduce flushing rates in already sheltered areas, which could be important for predicting system turbidity over time. More so than estuarine circulation, tidal amplitudes in particular were shown to be affected by both SAV presence and treatment type (Figure 7). Even in the main channel where depths and velocities are greatest, the M2 tidal constituent amplitude decreased by 8% for the SU and SV treatments and 11% for the DH and DP treatments; the same reduction in A M2 by each treatment was observed uniformly across all areas of the CBE. Estuarine-scale changes in tidal constituent amplitudes decreased nonlinearity A M4 /A M2 near the inlet and within SAV canopies, particularly so for the DH (and DP, not shown) treatments and less so for the SQ treatment ( Figure 8). Interestingly, the SQ treatment led to the greatest increase in area where the phase relationship 2ω M2 − ω M4 decreased towards flood-dominant values. The ability of discontinuous SAV habitats to collectively shift tidal dynamics in the CBE towards flood dominance is consistent with findings by Donatelli et al. [22], meaning that time-averaged hydro-and sediment transport fluxes are directed into the estuary and potentially lead to increased sediment retention.

Relative Role of SAV in the CBE
While it has been shown that enhanced bottom friction dissipates ebb current velocities (enhanced by streamflow) more than flood current velocities (reduced by streamflow) [22], we note that this assumption does not necessarily hold for the case of flexible SAV that may deflect more under stronger ebb velocities than flood velocities. Furthermore, although SAV habitats across the CBE were shown to modify tidal constituents and peak current velocities (Figure 9a,b), the presence and treatment of SAV are not shown to definitively drive large changes in the ratio of modeled observations of ebb and flood velocities (Figure 9c). Thus, neither the presence nor the treatment of SAV appear to play a significant role in controlling ebb and flood dominance in the CBE. The minor role of SAV on the CBE's peak current asymmetry contrasts with the strong influence by changes in estuarine bathymetry [55]. Eidam et al. [55] show that historical bathymetric changes (i.e., dredging of the main channel) have shifted the CBE from flood dominant to strongly ebb dominant. Finally, we note that the effect of SAV on the ratio of ebb and flood velocities is not expected to change under average seasonal streamflow variability. In JJA when streamflow is much lower, tidal duration asymmetry and peak velocity asymmetry are expected to be less [16]. If current velocities are less skewed, SAV deflection and, thus, SAV velocity attenuation, is likely to be fairly uniform under ebbing and flooding tidal periods. As under DJF streamflow conditions, dynamic treatments of SAV are likely to have a similar effect on peak current asymmetry as static treatments.
Despite the lack of significant change in peak current asymmetry by SAV, accretion of sediment within and on the shoreward edges of sheltered canopies (Figures 5 and 6) may still result in net accumulation of material within the estuary even though sediment infilling by tides and mean currents is not expected. Treatment type (static vs. dynamic) may also influence the magnitude and spatial distribution of sediment retention, as static treatments are shown to enhance flow routing between canopies but more strongly reduce velocities within canopies (Figure 4).

Applicability to Other Systems
Seagrasses are found in sheltered, shallow coastal waters across the globe at all but polar latitudes [108,109]. Thus, our investigation into the large-scale control of SAV on hydrodynamics has widespread relevance for estuaries with similar levels of SAV coverage. In the CBE sensitivity study, each of our modeled canopies was between 0.1 and 1 km 2 , totaling nearly 8 km 2 or 17% of the CBE's total area below MSL. SAV presence may have a less significant influence on circulation and tidal dynamics in estuaries and embayments with less extensive coverage. From a management perspective, the restoration or loss of a singular canopy may not have significant impacts on regional dynamics, but largescale die off events or habitat migration are more likely to impact tidal dynamics as well as the spatial distribution and frequency of sediment mobility. From a modeling perspective, a temporally varying dynamic friction routine may not be needed to capture the hydrodynamic impact of smaller habitats on large systems, such as in estuaries where total SAV habitat coverage is relatively small; values of n SAV may be nearly uniform across small, cohesive SAV canopies. However, a gradient in deflected stem height becomes more evident for larger SAV meadows due to the cascading effect of velocity and wave height attenuation with increasing distance into the meadow. Such spatial variation is also important to capture when modeling the localized hydrodynamics and sediment transport in and around vegetation canopies. As canopies and habitats become larger and/or increase in coverage over the hydrodynamic area of interest, the need for adaptive bottom friction treatments becomes greater.
Although not investigated in this study, variability in plant morphology should also be considered in site-specific parameterizations of bottom roughness by SAV, when known. In particular, SAV plant stem density has been shown to have the greatest control on its dissipative capacity [110]. Although variation in plant height l v for flexible species has little impact on its effective height l e and resulting dissipative capacity [110], increased plant stiffness (E and I) increases the sensitivity to l v uncertainty (Section 1.2). As E and I increase, the Cauchy number approaches zero, meaning plant deflection becomes negligible and stems remain static. In this study, modeled time-averaged values of l e are shown to be approximately half of l v throughout most of the estuary, excluding the near-inlet patches where time averaged l e /l v is shown to be closer to 0.2 (co-located with KE w /KE wc ∼ 1, Figure 3d). In very sheltered reaches of the estuary where U wc,rms ∼ 0 (Figure 3c), the friction routine is less sensitive to plant rigidity because l e is approximately equal to l v for the duration of the simulation. Using a bottom friction value representative of rigid vegetation would alter, but not eliminate, the time-varying nature of vegetationinduced bottom friction n because the formulation presented here is a function of plant submergence l e /h(t) [77]; for the same reason, values of n should also vary spatially for domains with spatially varying bathymetry. Thus, some form of adaptive friction routine should still be employed.

Conclusions
SAV species are ecologically valuable NNBFs in estuaries worldwide. Understanding how discontinuous SAV habitats control fluid and sediment transport pathways can be a critical component to planning and siting SAV restoration efforts. As numerical modelers move towards more physically representative parameterizations of SAV, we show that allowing for plant deflection is essential to realistically estimating the effect of flexible, buoyant SAV species on velocities within and adjacent to SAV habitats. Furthermore, estimating plant motion should account for wave-induced mean flow as well as tidally driven currents so as not to overestimate SAV drag forces during slack tide. In the Coos Bay Estuary investigated in this study, neglecting SAV deflection leads to over-attenuation of flow velocities within canopies and enhanced flow velocities in the unvegetated areas adjacent to the canopies. In the deep unvegetated main channel of the CBE, intertidal canopies of flexible SAV are also shown to dissipate tidal constituent energy more than, but reduce peak ebb and flood current velocities less than, static SAV treatments. This work demonstrates a valuable framework that can be utilized in any coupled wave-circulation numerical model to better account for sub-grid SAV canopy dynamics.  Acknowledgments: In addition to support by funding agencies, E.R.H. would like to acknowledge the ADCIRC Users Community for providing invaluable insight in support of model development.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
Symbolic key for variables and modeling constants used in this experiment. Wave friction coefficient m 2 /s 3