Large Eddy Simulation of Unidirectional and Wave Flows through Vegetation

The experimental conditions were chosen similar to the laboratory experiments of Liu et al. (2008), where the bulk velocity ($U_b$) in the channel was specified as $0.33\mathrm{m}/\mathrm{s}$. An initial ramp up time of 18 flow-through periods ($t_f=l/U_b$) were allowed before data collection started. The choice of this period was based on the fact that the mean and turbulence velocity profiles were devoid of any initial effects after 18 flow-through cycles. The mean velocity results became independent of the averaging period after 45 flow-through cycles beyond the initial ramp time, whereas the mean turbulence quantities became invariant of the averaging period after 50 flow-through periods beyond the initial ramp time. The independence of the results on the averaging period was quantified by comparing the difference in the RMS deviation (RMSD) with respect to the experimental data for the vertical profiles at the six chosen locations for each incremental averaging period with a step size of 5. When the difference in the RMSD between two profiles obtained from two successive averaging periods became less than 1%, the latest averaging period was chosen as being sufficient for data analysis purposes. Thus, the data presented in this paper were obtained by averaging over 55 flow-through periods after the initial 18 flow-through periods for flow ramping. Temporal data were collected at 333 Hz, which was sufficient to trap the tail of the Kolmogorov $-5/3$ spectrum.

In Fig. 3 the first and third row of panels from the top show the vertical profile of the variation of mean horizontal ($U$) and vertical ($W$) velocities, nondimensionalized by the bulk velocity ($U_b$) at each of the six locations shown previously in Fig. 1(b). The second and fourth rows show the vertical profiles of the turbulence component of the velocities in the two directions nondimensionalized by the shear stress defined as $U_{*}=\sqrt{ghS}$, where $h=0.114\mathrm{m}$ was the height of the domain as well as the depth of water; and $S=0.003$ was the slope of the channel in the original laboratory experiment. The turbulent velocities were calculated as RMS of the observed velocitieswhere $U_i$ and $W_i$ = resolved instantaneous horizontal and vertical velocity values.

Generally, it is observed that the mean horizontal velocities at all six locations match well with the experimental data and also compare favorably with the numerical results of Stoesser et al. (2009), considering the fact that the latter used a mesh with about twice the resolution as the present simulation. The streamwise velocities within the canopy layer are considerably lower than those in the fluid layer above the canopy and follow the general trend observed in flow over submerged vegetated canopies. The fluid on the top of the canopy accelerates and produces the classical hyperbolic tangent profile with the inflection point close to the top of the canopy. All six positions show the same global hyperbolic type profile, which suggests the spanwise homogeneity of this phenomenon, with the locations just behind ($V1$) and before ($V2$) the stem showing large gradients of horizontal velocity near the top of the canopy. This indicates the effect of the recirculation zone, which can be defined as the sheltered region within the wake of a cylinder, in the dominant flow direction, within which the centerline mean streamwise velocities are negative. It is worthwhile to note that the mean horizontal velocity within the recirculation zone is actually better predicted in this model than in Stoesser et al. (2009) who have noted that the accuracy of prediction of this region is strongly dependent on the choice of the subgrid-scale (SGS) model, LES data averaging period, or measurement period in the original experiment. In initial sensitivity runs it was found that the accurate prediction of the mean and turbulent quantities in the recirculation zone are much more sensitive to the near-wall mesh than the averaging period or the choice of the SGS model, assuming of course the averaging period is large enough so statistically steady-state turbulence profiles have been reached within a 1–2% error. Hence, the model consistently underpredicts the gradients of mean and turbulence quantities in the velocity profiles within the recirculation zone and even more so at vertical $V1$. An important goal of this research is to use LES simulations, with mesh sizes and temporal resolutions, which can successfully simulate experiments conducted in typical wave flumes found in engineering laboratories, within a reasonable level of accuracy. These experiments may either be in field scales or prototype scales of larger coastal wetland canopies mimicking the height, density, and diameter of the rigid portion of the stems. As a secondary condition, the total size of the problem in these simulations should be such that they can be run in conventional HPC clusters with relatively moderate computational effort. Hence, the near-wall resolution, which yields acceptable results for mean quantities in the majority of the locations, is considered satisfactory.

The streamwise turbulence velocity profiles (second row from the top in Fig. 3) agree very well with the experimental data, except at location $V1$, because of the previously explained reason. The peak of the intensities occurs at or near the stem tops with the peak lying a little below the tops at the free locations ($V2$ and $V5$), suggesting that streamwise vortices advected from preceding stems are pushed down when they reach the sideways free locations. It is also noted that the LES model slightly underpredicts the peaks at those two locations as was observed by Stoesser et al. (2009). The present model predictions are found to be better than those by Stoesser et al. (2009) in trapping the vertical gradients of the turbulent horizontal velocities at the free locations ($V2$ and $V5$). Also, in the above canopy region the turbulence streamwise velocities are better predicted consistently, particularly near the full-slip wall boundary at the top. The difference is caused by the full-slip boundary condition set at the top versus the symmetry condition as imposed by Stoesser et al. (2009), who instead explained the overprediction of turbulence intensities near the rigid lid to coarse resolution and to large, roughness scale turbulence elements originating at the actual free surface, which is treated as a rigid lid for the experiment presented in this paper. Simulations run with symmetry condition in this case produced similar overshoots caused by reflection of turbulent vortices back into the domain, and it is recommended that a full-slip top wall is a much better boundary condition than symmetry for setups of this type.

The vertical profiles for the vertical velocity (third row from top in Fig. 3) agree reasonably well with the experiment except within the recirculation regions ($V1$ and $V6$) in which the model does not predict the peaks very well. The vertical velocities also seem to be a little overpredicted both within and above the canopy. The LES model captures the peaks of the profile very well for the sheltered zones just outside the recirculation region ($V3$ and $V4$), but it seems to miss those in the free zones ($V2$ and $V5$) and is attributed to lower resolution at these locations.

Velocity profiles for the turbulent component of the vertical velocity (fourth row from top Fig. 3) are found to agree well consistently for all the locations, except for the recirculation zones ($V1$ and $V6$). The underprediction of the turbulent quantities within the recirculation zones is caused by a grid resolution lower than that used by Stoesser et al. (2009). The larger cell sizes within this zone are unable to trap the eddies with the peak turbulent energy, resulting in an under-prediction of the turbulent velocities. The peak of the profiles occur just below the canopy top for the free zones ($V2$ and $V5$) and at the canopy top for the sheltered zones ($V1$, $V3$, $V4$, and $V6$) and compares well with the results of Stoesser et al. (2009).

In addition to sampling velocity profiles, the depth-integrated forces on each cylinder were computed as a sum of pressure ($F_p$) and viscous ($F_v$) forces as $F_{\text{tot}}=F_p+F_v$where ${\tau }_{\mathrm{eff}}=(\tau +{\tau }_{\mathrm{sgs}})dU/dn{|}_{\mathrm{cyl}\text{-}\mathrm{wall}}$ is the effective shear stress on the cylinder wall with the gradient calculated with the wall parallel velocity. Here $r$ and $\varphi$ are the radius of the cylinder and the angle made by the infinitesimally small element located along the cylinder surface on which the forces act, with the axis of the cylinder in the $X$, $Y$-plane. The force data were sampled at 50 Hz and averaged over 55 flow-through periods to yield a representative steady-state force. The time-averaged and RMS forces for each of the cylinders is shown in Table 1. The mean forces showed little deviation from cylinder to cylinder with the percentage coefficient of variation ($\sigma /\mu \times 100$), where $\sigma$ and $\mu$ are, respectively, the standard deviation and mean of the forces on all cylinders, being about 0.6%, whereas the same for the RMS values was about 10%. The canopy-averaged mean force, $\overline{F_{\text{tot}}}=1/N{\sum }_{i=1}^N(\overline{F_{pi}}+\overline{F_{vi}})$, was used along with the mean velocity of the cells immediately at the top of the canopy ($\overline{U_{h_v}}$) to compute the canopy representative drag coefficient (Tanino et al. 2005) given by $C_d=\overline{F_{\text{tot}}}/0.5\rho {\overline{U_{h_v}}}^{2}D{h}_v$, where $h_v=12D$. The calculated values for $\overline{F_{\text{tot}}}$, $\overline{U_{h_v}}$, and $C_d$ from the model results were $16.77\times {10}^{-3}\mathrm{N}$, $43.2\mathrm{cm}/\mathrm{s}$, and 0.371, respectively. This value of $C_d$ has been compared with the empirically derived coefficient from Ghisalberti and Nepf (2004). The vertical variation of the ratio $\zeta =C_d(z)/C_{dA}$ is given in Ghisalberti and Nepf (2004) aswhere $C_{dA}=(1+10{\mathsf{R}}_{hv,D}^{-2/3})/1.16[1.16-9.31(ad)+38.6{(ad)}^2-59.8{(ad)}^3]=0.609$ is the empirical expression for the bulk drag coefficient of a random, emergent array of cylinders as obtained by Nepf (1999), and ${\mathsf{R}}_{h_v,D}=\overline{U_{h_v}}D/\nu =\mathrm{2,743}$ is the Reynolds number using the velocity at the canopy top ($\overline{U_{h_v}}$). The value of $ad$ for the present case is 0.071.

To derive the value of the depth-averaged drag coefficient, the expression in Eq. (8) is integrated as follows:which gives a depth-averaged drag coefficient of $C_{d,\mathrm{Nepf}}=0.386$. This was found to be in close agreement with the model predicted value of $C_d=0.371$, with a percentage difference of 3.9%. This analysis, in addition to providing an estimate of range of bulk drag forces, also validates numerically the method of Ghisalberti and Nepf (2004) for drag force estimation in submerged canopy flows and lends confidence to the present model, which can be used to supplement drag coefficient laboratory data sets of both submerged and emergent vegetation canopies in unidirectional flow in the future.

Six different wave conditions were simulated, as shown in Table 2, with each case run for an array of cylinders, an isolated cylinder, and a flume with no cylinders. The linear dimensions of the flume were the same (0.25 m wide, 0.75 m high, and $1.8+3L\mathrm{m}$ long, where $L$ = wavelength of the particular case) for a given case for the three setups, with the same refinement regions (near the free surface and the bottom boundary) as described earlier. Table 2 also shows the problem sizes in terms of the number of cells for each experiment. The degree of wave nonlinearity is represented by the Stokes’ parameter ($ϵ=kH/2$) and Ursell number ($Ur=H/h/{(h/L)}^2=H{L}^2/h^3$), where $H$ = wave height at the wavemaker; $h$ = uniform water depth; and $L$ = wavelength calculated by the dispersion relationship. For a given $h/L$ the degree of nonlinearity increases with an increase in $Ur$, whereas $ϵ$ gives an absolute estimate of the degree of nonlinearity among different $h/L$ cases, with higher values indicating greater nonlinearity. Velocity data were sampled at the four locations marked in Fig. 1(d) for both the array and isolated cylinder experiments. In addition to the isolated cylinder runs, velocity data were also sampled in the central vertical plane stretching from $-25$ to $25D$ with the cylinder placed at the origin. This was done to quantify the relative influence of the presence of the cylinder in modifying the wave field in either direction and compare it with the array cylinder velocity fields. Pressure ($p$) and velocities were collected for all the cells around the cylinders. Waves were started with a still-water condition in the flume, and it took approximately 18 wave periods from the start of the wavemaker to eliminate the effect of the initial conditions. The data were sampled at a frequency of $50/T\mathrm{Hz}$, where $T$ is the wave period of the wave conditions simulated over a period of 35 complete wave periods after the initial spin-up time of 20 wave periods. The choice of the averaging period was made because about 28–32 wave cycles were required for averaging before the phase-averaged force coefficients for the central cylinder showed less than 1% variation for a subsequent increase in number of averaging cycles. Thus, a total of 55 waves were introduced into the domain, and averaging was done over the last 35 wave periods. Because the goal in this paper is to understand the hydrodynamics and resultant forces associated with the mean wave energy, the data were analyzed and frequencies larger than the second-order peak of the surface elevation spectrum were filtered out using a fast Fourier transform algorithm in MATLAB.

To calculate the drag ($C_d$) and inertia ($C_m$) coefficients, Morrison’s equation (Sumer and Fredsøe 1997) as shown later was fitted to the sampled numerical phase-averaged total force ($F_{\text{tot}}$), calculated the same way as in Eq. (12), over one time period, using a nonlinear least-square fitting technique (López and García 1991)where $U(t)$ = phase-averaged velocity measured at the free surface, averaged at the sampling locations in Fig. 1. The least-square method minimizes the mean square error (${ϵ}^2$) between the measured and predicted forces to determine $C_d$ and $C_m$ and is expressed aswhere subscripts $n$ and $p$, respectively, represent the sampled (numerical) and predicted (using Morrison’s equation) values of the total force as in Eq. (15); and $N$ = total number of discrete data points. Using $\partial {ϵ}^2/\partial C_d=0$ and $\partial {ϵ}^2/\partial C_m=0$, the following two equations with two unknowns $C_d$ and $C_m$ are obtained, which were then solved using MATLABwhere $\dot{U}(t)$ = first derivative of the phase-averaged velocity at the free surface; and $f_d=\frac{1}{2}\rho C_{d}D$ and $f_m=\frac{\pi }{4}\rho C_m{D}^2$. The summation limits have been omitted for ease of reading and span from 1 to $N$, where $N$ is the total number of discrete intervals within a time period. $N=300$ was selected for all runs.

Fig. 4 shows the comparison of the variation of simulated drag and inertia coefficients with KC numbers ($\mathrm{KC}=U_{\max }T/D$, where $U_{\max }$ is the maximum orbital velocity from linear wave theory at the SWL) for the array and isolated cylinder cases with the experimental values along with the $\beta =\mathsf{R}/\mathrm{KC}$ numbers for each case. The agreement with the experimental data is good and shows that the model is able to simulate the fitted force coefficients accurately. Both $C_d$ and $C_m$ decrease with the KC number for the range of values tested for the $\beta =80$ case, whereas for the smaller $\beta =30$ the coefficients remain almost constant with KC showing only a reduction at $\mathrm{KC}=143$. The isolated cylinder shows marginally lower force coefficients than the corresponding array cylinder, whereas there is little overall difference between the array and isolated cylinder results in the simulations, which supports the observations of Hayashi and Chaplin (2012) who concluded that the spacing between the arrays ($10D$) was likely large for any significant vortex interactions. This is an interesting departure from the unidirectional flow results that show at similar spacings vortex interactions are quite active and influence the shedding dynamics of the cylinders.

Fig. 5(a) shows the time histories (about four wave cycles) of forces acting on the cylinders located along the centerline [Fig. 5(b) for exact cylinder locations] in the direction of wave propagation for Cases 1, 2, 5, and 6. Results for Case 3 are similar in trend to Cases 1 and 2, whereas those of Case 4 are similar in trend to Cases 5 and 6 and are not shown here. Although the RMS force magnitudes show little overall variation down the flume, the higher KC number cases show a greater variation, with the cylinders in the front of the flume (3 and 23) experiencing higher forces than those toward the back (83 and 103). There is a marginal increase in the force experienced by the central cylinder (53), which may be caused by a weak reflection from the back of the channel causing a partial standing wave antinode at this point. The isolated cylinder consistently experiences higher forces than the array cylinder at the same location. The column of panels on the right shows the force histories averaged over a single phase and time lagged over a single time window. It also clearly shows the variability in the force histories as well as the effect of the nonlinearity of the wave condition, particularly at higher KC numbers in which a clear inflection point is visible below the SWL line. The crest is more peaked than the trough at large KC numbers with the force direction showing an asymmetry with the wave phase.

The top row of panels in Fig. 6 shows the time series of pressure ($F_p$), viscous horizontal ($F_{v,x}$), viscous vertical ($F_{v,z}$), and the total force ($F_{\text{tot}}=F_p+F_{v,x}$) over about four wave cycles for Cases 1 and 6. The bottom row of panels shows the time series of the corresponding moments about the stem base. The surface elevation ($\eta$) is also plotted on a second $y$-axis on the right. The viscous forces are almost in phase with the surface elevation and, consequently, lead the pressure force by almost 90° at the lower KC number. This phase difference is, however, not so pronounced at the higher KC number. In terms of force magnitudes, it is found that at low KC numbers the viscous horizontal force contributes significantly to the total force, whereas at higher KC the contribution is negligible. The viscous vertical force is also higher than the viscous horizontal force at $\mathrm{KC}=9.1$. The viscous forces appear to be positive mostly with sharp peaks and shallow troughs. The phase difference between the pressure force and the surface elevation is near 90° at $\mathrm{KC}=9.1$, suggesting the contribution of the inertia effect, whereas the forces are almost in phase with the surface elevation at the higher KC. The addition of the viscous force to the pressure force decreases the phase difference somewhat between the total force and the surface elevation. The moments follow a similar trend to the forces except a change in scale as expected.

The upper row of panels in Fig. 7 shows the fitted drag and inertia forces according to Morrison’s equation at two representative KC numbers. The drag force and inertia force are 90° out of phase with each other because of the nature of Morrison’s equation with the inertia force leading the drag force. It is found that at lower $\mathrm{KC}=9.1$, the inertia force is significant in magnitude and is even greater than the drag force, whereas at the higher KC the drag force is much larger. This causes a temporal spreading of the total force for the lower KC number and results in a phase difference with the surface elevation. At the higher KC number in which the drag force is dominating and is almost in phase with the surface elevation, the total force follows a similar trend. The lower row of panels in Fig. 7 shows the time-lagged forces. The time lagging was done by calculating the correlation between the surface elevation time series and the force’s time series and finding out the lead or lag and shifting the force’s time series accordingly. It is seen that the effect of time lagging causes the predicted total forces to be almost double that of the measured total forces for the lower KC case, whereas the total force increases only slightly for the higher KC case. In most large-scale models it is a convenient practice to only consider time- and depth-averaged drag forces based on some a priori assumed drag coefficient value and to neglect the effect of the inertia forces for calculating the total forces acting on model plants because inertia forces are conservative (for linear wave assumption only). Based on the findings here it appears that although it may be fine to neglect the inertia effect at high KC numbers because they are too small, the same cannot be said at low KC numbers in which an overestimation of total forces might result by assuming drag contribution only.

In Fig. 8, the vertical variation of forces [Fig. 8(a)] and vertical variation of corresponding moments [Fig. 8(b)] acting on the central cylinder (53) in Fig. 5(b) for four representative cases (Cases 1, 2, 5, and 6) are shown. The crest or high-water level (HWL), SWL, and trough or low-water level (LWL) are marked by black dashed lines as well as the central zero line passing through the zero mark on the $X$-axis. For the force calculations, equations similar to Eq. (12) were used except the depth-integrating operation. Calculation of forces around the air-water interface can be difficult in VOF models because of the absence of a well-defined interface, thus to only consider fluid-filled regions, pressure ($p^{*}$), and effective viscosity (${\tau }_{\mathrm{eff}}$) values from cells with void fractions equal to or greater than 0.5 were used in the equations. Maximum and minimum forces and corresponding moments about the base of the cylinder on the flume bed signify the maximum positive and maximum negative forces acting on the cylinder over the entire wave record, whereas the RMS and mean quantities have their usual meaning. For all cases the peaks of the maximum (positive) force and moment occur just below the HWL but above the SWL. A similar phenomenon was observed by Torum (1989) who explained it by the analogy of unidirectional flow with free surface in which the pileup and increase of water level on the upstream side and the additional drawdown on the downstream face causes the pressure difference. As will be explained later in this section, a similar gradient in the mean velocity (and hence pressure) exists between the SWL and the HWL on either side of the cylinder in the wave direction due to the net positive velocity associated with the shoreward motion of the wave crest. The rising limb of the wave causes a recirculation in the horizontal plane, which spans vertically throughout between the wave crest and trough, whereas the falling limb of the wave has a recirculation primarily only below the SWL. This is also confirmed by the evidence that the minimum force (maximum of negative forces) almost vanishes between the SWL and HWL, suggesting that the recirculation region on the rising limb of the wave is stronger than the one developed during the falling limb. The peak of the RMS and minimum (maximum negative) forces and moments values occur just below the LWL. The mean values show a slight negative bias below the LWL but are entirely positive between the LWL and HWL with a peak lying at almost the same level as the maximum force. A local peak in the RMS value at that level is also noticed. A clear distinction can be seen between the higher and lower KC number cases. The lower KC number cases show a pronounced hyperbolic profile, whereas for the higher KC numbers the profiles are flatter in nature. Also interesting is the distinct boundary layer effect on the force profiles. The force profiles show a decrease within the wave bottom boundary layer, and this decrease is seen within a distance of 5–6 cm from the bed for the lower KC numbers to 2–3 cm for the higher KC numbers.

To analyze the hydrodynamics in the water column and investigate the spatial influence of the cylinder in the wave flow in the flume, the vertical and horizontal velocities at several locations in the central plane were compared with the array and isolated cases. Figs. 9 and 10 show the profiles of horizontal ($U_{\mathrm{rms}}$) and vertical ($W_{\mathrm{rms}}$) RMS of the wave velocities at 10 locations (five on each side of the cylinder, in which distances are made nondimensionalized by the cylinder diameter) including the four locations marked in Fig. 1(d) for the array case for Wave Cases 5 ($kH/2=0.027$) and 6 ($kH/2=0.054$). Here, for points between the crest and trough the temporal RMS operation considers velocity values from points below the water surface only. The results from the other cases show similar conclusions and are not shown here. The dotted dark black line shows the theoretical velocity profile at the center of the flume from the wave flume-only case (no cylinder) calculated using the inviscid Stokes’ second-order theory (Svendsen 2006). To calculate the Stokes’ theoretical profiles, the observed water surface elevation time series for each wave case was fitted for some unknown $x$ and $H$ against the Stokes’ second-order equation for surface elevation asand the best fit $x$ and $H$ were substituted in the equations for horizontal ($u_{\mathrm{stokes}}$) and vertical ($w_{\mathrm{stokes}}$) orbital velocities as

It was found that the horizontal and vertical velocity profiles at far-field locations $|X|=10D$ and $|X|=8.67D$ for the isolated cylinder case almost match the numerical results for the no-cylinder case, meaning that there is very little variation of the global wave flow field at these distances on either side of the cylinder. The array cylinder results (shown by green line) show that the velocities are slightly less than the isolated and no-cylinder cases because of the small damping effect of the cylinders in the front half of the canopy. Within $|X|\le 2D$ significant departure of the wave orbital horizontal RMS velocities is seen (Fig. 9). Below the wave trough or LWL the recirculation and turbulence action on either side of the cylinder decreases the wave excursions significantly compared with the no-cylinder case. The higher wave height case (Case 6) shows larger deviation caused by larger turbulence loss and thus greater wave orbital energy damping by the stem. Velocity profiles for both array and isolated cylinders in the region between the wave crest and the wave trough differ significantly on either side of the cylinder within the recirculation zone ($|X|\le 2D$). This apparent anomaly is explained further in Fig. 11, which shows the mean horizontal velocities at the same locations as in Fig. 9 for the two cases, with the temporal mean velocity calculated between the trough and the crest considering points below the water surface only. Below the wave trough the mean flow is almost zero for the far-field locations ($|X|=10D$ and $|X|=8.67D$), whereas within $X|\le 2D$, below the wave trough, the profiles show a net positive current, and behind the cylinder a net negative current caused by the recirculation region exists on either side of the cylinder. Between the wave trough and the crest, the far-field profiles show a net positive current caused by the forward motion of the wave crest that is entirely positive, which is expected for progressive waves. Near the cylinder, within the recirculation zone, is observed a decrease of this positive current behind the cylinder in the wave direction caused by the effect of the surface vortex (and hence negative velocities) created during the forward thrust of the crest behind the cylinder in the wave direction. This also explains the existence of a free surface velocity gradient on either side of the cylinder between the trough and the crest as was mentioned earlier. The RMS horizontal velocity profiles (Fig. 9) likewise show a reduction behind the cylinder between the trough and the crest compared with the points in front of the cylinder in the wave direction. The reduction in the free surface horizontal wave velocity seems to increase with wave nonlinearity (Case 5 versus Case 6), because the greater the nonlinearity the greater is the Stokes’ drift associated with the wave orbital excursion within the trough and the crest. In other words, because the time-averaged horizontal flow between the trough and the crest behaves as a uniform flow, with a vertical variation, it creates a mean damping effect similar to that observed in uniform flow around a cylinder, whereas the flow below the trough, having predominantly wave orbital signatures, is damped similar to wave damping in submerged vegetation. This separation of the hydrodynamic regime for an emergent vegetation problem has a potentially important influence on simplified quadratic law–based formulations commonly used in large-scale models and can be a topic of future research.

Vertical wave velocity RMS profiles (Fig. 10) show similar trends at far-field locations ($X=10D$ and $X=8.67D$) for array, isolated, and no-cylinder cases. Within the recirculation zone ($X\le 2D$) the profiles are similar, except near the free surface. In front of the cylinder, the obstruction by the cylinder causes run up on the seaward face during the forward motion of the wave crest followed by a localized recirculation vortex created below the SWL during the flow-reversal phase and impinging into the SWL. This causes a local increase in vertical velocity excursions. Behind the cylinder, because there is no impinging below the SWL of the surface vortex, the recirculation exists only when the crest is in shoreward motion causing lesser vertical velocity excursion amplitudes compared with that on the front face.

The previously mentioned free surface phenomenon is visualized in Fig. 12(a), in which the left panel represents Wave Case 5 and the right panel represents Wave Case 6. The figure shows from top to bottom, respectively, the mean horizontal, mean vertical, RMS horizontal, and RMS vertical velocities over 20 wave periods presented with the color bar on the right. The crest, mean, and trough locations are also indicated in the same figure. As explained in the line diagrams (Figs. 9–11), a net positive horizontal velocity flux above the SWL (top panels) is observed, which increases with increasing nonlinearity. Fig. 12(b) (bottom two rows) shows conceptually the wave height–integrated horizontal recirculation zones around the cylinder below and above the SWL during the rising limb (shoreward motion of wave crest) and falling limb (seaward motion of wave crest) of the free surface (locations shown relative to the cylinder in the top row). Thus, over a complete cycle, the excess flux above the SWL causes a uniform-flow type phenomenon with a pronounced recirculation zone behind the cylinder in the direction of wave propagation, particularly between the crest and the SWL, which makes the depth-integrated recirculation zone between the trough and the crest stronger in the shoreward direction of the cylinder. The shear effect of the flux also extends downward below the SWL. The mean vertical velocities in Fig. 12(a) show an asymmetrical signature with a net negative zone between the crest and the trough behind the cylinder and one between the trough and the SWL in front of it. This is caused by the impinging action of the surface vortex into the SWL during the flow reversal in the wave crest. The RMS horizontal velocity plots indicate the bulk of the free surface energy is damped between the SWL and wave crest, whereas below the trough the damping occurs on both sides with similar magnitudes similar to wave damping around submerged bluff bodies.

Vortex growth, shedding, and subsequent force variations near the free surface from an emergent surface-piercing cylinder has been a topic of active research (Kjeldsen et al. 2011; Torum 1989; Quetey et al. 2004; Barlas 2012). However, detailed visualization of the vortex shedding near the surface for an emergent cylinder has not been previously reported in the literature. Fig. 13 presents the 3D vortex visualizations from the isolated cylinder run corresponding to Case 5, as in Table 2, as a representative of the free surface damping phenomenon explained previously. Vortex cores are represented by the $Q$-criterion (Hunt et al. 1988), and isosurfaces are drawn at constant $Q=10$ with coloring done by vorticity (red-green-blue scale) about the $z$-axis with red representing a clockwise rotating vortex (positive) when viewed in the $z$-direction from the origin, and blue represents a counterclockwise (negative) vortex. The $Q$-criterion is mathematically represented aswhere $\Vert \mathbf{\Omega }\Vert =tr{[\mathbf{\Omega }{\mathbf{\Omega }}^t]}^{1/2}$; and $\Vert \mathbf{S}\Vert =tr{[\mathbf{S}{\mathbf{S}}^2]}^{1/2}$ = asymmetric and symmetric parts of $\nabla \overline{\mathbf{u}}$ defined as

A vortex region is said to exist in a flow field if $Q>0$ in the cells comprising that region and may be thought of as local regions having a greater rotation rate than strain rate. Figs. 13(a–h) show the instantaneous vortex isosurfaces generated near the free surface from the cylinder at eight equal intervals of time along a representative wave period ($T$). The figures have been magnified to show the region near the vicinity of the free surface with black lines marking the crest elevation or HWL, the mean water level or SWL, and trough elevation or LWL. At $t=T/8$ [Fig. 13(a)], the free surface is accelerating and descending downward, as shown by the arrow direction. Surface Vortex A, having its axis almost parallel to the SWL and forming a circular arc around the cylinder, is the remnant of the vortex shed in the preceding cycle, whereas two counterrotating vortex pairs (B-B’) that also shed in the preceding cycle are seen. A violent free surface also exists near the cylinder-water interface as well as strong isolated vortex streaks formed just below the free surface. At $t=T/4$ [Fig. 13(b)], as the downward acceleration of the free surface has reached a peak, Vortex A has formed almost a helical shape, is impinging below the free surface, and appears to be dissipating. This vortex aligns itself mostly between the SWL and somewhat below the LWL and is responsible for the increase in vertical velocity excursions in this zone. The counterrotating pair has weakened and is ready to be flipped over the cylinder. The rotation within the boundary layer also has changed direction with a weak vorticity structure developing. At $t=3T/8$ [Fig. 13(c)], as the downward motion of the free surface decelerates Vortex A is seen to be dissipating further and moving downward in a helical pattern. B-B’ vortices have flipped sides at this point and also have been pushed down with the top part almost dissipating under the strong vertical velocities. At $t=T/2$ [Fig. 13(d)], as the wave is about to change direction, Vortex A has almost dissipated with only weak helical signatures surrounding the cylinder. B-B’ has completely dissipated from the window at this stage. At $t=5T/8$ [Fig. 13(e)], as the free surface accelerates upward, the counterrotating vortex pair (C-C’) that was developing within the boundary layer since $t=T/4$ has now shed completely. It is worth noticing that the vortex is almost gone except for a few streaks near the free surface. The rise of the free surface actually causes a stretching of the counterrotating pair and they reach very near to the free surface. No prominent surface vortices are visible, unlike what was seen at $t=T/8$. At $t=3T/2$ [Fig. 13(f)], the free surface has reached the highest peak of upward acceleration and causes significant stretching and thinning of the C-C’ pair. Surface Vortex A has started shedding and is present only as a few streaks. The recirculation of this zone associated with the rising limb gives rise to the negative velocity zone behind the cylinder in the wave direction and will cause dampening of the surface energy in the subsequent two-phase instants. C-C’ is about to flip over the cylinder, and new counterrotating vortex signatures can be seen developing within the boundary layer, which will eventually become the B-B’ pair in the subsequent cycle. At $t=7T/8$ [Fig. 13(g)], the free surface rise has decelerated the detachment of the C-C’ vortex pair, and formation of a new A vortex is seen. Stronger vortices are formed under the free surface on the lee side of the cylinder during the rising limb of the wave than during the falling limb. This causes a stronger recirculation region in the wave direction and, consequently, higher positive forces on the cylinder under the SWL during the rising cycle of the wave than during the falling limb. It also explains the asymmetrical surface velocities observed previously on either side of the cylinder. Finally, at $t=T$ [Fig. 13(h)] C-C’ is completely dissipated while Vortex A aligns itself in a circular arc around the cylinder near the free surface.