Flow Structure in an Artificial Seagrass Meadow in Combined Wave-Current Conditions

Experiments were conducted in a laboratory flume using an artificial seagrass meadow, modeled after Zostera marina, to examine the impact of waves on the vertical structure of time-averaged current, Reynolds stress, and turbulent kinetic energy (TKE) under combined wave-current conditions. With the addition of smaller waves, defined by a ratio of wave velocity to current velocity Uw/Uc < 2.5, the time-averaged velocity peaked above the meadow, which was similar to pure current conditions. When Uw/Uc > 2.5, the presence of waves caused the time-averaged velocity to peak near the top of the meadow. For Uw/Uc > 1 the presence of waves reduced the magnitude of peak Reynolds stress. For all conditions considered, the wake production of turbulence dominated the shear production of turbulence in the meadow. However, the wave velocity was less efficient than the current velocity in generating TKE in the meadow because the movement of the blades forced by the oscillatory fluid motion reduced the relative velocity between the blades and the wave. A modified hybrid model for wake production of TKE in a flexible canopy under combined wave-current conditions was proposed to account for the relative contributions of waves and currents. Wake production of TKE was dominated by waves when Uw/Uc > 1 and dominated by currents when Uw/Uc < 1. The models and observations proposed in this study contribute to an enhanced understanding of the relative influences of waves and currents on seagrass meadow flow structure in realistic combined wave-current conditions.


INTRODUCTION
Meadows of submerged aquatic vegetation, such as seagrass, can damp wave energy (e.g., Knutson et al., 1982;Fonseca and Cahalan, 1992), reduce erosion, and improve water quality (e.g., Ginsburg and Lowenstam, 1958;Ward et al., 1984;Moore, 2004). The reduction of current velocity within the meadow enhances the creation of near-bed habitat relative to unvegetated regions (e.g., Fonseca et al., 1982;Homziak et al., 1982). These ecosystem services depend on the interactions between the meadow, waves, and current. For example, hydrodynamic intensity impacts seagrass nutrient uptake (e.g., Lei and Nepf, 2016;Gillis et al., 2017) and overall seagrass survival (Peralta et al., 2006;Fonseca et al., 2007). Furthermore, the initiation of sediment transport within vegetated regions has been linked to exceeding thresholds of vegetation-generated turbulence in pure current (Yang et al., 2016) and pure wave (Tinoco and Coco, 2018;Tang et al., 2019;Zhang and Nepf, 2019) conditions. This paper describes mean and turbulent flow within submerged aquatic vegetation under combined wavecurrent conditions, which can provide a deeper understanding of meadow hydrodynamics and enable a more accurate assessment of ecosystem services.
Under unidirectional flow, vegetation drag reduces velocity within the meadow and redirects some flow over the top of the meadow, creating a shear layer with an inflection point in the vertical profile of time-averaged velocity (Brunet et al., 1994;Raupach et al., 1996;Ghisalberti and Nepf, 2006). The shear layer at the top of the meadow, as well as the wakes of individual shoots and leaves, can convert mean kinetic energy into turbulent kinetic energy (e.g., Nepf, 1999;Nepf and Vivoni, 2000;Tanino and Nepf, 2008). In pure wave conditions, stem wake turbulence is generated in proportion to the wave velocity squared (Zhang et al., 2018;Tang et al., 2019;Zhang and Nepf, 2019), but shear layer turbulence is only generated for long waves for which the wave excursion significantly exceeds the drag length scale of the canopy (Ghisalberti and Schlosser, 2013).
Turbulence has been studied in pure current and pure wave conditions separately, but waves and currents can coexist in seagrass meadows (e.g., Koch et al., 2006). Multiple studies have considered turbulence in combined wave-current conditions in rigid and flexible meadows of model vegetation. Both Lou et al. (2018) and Chen et al. (2020) observed that the addition of waves increased turbulent kinetic energy (TKE) within dense arrays of rigid cylinders relative to pure currents. In a flexible meadow, the individual plants may deflect under currents and sway under waves, changing the structure of time-averaged and turbulent velocities. For example, in pure current conditions, oscillations in the meadow height (known as monami) decrease the sharpness of the drag interface at the top of the canopy and the magnitude of peak Reynolds stress, relative to stationary canopies (Ghisalberti and Nepf, 2006). The presence of waves can cause a mean pronation in the direction of wave propagation (e.g., Zhang et al., 2018). Paul and Gillis (2015) subjected a transplanted meadow of Zostera noltei to both pure current and combined wave-current conditions and observed no differences in timeaveraged velocity, but observed a reduction in TKE in combined wave-current conditions relative to corresponding pure current conditions, which is opposite to the trend observed for rigid models (e.g., Lou et al., 2018;Chen et al., 2020). In a field study on a floodplain of the Yangtze River and for multiple species of flexible vegetation, Zhang et al. (2021) observed little impact of waves on the magnitude of TKE associated with the current, but observed that the presence of waves increased the shear layer penetration depth, δ e , relative to that predicted for pure current , and they attributed this to the waving of the flexible leaves.
The primary goal of this study was to observe the influence of wave amplitude on time-mean flow structure, Reynolds stress, and TKE in combined wave-current conditions within a flexible meadow, including an evaluation of TKE prediction in combined wave-current conditions. Understanding how the combination of waves and currents impacts flow structure and turbulence in a seagrass meadow can improve the description of seagrass meadow hydrodynamics and the ecosystem services those conditions facilitate.

THEORETICAL BACKGROUND FOR PREDICTING VELOCITY AND TURBULENCE IN MEADOW
Several studies have described unidirectional flow over a rigid submerged meadow using a two-layer model: an overflow layer, with depth-and time-averaged velocity U 2 , and a canopy layer of height h, with depth-and time-averaged velocity U 1 (e.g., Huthoff et al., 2007;Cheng, 2011;Chen et al., 2013; see Figure 1). This was modified for a flexible meadow by Lei and Nepf (2021), who incorporated the impact of plant reconfiguration by using the deflected canopy height, h d , and the effective blade length, l e , defined in Luhar andNepf (2011, 2013) as the length of rigid blade that provides the equivalent drag to a flexible blade of length l. Specifically, for a depth-averaged velocity U c in total water depth D, the in-canopy velocity is in which φ is the solid volume fraction and a v is the frontal area per canopy volume. C D is the canopy drag coefficient (see Supplementary Section 1 for a dictionary of symbols). C is a coefficient representing the efficiency of turbulent momentum transfer between the canopy and overflow layers: in which K c is an empirical coefficient. For rigid canopies, Chen et al. (2013) found K c = 0.07 ± 0.02. The penetration length scale δ e describes the vertical distance into the canopy over which turbulent momentum flux is significant, and it is defined as the distance from the top of the canopy at which the Reynolds stress decays to 10% of the peak magnitude. The penetration length depends on the canopy height and drag length scale (C D a v ) −1 , specifically δ e = min 0.23 ± 0.06 Nepf and Vivoni, 2000;Konings et al., 2012). Shear is highest near the top of the canopy, such that the Reynolds shear stress u w reaches a peak magnitude at the canopy height (Figure 1), with u and w the velocity fluctuations in the streamwise and vertical directions, respectively, and the overbar denoting the timeaveraging operation. The peak magnitude of Reynolds stress, assumed to occur at the top of the canopy, scales with the velocity difference between the overflow and canopy layers (Konings et al., 2012;Chen et al., 2013): in which ρ is the density of water.
FIGURE 1 | Simplified sketch of time-averaged velocity profiles (denoted by arrows) and Reynolds stress u w profile in a submerged canopy (with canopy elements represented by vertical gray rectangles). The dotted vertical lines within and above the canopy in the time-averaged velocity profiles denote the spatial vertically averaged time-averaged velocity within and above the canopy (U 1 in the canopy and U 2 above the canopy).
In the presence of waves, Schaefer and Nepf (in review) proposed an extension to Equation 1 to account for the waveinduced current in the meadow: (4) The maximum wave-induced current U max is (Luhar, 2021) in which U w is the wave velocity amplitude, k is the wave number, and σ is the wave angular frequency. The addition in Equation 4 assumes that the waves do not modify the momentum transfer between the meadow and overflow layers (i.e., C is effectively unchanged), and that the generation of wave-induced current is not altered by presence of an imposed current. In addition to predicting the reduction of velocity, it is also important to predict the intensity of turbulence generated in the meadow. For unidirectional flow, a submerged meadow contributes two sources of turbulence production: shear production P s , associated with the shear at the top of the meadow, and stem production P w , associated with the wakes of individual stems and leaves. The shear production is in which z is the elevation above the bed and u is the timeaveraged velocity. For a stem diameter d and stem density m s (number of stems per unit bed area), the wake production is (e.g., Nepf, 2012).
Averaged over the meadow height, Equation 7 is approximated as Angle brackets indicate a canopy vertical spatial average. A prediction of TKE, k t , is possible by equating the rates of turbulence production to the rate of turbulence dissipation, ε ∼ k t 3/2 / l t , with l t the turbulence integral length scale (e.g., Tennekes and Lumley, 1972). Tanino and Nepf (2008) considered this model for an array of rigid emergent cylinders, for which the shear production can be neglected. Zhang et al. (2018) extended the Tanino and Nepf (2008) model to a flexible seagrass meadow with blades of width w b and blade density m b (number of blades per unit bed area). For shoot spacing S > 1.8w b , l t ∼ w b , such that P w ∼ ε, from which The scale factor δ = 1.1 for unidirectional flow in an emergent rigid meadow, as shown in Tanino and Nepf (2008). Equation 9 with δ = 1.1 has also been validated for pure wave conditions in rigid canopies, but with the velocity scale replaced with the root-mean-squared wave velocity, U w,RMS (Tang et al., 2019) or the wave velocity amplitude U w (Tinoco and Coco, 2018;Chen et al., 2020). However, Zhang et al. (2018) found that the scale coefficient δ was reduced for a submerged flexible artificial seagrass meadow, and this was attributed to the motion of the blades reducing the relative velocity between the waves and the blade.

MATERIALS AND METHODS
Experiments were performed in a 24 m long, 38 cm wide, and 60 cm deep laboratory flume (Figure 2). A piston-type paddle wavemaker controlled by a Syscomp WGM-101 waveform signal generator (refer to Appendix A in Luhar, 2012 for details) produced waves of period T = 2 s and varying amplitude. To reduce wave reflection, a 1:5 sloped aluminum ramp covered with rubberized coconut fiber was placed at the downstream end of the flume. The upstream edge of the ramp was lifted to allow for the passage of current. The ramp reduced wave reflection to below 11%, based on the analysis described in Goda and Suzuki (1977). Currents of varying speeds were recirculated through an inlet pipe (8 cm inner diameter) 0.8 m downstream of the wavemaker and an outlet drain downstream of the ramp.

Construction of Artificial Seagrass Shoot and Meadow
A seagrass shoot model was constructed to be geometrically and dynamically similar to Zostera marina (e.g., Ghisalberti and Nepf, 2002). Each model shoot had six blades laser-cut from lowdensity polyethylene (LDPE). Each blade was 13 cm long, 0.3 cm wide, and 0.1 mm thick. The blades were wrapped around the upper half of a 0.6 cm diameter, 1.3 cm long cylindrical dowel using tape, which increased the diameter by 0.1 cm. The erect height of each shoot was h = 13.6 cm.
The shoots were inserted in a staggered pattern into predrilled holes in polyvinyl chloride boards such that the rigid dowels, representing the seagrass sheaths, extended 0.6 cm above the baseboards. The shoot density m s = 950 shoots per square meter, and the blade density m b = 5,700 blades per square meter. A bare baseboard (1.2 m in length) was placed at the upstream and downstream ends of the 6.1 m long meadow. Experiments were performed with mean water depths D = 27 cm and 45 cm above the baseboards. For each water depth, four pure wave cases, four pure current cases, and sixteen combined wavecurrent cases were considered (see Supplementary Section 2). Turbulence generated at the current inlet slightly modified the wave amplitude relative to pure wave conditions with the same wavemaker setting.

Velocity Measurements, Wave Gauge Measurements, and Meadow Imaging
The velocities in the streamwise (x), lateral (y), and vertical (z) directions were defined as (u, v, w), respectively. Instantaneous velocities were measured with a Nortek Vectrino three-dimensional acoustic velocimeter, which was centered longitudinally between successive shoot rows and laterally between adjacent shoots, at the spanwise center of the flume. This position within a staggered array has been shown to offer accurate estimates of the laterally averaged velocity and turbulence in current and wave conditions (see Figure 2 in Chen et al., 2013 andFigure 2 in Zhang et al., 2018). To minimize interference with the measurement volume, blades were removed from shoots within a 10 cm diameter circle around the probe (see Luhar et al., 2010). Measurements were taken along a vertical transect in 1 cm intervals above the bed at the longitudinal center of the meadow (x = 2.8 m), for 240 s at 200 Hz for each case.
Velocity records were cleaned using the acceleration thresholding methods described in Goring and Nikora (2002). Each component of velocity was decomposed into the summation of the time-averaged (denoted by an overbar), phase-averaged (denoted by a tilde), and turbulent (denoted by prime) velocity. For example, for the streamwise component, in which t is time. The number of samples n s in each wave period was estimated through autocorrelation. The phase-averaged velocityŨ (θ) was calculated for each of the n s = 403 phase bins, in which θ is the wave phase. The time-averaged velocity was calculated as Frontiers in Marine Science | www.frontiersin.org The depth-averaged imposed current U c was defined by integrating U over depth upstream of the meadow. The timeaveraged TKE was calculated as (12) Water surface displacement η was measured with a wave gauge (1,000 Hz for 90 s) at the longitudinal center of the meadow, which provided both the wave amplitude and confirmed stationary wave conditions. The wave-to-wave variation in wave amplitude was less than 1% of the phase-averaged wave amplitude, determined using 1 2π 2π 0 η RMS (θ) dθ/η w,RMS , as described in Zhang et al. (2018). This indicated that the variation in wave form made negligible impact on the estimation of the turbulence within the phase-averaged method (Equation 12).
The in-canopy velocity U 1 for each case was estimated as a vertical average of the time-mean velocity between the bed and the time-mean deflected canopy height, h d . Mean deflected heights were estimated from digital videos of six blades painted black at the longitudinal center of the meadow collected using a Nikon D7500 camera. The videos were viewed and processed using MATLAB VideoReader, readFrame, and Sobel edge detection functions. A red-greenblue pixel identification algorithm was used to identify the painted blades. The vertical positions of the highest and lowest part of each blade were recorded for all cases under the wave crests and troughs. The mean deflected canopy height h d was defined as the average across all six blades. The characteristic wave velocity amplitude U w for each case was defined from Stokes second-order wave theory under the crest (LeMéhauté, 1976;Dean and Dalrymple, 1984) at the undeflected canopy height (z = h) using the measured wave amplitude at x = 2.8 m, the location of mid-meadow velocity measurements.

Drag Coefficient
The drag coefficient C D for the flat rectangular blades was assumed to be 1.95 (as in Zhang et al., 2018). Zhang et al. (2018) found that using a variable C D to account for the impact of different hydrodynamic conditions on blade reconfiguration did not improve TKE predictions, as the drag coefficient used in TKE predictions represents the drag along the vertical part of the plant. Furthermore, although Sarpkaya and Storm (1985) observed that the addition of a current to oscillatory flow could alter the drag coefficient, the impact depended on the Keulegan-Carpenter number KC = U w T/w b (Keulegan and Carpenter, 1958). They observed that the drag coefficients of cylinders in combined wave-current conditions converged for KC > 15, and the impact of current on drag coefficients became negligible for KC > 30. In this study all KC were greater than 24. Based on this, and for simplicity, in this study the drag coefficient in pure current, pure wave, and wave-current conditions was assumed to be 1.95.

Time-Averaged Velocity Profiles
Selected profiles of time-averaged velocity were used to illustrate the flow behavior (Figure 3), including the strongest imposed current in both depths. The full set of profiles is given in Supplementary Section 3. See Supplementary Section 2 for details of all conditions. For pure current cases (red triangles in Figure 3), the time-averaged velocity profiles included a mixinglayer at or just above the meadow interface and a peak velocity well above the meadow (Figures 3A,D). Wave-current cases (circles in Figure 3) with U w /U c < 2.5 retained this mixinglayer time-mean structure, but the center of the shear layer moved toward the bed as either the imposed current or wave velocity was increased, both of which were associated with a reduction in the canopy deflected height (horizontal dashed lines in Figure 3; see Supplementary Section 4).
For most cases with U w /U c > 2.5, the time-mean velocity peaked at the top of the canopy, reflecting the contribution from the wave-induced current (discussed in Section "Theoretical Background for Predicting Velocity and Turbulence in Meadow" and predicted with Equation 5), so that these conditions were considered to be wave-dominated. The transition at approximately U w /U c = 2.5 is supported by Figure 4, which illustrates the distribution of cases for which the time-averaged velocity peaked at the top of the canopy (wave-dominated, triangles) or for which the time-averaged velocity peaked above the canopy (circles). The transition was not a function of current Reynolds number (U c D/ν, in which ν is the kinematic viscosity of water). Finally, for the largest waves (black circles in Figure 3), the time-averaged current was reduced near the water surface. This was attributed to a wave-induced Eulerian drift, which has been observed in several previous studies in laboratory water channels and in the field (Nepf et al., 1995;Gjøsund, 2003;Smith, 2006;Monismith et al., 2007).

Reynolds Stress and Turbulent Kinetic Energy Profiles
For pure current (red triangles in Figures 3B,E) and wavecurrent (circles in Figures 3B,E) cases with a clear mixing layer flow structure, the magnitude of the Reynolds stress was approximately zero near the bed, increased to a peak near the mean deflected meadow height (dashed horizontal lines in Figure 3), and then decreased with distance above the meadow. As the wave amplitude increased, the peak Reynolds stress magnitude decreased (Figures 3B,E). A peak in TKE appeared near the top of the canopy (Figures 3C,F), coincident with the peak in Reynolds stress, which was consistent with shear production (Equation 6). As the wave amplitude increased, this TKE peak was diminished, consistent with the decrease in Reynolds stress peak magnitude. The opposite trend was observed within the canopy. Specifically, as wave amplitude increased, TKE within the canopy increased (see the progression of light gray to darker gray to black symbols below the dashed lines in Figures 3C,F), reflecting the contribution of wake production, which increased with the addition of waves. For In panels (D-F), U c = 6.8 cm s −1 and U w = 2.9-3.9 cm s −1 (lightest gray), U w = 8.1-8.8 cm s −1 (second-lightest gray), U w = 10.1-11.1 cm s −1 (second-darkest gray), and U w = 15.4-19.0 cm s −1 (black) (turbulence at the current inlet pipe slightly modified wave amplitudes for higher currents relative to pure wave conditions for the same programmed wave). Horizontal dashed lines denote measured h d , with matching colors (dashed-dotted lines correspond to pure current). Profiles could not be extended over the entire depth due to limitations of the instrument and the presence of waves. Profiles for cases with depth D = 45 cm are shortened to 30 cm to show detail within and above the canopy.
FIGURE 4 | Wave-to-current velocity ratio (U w /U c ) vs. current Reynolds number (U c D/ν). Circles denote cases for which the time-averaged velocity peaked above the canopy, similar to pure current conditions. Triangles denote cases for which the time-averaged velocity peaked at the top of the canopy, indicating that the wave-induced mean current was significant within the meadow.
larger waves, TKE in the upper water column increased toward the water surface (e.g., black circles in Figure 3F). The nearsurface enhancement in TKE was attributed to a secondary circulation. Specifically, the interaction between progressive waves and vertical vorticity, here associated with the sidewall boundary layers, generated a secondary circulation by Langmuir instability, with downwelling at the channel center (Nepf and Monismith, 1991). This circulation elevates nearsurface turbulence (Nepf et al., 1995).

Prediction of Magnitude of Peak Reynolds Stress
The addition of waves had little impact on the magnitude of peak Reynolds stress for U w /U c < 1, but for U w /U c > 1 the peak Reynolds stress magnitude decreased with increasing wave amplitude, compared to the corresponding pure current conditions (Figure 5A; also see the Reynolds stress profiles in Supplementary Figures 3.1, 3.2, noting the change in magnitude of the maximum Reynolds stress with increasing wave amplitude for each current). The reduction in Reynolds stress was predominantly associated with a decrease in velocity shear, rather than a change in the efficiency in turbulent momentum transport. This distinction was illustrated by considering the idealized two-layer model for canopy flow, which defines the Reynolds stress at the top of the meadow, τ h (max) , in terms of the layer-averaged time-averaged velocity within, U 1 , and above, U 2 , the meadow (Equation 3). When the momentum coefficient C (=0.014-0.020) was estimated using Equation 2, then Equation 3 predicted the maximum magnitude of Reynolds stress within uncertainty ( Figure 5B) for both pure current cases and the wave-current cases that exhibited shear layer behavior, indicating that the momentum coefficient was not significantly changed by the waves.
However, note that the best fit line through the combined wave-current cases (black dashed line fitted to circles in Figure  5B, u w max predicted = (0.80 ± 0.09) u w max measured ) fell below that for the pure current cases (red dashed line fitted to triangles in Figure 5B, u w max predicted = (1.11 ± 0.07) u w max measured ), which suggested that Equation 2 underpredicted C by 20-30% for combined wave-current cases. This trend suggested that the presence of waves augmented the vertical turbulent transfer of momentum. However, this conclusion is only speculative, because in fact Equation 3 produced predictions consistent with the measured stress. The trends in Figure 5 can be explained as follows. The addition of waves reduced h d , relative to pure current conditions, which increased the depth of overflow (D − h d ), which in turn decreased U 2 . The reduction in U 2 − U 1 led to a decrease in Reynolds stress. Because the efficiency of turbulent momentum change was not significantly altered (C was the same within uncertainty), the decrease in Reynolds stress was consistent with the decrease in velocity gradient ( Figure 5B). It is important to note that if C is not impacted by waves (as suggested here), then two-layer models developed for pure current conditions can be used to predict the in-canopy velocity in combined wave-current conditions, using the deflected height (i.e., Equation 4), which is supported by analysis in Schaefer and Nepf (in review). This is a useful result for modeling flows through submerged meadows.

Canopy Turbulent Kinetic Energy Measurements and Predictions
The canopy-averaged shear production P s of turbulence was estimated by averaging Equation 6 over the deflected canopy height. The canopy-averaged wake production P w of turbulence was estimated from Equation 7 using the in-canopy depthaveraged, time-averaged velocity, U 1 . This represents a lower bound on wake production, because u(z) 3 > U 1 3 , and contributions from wave velocity were neglected. Even using an underestimate for P w , the ratio P w / P s was greater than five for all but three cases and greater than 10 for all but eight cases (see Supplementary Section 4). Given this, it was reasonable to neglect the shear-production and thus to expect a form of FIGURE 6 | Measured canopy-averaged TKE ( TKE ) vs. the squared sum of the wave velocity and the in-canopy depth-averaged, time-averaged velocity ((U w + U 1 ) 2 ). Red and green dashed lines denote best fit lines for pure current and pure wave cases, respectively. For clarity among the cases with smaller TKE, the highest TKE cases were excluded (see Figure 7). Error bars indicate uncertainties assessed in the canopy-averaged TKE (propagating uncertainties using half of the range from each measurement when the records were split in half).
Equation 8 to predict TKE within the canopy. However, we must determine the appropriate velocity scale for TKE prediction.
For combined wave-current conditions, Chen et al. (2020) proposed that stem-generated turbulence scaled with the sum of the wave velocity and imposed depth-averaged current velocity, i.e., with the maximum velocity in the wave period, U max . Within a rigid submerged meadow, they observed a linear relationship between TKE and U 2 max . Here, we considered a similar relationship in a flexible canopy, using the sum of the wave velocity U w and in-canopy depth-averaged, time-averaged velocity U 1 . Recall that in the presence of waves, a wave-induced current will augment U 1 above the imposed current velocity (see Equation 4). Specifically, we considered the dependence of TKE and (U w + U 1 ) 2 , as shown in Figure 6.
First, note that there was a stronger dependence between TKE and velocity for pure current (red triangles in Figure 6) than for pure waves (green diamonds in Figure 6). Specifically, using the velocity scale (U w + U 1 ) in Equation 9, the scale factor for pure current, δ pc = 1.8 ± 0.3, was six times larger than for pure wave (δ pw = 0.30 ± 0.04). This makes physical sense, because flexible blades can move with the wave orbital velocity, which reduces the relative velocity, the drag, and the wake turbulence production Zhang et al. (2018). In pure current, the flexibility can allow plants to become more streamlined, but it cannot reduce the relative velocity. Zhang et al. (2018) also observed a smaller scale factor in Equation 9 for flexible blades under pure wave conditions. Specifically, using velocity scale U w,RMS in Equation 9 they found δ pw = 0.44, which converts to δ pw = 0.22 for the velocity scale U w , as used here. Note that the pure current scale factor δ pc was larger than the scale factor 1.1 found for rigid emergent cylinders (Tanino and Nepf, 2008), which may be explained by the difference in morphology between emergent cylinders and submerged model plants. First,Equation 9 assumes that the integral length scale is equal to the blade width, l t = w b , which may not be appropriate in the lower canopy, where the blades are bundled into a sheath of larger dimension (sheath diameter = 0.7 cm), or near the top of the canopy where larger shear layer vortices are present (Poggi et al., 2004;Zhang et al., 2020). Second, Equation 9 uses the canopy-averaged velocity U 1 .
FIGURE 7 | Measured canopy-averaged TKE ( TKE ) vs. the proposed hybrid model (δ 2 pc k t,pc + δ 2 pw k t,pw , Equation 13) with (A) axes scaled to highlight that the majority of cases collapse reasonably well by the hybrid model and (B) showing all cases. The solid black lines denote one-to-one lines. In panel (B), asterisks denote cases involving the highest waves, but atypically low TKE, as discussed in Section "Blade Reconfiguration and Wake Production" in the text. Error bars indicate uncertainties assessed in canopy-averaged TKE (as in Figure 7). However, for a submerged meadow the velocity varies over the canopy height, and U 1 2 will underestimate u(z) 2 , and this must be offset by a larger scale coefficient.
Next, consider the combined wave-current cases (gray to black circles in Figure 6), which predominantly fell in between the dependences observed for pure current (red triangles) and pure wave (green diamonds). This suggested that a simple hybrid wake production model might collapse all of the cases. Specifically, we proposed the following modification to Equation 9 for wavecurrent conditions k t,wc : Indeed, this hybrid model collapsed most of the data (Figure 7).
As an alternative to the hybrid model, it would be useful to evaluate when the canopy turbulence is dominated by either the current or the waves, allowing for prediction of k t using only the wave or current velocity. Consider the ratio of canopy-averaged turbulence in wave-current conditions TKE wc normalized by corresponding TKE in pure current TKE pc and pure wave FIGURE 9 | Conceptual sketch of the side view of a model seagrass plant showing blade motion under a wave crest (black curves) and wave trough (gray curves), corresponding to the strongest wave forcing (e.g., black circles in Figure 7). The rectangle represents the sheath. TKE pw vs. U w /U c , as shown in Figure 8. For U w /U c < 1, increasing U w /U c had minimal impact on canopy turbulence relative to pure current conditions, compared to the increase in TKE observed with increasing U w /U c for U w /U c > 1 ( Figure 8A). Alternatively, when current was added to waves, TKE wc / TKE pw was greater than 1 for U w /U c < 1, but converged to 1 for U w /U c > 1 (Figure 8B). Together, these trends suggest that for U w /U c > 1, a good prediction of TKE can be made using just the wave velocity. Conversely, for U w /U c < 1, a good prediction of TKE can be made using just the current velocity.

DISCUSSION
Laboratory experiments with an artificial seagrass meadow described the impact of waves on the time-averaged velocity, Reynolds stress, and TKE within a flexible submerged meadow. The study considered a range of wave and current conditions for a single meadow density (950 shoots m −2 , 5,700 blades m −2 ). One major result of this study is that for a flexible meadow, wave velocity is less efficient in generating TKE (smaller scale coefficient), because wave-induced motion of the blades reduces the relative velocity between the blade and the waves. The resultant proposed hybrid model (Equation 13) is valid when wake production dominates shear production. Overall, this study provided important validation of models that predict Reynolds stress and turbulence within a submerged meadow in realistic conditions of combined waves and currents. These models provide a quantitative framework to predict ecosystem services facilitated by meadow-mediated hydrodynamic conditions.

Impact of Waves on Mean Deflected Height
The observed reduction in h d in the presence of waves contrasted with observations in Paul and Gillis (2015), who observed that for combined wave and current conditions, the presence of a wave did not appear to impact the mean deflected canopy height. The discrepancy might be related to the fact that Paul and Gillis (2015) used the tip of the blade to estimate the canopy height, while the present study considered the highest position along the length of the blade. The pronation of blades in the direction of wave propagation, resulting in a decrease in mean deflected height, has been observed in previous studies and attributed to two mechanisms. First, observations showed that the vertical component of the wave orbital velocity induced asymmetric blade motion and deflection in the direction of wave propagation (Döbken, 2015). Second, the wave-induced current generated in the meadow can pronate the individual blades (Zhang et al., 2018).

Reynolds Stress at the Top of the Canopy
The direct contribution of waves to shear and turbulent momentum exchange (i.e., Reynolds stress) at the top of the canopy is limited to long period waves, as described in Ghisalberti and Schlosser (2013). They defined a wave Reynolds number, Re w = 2U w d/πν, and a Keulegan-Carpenter number KC s as the ratio of wave period to the time scale of vortex formation. Specifically, KC s = U w T/L D , with L D = 2 (C D a v ) −1 (1 − φ) (e.g., Chen et al., 2013). Pure waves can generate vortices at the top of a canopy only when both KC s > 5 and Reynolds number Re w > 1, 000. For the meadow in the present study, KC s < 2.2 and Re w < 600, indicating that the waves did not contribute to Reynolds stress or shear production at the top of the canopy, consistent with the fact that C was not impacted by the addition of waves.
In addition, note that the peak Reynolds stress did not always occur at the mean deflected height, but for most cases occurred within a few cm above the mean deflected height. This was in contrast to measurements in unidirectional flow made for a flexible meadow of lower density (230 shoot m −2 ), for which the peak Reynolds stress was consistently at the canopy interface (Ghisalberti and Nepf, 2006). The shoot density in the present study was much higher at 950 shoot m −2 , which may have reduced the shear penetration into the canopy, pushing the shear layer and peak Reynolds stress slightly above the canopy interface. Considering the measured maximum deflected height instead of the measured mean deflected height, as in Ghisalberti and Nepf (2006), did not explain the slightly shifted position of the Reynolds stress peak.

Blade Reconfiguration and Wake Production
For the three wave-current cases marked with asterisks in Figure 7, the canopy-averaged turbulence was noticeably lower than other wave-current cases of the same programmed wave amplitude and depth. Cases of the highest wave amplitudes had the smallest deflected heights (see Supplementary Section 4), with many blades touching the bed during the wave period (see the conceptual sketch in Figure 9). It was possible that the blades, which were in-line with the sheath, restricted flow around the sheath, which reduced or eliminated the shedding of vortices from the sheath for the marked cases. This is a potentially important observation, as it illustrates how wake turbulence can be eliminated when the meadow pronation due to reconfiguration becomes extreme. While this description does not fully explain the differences in canopy TKE among all of the strongest wave cases, it was likely a contributing factor.

Relative Strengths of Wave and Current Velocities
For clarity, we note that the aforementioned threshold of approximately U w /U c = 2.5 described changes in the timemean velocity profile. Specifically, the time-mean velocity peaked above the canopy when U w /U c < 2.5, and peaked at the top of the canopy when U w /U c > 2.5. Meanwhile, shifts in canopy turbulence were observed at a lower threshold U w /U c = 1. Combining these, when 1 < U w /U c < 2.5, the wave velocity dominated the production of TKE in the canopy, but did not significantly modify the current-induced time-mean flow structure behavior. We note that these transitions may have some dependence on canopy density, which was not explored in the present study. Finally, seagrass meadows exhibit a wide range of density (140-30,000 blades per bed area; see Table 3 in Luhar et al., 2010), both more and less dense than the model canopy considered in this study (5,700 blades m −2 ). The assumption that wake production dominates shear production, which was validated in the present study, likely extends to denser meadows, but may not be valid in sparser meadows.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HN and RS conceptualized the project and contributed to writing the draft. RS conducted the experiments and data processing. HN provided the laboratory space and equipment. Both authors contributed to the article and approved the submitted version.