Modeling complex flow structures and drag around a submerged plant of varied posture

Although vegetation is present in many rivers, the bulk of past work concerned with modeling the influence of vegetation on flow has considered vegetation to be morphologically simple and has generally neglected the complexity of natural plants. Here we report on a combined flume and numerical model experiment which incorporates time‐averaged plant posture, collected through terrestrial laser scanning, into a computational fluid dynamics model to predict flow around a submerged riparian plant. For three depth‐limited flow conditions (Reynolds number = 65,000–110,000), plant dynamics were recorded through high‐definition video imagery, and the numerical model was validated against flow velocities collected with an acoustic Doppler velocimeter. The plant morphology shows an 18% reduction in plant height and a 14% increase in plant length, compressing and reducing the volumetric canopy morphology as the Reynolds number increases. Plant shear layer turbulence is dominated by Kelvin‐Helmholtz type vortices generated through shear instability, the frequency of which is estimated to be between 0.20 and 0.30 Hz, increasing with Reynolds number. These results demonstrate the significant effect that the complex morphology of natural plants has on in‐stream drag, and allow a physically determined, species‐dependent drag coefficient to be calculated. Given the importance of vegetation in river corridor management, the approach developed here demonstrates the necessity to account for plant motion when calculating vegetative resistance.


Introduction
Vegetation is abundant in rivers and has a significant influence on the hydraulic functioning of these systems. Over multiple spatial and temporal scales, vegetation can influence mean and turbulent flow fields [Nepf, 2012a], flow conveyance [J€ arvel€ a, 2002], sediment dynamics [Sand-Jensen, 2003], and river form and morphodynamics [Gurnell, 2014]. As such, vegetation is a key component in controlling the hydrodynamics of aquatic ecosystems [Nikora, 2010].
The approach applied here assumes that vegetation is a dynamically moving porous blockage . Its porosity comes from the fact that within a controlled volume, flow is able to pass through individual plant elements (stems, leaves, etc.). Previously, vegetation has been represented as a porous medium with a permeability [Papke and Battiato, 2013;Battiato and Rubol, 2014;Rubol et al., 2016]. However, the volume is dynamic in that it interacts with flow forcing. As the flow strength increases, the plant will reconfigure and as such the plant volume will decrease, thus pushing the stems closer together and decreasing the porosity. This paper reports on a combined flume and numerical model study of depth-limited flow around an isolated, submerged riparian plant at three different flow Reynolds (Re) numbers. The influence of Re number was investigated by increasing the mean flow velocity whilst keeping flow depth constant. Initially, the unstressed (no flow forcing) three-dimensional morphology of the plant was captured using terrestrial laser scanning (TLS) following the methodology of Boothroyd et al. [2016a]. In the flume, flow velocities were quantified through acoustic Doppler velocimeter (aDv) measured velocity profiles, while plant motion was recorded and quantified through high-definition video imagery. This provided the necessary boundary conditions and validation data to incorporate the plant (unstressed and stressed) into a computational fluid dynamics (CFD) model by application of a mass flux scaling algorithm [Hardy et al., 2005].
This study enabled an examination of the feedbacks between flow and vegetation dynamics allowing us to: (i) gain an understanding of time-dynamic and time-averaged plant motion characteristics; (ii) quantify the three-dimensional mean and turbulent flow through and around an isolated, submerged riparian plant; (iii) evaluate how well the numerical model can predict measured flow; and (iv) quantify the influence of timeaveraged plant posture on vegetative resistance, through the calculation of drag forces and drag coefficients, the Vogel exponent, and Manning's n.

Plant Characteristics
The plant used in these experiments was the waxy leaved, evergreen shrub Hebe odora (Figure 1a SUBMERGED PLANT COMPLEX FLOW STRUCTURES Ecologically, Cockayne [1958] reports H. odora to be distributed across many river beds or near to streams flowing through tussock grass land or fell field sites. The species is particularly widespread throughout New Zealand, especially on flushed ground and stream banks [Wardle, 1991]. Practically, the shrub had a measured height of 0.22 m and a diameter of 0.20 m, therefore enabling complete submergence in the flume. The open framework, and internal structure of stems/leaves that were not especially dense, allowed the laser to penetrate into the plant interior and fully quantify the plant morphology in the unstressed state. Furthermore, the ability not to deteriorate in a laboratory environment for the duration of the measurement period (48 h) and the stable root ball which could be firmly attached to the base of the flume were further reasons for selection of this particular species.

Terrestrial Laser Scanning (TLS) and Voxelization for Capturing Plant Morphology
A three-dimensional representation of the plant morphology under unstressed conditions (no flow forcing) was acquired using terrestrial laser scanning (TLS), following the approach previously reported in Boothroyd et al. [2016a]. A RIEGL VZ-4000 scanner was used to scan the plant. The scanner has a beam divergence of 0.15 mrad, a field of view 608 3 3608 and an effective measurement rate of up to 222 000 measurements per second. Scans were collected at a distance of 5 m, with p and h increments set to 0.0158, controlling the horizontal and vertical alignment, respectively. At this scanning distance, the mean distance between neighboring points in the registered point cloud was 0.0024 m. The plant was scanned from four opposing perspectives to minimize occlusion [Moorthy et al., 2008].
The registered point cloud was postprocessed and voxelized following the workflow of Jalonen et al. [2015].
The neighborhood based approach of Rusu et al. [2008] to classify outliers and remove erroneous point returns was applied, where n 5 100 and r 5 1 (see Boothroyd et al. [2016a] for a full discussion). To simplify and reduce the number of postprocessed points, while maintaining the plant morphology, a voxelization procedure was applied (Figure 1b). A regular voxel size of 0.005 m adequately resolves plant morphologies in comparable plants [Boothroyd et al., 2016b]. Furthermore, the voxel size is more than double the scan resolution, so ensured an adequate representation of the internal structure of the plant (Figure 1c). From the voxelized plant representation, the areal and volumetric distributions are calculated along the x, y, and z axes.

Flume Experiments
Flume experiments were undertaken in an Armfield S6 MKII glass-sided tilting recirculatory flume with a 0.3 m wide cross section and a length of 12 m. The working section was the central 2 m, unaffected by inlet or outlet flows. The slope of the flume was adjusted to achieve constant local flow depth (0.29 m) in the central working section (with uniform flow conditions). Experimental conditions are shown in Table 1. We acknowledge that working at these width to depth ratios, there is the potential for wall-induced secondary circulation although secondary flow of this nature is typically less than a few percent of the average flow velocity [Colombini, 1993]. Furthermore, we acknowledge blockage ratio effects, as the plant width is of the same order as the flume width, and this has previously been shown to influence drag acting on submerged macrophytes [Cooper et al., 2007].  [Sontek, 1999]. Inflow measurements were made for the entire 2.5 h duration of the experimental runs.
Downstream of the plant, a Nortek Vectrino-II Profiler aDv was mounted on a moveable carriage and used to collect downstream velocity profiles, each composed of six 0.035 m high subprofiles centered at 0.5 Y/w. A streamwise spacing of 0.125 m was selected for the first five profiles, and 0.25 m thereafter so that the wake zone was fully measured, giving six profiles in total. The measurement volume of the Vectrino-II is displaced 0.04 m below the transmit transducer, meaning that sampling was limited to within 0.06 m of the water surface; therefore, velocity profiles capture 0.7 Z/h. The Vectrino-II has a reported accuracy of 0.001 m s 21 6 0.5% of the measured velocity [Nortek, 2012]. Cumulative variance associated with different time-averaging windows was used to test for stationarity in the velocity signal, with a 120 s record length sufficient for systematic convergence [Sukhodolov and Rhoads, 2001].
Postprocessing of the aDv data followed the recommendation of Thomas and McLelland [2015]. Initially, two-dimensional phase unwrapping or dealiasing was undertaken applying the two-step noncontinuous quality-guided path (TSNCQUAL) algorithm [Parkhurst et al., 2011]. Next, the phase-space filter of Wahl [2003] was applied to reduce the number of spikes resulting from Doppler noise or phase difference ambiguities [McLelland and Nicholas, 2000]. Finally, a noise analysis threshold of Zedel and Hay [2010] was applied, with the threshold set as a confidence interval of 99.9%. The time-averaged velocities presented herein were computed using only the data retained after dealiasing and filtering.
A weakness of the Vectrino-II aDv is that the signal-to-noise ratio (SNR) is not constant over the 0.035 m profiling distance, and in some cases this can lead to discontinuous mean velocity profiles when subprofiles are placed on top of one another and stacked [MacVicar et al., 2014]. Consequently, Brand et al. [2016] show mean velocities from the upper 0.025 m of subprofiles to be the most reliable. Sensitivity to the distance away from the transducer is heightened when analyzing turbulence statistics [MacVicar et al., 2014;Brand et al., 2016]. In this study, subprofiles were overlapped by 10%, and the extremities of overlapping subprofiles averaged to reduce the artifacts introduced at the periphery of the measurement footprint. This minimized any discontinuous effects on the mean velocity, thus maintaining the spatial coverage necessary to resolve fine-scale flow structures and enable comparisons with high-resolution CFD predictions.

Measuring Plant Motion Dynamics
To monitor the time-dynamic and time-averaged plant motion characteristics, a high-definition video camera with a 1440 3 1080 pixel resolution was fixed perpendicular to the plant on the outside of the flume. The video camera was focused on the plant, and a photogrammetric target attached to the outside of the flume for spatial scale. Video was recorded at 25 Hz, which in postprocessing far exceeded the resolution required to capture the plant motion (see section 3.3).
Individual image frames from the video were extracted and corrected for distortion [Wackrow et al., 2015].
To relate the size of pixels in the image space to measured distances on the photogrammetric target, the image-scale factor (ISF) was quantified [Wackrow et al., 2015]. This was completed by measuring the pixel distance on the photogrammetric target in the undistorted image space, and comparing this to the real distance in the object space. As such, pixels in the undistorted image frames are related to measured distances, enabling characteristics such as plant height and length to be determined with 60.005 m measurement error.
To distinguish between plant versus water pixels, frames were converted into a binary plant image, by applying the Otsu [1979] image classification algorithm ( Figure 2). Isolating the blue band of the RGB image space, a standard image processing technique was used to reduce the gray-scale image into a binary image, by calculating the global threshold using the method of Otsu [1979]. Following Hardy et al. [2011a], the binary images were postprocessed to remove any isolated, spur, or H-connected pixels, with disconnected areas containing less than 50 pixels assumed to be flume debris and removed. Similarly, holes in the binary image were filled, with the final binary image for each frame detailing the outer extent of the plant (0 5 unoccupied, 1 5 occupied, Figure 2). To investigate the time-dynamic plant motion characteristics at the finest spatial scale, tracks of the centroids from isolated plant tips are monitored. Plant tips are selected at three distinct locations across the plant (1 5 uppermost front, 2 5 uppermost back, 3 5 middle, Figure 2b), and displayed as time-dynamic tracks that represent two-dimensional motion in the downstream and vertical ( Figure 3).
The time-averaged plant posture is quantified by calculating the mean position of the plant from the binary images over the cycle of plant motion. This is displayed as a relative probability, where a probability value of 1 indicates that the pixel was constantly occupied, and a probability of 0 indicates no occupancy ( Figure 4a). To highlight zones of time-averaged plant motion, a transition frequency matrix was  constructed using the sequence of binary values in each pixel. By concentrating on transitions from 0 to 1, or from 1 to 0, it is possible to detect the zones which were intermittently occupied, and therefore represent motion. This is displayed as fraction of the binary time series where transitions occur (Figure 4b), with a greater fraction indicating comparably more plant motion.

Unstressed and Stressed Plant Representations
The unstressed plant representation is discretized directly from the voxelized point cloud, and herein referred to as the unstressed model. For the stressed plant representation, the time-averaged plant postures are used as boundary conditions to inform the discretization of the plant for the stressed model. For this, the unstressed voxelized representation is equally sliced at 0.02 m intervals along the vertical axis, and these slices are incrementally translated to match the time-averaged plant posture in the model coordinate system under increasing Re. No translation occurs in the spanwise direction, with posture only shifted in the vertical and streamwise directions. The approach ensures that the number of voxels remains constant between unstressed and stressed plant representations, but the porosity and void spaces are reduced as the plant is vertically compressed, altering the volumetric canopy morphology. We acknowledge that plant reconfiguration is a three-dimensional process and that the spanwise thickness of the plant is likely to change under flow, further altering the volumetric canopy morphology, and reducing the porosity.

Incorporation of Unstressed and Stressed Plant Representations Into a CFD Model
The CFD model solves the full three-dimensional Navier-Stokes equations discretized with a finite-volume method. Simulations were run by applying a two-equation j-e turbulence model, modified using Renormalization Group Theory (RNG) [Yakhot and Orszag, 1986]. This model has been shown to outperform the standard j-e turbulence model in regions of high strain, flow separation, and reattachment [Yakhot and Orszag, 1986;Lien and Leschziner, 1994;Hodskinson and Ferguson, 1998;Bradbrook et al., 2000] and has been widely adopted in geomorphological CFD applications [Marjoribanks et al., 2016]. It is applied here due to the large degree of fluid strain that is associated with the shearing of flow around the plant. Furthermore, no modification has been applied to the turbulence model for either bed and domain sidewalls, which were treated as a no-slip boundary, and the nonequilibrium wall function is applied which assumes local equilibrium of turbulence (y 1 5 47) [Launder and Spalding, 1974].
The differencing scheme used was hybrid-upwind, where upwind differences are used in high-convection areas (Peclet number > 2) and central differences are used where diffusion dominates (Peclet number < 2). The pressure and momentum equations are coupled by applying the SIMPLEST algorithm [Patankar and Spalding, 1972], where the velocity field was solved using the momentum equation followed by a pressure correction to solve the mass equation, ensuring a divergence-free velocity field. Weak linear relaxation was used for the pressure correction, while false time step relaxation was used for the other variables. The convergence criterion was set such that the residuals of mass and momentum flux were reduced to 0.1% of the inlet flux.
The model domain was designed to replicate the working section of the flume; 400 cells long, 60 cells wide, and 59 cells high (1,416,000 grid cells) created with a spatial resolution of 0.005 m. Inlet velocity was set to match the bulk inlet boundary conditions measured using the inlet aDv (0.22, 0.30, and 0.37 m s 21 , Table 1), with a turbulent intensity of 5%. The outlet was defined with a fixed-pressure boundary condition where mass is allowed to enter and leave the domain.
The unstressed and stressed plant representations are read into the numerical model using a mass flux scaling algorithm (MFSA) [Lane et al., , 2004Hardy et al., 2005] representing the plant as a numerical porosity [Marjoribanks et al., 2014c;Boothroyd et al., 2016a;Marjoribanks et al., 2016]. The approach we apply here has previously been reported by Marjoribanks et al. [2016] where the MFSA accounts for the blockage created by the plant. No additional drag force term is needed because the plant is explicitly represented through a grid-scale blockage, so wall conditions are automatically set at the edge of each blocked cell. As such, the finite volume continuity equation takes the form ; (1) Water Resources Research where f is the variable of interest (u,v), p represents the value at the cell center, the index k represents the values at neighbouring cell centers and the previous time step, and S is the linear source coefficient. The neighbor links (a k ) have the form where A k is the cell-face area, / is the cell-face porosity, q is the fluid density, u is the local velocity perpendicular to the face, D is the diffusion term, and T is the transient term. To introduce the MFSA, the value of / is calculated for each face according to the presence of vegetation. This information comes straight from the voxelization procedure and follows a binary occupied/unoccupied porosity, as the 0.005 m voxel size is equal to that of the 0.005 m grid cell size. The numerical grid was defined with vertices that are exactly collocated with the voxelized representation, meaning that the voxelized plant maps directly onto the grid cells (having equal density), and therefore a binary blocked/unblocked porosity treatment follows [Lane et al., 2004]. This means that all permeability through the plant is explicitly represented through the gridscale blockage. Were the numerical porosity nonbinary, the drag force implemented as a momentum sink term would be required to account for permeability. The plant, represented as the voxelized blockage, is incorporated 0.35 m downstream from the inlet (0.175 X/l), and centered in the domain (0.5 Y/w). The computational domain is shown in Figure 5.

Evaluation of Numerical Predictions Against Measured Data
The general agreement of velocity profiles between measured and modeled data are first quantified. Velocity profiles were selected to cover the entire downstream range of wake separation and reattachment. Each of the measured velocity profiles consists of 195 points at a 0.001 m spatial resolution, therefore covering 0.7 Z/h due to sampling constraints near the surface. For the modeled velocity profiles, the 59 points at a 0.005 m spatial resolution were linearly interpolated to the 0.001 m spatial resolution of the measured data, meaning that point densities between the data sets were equal. To quantify the general level of agreement, it is necessary to consider both the degree of scatter and the degree of bias in the data. This was analyzed using the correlation coefficient (r) and the slope (b) obtained from reduced major axis (RMA) regression. RMA was selected to reflect the expected uncertainty on both the dependent and independent variables , as both velocity data has uncertainty, and this is especially relevant given the known data quality issues surrounding the Nortek Vectrino-II Profiler aDv. Following this, we quantify differences between measured and modeled profiles through the root-mean-square error (RMSE) and mean absolute error (MAE).
To assess in detail the agreement between measured and modeled data, a shape-based similarity statistic is used to calculate the visual distance (d V ) between velocity profiles. The d V statistic was first proposed by Marron and Tsybakov [1995] in the context of qualitative smoothing, with the distance measure analogous to that detected by the naked eye. The closer the d V value to zero, the greater the similarity between velocity profiles. The d V statistic therefore accounts for the minimum Euclidean distances between points on curves/profiles in both the horizontal and vertical directions, following the removal of scaling (see Minas et al. [2011] for overview). Conventional distance based methods, such as the area between curves, only compute the distance in one direction, and therefore do not truly capture the notion of shape [Minas et al., 2011]. The plant tips therefore begin to follow a more steeply inclined trajectory with increasing Re. Divergence from a 1:1 ratio between downstream and vertical movement extents, and the more pronounced variation about the motion tracks at higher Re, imply a transition in the oscillation regime as the plant reacts to the increased flow. These different movement extents and characteristics of motion result from differences in the internal plant structure, specifically the branch structure, which leads to different exposures to the flow and implies different natural frequencies of movement. Although these time-dynamic plant motions are only highlighted for three isolated plant tips, differences in motions are demonstrated across the extent of the plant body, meaning that over the whole plant, the potential range of movement is likely to be large.
Changes in the time-averaged plant posture with increasing Re are shown in Figure 4, detailing the mean position of the plant (Figure 4a), and the zones where greatest plant motion is detected (Figure 4b). The plant is clearly deflected and vertically compressed when stressed, with the plant positioned lower in the flow depth and slightly shifted further downstream. Compared to the unstressed state, the mean plant height is lowered by 6.4, 11.2, and 17.7% with increasing Re. In contrast, the mean plant length increased by 4.8, 9.8, and 14.4%. This demonstrates that as the time-averaged plant posture is shifted under the constant flow, the plant is vertically compressed and reduces the volume. In doing so, this causes a reduction of the volumetric canopy morphology, and a reduction in the plant porosity, causing more flow to pass around the blockage. The postural changes influence the lead angle at the plant front and lee angle at the plant back, both of which are measured from the upright ( Table 1). The lead angle is consistently greater than the lee angle (>3.58) throughout the Re range, and as Re increases from 65,000 to 110,000, both lead and lee angles are shown to more than double, again indicating a Re dependence.
Time-averaged plant motion is highlighted through the fraction of time occupied by transitions in the binary time series image sequence ( Figure 4b). Regions of greatest motion are shown by thicker zones around the plant extent, and tend to be positioned on the upper and leeward sides of the plant body. As Re increases from 65,000 to 89,000, the area of the transition zones increases by 12.7%; with a further 1.4% increase in area as Re increases from 89,000 to 110,000. Due to biomechanical constraints, it is expected that there will be an upper limit as to how much the plant can reconfigure, and because of this, the amount of change will be initially large but then decrease with increasing Re.
Traditionally, the spectrum of vegetation canopy motion in response to increasing flow speeds can be categorized into four distinct regimes; erect, gently swaying, strong coherent swaying (monami, or vortexinduced vibrations for the case of a single plant), and prone [Nepf and Vivoni, 2000]. Here, however, the time-dynamic and time-averaged plant motions vary across the single plant body. The time-averaged plant motion indicates that the plant has shifted its general posture, associated with static reconfiguration and streamlining to the mean flow [Sand-Jensen, 2003;Siniscalchi and Nikora, 2013]. Dynamic reconfiguration is related to smaller-scale oscillations, associated with the instantaneous flow and correlated with drag fluctuations and upstream turbulence [Siniscalchi and Nikora, 2013]. For a single plant, a combination of the four motion regimes can therefore occur simultaneously at a constant flow speed. This is because of the range of biomechanical properties in the plant, and the influence of the plant on the local flow conditions [Hurd, 2000]. There will be different stem widths, lengths, thicknesses, and flexural rigidities throughout the plant, and therefore the drag force will affect each component differently, resulting in a range of responses to the flow. This is not a random process with respect to time and space. Some parts of the plant will respond to the flow first, while other parts of the plant will take longer to readjust and reconfigure. The ability to reconfigure to the changing flow stress will therefore vary over the plant body.

Velocity Profiles
To investigate the general agreement between the measured and modeled time-averaged velocity profiles, results from the RMA regression are shown (stressed model in Figure 6; correlation coefficient (r) and slope (b) of unstressed and stressed models in Table 2 Lane et al., 2004;Hardy et al., 2005Hardy et al., , 2011bSandbach et al., 2012], and therefore demonstrate a very good general agreement between the measured and modeled data. To further quantify the differences between the measured and modeled data, the root-mean-square error (RMSE) and mean absolute error (MAE) are calculated, and the mean average of these errors over the six profiles displayed (Table 3). Although both the unstressed and stressed models represent the measured values well, with maximum RMSE of 0.072, 0.020, and 0.017 m s 21 , and MAE of 0.059, 0.016, and 0.015 m s 21 for the u, v, and w velocities; errors are consistently lower for the u velocity in the stressed CFD model. For v and w velocity, errors remain very small (<0.02 m s 21 ) for both the unstressed and stressed CFD models. Relating this back to the general agreement of velocity profiles from the RMA regression ( Figure 6 and Table 2), these very small errors are associated with low-magnitude velocity values (<0.05 m s 21 ), implying that minor discrepancies between measured and modeled data result in relatively large differences in the correlation coefficients.
To understand the specifics of the u velocity profiles, Figure 7 shows profiles of measurements and both unstressed and stressed models in the downstream region 0.25-0.625 X/l. An obvious agreement in the overall profile shape between the measured and modeled data is visible. In the measured and both of the modeled profiles, a zone of reduced velocity in the plant wake exists in the range 0.25-0.375 X/l. For the unstressed model, velocities in this wake zone appear more uniform and gradually transition from low to high, whereas for the stressed model more velocity fluctuations are shown (Figure 7c at 0.25 X/l). Velocity minima in the low-velocity zone are better represented by the stressed model, with a maximum difference in minima of 60.037 m s 21 against the measured profiles over the entire Re range, compared to 60.049 m s 21 with the unstressed model. Beneath 0.2 Z/h, a zone of flow acceleration associated with subcanopy flow in the near-bed region is measured and modeled, although the model consistently fails to capture the magnitude in the velocity peak. With increasing distance downstream, velocity profiles show signs of recovery, reverting toward a fully developed profile.
To quantify the shape-based similarity between these profiles, the visual distance statistic is calculated and displayed (Figure 7). The closer the d V statistic to zero, the greater the similarity in shape between the measured and modeled profiles [Marron and Tsybakov, 1995;Minas et al., 2011]. Proximal to the plant (0.25-0.3125 X/l) where the velocity profiles are most heterogeneous, when the visual distance statistics are  5 1.06), to the six measured profiles. Velocity profile shape is therefore better modeled closer to the plant. Overall, given the closer general agreement and similarity in profile shapes of the stressed CFD model for the u velocity, and the smaller quantified errors of the v and w velocity, we argue that the stressed model is more representative and more able to predict threedimensional mean flow than the unstressed model.
Next, we examine the spatial patterns of the u velocity difference between the measured and stressed modeled data (aDv-CFD), and the cumulative difference over the six profiles, to investigate areas of model under Further downstream as flow recovers, the magnitude of this velocity difference becomes smaller. Over the vertical extent, measured u velocities tend to be greater than modeled u velocities in the range 0-0.3 Z/h, suggestive of model under-prediction in the subcanopy region. This effect is most prominent closest to the plant, before flow recovers. This model under-prediction in the subcanopy region is present for the entire Re range. Above this zone, model overprediction is noted in the plant wake region. These differences between the measured and stressed modeled velocities, and specifically their vertical dependence, explain the systematic variation shown in Figure 6. The structural difference is caused by the transition from zones of model under-prediction in the subcanopy region, to zones of model overprediction in the wake region. This is attributed to an under-representation of the plant blockage in the model, meaning that the plant representations are more porous than the plant in the flume. This could be caused by not representing the spanwise deformation of the plant during reconfiguration, which would likely further reduce the porosity of the plant, and is a limitation of the current approach applied here. This would help account for the limited ability to capture the magnitude of velocity peaks in the subcanopy region, and overprediction of velocity in the wake region, whereby penetration of flow through the interior of the plant is responsible for modifying the velocity immediately behind the plant blockage (Figures 7 and 8). If the plant porosity was further reduced, the magnitude of subcanopy flow would increase, with reduced velocity expected in the wake region.

Streamwise Velocity Field and Wake Structure
To further investigate the streamwise velocity field and wake structure, measurements of u velocity have been linearly interpolated to allow for the comparison of the region 0.25-0.625 X/l against the same region in the stressed model. Time-averaged u velocity was normalized by the inlet velocity (0.22, 0.30, and 0.37 m s 21 , Figure 9), with the difference then calculated (aDv-CFD). Vertical slices of the u velocity flow field along the midline of the model domain were compared (0.5 Y/w). A number of similarities between the u velocity fields are clear, most notably in the length and shape characteristics for each of the plant wakes. The calculated differences (Figure 9) mirror the zones of model under/overprediction identified previously, associated with an under-representation of the plant blockage in the model.
Flow separation and reattachment results in the formation of a low-velocity wake zone behind the plant blockage; here defined as <0.5 of the streamwise velocity normalized by inlet velocity (Figure 10). A number of similarities are evident between the measured and stressed modeled plant wakes. For the stressed modeled data, the wake shape remains approximately constant between Re 5 65,000-89,000, forming a dual-layered structure, with the larger upper wake inclined slightly upward. At Re 5 110,000, a singlelayered, thicker wake exists, and this is positioned lower in the flow depth. In each case the wake markedly thins in the downstream direction. The modeled wake thickness has decreased by 41.5, 59.2, and 41.2% between 0.25 and 0.4 X/l for each Re. Also, the modeled wake length decreases from 3.9 plant lengths at Re 5 65,000, to 3.3 plant lengths at Re 5 110,000.
In the measured case, a single-layered wake maintains an almost consistent shape throughout the entire Re range, although ends more abruptly than in the modeled case. The measured wake thickness has decreased by 38.7, 44.1, and 40.2% between 0.25 and 0.4 X/l for each Re, and therefore demonstrates a similar thinning pattern to the modeled wake. The measured wake length is similar to the modeled wake length, showing a decrease with increasing Re, from 3.9 plant lengths at Re 5 65,000, to 3.1 plant lengths at Re 5 110,000. The trailing edge of the wake, which indicates the downstream wake limit, is progressively shifted toward lower values of Z/h with increasing Re, and is expected to be associated with the vertical compression of the plant. Combining Figures 9 and 10, a rapid gradation in velocity is modeled to occur between the low-velocity wake zone and the free-stream zone, with this velocity discontinuity indicative of shear layer formation and the presence of Kelvin-Helmholtz instabilities, the frequency of which we calculated to be 0.20, 0.23, and 0.30 Hz for the different Re [Ho and Huerre, 1984;Ghisalberti and Nepf, 2002].

Turbulent Flow Structures
We have shown that the modeling system is capable of predicting the measured time-averaged flow conditions, and therefore extend our analysis to investigate the turbulent structures associated with the plant. The horizontal distribution of uv vorticity in a plane at 0.5 Z/h is shown in Figure 11, with twodimensional uv flow lines overlain, to show the modeled flow pathways through and around the plant. In each Re case, the shedding of clockwise and anticlockwise rotating regions of high vorticity occurs, with the development of a large region of high vorticity at the outer edge of the plant boundary, whose position correspond with deflections to flow lines. It is acknowledged that vorticity alone is unable to distinguish vortices from the strain field, which is expected to be high in the region due to the velocity gradient present. Here it is suggested that this plant shear layer turbulence is dominated by Kelvin-Helmholtz and G€ ortler-type vortices generated through shear instability [Ghisalberti and Nepf, 2002]. Within the plant body, smaller regions of high vorticity are distributed in the central region of the plant, and the flow lines distorted. This illustrates the internal flow dynamics and the forcing of flow through conduits in the plant morphology. As Re increases, and the volumetric canopy morphology is vertically compressed and its porosity reduced, both the magnitude of the uv vorticity and the length of the highvorticity regions at the outer plant boundary increase. To exemplify this, by taking an arbitrarily defined value of uv vorticities as 6>1.5 Hz, the maximum length of the high-vorticity region increased by 41.2% over the entire Re range. Changes in the volumetric canopy morphology therefore influence the velocity gradient and regions of vorticity present.
An impression of the three-dimensional turbulent structures forming around the plant blockage is visualized by plotting isosurfaces of the Q criterion [Hunt et al., 1988] (Figure 12). For a physical interpretation of the Q criterion, a vortex is assumed to be present where the vorticity tensor exceeds the rate of the strain tensor, and a local pressure minima exists, forming the vortex envelope. With increasing Re, it is shown that the vortices extend further downstream, with the maximum vortex length almost doubling over the modeled Re range. Vorticity magnitude increases with Re, with higher magnitudes distributed toward the upper region of the plant. Between Re 5 65,000-89,000, the vortex appears to be stretched but retains a similar form, although vorticity magnitude increases. For Re 5 110,000, the vortex shape is modified, with a lengthening of the uppermost vortex tail, and this is suggested to be associated with the stronger shear layer instability formed between the low-velocity wake zones and the free-stream zone above. Again, this indicates an interplay between the volumetric canopy morphology and turbulent flow structure.

Implications for Vegetative Flow Resistance
These high-resolution process predictions of flow and vegetation dynamics are useful when considering vegetative flow resistance. Using the modeled pressure fields and integrating the upstream/downstream difference in pressure acting normal to the plant surface over the entire lateral extent of the blockage, drag forces (F d ) and drag coefficients (C d ) are calculated following the method described by Marjoribanks et al.  Table 4 we compare F d and C d between the unstressed and stressed models, to analyze how differences in plant posture influence drag.
For both plant representations, F d increases almost linearly with Re. However, the modeled F d is consistently lower when using the stressed model, by as much as 27.1%. Incorporation of the stressed plant representation, which captures the shift in time-averaged plant posture to the mean flow, results in a reduced drag force relative to the unstressed plant representation where the plant remains fully upright to flow. Similarly, the C d is lower in the stressed model, and this decreases with increasing Re. The more streamlined shape of the plant helps account for this reduction and emphasizes the dynamic nature of C d under changing plant posture. These results signify the sensitivity of F d and C d to changes in plant posture, and the importance of explicitly accounting for this when investigating vegetative resistance.
The Vogel exponent, w, quantifies the drag reduction by reconfiguration of a flexible body through a power law dependence with velocity [Vogel, 1984] F d / U 21w : Using this F d -U relationship, w has previously been found to typically range from 20.2 to 21.2 for natural vegetation [de Langre et al., 2012], with w relating to the flexibility of the plant . For a rigid plant, w 5 0, equation (3) returns the classical squared relation, while for a flexible plant, w 5 21, equation (3) returns a linear force-velocity relation [Cullen, 2005;Whittaker et al., 2013]. For the H. odora specimen investigated here, w 5 20.37, therefore suggesting the plant falls at the flexible end of the stiffer range. It is important to remember that the F d -U relationship is representative of only the relatively lowvelocity range investigated (<0.37 m s 21 ), and therefore could be modified at higher velocities. However, the marked plant deflection and considerable motion detected (Figures 3 and 4) suggest the plant has moved beyond the trans-flexing zone where deflection is negligible and into the flexing zone where the plant streamlines and frontal area reduces with velocity [Xavier et al., 2010]. Although providing an empirical relationship, the Vogel exponent is not dimensionally correct, and cannot be used to calculate the drag force and subsequent energy loss within vegetated rivers [Marjoribanks et al., 2014a]. However, it does allow for the broader quantification of flexibility in plants, as demonstrated here . Furthermore, using the Vogel exponent approach, a number of authors have characterized bulk vegetative resistance terms including parameterization of plant biomechanical and plant geometry components, with inclusion of separate foliage and stem contributions [V€ astil€ a and J€ arvel€ a, 2014], and species-specific drag coefficients [J€ arvel€ a, 2004; Aberle and J€ arvel€ a, 2013] to improve process representation in such equations [Marjoribanks et al., 2014a].
Here we link the Re-dependent drag coefficient to Manning's n. In submerged cases, the presence of a turbulent mixing layer between the vegetated low-velocity and free-stream zones adds complexity to derivations of vegetative resistance [Shucksmith et al., 2011], and therefore an equation applicable to submerged vegetation is required. For submerged grasses, Wilson and Horritt [2002] consider the momentum absorbing plant area for the calculation of Manning's n, following: where g is gravitational acceleration, R is the hydraulic radius, A is the momentum absorbing area or frontal area calculated from the TLS, and a is the cross-sectional flow area. Manning's n values of 0.086, 0.080,

Discussion
Here we have investigated time-dependent and time-averaged plant dynamics for a submerged, riparian plant under three different flow Re through a combined experimental and numerical modeling approach. Plant dynamics investigated through experimental measurements, demonstrated that for the timeaveraged plant motion there is a shift in the general plant posture as the plant reconfigures to the mean flow [Sand-Jensen, 2003;Siniscalchi and Nikora, 2013]. In doing so, the volumetric canopy morphology is vertically compressed and therefore the volume available for flow to pass through the plant (i.e., porosity of the plant) is reduced. When time-dynamic motion is analyzed, we show the oscillatory motion to be Re and location-dependent. This would increase the overall drag of the plant and increase the turbulence on the lee side of the plant. However, there is no detectable coherence of movement across the plant. It may be expected that the plant would demonstrate coherence of plant motions when coupled to strong oscillations in flow velocity as shown for large-scale canopy flow processes [Ghisalberti and Nepf, 2006;Okamoto and Nezu, 2009;Okamoto et al., 2016], and where instantaneous motion is closely related to the passage of large-scale eddies that interact with the entire plant [Siniscalchi and Nikora, 2012]. However, a natural plant is not a single homogenous unit and there are a range of stem widths, lengths, thicknesses, and flexural rigidities throughout, thereby affecting plant motion [Hurd, 2000]. This motion will be further affected by the internal plant structure influencing the exposure to flow, and the flow passing through the plant, which is discussed below. This highlights the difficulty in predicting motion dynamics in riparian plants, and is further exemplified by Weissteiner et al. [2015] who show that structural properties of riparian vegetation, such as the vertical distribution of stem area, can control the extent of plant compression during reconfiguration.
The plant morphology captured for the unstressed state through TLS, and discretized from high-resolution video imagery for the stressed state, enabled the time-averaged plant posture to be used as boundary conditions in a CFD model. For the riparian species investigated here, TLS has captured the static plant posture in air, and when submerged in the flume this structure remained similar due to the flexural rigidity supporting the plant. For macrophytes where the flexural rigidity is lower, however, once submerged in the flume arguably the plant posture would be different due to buoyancy effects influencing the plant structure. As such, we justify the use of a static model of plant morphology to model the stressed posture by altering the volumetric canopy morphology to account for the reduction in porosity. The modeling system was shown to correctly predict flow measurements by comparison with measured data, although the model slightly under-predicted flow in the subcanopy region and very slightly overpredicted flow in the plant wake region. Explicit consideration of the spanwise deformation to the plant under three-dimensional reconfiguration, and the effect this has on plant porosity, would likely improve predictions. In the model we replicate a region of subcanopy flow acceleration beneath the bulk of the plant blockage, which interacts and inclines upward the low-velocity wake, resulting in a complex and three-dimensional wake structure ( Figure  10). Similar observations have been made in airflow around a fir tree windbreak and have revealed largescale recirculation caused by an upwelling vortex immediately behind the tree, sustained by this subcanopy acceleration [Lee and Lee, 2012], while in river applications, Freeman et al. [2000] reported flow diversion and acceleration beneath the plant canopy. B€ olscher et al. [2005] quantified this phenomenon in the field with velocity profiles. In the plant wake region, the modeled wake shape remains approximately constant for the low and mid Re, forming a dual-layered structure, with the larger upper wake inclined slightly upward and markedly thinning in the downstream direction. For the highest Re, a thicker single-layered wake is positioned lower in the flow depth. The measured and modeled wake lengths show very good agreement, decreasing in length with increasing Re.

10.1002/2016WR020186
In the plant wake, regions of high vorticity develop at the outer edge of the plant boundary. Within the plant body, smaller regions of high vorticity in the central region of the plant illustrate the internal flow dynamics and the forcing of flow through conduits in the plant morphology. As Re increases, and the volumetric canopy morphology is reduced, both the magnitude of the uv vorticity and the length of the highvorticity region from the outer plant boundary increase. This indicates an association between vortex length and the strength of the shear instability and demonstrates the interplay between volumetric canopy morphology, changes in plant porosity, and turbulent flow structures.
Plant shear layer turbulence is dominated by Kelvin-Helmholtz and G€ ortler-type vortices generated through shear instability [Ghisalberti and Nepf, 2002]. However, the plant is not a classic bluff body [Schnauder et al., 2007], instead tending toward porous media flows [Yagci et al., 2010], with penetration of fluid through the canopy resembling ''bleed-flow'' [Raine and Stevenson, 1977]. This aligns with our observation of flow through conduits in the internal structure of the plant. As the flow strength increases and the plant reconfigures to a dynamic equilibrium position with flow, the volumetric canopy morphology decreases and the porosity of the plant decreases. Similarities are drawn with experiments around rigid arrays of cylinders to resemble a porous cylinder, where Zong and Nepf [2012] report a weakening of the vortex street as the solid volume fraction decreases (and so porosity increases). Here we show that as porosity is reduced, more flow is forced around the plant, and with increasing Re vortices migrate further downstream as the shear instability grows stronger (Figures 11 and 12).
These changes in flow dynamics have direct implications for the pressure regime, thus modifying the drag response. We quantify a 13-27% drag reduction when comparing the stressed plant posture to the unstressed plant posture. This highlights the significance of plant postural changes in reducing flow separation through streamlining and reconfiguration [Vollsinger et al., 2005], leading to a subsequent reduction in form drag [Nikora, 2010] which accounts for 90-97% of the total drag for turbulent flows [Lilly, 1967;Vogel, 1994;Stoesser et al., 2010]. We acknowledge the limitation that only form drag is directly modeled in the current approach, although for riparian plant species, form drag has been shown to dominate over viscous drag [Nikora, 2010; V€ astil€ a and J€ arvel€ a, 2014]. This process understanding can be extended to consider drag coefficients which are generally poorly understood for the complex geometries that characterize natural vegetation [Marjoribanks et al., 2014a] and vary as a function of both vegetation density and stem Reynolds number [Tanino and Nepf, 2008]. The wide range of plant morphologies associated with natural vegetation, and the effects of their foliage contribution, will further act to modify the drag coefficient [Wilson et al., 2003;Boothroyd et al., 2016a]. However, practical applications dealing with vegetation often assume a predefined, constant drag coefficient of 1.0 or 1.2 . In the example here, we show that the drag coefficient varies from 1.45 to 1.34 for the H. odora specimen. This has direct implications for bulk vegetative resistance, and demonstrates the necessity to account for plant motion and morphology when defining drag coefficients, instead of relying on predefined parameterizations.
In order to demonstrate this uncertainty, calculated Manning's n values between 0.086 and 0.078 were higher than those of established texts [e.g., Chow, 1959], although we acknowledge the differences in spatial scale to which these values apply e.g. reach scale versus flume. Such improved parameterizations could be applied to conveyance estimators either applying a one-dimensional hydrodynamic model or using a Conveyance Estimation System [Wallingford, 2004]. These are typically based upon the Manning's equation which parameterizes and/or calibrates all frictional resistance and is applied to produce the correct relationship between flow and water level. The same is true for industry-standard two-dimensional hydraulic models, which are shown to be particularly sensitive to the floodplain roughness values selected [Straatsma and Baptist, 2008]. For example, Abu-Aly et al.
[2014] report a maximum Manning's n value of 0.343 for woody vegetation patches (which is far higher than established texts) and find that application of spatially distributed vegetation roughness values causes a 7.4% increase in mean depth, and 17.5% decrease in mean velocity relative to an unvegetated roughness scenario. Changes in the Manning's n values associated with the vegetation component of resistance will therefore alter hydraulic model outputs. In practice, the roughness term in these models are used as a calibration parameter [e.g., Mason et al., 2003], and may no longer reflect the intended process representation. The Manning's n value therefore becomes effective as it represents several processes that contribute to energy loss (e.g., momentum loss, dispersion associated with secondary circulation, and diffusion). Applications of these models therefore tend to be run for a matrix of effective Manning's n values over a parameter range to ensure an optimum value is identified, meaning the Water Resources Research

10.1002/2016WR020186
values used may differ significantly from their measured or estimated values [Lane, 2005]. However, by deriving Manning's n for a single species following the approach used in this paper, which considers the momentum loss, dispersion associated with secondary circulation, and diffusion, the Manning's n value has greater physical basis.
An explicit consideration of plant postural changes is shown to be essential in the prediction of threedimensional flow and wake structures. Previous studies have shown how plant characteristics such as the plant morphology and foliation state can influence velocity distributions and drag forces [Gross, 1987;Freeman et al., 2000;Wilson et al., 2008;Endalew et al., 2011]. Here we show how changes in the time-averaged plant posture, influencing the volumetric canopy morphology and porosity, can result in similar effects. These postural changes are relevant to a number of applications concerning vegetation-flow interactions. This is highly appropriate to patch-scale studies of riparian vegetation roughness at the margins of active channels [e.g., Manners et al., 2013] where a simplified beam elasticity equation was used to quantify deflections of stem height under flow. Whether such an equation accounts for the plant dynamics remains unclear, especially as we highlight plant motions to vary spatially over the plant body across different scales, and have a Re dependence. We are therefore currently developing the static plant representation through the approach of Marjoribanks et al. [2015Marjoribanks et al. [ , 2016 where an Euler beam equation is used to dynamically reconfigure the plant or patch through hydrodynamic forcing.

Conclusions
Results presented here successfully demonstrate that by incorporating the time-averaged plant posture into a high-resolution CFD model, we are able to predict three-dimensional mean and turbulent flow around a submerged riparian plant. Key findings include: 1. Under hydrodynamic loading, plant motion can be separated into time-dynamic and time-averaged components. Time-dynamic plant motions were investigated by tracking the motion of plant tips, showing a transition from horizontally dominated to vertically dominated movement and increased movement extents as Re increases. Plant tip motion varies across the plant body in reaction to the local flow and shows a dependence with Re. We suggest this is because stem widths, lengths, thicknesses, and flexural rigidities vary throughout the plant and therefore the drag force will affect each component differently, meaning there will be a range of responses to the flow. Some parts of the plant will respond to the flow first, while other parts of the plant will take longer to readjust and reconfigure. This results from differences in the exposure of the plant to the flow, caused by the internal plant structure. Time-averaged plant motion is associated with shifts in the general plant posture which reconfigures to the mean flow, resulting in an 18% reduction in plant height, a 14% increase in plant length, and a doubling of the lead and lee angles of the plant body. These motions are responsible for vertically compressing the plant in the flow, thereby reducing the volumetric canopy morphology and plant porosity. 2. Vertical velocity profiles illustrate a zone of subcanopy flow acceleration which interacts with the lowvelocity wake region behind the plant, and helps generate an upwardly inclined wake structure that thins in the downstream direction. With the wake defined as <0.5 of u velocity normalized by inlet velocity, with increasing Re the wake length decreases from 3.9 plant lengths to 3.1 plant lengths. We suggest that the plant shear layer turbulence is dominated by Kelvin-Helmholtz and G€ ortler-type vortices generated through shear instability [Ghisalberti and Nepf, 2002], the frequency of which we estimate to be 0.20, 0.23, and 0.30 Hz at Re 5 65,000, 89,000, and 110,000. 3. The time-averaged plant posture has been used as boundary conditions to discretize the plant in a highresolution CFD model, with the modeling system shown to correctly predict flow measurements. General agreement is quantified through RMA regression, with u velocity correlation coefficients >0.8. Furthermore, the visual distance statistic demonstrates the similarity between aDv-measured and CFD-modeled velocity profile shapes. Following validation, analysis is extended to investigate the turbulent flow structures associated with the plant, where it is shown that the modeled vortices migrate further downstream as the shear instability grows stronger. 4. Using this process understanding of flow and vegetation dynamics to consider vegetative resistance, we demonstrate the drag force experienced by the stressed plant representation is up to 27% lower than that of the unstressed plant representation, highlighting the important role of reconfiguration as a form Water Resources Research 10.1002/2016WR020186 drag reduction mechanism [Nikora, 2010]. We quantify the F d -U relationship, and calculate a Vogel exponent of 20.37, indicative of a plant at the flexible end of the stiffer range. Similarly, we show the dynamic nature of the drag coefficient under increasing Re, and how this impacts upon Manning's n values, falling from 0.086 to 0.078. These values remain considerably higher than traditional bulk vegetative resistance terms for comparable vegetation, selected from classical look-up tables applied to channels/floodplains.