3.1. Instantaneous flow features within the roughness elements

To assess the impact of the roughness elements on the near-wall flow features, figure 3 compares instantaneous horizontal velocity contours for the four cases discussed in this manuscript. Owing to the small non-dimensional roughness $\bar {k}_s^{+}=u_*\bar {k}_s/\nu = 6$ (formally, roughness is hydraulically smooth when $\bar {k}_s^{+} \le 4$ (Jiménez Reference Jiménez2004)), there are no obvious separation regions for the steady flow case with bumps. Furthermore, despite the stronger velocity fluctuations for the wave-current case with bumps (compare the flat and bumpy wave-current cases at $\omega t = 8{\rm \pi} /5$ in figure 3), based on the small value of $KC=u_*/(\omega \bar {k}_s)=1.64$, there are very weak flow separation features when bumps are added to the wave-current case (at $\omega t = 3{\rm \pi} /5$). This is consistent with Nielsen (Reference Nielsen1992), who points out that the relative importance of form drag (or flow separation induced pressure drag) decreases with decreasing Keulegan–Carpenter number. Furthermore, as discussed in § 3.5, the contributions due to the dispersive stresses that are associated with flow separation are minimal compared with the contributions from other flow features that are observed. Rather than exhibiting flow separation which would imply an increase in the form drag by the roughness elements, the flow in the roughness elements is dominated by viscous effects which essentially decelerate the flow, much like the flow in a wave-current, rough-wall boundary layer that resembles a canopy flow (Egan et al. Reference Egan, Cowherd, Fringer and Monismith2019).

Figure 3. Contours of instantaneous streamwise velocity ($U_1/u_*$) for the four cases discussed in this manuscript along the streamwise transect at $x_2^{+}= 330$. Panels (a–d) correspond to case C350F, (e–h) correspond to case C350B, (i–l) correspond to case WC350F and (m–p) correspond to case WC350B. For the steady flow cases (i.e. C350F and C350B), the four panels correspond to different instantaneous turnover times $t^{+}=u_* t/H$, while the four panels for the wave-currents cases (i.e. WC350F and WC350B) correspond to different wave phases $\omega t$ relative to time $t^{+}=20$. The magenta contours indicate $U_1/u_*=0$, while the red arrows in (i–l) indicate the mean wave velocity $\sin (\omega t)$ for both wave-current cases. The roughness elements are indicated by the red regions in the figures. Transects at other locations and times exhibit similar behaviour lacking obvious flow separation features.

The viscous effects in the bumps merely alter the behaviour of the Stokes boundary layer in the oscillatory flow. During the acceleration phase of the wave cycle (${\rm \pi} /10 \leq \omega t \leq 3{\rm \pi} /5$), although the boundary layer thickness is substantially reduced for both wave-current cases when compared with the steady cases, it is thicker for the wave-current case with bumps. Similarly, during the deceleration phase of the wave cycle (${\rm \pi} \leq \omega t \leq 8{\rm \pi} /5$), there is a subsequent increase in the boundary layer thickness for both wave-current cases, although again it is thicker in the presence of bumps. This behaviour of the boundary layer flow that is out of phase with the driving pressure gradient is consistent with the modified Stokes boundary layer depicted in figure 10. Despite the weak flow separation induced by the roughness features, they induce stronger flow variability and velocity fluctuations which produce stronger velocity shear and Reynolds stress, thus increasing the bottom drag coefficient as discussed in what follows.

3.2. Measure of convergence for turbulent statistics

The steady flat and bumpy wall cases C350F and C350B are simulated for a total of 40 eddy turnover times ($T_\epsilon$). The first $10 T_\epsilon$ involve an initial transient after which the flow requires approximately $20T_\epsilon$ for the turbulent kinetic energy (TKE) to equilibrate, beyond which the turbulence statistics are computed with $T_{avg} = 10T_\epsilon$. The initial conditions for these cases are given by

(3.1)\begin{gather} u_1(x_1,x_2,x_3) = 2 U_c \left( 1 - x_3/H \right) + \alpha_{TI} R(x_1,x_2,x_3), \end{gather}

(3.2)\begin{gather} u_2(x_1,x_2,x_3) = \alpha_{TI} R(x_1,x_2,x_3), \end{gather}

(3.3)\begin{gather} u_3(x_1,x_2,x_3) = \alpha_{TI} R(x_1,x_2,x_3), \end{gather}

where $U_c$ is the mean of the velocity profile given composed of the viscous sublayer and log law,

(3.4)\begin{equation} U_c = \frac{1}{H} \left\{ \int_{0}^{11.6} x_3^{+} \,{\rm d}x_3^{+} + \int_{11.6}^{\textit{Re}_{*}} \left[ \frac{u_{*}}{\kappa} \ln \left( x_3^{+} \right) + 5.2 \right] \,{\rm d}x_3^{+} \right\}, \end{equation}

which ensures that the volume- and time-averaged velocity is given approximately by $U_c$. The magnitude of the initial perturbations is $\alpha _{TI} = 0.05$, $R(x_1,x_2,x_3)$ is a uniformly distributed random number in the range $[-1,1]$, and $\kappa = 0.4$ is the von Kármán constant. Figure 4 shows the time evolution of the friction Reynolds number which converges after roughly $25T_\epsilon$, after which the deviations of the moving average over one eddy turnover from the target value are less than $0.5\,\%$. The flow transitions to a turbulent state at $3T_\epsilon$, subsequently leading to the evolution of streamwise turbulent structures that require approximately $20 T_\epsilon$ to reach the target levels dictated by the driving pressure gradient. The linear stress profile in figure 5 is used as an indicator for the level of convergence and validation of the results attained in the DNS. The corresponding streamwise velocity profile also follows the linear and log-law analytic predictions when averaged over $10T_{\epsilon }$ following the initial transient of $30T_{\epsilon }$, as shown in figure 6.

Figure 4. Time evolution of $Re_*$ for case C350F.

Figure 5. Time- and planform-averaged viscous stress, Reynolds stress and total stress profiles for case C350F. The time averaging is carried out over $10T_{\epsilon }$ for $t/T_{\epsilon } \geq 30$.

Figure 6. Comparison of time- and planform-averaged, streamwise velocity profile for case C350F against the linear and log-law analytic expressions. The time averaging is carried out over $10T_{\epsilon }$ for $t/T_{\epsilon } \geq 30$.

Once the steady cases reach equilibrium after $30T_{\epsilon }$, the three-dimensional velocity fields are used as initial conditions for the wave-current cases. The convergence of time-averaged energetics in the wave-current cases can be estimated by observing the departure of the time rate of change of TKE from zero. This convergence criterion for the flat wall case can be formulated using the volume-integrated TKE balance equation given by

(3.5)\begin{equation} \frac{\partial \langle k\rangle_v }{\partial t} ={-}\langle P_k \rangle_v + \langle \epsilon_k \rangle_v, \end{equation}

where $\partial \langle k\rangle _v / \partial t$ is the time rate of change of volume-integrated TKE, $\langle P_k \rangle _v$ is the volume-integrated TKE production by mean shear and $\langle \epsilon _k \rangle _v$ is the volume-integrated dissipation of TKE. Time-averaging equation (3.5) over an integer number of wave periods implies that the time-averaged TKE production is balanced by the time-averaged dissipation of TKE. Thus, we can conclude that the wave-current system reaches equilibrium once $\partial \langle k\rangle _v / \partial t = 0$ in a time-averaged sense. Based on the definitions presented above, the flat wall, wave-current case appears to converge after around $40T_w$, as seen in figure 7. This suggests that the flow in the flat wall, wave-current case requires an additional $6T_\epsilon \sim 40T_w$ after the introduction of the oscillatory pressure gradient to reach equilibrium, as defined earlier. As for the bumpy wall, wave-current case, the system reaches equilibrium within $20T_w$ or $3T_\epsilon$, which is consistent with the elevated bottom drag and dissipation for flows over bumpy walls. The convergence of flow statistics has been verified by integrating the governing equations for an additional 200 wave periods and it was found that the statistics change by no more than $1\,\%$. In what follows, the turbulent statistics are computed with $T_{avg} = 100T_w$ after an initial transient of 100 wave periods to avoid any residual transitional effects.

Figure 7. Time evolution of the time rate of change of the volume-integrated TKE. (a) Case WC350F and (b) case WC350B. The blue dash-dot lines correspond to the instantaneous value, the solid black lines correspond to the cumulative mean starting from $t/T_w=0$ until the instantaneous value, while the red dashed lines correspond to the mean over the entire time series.

3.3. Mean and wave-driven velocity profiles

Figure 8 shows that the flat wall cases C350F and WC350F follow the log law, with minimal differences observed between the two cases. Close to the wall ($x_3^{+} < 6$), the two cases are identical, while minor differences can be observed outside the buffer layer ($x_3^{+} > 30$). These results are consistent with the findings of Lodahl et al. (Reference Lodahl, Sumer and Fredsøe1998), Manna et al. (Reference Manna, Vacca and Verzicco2012) and Nelson & Fringer (Reference Nelson and Fringer2018) that the addition of a laminar wave to a turbulent current over flat walls in the current-dominated flow regime does not alter the mean velocity profile significantly. The bumpy wall log law is given by (Raupach et al. Reference Raupach, Antonia and Rajagopalan1991)

(3.6)\begin{equation} \frac{\langle \bar{u}_1 \rangle}{u_*} = \frac{1}{\kappa} \ln \left( \frac{x_3 - \bar{k}_s}{z_0} \right), \end{equation}

where $\bar {k}_s$ is known a priori (see figure 2) and $z_0 = \bar {k}_s/ \alpha _k$, where $\alpha _k$ is the regression factor used to fit the velocity profiles. The values for $u_*$ and $\kappa$ are constant, as discussed earlier. Regressing for $\alpha _k$ for cases C350B and WC350B yields $\alpha _k = 26$ and $\alpha _k = 19$, respectively. For conventional sand-grain type roughness, Nikuradse (Reference Nikuradse1932) proposed $\alpha _k = 30$, where smaller values of $\alpha _k$ imply a larger effective roughness $z_0$. It is also crucial to note that the sand-grain type roughness discussed in Nikuradse (Reference Nikuradse1932) has a $C_o = 1.0$. As seen in figure 8, the mean velocity profile for WC350B shifts further away from the flat wall log law when compared with C350B. Unlike the flat wall cases, the bumpy wall cases are expected to show enhanced drag coefficients when compared with the flat wall cases. This drag coefficient enhancement is a function of the Corey shape factor presented in (2.6). As suggested by Julien (Reference Julien2010), with decreasing $C_o$ the drag coefficient should increase. Therefore, case C350B with a $C_o = 0.6$ exhibits a larger drag coefficient when compared with the sand-grain type roughness elements with a $C_o = 1.0$ (Raupach & Thom Reference Raupach and Thom1981; Ghodke & Apte Reference Ghodke and Apte2017). The smaller value of $\alpha _k$ observed for case WC350B suggests that despite the laminar nature of the wave, it greatly affects the mean flow for bumpy wall cases.

Figure 8. Time- and planform-averaged velocity profiles for the flat and bumpy wall cases.

Table 2. Comparison of drag coefficient computed in the present study to wave-current experiments by Fredsøe (Reference Fredsøe1984), Myrhaug & Slaattelid (Reference Myrhaug and Slaattelid1990) and Huynh-Thanh & Temperville (Reference Huynh-Thanh and Temperville1991), in addition to the steady current flat wall DNS cases by Moser et al. (Reference Moser, Kim and Mansour1999) and del Álamo & Jiménez (Reference del Álamo and Jiménez2003). Wave-current experimental data adapted from Soulsby et al. (Reference Soulsby, Hamm, Klopman, Myrhaug, Simons and Thomas1993).

A quantitative measure of the effects of waves is given through a drag coefficient defined as

(3.7)\begin{equation} C_d^{*} = \frac{u_*^{2}}{\langle \bar{u}_1 \rangle_v^{2}}. \end{equation}

Obtaining vertical profiles of the velocity in experimental or in situ studies of wave-current boundary layers may not always allow computing the drag coefficient using (3.7). In this case, the drag coefficient can be defined by assuming a log-law velocity profile, which gives

(3.8)\begin{equation} C_d = \left[ \frac{\kappa}{\ln(H/z_0) - 1} \right]^{2}, \end{equation}

where $z_0 = \bar {k}_s/\alpha _k$ for bumpy walls and $z_0 = \nu /(9u_*)$ for flat walls. As shown in table 2, $C_d$ and $C_d^{*}$ for the flat wall case C350F are comparable to previous numerical studies (Moser, Kim & Mansour Reference Moser, Kim and Mansour1999). Case WC350F, as expected, does not show any substantial changes in the drag coefficient. However, there is a significant increase in the drag coefficient with the addition of bumps when compared with the flat wall case. Case WC350B surprisingly shows enhanced drag coefficient when compared with case C350B, despite the laminar nature of the waves. Comparison of $C_d$ for case WC350B against the experimental studies by Fredsøe (Reference Fredsøe1984), Myrhaug & Slaattelid (Reference Myrhaug and Slaattelid1990) and Huynh-Thanh & Temperville (Reference Huynh-Thanh and Temperville1991) with similar $z_0/H$ values supports the validity of the numerical results. However, we note that the flow conditions in the experimental studies listed in table 2 correspond to wave-dominated flow conditions since $U_b/\langle \bar {u}_1 \rangle _v > 1$ for all cases. These results suggest that the principal conclusion of drag coefficient enhancement under the combined action of waves and currents over rough walls from the wave-current turbulence model developed by Grant & Madsen (Reference Grant and Madsen1979) hold even for weak wave-current interactions in the current-dominated flow regimes. This is despite the fact that the wave-current model proposed by Grant & Madsen (Reference Grant and Madsen1979) does not apply in the hydraulically smooth wall and current-dominated flow conditions. Using the Grant & Madsen (Reference Grant and Madsen1979) wave-current friction factor formulation, it can be shown that the lower limit for its validity is $U_b/ \langle \bar {u}_1 \rangle _v = 0.59$, while $U_b/ \langle \bar {u}_1 \rangle _v = 0.38$ for case WC350B. These results suggest that $C_d^{*}$ is a good predictor of the drag coefficient as long as $z_0/H$ is accurately estimated as a function on ${\textit {Re}}_*$, ${\textit {Re}}_w$, $H/\bar {k}_s$ and $u_*/(\omega \bar {k}_s)$. However, the central challenge is to accurately estimate $z_0/H$ as a function of the relevant non-dimensional parameters as discussed in Grant & Madsen (Reference Grant and Madsen1979). Figure 9 supports this observation, and suggests that the differences observed between the drag coefficient presented in (3.7) and (3.8) can be attributed to the mismatch of $\langle \bar {u}_1 \rangle _v(x_3 = H)$ and the log-law prediction. This mismatch is expected to increase with increasing ${\textit {Re}}_*$, due to the presence of a wake region. Wave-driven turbulence, stratification and inclusion of a buffer region may also affect this comparison (Egan et al. Reference Egan, Manning, Chang, Fringer and Monismith2020).

Figure 9. Comparison of the drag coefficient computed in the present study against experimental, numerical and analytic expressions.

Figure 10 shows the differences between the phase-averaged wave velocity for the flat and bumpy wall cases. The analytical solution proposed by Stokes (Reference Stokes1851) compares well with the DNS after regressing to determine the wall height. The DNS predicts slightly enhanced peak wave velocities when compared with the analytic solution, although the solutions agree well far from the wall. These results suggest that the wave and the turbulent flow fields are decoupled and any correlations between these two flow components will be small for the current-dominated flow system. This also suggests one-way coupling in which the waves affect the turbulence but the turbulence does not affect the waves (Manna et al. Reference Manna, Vacca and Verzicco2012).

Figure 10. Phase- and planform-averaged wave velocity for the bumpy wall case compared with the Stokes wave solution. The black solid line corresponds to the WC350B case while the red dashed line corresponds to the Stokes wave solution. The black dash-dot line denotes the roughness crest level ($k_c$). The vertical coordinate is normalised by the Stokes boundary layer thickness $\delta _w = \sqrt {2\nu /\omega }$.

3.4. Time- and planform-averaged stress profiles

The time- and planform-averaged momentum balance for the wave-current channel flow system above the roughness crest is given by

(3.9)\begin{equation} 0 = \varPi_c + \frac{\partial \langle \bar{\tau} \rangle}{\partial x_3}, \end{equation}

where the oscillatory forcing vanishes due to time averaging, $\langle \bar {\tau } \rangle$ represents the time-averaged and planform-averaged total stress, which is the sum of the viscous and the Reynolds stress, viz. $\langle \bar {\tau } \rangle = \nu ({\partial \langle \bar {u}_1 \rangle }/{\partial x_3}) - \langle \overline {u_1^{\prime }u_3^{\prime }} \rangle$. Integrating (3.9) from some height $x_3$ to the channel height $x_3=H$ (where $\langle \bar {\tau } \rangle = 0$) and substituting $\varPi _c = u_*^{2}/H$ gives

(3.10)\begin{equation} \frac{\langle \bar{\tau} \rangle}{u_*^{2}} = 1 - \frac{x_3}{H}. \end{equation}

This result suggests that the time- and planform-averaged stress profile above the roughness elements for the wave-current bumpy wall case follows the traditional linear stress profile. Figure 11 shows the time- and planform-averaged stress profiles for the four cases discussed earlier. Above the roughness crest, the total stress profiles collapse onto the analytic solution presented in (3.10). For the bumpy wall cases below the roughness crest, the immersed boundary force (not shown in figure 11) accounts for the deficit stress. This immersed boundary force cannot be easily computed due to the irregular nature of the roughness elements, although it has been validated in previous studies (Scotti Reference Scotti2006; Yuan & Piomelli Reference Yuan and Piomelli2014). There are no significant differences between cases C350F and WC350F. For the bumpy wall cases, there is appreciable attenuation of the streamwise root-mean-square (r.m.s.) velocity $( u_{i,{rms}} = \sqrt {\langle \overline {u_i^{\prime 2}}\rangle } )$ for the wave-current case when compared with its steady current counterpart. This attenuation is limited to a region above the roughness crest throughout the effective buffer layer region ($5 \lesssim x_3^{+} -\bar {k}_s^{+} \lesssim 30$). Below the roughness crest, the wave-current case shows elevated streamwise r.m.s. velocities compared with the other cases. Similar trends are observed for case C350B, suggesting that compared with the flat wall, the roughness elements act as a momentum sink that is stronger for case WC350B when compared with case C350B. Above the roughness elements, the log-law region shows the opposite trend for the streamwise r.m.s. velocity, with case WC350B showing larger values compared with the other cases. For the spanwise and vertical r.m.s. velocity profiles, the general trend compared with that of the streamwise r.m.s. velocity profile is opposite, suggesting that there are exchanges between the streamwise, spanwise and vertical momentum fluctuations.

Figure 11. (a) Comparison of time-averaged and planform-averaged viscous $( \nu ({\partial \langle \bar {u}_1 \rangle }/{\partial x_3}))$, Reynolds $( - \langle \overline {u_1^{\prime } u_3^{\prime }} \rangle )$ and total stress $( \langle \bar {\tau } \rangle )$. Red circles represent case C350F, black squares represent case C350B, green asterisks represent case WC350F and magenta diamonds represent case WC350B. Solid lines indicate total stress while the blue dashed line represents (3.10). The dashed horizontal line marks the location of the top of the roughness elements ($k_c^{+}$) on both the panels.

The vertical and wave-phase variation of the TKE $( \langle \widetilde {u_i^{\prime } u_i^{\prime }} \rangle )$ and Reynolds stress $( -\langle \widetilde {u_1^{\prime } u_3^{\prime }} \rangle )$ are shown in figure 11. Both wave-current cases show similar wave-phase dependence with varying degrees of attenuation as a function of wave phase and distance away from the wall. Consistent with the time- and planform-averaged r.m.s. profiles in figure 11, case WC350B shows slightly lower levels of TKE throughout the wave cycle when compared with case WC350F. The TKE peaks during the deceleration $(({{\rm d}}/{{\rm d}t}) U_b \sin (\omega t) < 0 )$ portion of the wave cycle similar to the experimental studies by Hino, Sawamoto & Takasu (Reference Hino, Sawamoto and Takasu1976) and Fishier & Brodkey (Reference Fishier and Brodkey1991). As detailed by Bae & Lee (Reference Bae and Lee2021), these enhanced TKE levels are indicative of bursting events that transport the TKE away from the wall, as bursting is associated with strong turbulent ejection events (second quadrant events in the $u_3^{\prime }$ vs $u_1^{\prime }$ plot). Figure 12 shows the difference between the second quadrant ejection events that are characteristic of the streak lifting or breakdown events. Here $A_1 A_3 = (u_1^{\prime } | u_1^{\prime } < 0) (u_3^{\prime } | u_3^{\prime } > 0)$ is the conditional product that is time averaged and planform averaged for the two wave-current cases. The bumpy wall case shows elevated $\langle \overline {A_1 A_3} \rangle$ magnitudes when compared with the flat wall case, and these events are concentrated during the deceleration portion of the wave cycle (not shown).

Figure 12. Time- and planform-averaged conditional product of the streamwise and vertical turbulent velocity component ($u_1^{\prime }u_3^{\prime }$) comparison for the flat and bumpy wall, wave-current cases.

Once the TKE within the roughness elements is enhanced at the expense of the TKE above the roughness elements, it becomes available for conversion to the Reynolds stress. As seen in figure 13, the Reynolds stress is larger during the acceleration $(({{\rm d}}/{{\rm d}t}) U_b \sin (\omega t) > 0 )$ part of the wave cycle. This suggests that the enhancements in TKE as a result of the turbulent bursting events in the deceleration wave cycle are transported towards the wall during the accelerating part of the wave cycle. While the TKE enhancement occurs in the effective buffer region (defined earlier), the Reynolds stress enhancements appear to occur mainly outside the buffer layer. This phase evolution supports the life cycle of buffer layer streaks proposed by Bae & Lee (Reference Bae and Lee2021), where the turbulent streaks are generated within the buffer layer and transported away from the wall during the ejection events, followed by sweep events that transport the associated TKE towards the wall. However, as Bae & Lee (Reference Bae and Lee2021) rightly note, this process can only provide a diagnostic explanation for the observations presented in figure 13.

Figure 13. Phase- and planform-averaged TKE (a,b) and Reynolds stress (c,d) as a function of height and wave phase. Panels (e,f) show the wave orbital velocity ($U_b \sin (\omega t$)) which is out of phase with the oscillatory pressure gradient ($U_b \omega \cos (\omega t$)) on the right-hand side of (2.5). Panels (a,c,e) are the flat wall case WC350F and (b,d,f) are the bumpy wall case WC350B. The horizontal dash-dot line marks the roughness crest level ($k_c^{+}$) while the horizontal dashed line marks twice the turbulent Stokes length offset above the roughness crest ($2l_t^{+} + k_c^{+}$).

Interactions between the TKE and Reynolds stress can be analysed with the structure parameter defined as

(3.11)\begin{equation} \langle \tilde{M}_r \rangle ={-}\frac{\langle \widetilde{u_1^{\prime} u_3^{\prime}} \rangle}{\langle \widetilde{u_i^{\prime} u_i^{\prime}} \rangle}, \end{equation}

which is a measure of the efficiency with which the turbulent eddies extract the Reynolds stress from the available TKE (Scotti & Piomelli Reference Scotti and Piomelli2001). Figure 14 shows the phase- and planform-averaged structure parameter comparison for the two wave-current cases (i.e. WC350F and WC350B) normalised by case C350F. Case WC350F shown in figure 14(a) exhibits elevated values of $\langle \tilde {M}_r \rangle$ during the acceleration part of the wave cycle, followed by lower values of $\langle \tilde {M}_r \rangle$ during the deceleration part of the wave cycle. These wave-phase variations occur in a region close to the wall and below $2 l_t^{+}$ (defined in (2.16)). The enhancement that occurs in $\langle \tilde {M}_r \rangle$ close to the wall during the acceleration wave phase is not equal in magnitude to that occurring during the deceleration wave phase despite the sinusoidal nature of the wave forcing in (2.5). Far from the wall, $\langle \tilde {M}_r \rangle$ does not vary with the wave phase as it does closer to the wall. As lower values of $\langle \tilde {M}_r \rangle$ are indicative of non-equilibrium flows (Scotti & Piomelli Reference Scotti and Piomelli2001), case WC350F exhibits a near-wall region $x_3^{+} < 2 l_t^{+}$ where the flow is in relative non-equilibrium compared with the outer region $x_3^{+} > 2 l_t^{+}$. Since the flow is not in equilibrium within the inner region ($x_3^{+} < 2 l_t^{+}$), it affects the nonlinear production cycle (Jiménez & Moin Reference Jiménez and Moin1991) as the inner region includes the location where peak production occurs ($x_3^{+} \approx 11.8$).

Figure 14. Wave-phase variations of the normalised structure parameter $( \langle \tilde {M}_r \rangle / \langle \tilde {M}_r^{C350F} \rangle )$ comparing the two wave-current cases (a) WC350F and (b) WC350B. The horizontal dashed line marks the roughness crest level, while the horizontal dash-dot line marks twice the turbulent Stokes length offset above the roughness crest $(2l_t^{+} + k_c^{+} )$.

As for case WC350B, higher values of $\langle \tilde {M}_r \rangle$ are generated during the acceleration portion of the wave cycle primarily below the roughness crests and transported above them during the deceleration portion of the wave cycle. Most of the $\langle \tilde {M}_r \rangle$ enhancements compared with case C350F occur below $2l_t^{+}$, suggesting that the effect of waves is confined to a region close to the wall. Similar to case WC350F, case WC350B also exhibits asymmetric $\langle \tilde {M}_r \rangle$ modifications as a function of wave phase within the roughness region ($x_3^{+} < k_c^{+}$), with strongly enhanced values during the acceleration phase and strongly attenuated values during the deceleration phase. Above the roughness region and below twice the turbulent Stokes length scale ($k_c^{+} \leq x_3^{+} < 2 l_t^{+}$), minor modifications of $\langle \tilde {M}_r \rangle$ can be observed with respect to the wave phase and increasing $x_3^{+}$. Further away from the wall ($x_3^{+} > 2l_t^{+}$), $\langle \tilde {M}_r \rangle$ does not change significantly as a function of the wave phase and $x_3^{+}$, similar to case WC350F.

The time- and planform-averaged vertical profiles of $\langle \bar {M}_r \rangle$ for the four cases are shown in figure 15. It is evident that the time-averaged effects for cases WC350F and WC350B are significantly different. Cases C350F and WC350F show similar vertical $\langle \bar {M}_r \rangle$ profiles, while case WC350B shows enhanced $\langle \bar {M}_r \rangle$ values below the roughness crests and twice the turbulent Stokes length scale. This further illustrates the utility of $2 l_t^{+}$ as a length scale that predicts the height below which the effects of waves is dominant in a wave-current boundary layer type flow. In the outer layer ($x_3^{+} > 2l_t^{+}$), cases WC350F, C350B and WC350B show reduced values of $\langle \bar {M}_r \rangle$ when compared with case C350F. This suggests that these cases do not effectively convert the available TKE to Reynolds stress in the outer region (or the log-law region) consistent with the observations of Scotti & Piomelli (Reference Scotti and Piomelli2001). Comparison of the vertically integrated $\langle \tilde {M}_r \rangle$ profile as a function of wave phase within the three regions defined in (2.17)–(2.19) is shown in figure 16. Within region I, the vertically integrated structure parameter for case WC350F is sinusoidal and in phase with the wave velocity, as seen in figure 16. Case WC350B, on the other hand, shows a prominent peak during the acceleration phase, but does not show an equally strong negative peak during the deceleration phase of the wave cycle. This asymmetry leads to a higher $\langle \bar {M}_r \rangle$ for case WC350B when compared with case WC350F in region I. In region II a similar trend is observed when cases WC350F and WC350B are compared with each other. The phase dependence now is in sync with the driving pressure gradient instead of the wave velocity as observed in region I. One crucial difference between regions I and II is that the phase dependence for case WC350B in region II is sinusoidal. Finally, in region III the structure parameter is wave-phase independent, and case WC350B shows a slightly smaller magnitude when compared with case WC350F. The time- and planform-averaged, vertically integrated values of the structure parameter ($\langle \bar {M}_r \rangle _v$) shown in table 3 further support the net behaviour for the four cases under consideration. Both the bumpy wall cases C350B and WC350B show a net increase in the structure parameter within the roughness region (i.e. below the roughness crest $x_3^{+} < k_c^{+}$) at the expense of the other two regions. This is despite the relative increase in $\langle \bar {M}_r \rangle _v$ for cases C350B and WC350B when compared with C350F and WC350F. These results suggest that the effect of waves is confined to a region close to the wall $(x_3^{+} < 2l_t^{+})$, while far away from the wall $(x_3^{+} > 2l_t^{+})$ there are minor changes observed in the turbulence statistics. Additionally, the strong wave-phase dependence of $\langle \tilde {M}_r \rangle$ for case WC350B when compared with WC350F suggests that the turbulence is out of equilibrium for a portion of the cycle (i.e. acceleration phase) and equilibrates during the deceleration portion of the wave cycle.

Figure 15. Time- and planform-averaged comparison of $\langle \bar {M}_r \rangle$ for the four cases under consideration in this paper. The vertical dash-dot line represents the roughness crest level while the vertical dashed line represents twice the turbulent Stokes length offset above the roughness crest level ($2l_t^{+} + k_c^{+}$).

Figure 16. Comparison of the time- and planform-averaged, vertically integrated structure parameter over the three regions for the two wave-current cases as a function of wave phase.

Table 3. Time- and planform-averaged, vertically integrated structure parameter ($M_r$) comparison across three regions. The first region is below the top of the roughness elements ($k_c^{+}$), the second region is between the top of the roughness elements and twice the turbulent Stokes length scale ($l_t^{+}$) defined in (2.16), and the third region is above $2l_t^{+}$. The percentages next to the numerical values represent the contribution for the corresponding region towards the total value.

3.5. Phase-averaged TKE and Reynolds stress budgets

The detailed mechanisms governing the behaviour of the structure parameter can be understood by analysing the phase-averaged TKE and Reynolds stress budgets. The phase-averaged Reynolds stress budget for channel flow geometries over bumpy walls can be written as

(3.12)\begin{equation} \omega \frac{\partial \langle \widetilde{u_i^{\prime} u_k^{\prime}} \rangle}{\partial \theta } = P_s + P_w + P_r + T_k + T_r + T_p + T_v + S_k - \epsilon, \end{equation}

where the production by mean shear is

(3.13)\begin{equation} P_s ={-} \langle \widetilde{u_i^{\prime} u_3^{\prime}} \rangle \frac{\partial \langle \bar{u}_k \rangle}{\partial x_3} - \langle \widetilde{u_k^{\prime} u_3^{\prime}} \rangle \frac{\partial \langle \bar{u}_i \rangle}{\partial x_3}, \end{equation}

the production by wave shear is

(3.14)\begin{equation} P_w ={-} \langle \widetilde{u_i^{\prime} u_3^{\prime}} \rangle \frac{\partial \langle \tilde{u}_{w,k} \rangle}{\partial x_3} - \langle \widetilde{u_k^{\prime} u_3^{\prime}} \rangle \frac{\partial \langle \tilde{u}_{w,i} \rangle}{\partial x_3}, \end{equation}

the production by roughness-induced fluctuations on mean dispersive shear is

(3.15)\begin{equation} P_r ={-} \left\langle \left( \widetilde{u_k^{\prime} u_j^{\prime}} \right)_r \frac{\partial \tilde{u}_{i,r}}{\partial x_j} \right\rangle - \left\langle \left( \widetilde{u_i^{\prime} u_j^{\prime}} \right)_r \frac{\partial \tilde{u}_{k,r}}{\partial x_j} \right\rangle, \end{equation}

the turbulent transport is

(3.16)\begin{equation} T_k ={-} \frac{\partial}{\partial x_3} \langle \widetilde{u_i^{\prime} u_k^{\prime} u_3^{\prime}} \rangle, \end{equation}

the dispersive transport is

(3.17)\begin{equation} T_r ={-} \frac{\partial}{\partial x_3} \langle \widetilde{u_i^{\prime} u_k^{\prime} \tilde{u}_{r,3}} \rangle, \end{equation}

the pressure transport is

(3.18)\begin{equation} T_p ={-} \frac{1}{\rho_0} \frac{\partial}{\partial x_3} \left[ \langle \widetilde{u_k^{\prime} p^{\prime}} \rangle + \langle \widetilde{ u_i^{\prime} p^{\prime}} \rangle \right], \end{equation}

the viscous-diffusion term is

(3.19)\begin{equation} T_v = \nu \frac{\partial^{2}}{\partial x_3^{2}} \langle \widetilde{u_i^{\prime} u_k^{\prime} } \rangle, \end{equation}

the pressure-strain rate correlation is

(3.20)\begin{equation} S_k = \frac{1}{\rho_0} \left[ \left\langle \widetilde{ p^{\prime} \frac{\partial u_i^{\prime}}{\partial x_k}} \right\rangle + \left\langle \widetilde{ p^{\prime} \frac{\partial u_k^{\prime}}{\partial x_i}} \right\rangle \right], \end{equation}

the dissipation is

(3.21)\begin{equation} \epsilon = 2\nu \left\langle \widetilde{\frac{\partial u_i^{\prime}}{\partial x_j} \frac{\partial u_k^{\prime}}{\partial x_j}} \right\rangle, \end{equation}

and $\theta$ is the wave phase. The dispersive terms are similar to those derived in the energetic budgets by Raupach et al. (Reference Raupach, Antonia and Rajagopalan1991), Mignot et al. (Reference Mignot, Barthelemy and Hurther2009) and Yuan & Piomelli (Reference Yuan and Piomelli2014). Note that eliminating the dispersive transport and dispersive production terms recovers the flat wall TKE and Reynolds stress budget.

3.5.1. Turbulent kinetic energy budget

As discussed in § 3.5, the time- and volume-averaged TKE equation can be obtained from (3.12) by setting $i = k$ to give a balance between shear production $\langle \bar {P}_s \rangle _v$ and dissipation $\langle \bar {\epsilon } \rangle _v$. For the flat wall, wave-current cases, Scotti & Piomelli (Reference Scotti and Piomelli2001) found that the evolution of TKE is primarily governed by the changes in $\langle \tilde {P}_s \rangle$ and $\langle \tilde {\epsilon } \rangle$. Thus, comparison of $\langle \tilde {P}_s \rangle$ and $\langle \tilde {\epsilon } \rangle$ across the four cases can help explain the $\langle \tilde {M}_r \rangle$ evolution discussed in the previous section. Figure 17(a) shows the time- and planform-averaged shear production over dissipation profiles (i.e. $\langle \bar {P}_s \rangle / \langle \bar {\epsilon } \rangle$) for the four cases under consideration. Case C350F has peak production at $x_3^{+} \approx 11.6$, as can be shown using analytic form of the shear production and dissipation for steady channel flows (Pope Reference Pope2000). Within the log-law region (i.e. $80 \lessapprox x_3^{+} \lessapprox 200$) where shear production is balanced by dissipation, this ratio has a value of unity, while further away from the wall, the ratio decreases due to the presence of a wake region near the surface (i.e. at $x_3^{+} = Re_*$). The flat wall, wave-current case WC350F closely follows case C350F except in the region where peak production occurs. The $\langle \bar {P}_s \rangle /\langle \bar {\epsilon } \rangle$ value shows slight attenuation in the peak production region when compared with case C350F. As seen in figure 17(b), this attenuation for case WC350F is a result of a minor enhancement of the dissipation when compared with case C350F. Further away from the wall, case WC350F follows case C350F identically. The bumpy wall cases C350B and WC350B show substantial differences in the shear production and dissipation profiles. Cases C350B and WC350B show significantly attenuated $\langle \bar {P}_s \rangle /\langle \bar {\epsilon } \rangle$ values below the location of peak production (i.e. $x_3^{+} \approx 11.8$). Far from the wall, both the cases follow their flat wall counterparts. These differences in the $\langle \bar {P}_s \rangle /\langle \bar {\epsilon } \rangle$ profiles are mainly attributed to the enhanced time-averaged dissipation for both the bumpy wall cases. As seen in figure 17(b), the location of the peak in magnitude of $\langle \bar {\epsilon } \rangle$ for the bumpy wall cases is just below the top of the roughness elements (i.e. $k_c^{+}$). The addition of waves in case WC350B accentuates the effect of the bumps, as indicated by the further decrease in the $\langle \bar {P}_s \rangle /\langle \bar {\epsilon } \rangle$ when compared with case C350B. In both bumpy wall cases, the presence of roughness elements shifts the peak production further away from the wall in addition to enhancing dissipation below the peak production region. The additional decrease in $\langle \bar {P}_s \rangle /\langle \bar {\epsilon } \rangle$ for case WC350B is due to a substantial increase in the time-averaged dissipation within the roughness region ($x_3^{+} \leq k_c^{+}$) in comparison to case C350B. These changes lead to a substantial reduction in the efficiency with which the available TKE is converted to Reynolds stress, as shown in § 3.4.

Figure 17. (a) Time- and planform-averaged shear production over dissipation comparison for the four cases. (b) Time- and planform-averaged profiles of shear production and dissipation normalised using $u_*^{4}/(\kappa \nu )$ for the four cases under consideration. The vertical dash-dotted line marks the top of the roughness elements ($k_c^{+}$) and the vertical dashed line marks $2l_t^{+} + k_c^{+}$.

Wave-phase variations of the ratio of shear production to dissipation across the three regions defined in (2.17)–(2.19) are shown in figure 18. The largest variations in $\langle \tilde {P}_s \rangle /\langle \tilde {\epsilon } \rangle$ as a function of wave phase are observed in region II followed by region I and region III, respectively. Starting with region III which is furthest from the wall, cases WC350B and WC350F are roughly the same order of magnitude with case WC350B showing slightly larger magnitude. As most of the production and dissipation is concentrated in regions I and II, the phase variations of $\langle \tilde {P}_s \rangle ^{(III)}/\langle \tilde {\epsilon } \rangle ^{(III)}$ are in phase with the wave velocity ($U_b \sin (\omega t$)) and peak at $\theta = {\rm \pi}/2$ showing nearly periodic behaviour as a function of the wave phase. The value of $\langle \tilde {P}_s \rangle ^{(III)}/\langle \tilde {\epsilon } \rangle ^{(III)}$ is less than one due to the presence of a wake region just above the log-law region close to the channel surface (at $x_3^{+} = Re_*$). Closer to the wall in regions I and II, $\langle \tilde {P}_s \rangle /\langle \tilde {\epsilon } \rangle$ is less sinusoidal, especially for the bumpy wall cases. The phase variations in regions I and II can be understood by comparing the phase variations of shear production and dissipation shown in figure 19(a,b,e,f). The phase variations in shear production are due to phase variations in the Reynolds stress ($\langle \widetilde {u_1^{\prime } u_3^{\prime }} \rangle$). Consequently, the elevated dissipation levels for case WC350B in regions I and II contribute to the suppression of $\langle \tilde {P}_s \rangle /\langle \tilde {\epsilon } \rangle$. Additionally, for case WC350F, the peak values in regions I and II are in phase with the wave velocity and occur at the same location. However, case WC350B leads case WC350F in region II by ${\rm \pi} /4$ while it lags case WC350F in region I, also by ${\rm \pi} /4$. One possible explanation for this behaviour is the phase variation of dissipation which is significantly different for the two wave-current cases. However, the exact mechanisms leading to this behaviour are beyond the scope of this paper.

Figure 18. Phase- and region-averaged production to dissipation ratio with panels (a), (b) and (c) corresponding to regions I, II and III, respectively. Note that each of the panels have different scales in the $y$-axis.

Figure 19. Phase variability of the dominant terms in the TKE budget. Panels (a–h) show a side-by-side comparison of cases WC350B and WC350F of the shear production (a,b), pressure transport (c,d), dissipation (e,f) and phase rate of change of TKE (g,h), respectively. For each of the pairs, panels (a,c,e,g) correspond to case WC350B while panels (b,d,f,h) correspond to case WC350F. All terms are normalised by $u_*^{4}/(\kappa \nu )$. The values indicated on the colourbar are clipped at the maximum and minimum values observed for case WC350F for comparison. The dash-dot line marks the top of the roughness elements ($k_c^{+}$) while the dashed line marks the turbulent Stokes length offset above the roughness elements (i.e. $2l_t^{+} + k_c^{+}$).

The phase variability of $\langle \tilde {P}_s \rangle / \langle \tilde {\epsilon } \rangle$ can be quantified by studying the phase variations of the dominant terms in the TKE budget, as shown in figure 19. As indicated in panels (a,b), case WC350B shows a smaller wave-cycle region over which the peak production occurs when compared with case WC350F. Additionally, shear production peaks slightly later in the wave cycle for case WC350B when compared with case WC350F. This behaviour confirms the slightly reduced shear production magnitudes observed in the time-averaged profiles shown in figure 17(b). As discussed earlier, these phase variations suggest that most of the changes observed in $\langle \tilde {P}_s \rangle / \langle \tilde {\epsilon } \rangle$ are dominated by the phase variations in dissipation. Figure 19(e,f) compares the dissipation for case WC350B and WC350F, clearly showing the large differences. During the acceleration portion of the wave cycle, dissipation occurs primarily just above the top of the roughness elements ($k_c^{+}$) for case WC350B. However, during the deceleration portion of the wave cycle, substantial enhancements in dissipation are observed for case WC350B below the top of the roughness elements when compared with case WC350F. In region II (i.e. the buffer layer) similar magnitudes of dissipation are observed for both cases, and far away from the wall (i.e. region III) they are identical. Panels (c,d) indicate significant differences in the pressure transport of TKE between cases WC350B and WC350F, implying that pressure transport plays a major role for case WC350B near the wall. Additionally, the changes observed in $\partial k / \partial \theta$ are confined to a near-wall region, as indicated in panels (g,h). These results suggest that most of the modifications in TKE occur in a region below the top of the roughness elements for case WC350B. Furthermore, the effects of waves are confined to a region below twice the turbulent Stokes length, while the region far from the wall behaves identically for cases WC350B and WC350F and appear to be dynamically decoupled from the near-wall region. These results support the two-layer model proposed by Townsend (Reference Townsend1976) for steady boundary layer flows, in which the region with the highest production (i.e. the buffer layer) is dynamically decoupled from the region far from the wall (i.e. the log-law region).

3.5.2. Vertical Reynolds stress budget

The vertical Reynolds stress budget can be derived by setting $i=1$ and $k=3$ in (3.12). Since the vertical Reynolds stress (hereafter Reynolds stress) contributes to the total shear stress for the channel geometry, studying the Reynolds stress as a function of wave phase and $x_3^{+}$ can provide a great deal of insight into the mean statistics observed in the previous sections. For hydraulically smooth, steady turbulent boundary layers, Yuan & Piomelli (Reference Yuan and Piomelli2014) showed that the relative magnitude of some of the Reynolds stress budget terms is small and, thus, similar to the TKE budget, not all terms have significant phase variability. Consequently, the phase variations of the Reynolds stress, turbulent transport, dispersive transport, viscous-diffusion and dissipation terms will not be presented in this section for the sake of brevity. We note that, for a higher roughness Reynolds number, it may not be justified to neglect some of the terms. After ignoring the terms listed above, the time- and planform-averaged Reynolds stress budget simplifies such that the pressure-strain rate correlations balance the production of Reynolds stress (Pope Reference Pope2000). As detailed by Lumley (Reference Lumley1975); Hao & Gorlé (Reference Hao and Gorlé2020), the role of the pressure-strain rate correlations is to disorganise the turbulent eddy structures resulting in decorrelation between the streamwise and the vertical turbulent flow components, i.e. they act as a sink term in the Reynolds stress budget. Another perspective on pressure-strain rate correlations is provided by Brasseur & Lee (Reference Brasseur and Lee1989), who interpret the role of the pressure-strain rate correlations as strong local inter-component transfer of energy associated with vortical structures. As seen in the r.m.s. velocity profiles in figure 11, the role of pressure-strain rate correlations is evident as the vertical and spanwise turbulent velocities are accentuated at the expense of the streamwise component for the bumpy wall cases C350B and WC350B.

As shown in figure 20, case WC350B shows enhanced pressure-strain rate $( \bar {S}_{1,3} )$ and production of Reynolds stress by mean shear $( \bar {P}_{1,3})$ when compared with case WC350F. The peak value of $\bar {S}_{1,3}$ and $\bar {P}_{1,3}$ occurs at the top of the roughness elements for case WC350B, while $\bar {S}_{1,3}$ and $\bar {P}_{1,3}$ both peak within the buffer layer region for case WC350F. Most of the changes observed occur in the vicinity of the wall for case WC350B when compared with case WC350F, while far from the wall the two wave-current cases have identical behaviour. The enhanced r.m.s. spanwise and vertical velocity profiles observed in figure 13(b) can be explained by the enhanced pressure-strain rate terms near the wall (not shown here) as observed by Huang, Wang & Fu (Reference Huang, Wang and Fu2021). Since $\langle \widetilde {u_3^{\prime } u_3^{\prime }}\rangle \partial \langle \tilde {u}_1 \rangle / \partial x_3$ is a source term in the evolution equation for $-\langle \widetilde {u_1^{\prime } u_3^{\prime }} \rangle$, the enhancements of turbulent flow components explains the elevated stress in the wave-current boundary layer.

Figure 20. Time- and planform-averaged comparison of production of Reynolds stress by mean shear $( \langle \bar {P}_{1,3} \rangle )$ and pressure-strain rate $( \langle \bar {S}_{1,3} \rangle )$ for cases WC350F and WC350B.

Figure 21(a,b) compares the production of Reynolds stress by mean shear and shows higher production during the deceleration portion of the wave cycle. As the TKE peaks during the deceleration portion of the wave cycle, the available $\langle \widetilde {u_3^{\prime } u_3^{\prime }} \rangle$ is relatively higher with $\partial \langle \bar {u}_1 \rangle / \partial x_3$ being constant. Thus, the net production by mean shear is highest during the deceleration portion of the wave cycle. Panels (e,f) compare the Reynolds stress production by wave shear for cases WC350B and WC350F. Opposing behaviour is observed for the wave shear based production, as the values during the acceleration portion of the wave cycle near the wall are negative, i.e. Reynolds stress production by wave shear is negative, while during the deceleration portion positive values are observed suggesting the opposite. It is clear to see that both production mechanisms are significantly enhanced for case WC350B when compared with case WC350F. Similar observations can be made for the pressure transport and pressure-strain rate terms. It is important to note that most of the differences observed are confined to a near-wall region for case WC350B.

Figure 21. Phase variability of the dominant terms in the Reynolds stress budget. Panels (a–h) show a side-by-side comparison of cases WC350B and WC350F of the shear production (a,b), pressure transport (c,d), production by wave shear (e,f) and pressure-strain rate (g,h), respectively. For each of the pairs, panels (a,c,e,g) correspond to case WC350B while panels (b,d,f,h) correspond to case WC350F. All terms are normalised by $u_*^{4}/(\kappa \nu )$. The values indicated on the colourbar are clipped at the maximum and minimum values observed for case WC350F for comparison. The dash-dot line marks the top of the roughness elements ($k_c^{+}$) while the dashed line marks the turbulent Stokes length offset above the roughness elements (i.e. $2l_t^{+} + k_c^{+}$).