1. Introduction
Estuarine bottom boundary layers are primarily driven by the combined action of tidal or wave-driven mean flows and oscillatory wave motions. Oscillatory wave motions provide the necessary bottom shear stress to erode the sediment bed, while the mean flows transport the sediment horizontally (Lacy & MacVean Reference Lacy and MacVean2016). Such bottom boundary layers often exhibit highly varying roughness characteristics both spatially and seasonally (Egan et al. Reference Egan, Cowherd, Fringer and Monismith2019), leading to non-trivial effects in the bottom boundary layer dynamics. Due to practical limitations, most experimental literature on wave-current boundary layers has focused on flow conditions where the strength of the oscillatory wave motions is greater than the mean flow and the wave orbital excursion is larger than the roughness height. Numerically, the challenges in adequately resolving roughness features along with the turbulent physics impose steep requirements on the grid resolution. In the presence of laminar waves, the computational grid only needs to resolve the turbulence generated by the mean flow component. However, as the oscillatory wave motion strengthens, turbulence generation occurs due to the mean flow shear and the instantaneous wave shear. This results in larger instantaneous flow Reynolds numbers, thus requiring finer grid resolution to resolve the turbulence. Accurate simulation of complex roughness features in a DNS framework imposes additional computational constrains. Consequently, most, if not all numerical studies investigating wave-current boundary layer flows have focused on understanding the dynamics for flat walls.
Typically, wave-current boundary layer flows are characterised by prescribing the strength of the mean flow, the strength of the oscillatory wave motion and the ratio of the oscillatory wave excursion to the bed roughness height. The corresponding non-dimensional parameters are (a) the mean flow friction Reynolds number 
$({\textit {Re}}_* = u_* H/\nu )$, where 
$u_*$ is the friction velocity, 
$H$ is the flow depth and 
$\nu$ is the kinematic viscosity of the fluid; (b) the wave Reynolds number 
$({\textit {Re}}_w = U_b^{2}/(\omega \nu ))$, where 
$U_b$ is the maximum wave orbital velocity and 
$\omega$ is the wave frequency; (c) the relative roughness 
$(A/\bar {k}_s)$, where 
$A = U_b/\omega$ is the wave orbital excursion length and 
$\bar {k}_s$ is the mean bed roughness height. As shown in figure 1(a), for a fixed relative roughness 
$A/\bar {k}_s$, purely oscillatory flow transitions from hydraulically smooth and laminar wave flow conditions to hydraulically rough and turbulent wave flow conditions with increasing 
${\textit {Re}}_w$. The flow may also transition to a hydraulically rough and turbulent wave flow state if the roughness height increases while 
${\textit {Re}}_w > 10^{4}$. For flat walls corresponding to the limit 
$A/\bar {k}_s \to \infty$, Lodahl, Sumer & Fredsøe (Reference Lodahl, Sumer and Fredsøe1998) studied the effect of oscillatory wave motion over a turbulent mean flow and found two distinct flow regimes depending on 
$Re_w$ and 
$Re_*$, as shown in figure 1(b). The first regime (green line in figure 1b) corresponds to a low enough mean flow friction Reynolds number such that, upon increasing the wave strength (increasing 
${\textit {Re}}_w$), the wave-current boundary layer becomes wave dominated (i.e. 
$U_c/U_b > 1$, where 
$U_c$ is the mean flow velocity) before the wave transitions to a turbulent state (i.e. before 
$Re_w = 1.5 \times 10^{5}$). This flow regime is characterised by a reduction in the bottom stress with the addition of waves between points I and II. Beyond point II, the wave becomes turbulent and the bottom stress increases monotonically with increasing wave strength. The second flow regime (red line in figure 1b) corresponds to a mean flow friction Reynolds number in which, upon increasing the wave strength, the wave-current boundary layer transitions to wave-dominated flow conditions after the critical value of 
${\textit {Re}}_w = 1.5 \times 10^{5}$. For this flow regime, the bottom stress remains constant and increases monotonically only after the flow becomes wave dominated (i.e. 
$U_b > U_c$). The findings of Lodahl et al. (Reference Lodahl, Sumer and Fredsøe1998) have since been validated numerically by Scotti & Piomelli (Reference Scotti and Piomelli2001), Manna, Vacca & Verzicco (Reference Manna, Vacca and Verzicco2012, Reference Manna, Vacca and Verzicco2015) and Nelson & Fringer (Reference Nelson and Fringer2018), to explain the underlying mechanisms leading to the non-monotonic bottom stress for the first regime.
(a) Flow regime classification for purely oscillatory wave motion over bumpy walls. Black solid lines demarcate the boundaries of different wave flow conditions for varying 
${\textit {Re}}_w$ and 
$A/\bar {k}_s$ (adapted from Lacy & MacVean Reference Lacy and MacVean2016). The red symbol marks case WC350B, one of the cases simulated in this paper as detailed in table 1. (b) Bottom stress (
$\tau$) for flat wall, wave-current boundary layer flows for two mean flow Reynolds numbers 
$Re_*^{1}$ and 
$Re_*^{2}$ (adapted from Lodahl et al. Reference Lodahl, Sumer and Fredsøe1998). Here 
$U_c$ is the mean flow velocity, 
$U_b$ is the wave orbital velocity and 
$\tau _c$ is the bottom stress without waves.
Wave-current boundary layer flows over bumpy walls (finite 
$A/\bar {k}_s$ in figure 1a) have been extensively investigated experimentally (Grant & Madsen Reference Grant and Madsen1979; Kemp & Simons Reference Kemp and Simons1982, Reference Kemp and Simons1983; Arnskov, Fredsøe & Sumer Reference Arnskov, Fredsøe and Sumer1993). The Grant & Madsen (Reference Grant and Madsen1979) wave-current model is the most widely accepted theory for wave-current boundary layer flows over rough walls. While the model holds for varying flow conditions, studies by Sleath (Reference Sleath1987) and recent in situ measurements in San Francisco Bay by Cowherd et al. (Reference Cowherd, Egan, Monismith and Fringer2021) have observed that the time invariance of the eddy-viscosity assumption does not hold for wave-current boundary layer flows. Cowherd et al. (Reference Cowherd, Egan, Monismith and Fringer2021) also found that the instantaneous boundary layer response assumed in Grant & Madsen (Reference Grant and Madsen1979) may not hold.
Although the flat wall, wave-current boundary layer drag reduction has been thoroughly investigated numerically (Scotti & Piomelli Reference Scotti and Piomelli2001; Manna et al. Reference Manna, Vacca and Verzicco2012, Reference Manna, Vacca and Verzicco2015; Nelson & Fringer Reference Nelson and Fringer2018), the energetics and mechanics of enhanced drag over bumpy walls have not been well studied for current-dominated flow conditions. The only exception is the work of Bhaganagar (Reference Bhaganagar2008) presenting first-order statistics for wave-current flows over egg-carton type roughness features in current-dominated flow conditions. In addition to the study by Bhaganagar (Reference Bhaganagar2008), numerical investigations of pulsative, i.e. wave-current boundary layer flows, have been carried out by Jelly et al. (Reference Jelly, Chin, Illingworth, Monty, Marusic and Ooi2020) over a cosine based roughness topography in the recent past. Jelly et al. (Reference Jelly, Chin, Illingworth, Monty, Marusic and Ooi2020) found that the contribution of pressure drag can be significant when compared with the skin friction drag during some portions of the wave cycle. They also found that outer-layer similarity proposed by Townsend (Reference Townsend1976) holds for such unsteady forcing conditions. Despite these numerical explorations and the wide range of aforementioned experimental studies, the mechanisms explaining enhanced drag over rough walls predicted by Grant & Madsen (Reference Grant and Madsen1979) are not well understood, particularly in current-dominated flow conditions. These deficiencies reinvigorate the need to understand wave-current boundary layers using numerical simulations.
Although weak wind waves interacting with a turbulent current over roughness elements are ubiquitous (Lacy & MacVean Reference Lacy and MacVean2016), numerical investigations of such systems are lacking in the literature. The present study aims to bridge this gap by using direct numerical simulation (DNS) to study the dynamics of a current-dominated, wave-current boundary layer over bumpy walls in hydraulically smooth (based on 
$A/\bar {k}_s$) flow conditions. In this paper a hydraulically smooth bed corresponds to one in which 
$\bar {k}_s/\delta _w \lesssim 4$, where 
$\delta _w$ is the wave boundary layer thickness (Lacy & MacVean Reference Lacy and MacVean2016). The numerical wave-current flume replicates a U-tube type experimental set-up that is commonly used to study wave-current boundary layers (Lodahl et al. Reference Lodahl, Sumer and Fredsøe1998; Yuan & Madsen Reference Yuan and Madsen2015).
2. Problem formulation
2.1. Governing equations and computational framework
We perform DNS of wave-current boundary layer flows over flat and bumpy walls in a channel flow configuration using the immersed boundary method (IBM) to simulate the bumps. The governing equations are given by
\begin{gather} \frac{\partial u_i}{\partial t} + \frac{\partial}{\partial x_j} (u_i u_j) ={-}\frac{1}{\rho_0} \frac{\partial p}{\partial x_i} + \nu \frac{\partial^{2} u_i}{\partial x_j \partial x_j} + U_b \omega \cos(\omega t) \delta_{i1} + \varPi_c \delta_{i1} + F_{{IBM}}, \end{gather}
\begin{gather} \frac{\partial u_i}{\partial x_i} = 0, \end{gather}
where 
$u_i$ is the velocity vector, 
$t$ is time, 
$x_j$ is the Cartesian coordinate vector, 
$\delta _{ij}$ is the Kronecker delta function, 
$\rho _0$ is the reference density of the fluid, 
$p$ is the pressure, 
$\nu$ is the kinematic viscosity, 
$U_b$ is the maximum wave orbital velocity, 
$\omega$ is the wave frequency, 
$\varPi _c$ is the constant pressure gradient driving the flow and 
$F_{{IBM}}$ is the immersed boundary force to represent the bumps (see below). Coordinate axes are aligned as 
$x_1$, 
$x_2$ and 
$x_3$ in the streamwise, spanwise and vertical directions, respectively. Periodic boundary conditions are applied in the streamwise and spanwise directions, while a no-slip boundary condition is applied at the bottom wall and a free-slip boundary condition is applied at the top wall to simulate open-channel like geometries. The boundary conditions at the top wall are given by
\begin{equation} u_3(x_3=H) = 0,\quad \frac{\partial u_i}{\partial x_3} (x_3=H) = 0,\quad \forall i \in {1,2}, \end{equation}
where 
$H$ is the channel height. Choosing 
$U_c$, 
$\bar {k}_s$ and 
$\nu$ as the repeating variables, the bottom stress has a functional dependence given by (in this paper, the bottom stress is assumed to have units of velocity squared)
\begin{equation} \frac{\tau}{U_c^{2}} = f \left( \frac{U_b}{U_c}, \frac{H}{\bar{k}_s}, \frac{U_c}{\omega \bar{k}_s}, \frac{U_c \bar{k}_s}{\nu} \right), \end{equation}
where the parameters on the right-hand side are the flow dominance parameter 
$( U_b/U_c )$, relative roughness (
$H/\bar {k}_s$), Keulegan–Carpenter number (
$KC = U_c/(\omega \bar {k}_s$)) and roughness Reynolds number (
$Re_k = U_c \bar {k}_s/ \nu$). The left-hand side of (2.4) is the drag coefficient. Note that the last two non-dimensional numbers can be represented using the friction velocity 
$u_*$ instead of 
$U_c$, where 
$u_*/(\omega \bar {k}_s)$ is the friction velocity based Keulegan–Carpenter number and 
$u_* \bar {k}_s / \nu$ is the roughness Reynolds number. It is important to note that this non-dimensional scaling is not unique for the governing equations described above. Choosing oscillatory motion based velocity (
$U_b$), length (
$A = U_b/\omega$) and time (
$1/\omega$) scales give
\begin{equation} \frac{\partial u_i^{*}}{\partial t^{*}} + \frac{\partial}{\partial x_j^{*}} (u_i^{*} u_j^{*}) ={-}\frac{\partial p^{*}}{\partial x_i^{*}} + \frac{1}{{\textit{Re}}_w} \frac{\partial^{2} u_i^{*}}{\partial x_j^{*} \partial x_j^{*}} + \left( \cos(t^{*}) + T^{*} D^{*} \right) \delta_{i1} + F^{*}_{{IBM}}, \end{equation}
where 
$({\cdot })^{*}$ denotes a non-dimensional quantity, 
${\textit {Re}}_w$ is the wave Reynolds number, 
$T^{*} = u_*/(H\omega )$ is the response time scale ratio of the wave component to the characteristic turbulent component and 
$D^{*} = u_*/U_b$ is the inverse of the flow dominance parameter as defined in (2.4) since 
$u_*$ and 
$U_c$ correspond to the steady component of the external forcing. By choosing a specific response time scale ratio such that ‘frozen’ turbulence exists, i.e. 
$\omega ^{+} > 0.04$ (Jelly et al. Reference Jelly, Chin, Illingworth, Monty, Marusic and Ooi2020), the flow system can be studied by only varying 
$U_b$.
The governing equations are solved with a staggered-grid, second-order accurate, finite-difference spatial discretisation (Orlandi Reference Orlandi2000; Moin & Verzicco Reference Moin and Verzicco2016). The fractional-step method with a third-order accurate Runge–Kutta time-advancing scheme is used to integrate the governing equations in time (Kim & Moin Reference Kim and Moin1985). The code has been validated in previous studies for turbulent channel flows (Lozano-Durán & Bae Reference Lozano-Durán and Bae2016, Reference Lozano-Durán and Bae2019). We implement a direct forcing IBM to include irregular bumps at the bottom wall as proposed by Scotti (Reference Scotti2006). It is important to note that computing surface integrals over the roughness elements is a non-trivial procedure due to the irregular nature of the bottom bathymetry. Consequently, computing forces over the roughness elements is not possible due to the direct forcing nature of the IBM. The principal utility of this approach is to model the effects of roughness without introducing additional control parameters. This computational approach to model the effect of roughness elements has been thoroughly validated by Yuan & Piomelli (Reference Yuan and Piomelli2014).
2.2. Computational grid and simulation parameters
The channel has dimensions 
$2{\rm \pi} H$, 
${\rm \pi} H$ and 
$H$ in the streamwise, spanwise and vertical directions, respectively. These channel dimensions are sufficiently large to correctly predict one-point statistics for 
${\textit {Re}}_{*} \leq 4200$ (Lozano-Durán & Jiménez Reference Lozano-Durán and Jiménez2014). The constant pressure gradient in (2.5) is prescribed as 
$\varPi _c = u_{*}^{2}/H$, where 
$u_* = 0.0035\,{\rm m}\,{\rm s}^{-1}$ and 
$H = 0.1$ m. The oscillatory pressure gradient in (2.5) is prescribed by fixing 
$U_b = 0.021\,{\rm m}\,{\rm s}^{-1}$ and 
$\omega = 2{\rm \pi} /T_w$, where 
$T_w = 5$ s is the wave period. The non-dimensional forcing frequency 
$\omega ^{+} = 0.1026 > 0.04$ indicates that the turbulence is expected to display asymptotic behaviour towards the ‘frozen’ state, as detailed in Jelly et al. (Reference Jelly, Chin, Illingworth, Monty, Marusic and Ooi2020). Based on these body forcing conditions, the flow is expected to be in the current-dominated regime 
$U_b/U_c \leq 1$. Additionally, the roughness conditions for the bumpy wall case correspond to hydraulically smooth flow conditions since the relative roughness is 
$\bar {k}_s/\delta _w = 1.34 < 4$, where 
$\delta _w$ is the wave boundary layer thickness, the friction Reynolds number is 
$Re_* = u_*H/\nu = 350$ and the wave Reynolds number is 
${\textit {Re}}_w = U_b^{2}/(\omega \nu ) = 351$. For the bumpy wall case, the channel is discretised with 
$512 \times 256 \times 128$ grid points, and for the flat wall case, the channel is discretised with 
$512\times 256 \times 92$ grid points. Uniform grid spacing is used over the bumpy surface (roughness elements) with a resolution of 
$\Delta x_3^{+} = 0.45$, where the plus unit indicates normalisation by wall units, i.e. 
$\Delta x_3^{+} = u_* x_3/ \nu$ . Above the roughness crest, the grid is stretched so that the maximum vertical grid spacing at 
$x_3 = H$ is 
$\Delta x_{3,max}^{+} = 6.6$. For the flat wall case, the first vertical grid cell has a height 
$\Delta x_3^{+} = 0.5$ and stretched until 
$\Delta x_{3}^{+} = 8.0$. Uniform grid spacing is used in the streamwise and spanwise directions with 
$\Delta x_1^{+} = \Delta x_2^{+} = 4.2$ for both the flat and bumpy wall cases. The mean roughness height for the bumpy wall cases is 
$\bar {k}_s^{+} = 6$, with the grid resolution comparable to Yuan & Piomelli (Reference Yuan and Piomelli2014). It is important to note that this grid resolution is sufficient to resolve the maximum friction velocity based on the superposition of the steady and oscillatory flows (Stokes Reference Stokes1851). Using the maximum instantaneous friction velocity, the resolution over the roughness elements is 
$\Delta x_3^{+} = 0.7$, 
$\Delta x_{3,max}^{+} = 11.63$ and 
$\Delta x_1^{+} = \Delta x_2^{+} = 7.27$, which is sufficient for resolving the requisite turbulent features. Note that the instantaneous friction velocity is twice as large as the mean friction velocity for the flat wall case. However, since the waves are laminar, they are not expected to generate associated turbulent structures for the flat wall case. Therefore, resolving the mean friction velocity is sufficient for the flat wall case. A time step size of 
$\Delta t^{+} = \Delta t/T_{\epsilon } \equiv 1.75 \times 10^{-4}$ (
$T_{\epsilon } = H/u_*$ is the eddy turn over time) is used for all cases based on ensuring a maximum Courant number of 0.4 for a total simulation time of 
$10^{3}$ s or 
$200$ wave periods for the wave-current cases. Simulations are run at the Texas advanced computing cluster on Stampede2 (Intel KNL) using 64 processors. On average, 6144 processor hours are required to simulate 
$10^{3}$ s of real time. The various flow simulations carried out are listed in table 1.
The various DNS cases carried out with C denoting the steady channel flow cases, WC denoting the wave-current case, 350 denotes the friction Reynolds number for the steady component, the letter F denotes the flat wall case, the letter B denotes the bumpy wall case and 
$N_w$ denotes the number of wave periods after an initial transient of 100 periods over which the statistics are gathered for the wave-current cases. The wave boundary layer thickness is defined as 
$\delta _w = \sqrt {2\nu /\omega }$. The wave-current cases are hydraulically smooth based on 
$\bar {k}_s/\delta _w < 4$ (Lacy & MacVean Reference Lacy and MacVean2016).
The rough wall at the bed is generated by placing an array of randomly oriented ellipsoids centred at 
$x_3 = -0.5 \bar {k}_s$, with their semi-axes lengths 
$k_{s,x_1} = \bar {k}_s$, 
$k_{s,x_2} = 1.4\bar {k}_s$ and 
$k_{s,x_3} = 2\bar {k}_s$, as originally proposed by Scotti (Reference Scotti2006). Using this algorithm, the value of 
$\bar {k}_s^{+}$ is known a priori as seen in figure 2(c). Shape characterisation can be achieved through the definition of the Corey shape factor (Corey Reference Corey1949)
\begin{equation} C_o = \frac{k_{s,x_1}}{\sqrt{k_{s,x_3} k_{s,x_2}}}. \end{equation}
The IBM algorithm results in a roughness function (or area fraction) 
$\psi _r(x_3)$ that is a function of the vertical coordinate axis, thus eliminating the need to include streamwise and spanwise separation length scales in the roughness function definition. Figure 2(a) shows a schematic of the channel with roughness elements at the bottom wall with the mean roughness height 
$\bar {k}_s^{+} = 32$, and panel (b) shows one such ellipsoidal roughness element represented on the computational grid. The blue region corresponds to the grid points within the solid, while the yellow shaded region corresponds to the grid points that are within the fluid. Panel (c) shows the roughness function 
$\psi _r(x_3)$ for cases C350B and WC350B. The roughness features employed in this paper have 
$C_o=0.6$. As discussed in § 3.3, the smaller value of 
$C_o$ when compared with sand-grain type roughness (spheres with 
$C_o=1$) leads to a larger drag coefficient due to the protruding nature of the ellipsoids into the boundary layer.
(a) Bumpy wall generated using an array of randomly oriented ellipsoidal elements for 
$\bar {k}_s^{+} = 32$. (b) Close-up view of the grid and roughness element for 
$\bar {k}_s^{+} = 32$. (c) Roughness function depicting the mean roughness height 
$\bar {k}_s^{+} = 6$ and the top of the roughness elements 
$k_c^{+}$ for the current-dominated case discussed in this study.
2.3. Flow velocity decomposition and averaging methods
We decompose the flow variables 
$f_i(x_1,x_2,x_3,t)$ into four components
\begin{equation} f_i(x_1,x_2,x_3,t) = \langle \bar{f}_i \rangle (x_3) + \tilde{f}_{r,i}(x_1,x_2,x_3,t) + \langle \tilde{f}_{w,i} \rangle (x_3,t) + f^{\prime}_i(x_1,x_2,x_3,t), \end{equation}
where the terms on the right-hand side are the double-averaged, dispersive, wave and turbulent flow components, respectively. This velocity decomposition is similar to Nikora et al. (Reference Nikora, McEwan, McLean, Coleman, Pokrajac and Walters2007) and Mignot, Barthelemy & Hurther (Reference Mignot, Barthelemy and Hurther2009) except for the wave component. In (2.7) the time average 
$\bar {{\cdot }}$, the phase average 
$\tilde {{\cdot }}$ and the planform average 
$\langle {\cdot } \rangle$ are given by
\begin{gather} \overline{f_i}(x_1,x_2,x_3) = \frac{1}{T_{avg}} \int_{t}^{t+T_{avg}} f_i(x_1,x_2,x_3,t^{\prime})\, {\rm d}t^{\prime}, \end{gather}
\begin{gather} \tilde{f}_i(x_1,x_2,x_3,t) = \frac{1}{N_w} \sum_{j=1}^{N_w} f_i(x_1,x_2,x_3,t+jT_w) - \overline{f_i}(x_1,x_2,x_3), \end{gather}
\begin{gather} \langle\, f_i \rangle (x_3,t) = \frac{1}{A_f(x_3)} \int_{A_f(x_3)} f(x_1,x_2,x_3,t) \,{\rm d}A , \end{gather}
where 
$T_{avg}$ is the time over which time averaging is carried out, 
$N_w$ is the number of waves over which the phase averaging is computed, and 
$A_f(x_3) = 2{\rm \pi} ^{2} H^{2} [ 1 - \psi _r(x_3) ]$ is the planform area occupied by the fluid which varies with height 
$x_3$ due to the bumps. Additionally, the vertical integral of the planform-averaged quantity gives the volume average
\begin{equation} \langle f_i \rangle_v (t) = \frac{1}{H} \int_0^{H} \langle f_i \rangle(x_3^{\prime},t) \,{\rm d}x_3^{\prime}, \end{equation}
and the cumulative mean at time 
$t_{cm}$ is given by
\begin{equation} \bar{f}_{i}^{cm} = \frac{1}{t_{cm}} \int_0^{t_{cm}} f_i(t)\,{\rm d}t. \end{equation}The planform-averaging equation (2.7) gives the vertical profile of the wave component
\begin{equation} \langle \tilde{u}_{w,i} \rangle (x_3,t) = \langle u_i \rangle(x_3,t) - \langle \bar{u}_i \rangle(x_3), \end{equation}and the phase-averaging equation (2.7) gives the dispersive component (Raupach, Antonia & Rajagopalan Reference Raupach, Antonia and Rajagopalan1991)
\begin{equation} \tilde{u}_{r,i}(x_1,x_2,x_3,t) = \tilde{u}_i(x_1,x_2,x_3,t) - \langle \bar{u}_i \rangle(x_3) - \langle \tilde{u}_{w,i} \rangle (x_3,t). \end{equation}Using the time-averaged velocity component along with the identities in (2.13) and (2.14), the turbulent velocity component can be isolated using the identity in (2.7) as
\begin{equation} u_i^{\prime} (x_1,x_2,x_3,t) = u_i(x_1,x_2,x_3,t) - \langle \bar{u}_i \rangle(x_3) - \langle \tilde{u}_{w,i} \rangle (x_3,t) - \tilde{u}_{r,i} (x_1,x_2,x_3,t). \end{equation}
Using the time- and phase-averaging definitions in (2.8) and (2.9), it can be shown that for any two flow variables 
$f$ and 
$g$, 
$\overline {f^{\prime }\tilde {g}} = \widetilde {f^{\prime } \tilde {g}} = 0$, while, 
$\widetilde {f^{\prime }} = \bar {\tilde {f}} =0$, 
$\widetilde {\tilde {f} \tilde {g}} = \tilde {f}\tilde {g} - \overline {\tilde {f}\tilde {g}}$ and 
$\overline {f^{\prime }} = 0$ (Hussain & Reynolds Reference Hussain and Reynolds1970). Velocity and pressure data are stored every 
$T_w/20$ for the wave-current cases to compute the statistics, and these data are used to study the phase variability.
The turbulent Stokes length (
$l_t^{+}$), defined as
\begin{equation} l_t^{+} = l_s^{+} \left[ \frac{\kappa l_s^{+}}{2} + \sqrt{1 + \left(\frac{\kappa l_s^{+}}{2}\right)^{2}} \,\right] = 9.8, \end{equation}
measures the height below which the waves affect the turbulence (Scotti & Piomelli Reference Scotti and Piomelli2001). In (2.16), 
$l_s^{+} = \sqrt {2u_{*}^{2}/\omega \nu } = u_*\delta _w/\nu = 4.4$ is the non-dimensional wave boundary layer thickness. Using the roughness crest (
$k_c^{+}$) and the turbulent Stokes length (
$l_t^{+}$), the channel depth 
$H$ can be divided into three distinct regions in the vertical. The vertically integrated flow quantity corresponding to each of these regions is given by
\begin{gather} f_i^{(I)} = \int_0^{k_c^{+}} f_i \,{\rm d}x_3^{+}, \end{gather}
\begin{gather} f_i^{(II)} = \int_{k_c^{+}}^{2l_t^{+}} f_i \,{\rm d}x_3^{+}, \end{gather}
\begin{gather} f_i^{(III)} = \int_{2l_t^{+}}^{Re_*} f_i \,{\rm d}x_3^{+}. \end{gather}3. Results and discussion
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).
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
\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}
\begin{gather} u_2(x_1,x_2,x_3) = \alpha_{TI} R(x_1,x_2,x_3), \end{gather}
\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,
\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.
Time evolution of 
$Re_*$ for case C350F.
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$.
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
\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.
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)
\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.
Time- and planform-averaged velocity profiles for the flat and bumpy wall cases.
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
\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
\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).
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).
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
\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
\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.
(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).
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.
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
\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$).
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.
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^{+}$).
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.
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
\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
\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
\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
\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
\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
\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
\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
\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
\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
\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.
(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.
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.
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.
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.
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^{+}$).
4. Conclusion
We studied wave-current boundary layer dynamics in current-dominated flow conditions over roughness elements for hydraulically smooth walls to understand the flow drag and energetics. It was found that, although there is negligible flow separation over the hydraulically smooth roughness elements, their presence leads to an increased drag coefficient for the wave-current boundary layer when compared with the steady boundary layer case with identical roughness elements. However, for identical wave forcing conditions, the flat wall case does not undergo any perceptible change in the mean velocity profile when compared with its steady boundary layer counterpart. Despite the three disparate time scale processes, i.e. mean flow, waves and turbulence, the three components are decoupled, suggesting that there is minimal interaction between the wave and the turbulent components for current-dominated flow regimes. However, there seems to exist a one-way coupling such that the waves modulate the turbulence, but not vice versa. Comparison of the drag coefficients suggests a drag enhancement of 3 
$\%$–11 
$\%$ (depending on the drag coefficient formulation used, i.e. (3.7) or (3.8)) when adding waves to the steady case with bumpy walls.
Comparison of the total stress profiles suggests that above the top of the roughness elements, the total stress profile is linear, as expected for the steady turbulent channel flows. Comparison of the r.m.s. velocity profiles suggests that the presence of roughness elements significantly alters the streamwise near-wall component, but shows slight enhancement away from the wall. As for the other two normal stresses, opposite trends are observed when compared with the steady and wave-current turbulent boundary layer flows over flat walls and steady turbulent boundary layer flows over bumpy walls. These observations suggest that the role of pressure transport and pressure-strain rate for the bumpy wall cases can be significant even for weak wave conditions.
The phase variations in the TKE and Reynolds stress provide a diagnostic process based explanation of the differences observed in the drag coefficient. During the acceleration portions of the wave cycle, the available TKE is more readily converted to Reynolds stress for the bumpy wall case when compared with the flat wall case. During the deceleration portion of the wave cycle, dissipation intensive turbulent bursting processes lead to a reduction in the Reynolds stress. Most of the alterations observed in the bumpy wall case occur below the top of the roughness elements and between the top of the roughness elements and the turbulent Stokes length.
Phase-dependent TKE budget analysis supports the diagnostic mechanism of drag enhancement observed with the variations of TKE and Reynolds stress. Observing the ratio of production and dissipation of TKE suggests that there is strong attenuation in the net production for the wave-current bumpy wall case when compared with the flat wall cases. Away from the wall, no significant differences were observed between the four cases discussed in this investigation. The TKE budget phase variations for the bumpy wall cases shows strong variations for most of the terms, while the flat wall, wave-current case shows minimal changes. The Reynolds stress budget, on the other hand, reveals strong changes in the pressure transport and pressure-strain rate terms for the wave-current bumpy wall case. However, the flat wall cases show minimal to no changes in the Reynolds stress budget.
To summarise, our results suggest that the addition of a weak wave over a steady turbulent current for flat walls does not significantly alter the flow. However, the addition of a weak wave over turbulent currents with bumpy walls acts to enhance the drag coefficient. These results have implications for the development and application of wave-current models for estuarine flow conditions. While simple eddy-viscosity wave-current models imply similar consequences for transitional and turbulent wave flow conditions (Grant & Madsen Reference Grant and Madsen1979), the results in this study suggest that such a drag enhancement is also observed for hydraulically smooth, current-dominated flow conditions in the presence of roughness elements. As the drag coefficient is used to model estuarine sediment transport, large-scale variability in the prediction of drag coefficients can drastically change the outcomes of simpler wave-current turbulence models even in current-dominated flow conditions like those discussed in this study.
Acknowledgements
We thank A. Lozáno-Duran for providing the computational core of the finite-difference code used in this study. We also thank P. Yi and M. Lindhart for the insightful discussions during the development of this manuscript. We would also like to thank the three anonymous reviewers for their comments and recommendations that helped improve the manuscript.
Funding
This research was supported by NSF award OCE 1736668. Simulations were conducted with supercomputer resources under XSEDE project CTS190063.
Declaration of interests
The authors report no conflict of interest.
Open access
