2.1. Experimental facility

The experimental facility is a newly-built OWT, named the wave–current–sediment facility (WCS), in the Hydraulic Engineering Laboratory of the Department of Civil and Environmental Engineering at the National University of Singapore. The facility has a 10 m long, 50 cm deep and 40 cm wide horizontal test channel and a powerful piston system which can generate a variety of periodic oscillations corresponding to full-scale flow conditions along the entire test channel. These oscillatory flows in the WCS will be simply referred to as ‘waves’ hereafter. A Börger EL1550 Rotary Lobe pump is connected to introduce collinear currents of up to $60~\text{cm}~\text{s}^{-1}$ average velocity in the test channel. The current direction can be easily reversed by reversing the pump’s rotation. A two-dimensional PIV system, supplied by the TSI Corporation, is used to obtain velocity measurements. For most tests in this study, the vertical resolution is approximately $0.4~\text{mm}/\text{grid}$ , which is sufficiently fine to reveal details of the boundary layer flows. Each measurement is obtained over $N$ wave periods, and the velocities measured at the same vertical level at $M$ different longitudinal positions, $x_{m}$ , are effectively homogeneous, so the Reynolds-averaged velocity profiles over a wave period are obtained by performing both spatial and phase averaging:

(2.1) $$\begin{eqnarray}\langle \hat{{\it\xi}}\rangle (z,t)=\frac{1}{MN}\mathop{\sum }_{m=1}^{M}\mathop{\sum }_{n=1}^{N}{\it\xi}(x_{m},z,t+(n-1)T),\quad 0<t<T\end{eqnarray}$$

where ${\it\xi}$ is either the horizontal or vertical component of flow velocity $(u,w)$ , $t$ is the time, and $(x_{m},z)$ are the horizontal and vertical coordinates. In the following we will, unless otherwise indicated, for simplicity use ${\it\xi}(z,t)$ to denote the double-averaged quantities. The reader is referred to Yuan & Madsen (Reference Yuan and Madsen2014) (YM14 hereafter) for more details on the WCS and the PIV system.

2.2. Bottom conditions

In this study, three bottom conditions, a smooth bottom and two fixed rough bottoms, are included. The smooth bottom is formed by smooth aluminium plates. One fixed rough bottom is created by gluing $3\text{M}^{\text{TM}}$ 710 Safety-WalkTM Slip-Resistant Coarse tapes (physical roughness scale of approximately 1 mm) onto the aluminium plates. This bottom will hereafter be referred to as the ‘sandpaper’ bottom. The other rough bottom consists of a mono-layer of 12.5 mm-diameter ceramic marbles glued onto aluminium plates.

Figure 1. Free-stream velocities of nonlinear waves: (a) Stokes wave; (b) forward-leaning wave.

The theoretical bottom location $z=0$ and bottom roughness $k_{b}$ were determined by YM14 for each bottom condition by logarithmic profile fitting the near-bottom velocity measurements for several pure current and pure sinusoidal wave tests. For the smooth bottom, $z=0$ is directly determined from the PIV images with only 0.1 mm uncertainty, and the effective bottom roughness formula for steady turbulent flows in smooth pipes by Nikuradse (Reference Nikuradse1932), $k_{b}=3.3{\it\nu}/u_{\ast }$ , is found also to be applicable for smooth turbulent oscillatory flows if the period-averaged magnitude of shear velocity $u_{\ast }$ is used as the characteristic shear velocity. For the sandpaper bottom, $z=0$ is found to be $0.6~\text{mm}\pm 0.1~\text{mm}$ below the mean crest level of bottom roughness elements and the Nikuradse equivalent sand grain roughness $k_{N}$ is $3.7~\text{mm}\pm 0.1~\text{mm}$ . For the ceramic-marble bottom, $z=0$ is $4.0~\text{mm}\pm 0.4~\text{mm}$ (roughly $1/3$ of the ceramic marbles’ diameter) below the top of the marbles and $k_{N}$ is $20~\text{mm}\pm 3~\text{mm}$ . For further details on the determination of the bottom conditions summarized above, including the log-profile fitting analysis, the reader is referred to YM14.

2.3. Flow conditions

Three periodic wave shapes, sinusoidal, Stokes and forward-leaning waves, are included in this study. The latter two are the sum of two harmonics:

(2.2) $$\begin{eqnarray}u_{\infty }(t)=U_{\infty ,1}\cos ({\it\omega}t)+\frac{U_{\infty ,1}}{4}\cos (2{\it\omega}t+{\it\varphi}_{\infty ,2}).\end{eqnarray}$$

The second-harmonic phase ${\it\varphi}_{\infty ,2}$ is 0° for Stokes waves and 90° for forward-leaning waves. As shown in figure 1, they represent the skewed and forward-leaning nature of near-bottom flows under nonlinear waves, respectively. For each wave shape, a variety of wave amplitudes and two wave periods (6.25 and 12.5 s) are considered, as summarized in table 1. Here the amplitudes are controlled by the first-harmonic displacement amplitude of the piston $s_{1}$ , so the measured $U_{\infty ,1}$ may deviate slightly ( $1{-}5\,\%$ ) from the nominal values listed in the fourth column of table 1 which are the cross-section averaged $U_{\infty ,1}$ based on $s_{1}$ . The PIV sampling frequency is 5.12 Hz for all tests, so the number of samplings per period is 32 for the short-period (6.25 s) tests and 64 for the long-period (12.5 s) tests. Each test is sampled for 32 wave periods, which was determined to produce reliable phase averaging. Current generation is specified by the pump’s working frequencies $f$ . Higher working frequency gives stronger current. Two currents for which the pump’s frequencies are 13 and 40 Hz are considered in this study. For various wave conditions the ratio of current bottom shear stress ${\it\tau}_{cb}$ to maximum wave bottom shear stress ${\it\tau}_{wm}$ varies between 8 % and 60 %. This range should be sufficient to cover most wave–current flows encountered in the coastal environment. Since the nonlinear waves have direction dependence, currents in both positive and negative directions are considered. Here the positive direction is the direction of positive wave velocity, e.g. the direction of the maximum velocity of Stokes waves.

Table 1. Target wave conditions ( $s_{1}$ : first-harmonic displacement amplitude of the piston; $U_{\infty ,1}$ : approximate amplitude of the first-harmonic free-stream velocity; $T$ : wave period; $Re=U_{\infty ,1}A_{\infty ,1}/{\it\nu}$ : Reynolds number based on first-harmonic free-stream velocity).

A three-part scheme, ‘wave_current_bottom’, is used to identify tests: ‘wave’ is the wave identifier chosen from the test IDs listed in table 1; ‘current’ is the current identifier in the form of ‘C’ plus a number indicating $f$ , e.g. C13 denotes a current with a 13 Hz pump frequency. For currents in the negative direction, a letter ‘r’ (for ‘reverse’) is added to the identifier, e.g. ‘C13r’. Finally ‘bottom’ is the bottom identifier which is chosen from ‘sm’ (the smooth bottom), ‘sa’ (the ‘sandpaper’ bottom) and ‘ce’ (the ceramic-marble bottom). For example, the test ST400a_C40r_ce is a combination of the Stokes wave given by the fifth row in table 1 and a 40 Hz current in the negative direction over the ceramic-marble bottom.

2.4. Evaluation of current generation

YM14 showed that the WCS can generate a desired wave motion with excellent accuracy, i.e. the amplitude of the oscillatory velocity is produced within 1 % of the target value, so here only the evaluation of current generation is presented.

Figure 2. Spatial-average velocity of 40 Hz test measured at 90 mm above the sandpaper bottom: (a) time series; (b) frequency spectrum.

To check the pump’s ability to generate steady discharges when operating alone, five preliminary current-only tests for which the pump’s working frequencies are 13, 26, 40, $-$ 26 and $-$ 40 Hz are performed over the sandpaper bottom (minus signs indicates reversed rotation). The current velocity profiles within 10 cm from the bottom are measured using the PIV system with a 5 Hz sampling frequency. Figure 2(a) shows the spatial-averaged velocity of the 40 Hz test measured at 90 mm above the bottom. The time series shows no visually detectable trend of variation over the entire test duration. The random fluctuation gives a standard deviation of $2.0~\text{cm}~\text{s}^{-1}$ which is less than 5 % of the mean value ( $44.3~\text{cm}~\text{s}^{-1}$ ). It is due to the residual turbulence after spatial averaging and the pulsating nature of the rotary lobe pump’s current generation. The frequency spectrum shown in figure 2(b) suggests that this random fluctuation is spread over the entire resolvable frequency range (0–2.5 Hz). For typical wave–current flows with over $1~\text{m}~\text{s}^{-1}$ wave velocity amplitudes and $0.1{-}0.5~\text{Hz}$ wave frequency, such noise (generally less than $0.2~\text{cm}~\text{s}^{-1}$ in amplitude) will add less than 1 % error to the measured wave velocity. Thus, it should not be a concern for either current generation or wave velocity measurements.

The total pump discharge should be proportional to the pump’s rotation frequency. Here the current velocity measured at $250~\text{mm}/2.72=90~\text{mm}$ above the bottom is used as a rough estimate of the cross-section average velocity by assuming that (a) the cross-section average velocity is close to the depth-average velocity of the bottom boundary layer, (b) the bottom boundary layer thickness is half of the channel working depth (250 mm) and (c) the current velocity profile within the bottom boundary layer is logarithmic. The measurements are plotted against the pump rotation frequency $f$ in figure 3. The signs of reversed currents are changed for easy comparison. Clearly, the measurements fall nicely on a fitted straight line giving a slope of $1.14~\text{cm}~\text{s}^{-1}~\text{Hz}$ with a $\pm 8\,\%$ $95\,\%$ -confidence interval:

(2.3) $$\begin{eqnarray}\bar{u}_{c}\;(\text{cm}~\text{s}^{-1})=1.14\,f\;(\text{Hz}).\end{eqnarray}$$

Thus, the average current velocity nominal magnitude is accurately determined from the pump rotation frequency.

Figure 3. Cross-section average current velocities under different pump rotation frequencies, showing current in the positive and negative directions and a linear fit).

For current generation, the most important requirement is its stability, i.e. the pump must be able to produce a steady discharge against the time-varying pressure produced by the piston oscillatory motion. This can be demonstrated if (a) no sizeable harmonics other than the intended harmonics exist in the frequency spectrum of free-stream velocity and (b) the generation of intended harmonics is as excellent as when waves are generated alone. To evaluate the extent to which our current-generation system meets these requirements, two preliminary tests of combined wave–current flows and a baseline pure wave test over the sandpaper bottom were performed. In the three tests, the piston produces a SP400a wave, sinusoidal oscillation with approximately $160~\text{cm}~\text{s}^{-1}$ first-harmonic free-stream velocity amplitude and 6.25 s period, while the pump rotation frequencies are 13 and 40 Hz for the two wave–current tests.

To evaluate our current-generation system performance, we should consider the flow velocity within the free-stream region of the wave boundary layer, so the spatial-averaged but not phase-averaged measurements at 17 cm above the bottom are selected in the following analysis. Figure 4 shows the frequency spectrum of the measurements for test SP400a_C40_sa. For clarity, the frequency is normalized by the wave frequency (0.16 Hz). The intended first-harmonic velocity and the mean (current) velocity are obviously much larger than any other harmonic. The residual noise (mostly less than $0.4~\text{cm}~\text{s}^{-1}$ ) is comparable to the noise in figure 2(b). Some integer-number harmonics, indicated by the circles, e.g. the second and third harmonics, are much larger than the residual noise but much smaller than the intended first harmonic. These higher harmonics are most likely produced by boundary layer processes related to the temporal variation of the turbulent eddy viscosity as discussed by YM14. Therefore, no unwanted harmonics in the frequency spectrum are produced by current generation. The first three harmonics of the measured free-stream velocity are shown in table 2. The free-stream velocity can also be inferred from the measured piston displacements based on continuity, and the results (denoted by ‘inferred’) are listed together with PIV measurements in table 2 for comparison. The excellent agreement among the three inferred first-harmonic amplitudes ( $U_{\infty ,1}$ ) suggests that the piston oscillatory movement is not affected by current generation. The $U_{\infty ,1}$ of wave–current tests deviate from the inferred values by only $0.2{-}0.3~\text{cm}~\text{s}^{-1}$ (0.2 % of $U_{\infty ,1}$ ), which is similar to wave-alone tests reported by YM14. The PIV measured second and third harmonics ( $U_{\infty ,2}$ and $U_{\infty ,3}$ ) are affected by boundary layer processes and therefore not in agreement with the inferred values, but they are not intended to exist and are sufficiently small to be considered negligible. These results, together with the fact that no unwanted harmonics are observed, demonstrate excellent stability of the current-generation system.

Figure 4. Frequency spectrum of spatial-averaged free-stream velocity of test SP400a_C40_sa (frequency is normalized by the wave frequency 0.16 Hz).

Table 2. Comparison of measured free-stream velocity to the free-stream velocity inferred from piston displacement measurements ( $U_{\infty ,1}$ , $U_{\infty ,2}$ and $U_{\infty ,3}$ are the amplitudes of the first three harmonics).

3.1. Sinusoidal-wave–current boundary layers

3.1.1. Velocity profiles

The basic wave–current interaction illustrated by the Grant & Madsen (Reference Grant and Madsen1979) model (GM hereafter) is ideal for describing the sinusoidal-wave–current flows in OWTs. In the GM model, the following two-layer time-invariant turbulent eddy viscosity ${\it\nu}_{t}$ is proposed:

(3.1) $$\begin{eqnarray}{\it\nu}_{t}=\left\{\begin{array}{@{}l@{}}{\it\kappa}u_{\ast m}z,\quad z\leqslant {\it\delta}_{cw},\\ {\it\kappa}u_{\ast c}z,\quad z>{\it\delta}_{cw},\end{array}\right.\end{eqnarray}$$

where $u_{\ast c}$ is the current shear velocity $\sqrt{{\it\tau}_{cb}/{\it\rho}}$ , $u_{\ast m}$ is the maximum shear velocity $\sqrt{{\it\tau}_{bm}/{\it\rho}}$ , ${\it\delta}_{cw}$ is a transition level where wave-produced turbulence vanishes and ${\it\kappa}$ is the von Kármán constant which is found to be between $0.38{-}0.42$ in various studies (we simply use 0.40 hereafter). Here ${\it\tau}_{cb}$ and ${\it\tau}_{bm}$ are the current and maximum combined wave–current bottom shear stresses, respectively. The vertical variation of the current shear stress ${\it\tau}_{c}$ is neglected by assuming the region of interest is the very near-bottom part of the entire current boundary layer, so the linearized boundary layer equation is solved analytically to give a two-log-profile structure for the current velocity profile:

(3.2) $$\begin{eqnarray}\bar{u}=\left\{\begin{array}{@{}l@{}}\displaystyle \frac{u_{\ast c}^{2}}{{\it\kappa}u_{\ast m}}\ln \left(\frac{z}{z_{0}}\right)=\frac{{\it\alpha}u_{\ast c}}{{\it\kappa}}\ln \left(\frac{z}{z_{0}}\right)=\frac{u_{\ast c,1}}{{\it\kappa}}\ln \left(\frac{z}{z_{0}}\right),\quad z\leqslant {\it\delta}_{cw},\\ \displaystyle \frac{u_{\ast c}}{{\it\kappa}}\ln \left(\frac{z}{z_{0a}}\right)=\frac{u_{\ast c,2}}{{\it\kappa}}\ln \left(\frac{z}{z_{0a}}\right),\quad z>{\it\delta}_{cw},\end{array}\right.\end{eqnarray}$$

where ${\it\alpha}=u_{\ast c}/u_{\ast m}$ and $z_{0a}$ defines a new roughness scale $k_{a}=30z_{0a}$ , the apparent roughness. If waves are present, $z_{0a}$ is always larger than $z_{0}=k_{b}/30$ , sometimes by several orders of magnitude, which means that the current velocity profile, except for the very near-bottom region, experiences a dramatically increased apparent bottom roughness. The only influence of currents on waves is a slight increase of the turbulent eddy viscosity due to ${\it\tau}_{cb}$ , so waves effectively do not feel the existence of currents. The model prediction suggests that the amplitude of the first-harmonic velocity $U_{1}$ follows a logarithmic profile controlled by $z_{0}$ in the very near-bottom region:

(3.3) $$\begin{eqnarray}U_{1}=\frac{u_{\ast wm}^{2}}{{\it\kappa}u_{\ast m}}\ln \left(\frac{z}{z_{0}}\right)=\frac{u_{\ast w,1}}{{\it\kappa}}\ln \left(\frac{z}{z_{0}}\right),\end{eqnarray}$$

where $u_{\ast wm}$ is the maximum wave shear velocity $\sqrt{{\it\tau}_{wm}/{\it\rho}}$ with ${\it\tau}_{wm}$ being the maximum wave bottom shear stress. The predictive abilities of the original GM model are not always satisfactory due to some oversimplifications, e.g. the proposed discontinuous two-layer structure is only conceptually correct, so it was further modified several times (e.g. Madsen Reference Madsen1994). Humbyrd (Reference Humbyrd2012) provided the latest and most consistent modification based on a continuous three-layer structure for ${\it\nu}_{t}$ (see appendix A for details). Therefore, we shall use it for quantitative comparisons with our measurements and hence name it the ‘improved GM model’.

Figure 5. Velocity profiles of test SP400a_C40_ce: (a) current velocity profile; (b) first-harmonic amplitude profile (dots: measurements; dashed lines: fitted logarithmic profiles; crosses: limits for selecting data points for log-profile fitting).

Figure 5 shows the measured current velocity profile and the amplitude profile of the first-harmonic velocity for a typical test over the ceramic-marble bottom, SP400a_C40_ce. This test has a sinusoidal wave of approximately $160~\text{cm}~\text{s}^{-1}$ free-stream velocity amplitude with 6.25 s period and a cross-section average current velocity of approximately $46~\text{cm}~\text{s}^{-1}$ . The current velocity profile clearly exhibits the two-log-profile structure suggested by the GM model. The lower current profile has a relatively larger slope, which indicates that it is scaled by a reduced shear velocity, i.e. ${\it\alpha}=u_{\ast c}/u_{\ast m}<1$ in (3.2). The amplitude profile of the first-harmonic velocity can be well represented by a straight line in the very near-bottom region, which is in agreement with (3.3).

The following data selection rules are applied to select measurements for log-profile fittings. According to Grant & Madsen (Reference Grant and Madsen1979), the first-harmonic velocity amplitude profile is logarithmic in the very near-bottom region. Therefore, to fit the first-harmonic velocity amplitude profile and the lower current velocity profile which is supposed to be within the wave boundary layer, an upper data selection limit is imposed:

(3.4) $$\begin{eqnarray}z/l<0.15\end{eqnarray}$$

where the characteristic wave boundary layer length scale $l$ is

(3.5) $$\begin{eqnarray}l=\frac{{\it\kappa}u_{\ast m}}{{\it\omega}}.\end{eqnarray}$$

The choice of the 0.15 as the upper limit in (3.4) allows at least $7{-}10$ data points for tests with thin wave boundary layers. To account for the laminar sublayer and the buffer layer (Jiménez Reference Jiménez2004), a lower data selection limit is imposed for flows over a smooth bottom:

(3.6) $$\begin{eqnarray}z>100{\it\nu}/u_{\ast m}.\end{eqnarray}$$

It is found that PIV measurements are invalid for flows over rough bottoms in the region very close to the bottom roughness elements due to strong local laser reflection, so we also apply an alternative lower limit which requires that the percentage of good PIV measurements at a certain vertical level must exceed 75 %. This usually gives a lower limit between 0.5 and 1.5 mm above the top of bottom roughness elements. It should be noted that the spatial inhomogeneity of flow velocity due to individual bottom roughness elements is negligible at such a lower limit, i.e. the relative streamwise variation of the first-harmonic velocity amplitude is less than 5 %. The validity of these data selection limits is demonstrated by YM14. For fitting the upper current velocity profile, based on the improved GM model the lower limit for data selection is defined as

(3.7) $$\begin{eqnarray}z>1.5{\it\delta}_{ct},\end{eqnarray}$$

Table 3. Log-profile fitting of the current velocity profile and the first-harmonic amplitude profile of Test SP400a_C40_ce ( $1-R^{2}$ : coefficient of determination; $u_{\ast }$ : fitted shear velocity; $\pm {\rm\Delta}u_{\ast }/u_{\ast }$ : relative 95 % confidence interval of $u_{\ast }$ ; $k_{b}$ : fitted bottom roughness; $r_{{\rm\Delta}k}$ : 95 % confidence factor of $k_{b}$ ).

where ${\it\delta}_{ct}$ is the beginning of the upper logarithmic current profile according to the improved GM model. It is given by (A 2)–(A 7) using experimental values of $A_{bm}/k_{b}$ and ${\it\alpha}=u_{\ast c}/u_{\ast m}$ . Here $u_{\ast c}$ and $u_{\ast m}$ can be experimentally determined in advance by analysing the instantaneous bottom shear stress inferred from instantaneous velocity profiles (see § 3.1.2). The safety factor of 1.5 in (3.7) is applied to account for the uncertainty in ${\it\delta}_{ct}$ . Since the current in the WCS is driven by a depth-invariant mean pressure gradient, ${\it\tau}_{c}$ should decay approximately linearly from the bottom to the edge of the bottom current boundary layer $z={\it\delta}_{c}$ , resulting in the following upper current velocity profile (with (3.1) for ${\it\nu}_{t}$ ):

(3.8) $$\begin{eqnarray}\bar{u}=\frac{u_{\ast c}}{{\it\kappa}}\ln \left(\frac{z}{z_{0a}}\right)-\frac{u_{\ast c}z}{{\it\kappa}{\it\delta}_{c}}.\end{eqnarray}$$

For all tests in this study, the apparent roughness is in general between 1 and 20 cm, and current boundary layer thickness, ${\it\delta}_{c}$ , can be roughly estimated to be 250 mm, i.e. half of the test channel depth. With these representative values, we estimate that the second term in (3.8) is no more than $7{-}15\,\%$ of the first term below $z=100$  mm, indicating that the measured current velocity profile can be considered in agreement with the GM model’s prediction. Thus, $z=100$  mm is applied as the upper limit for data selection when analysing current velocity profiles in the outer region.

The results of log-profile fittings for the representative test are shown in table 3. The quality of log-profile fitting is quantified by the coefficient of determination $R^{2}$ ( $R^{2}=1$ indicating a perfect fit). The confidence level for $u_{\ast }$ is given by a normalized 95 % confidence interval ${\rm\Delta}u_{\ast }/u_{\ast }$ , and the confidence level for $k_{b}$ is indicated by a 95 % confidence factor $r_{{\rm\Delta}k}\geqslant 1$ , i.e. the true $k_{b}$ is $95\,\%$ likely to fall between $k_{b}/r_{{\rm\Delta}k}$ and $k_{b}\,r_{{\rm\Delta}k}$ . The lower current velocity profile has only 8 data points for log-profile fitting, so it gives the largest value of $1-R^{2}$ (of the order of $O(10^{-3})$ ), and consequently the largest ${\rm\Delta}u_{\ast }/u_{\ast }$ (5 %) and $r_{{\rm\Delta}k}$ (1.14). These confidence intervals are small enough to be considered acceptable. For collinear wave–current flows, the maximum combined bottom shear stress ${\it\tau}_{bm}$ is the sum of current bottom shear stress ${\it\tau}_{cb}$ and the maximum wave bottom shear stress ${\it\tau}_{wm}$ , so the following relationship between shear velocities is obtained:

(3.9) $$\begin{eqnarray}u_{\ast c}^{2}+u_{\ast wm}^{2}=u_{\ast m}^{2}.\end{eqnarray}$$

According to (3.2), (3.3) and (3.9), adding the shear velocity inferred from the lower current profile $u_{\ast c,1}=u_{\ast c}^{2}/u_{\ast m}=3.3~\text{cm}~\text{s}^{-1}$ and the shear velocity inferred from the first-harmonic velocity amplitude profile $u_{\ast w,1}=u_{\ast wm}^{2}/u_{\ast m}=17.9~\text{cm}~\text{s}^{-1}$ gives the maximum shear velocity $u_{\ast m}=3.3~\text{cm}~\text{s}^{-1}+17.9~\text{cm}~\text{s}^{-1}=21.2~\text{cm}~\text{s}^{-1}$ . Invoking $u_{\ast c,1}=u_{\ast c}^{2}/u_{\ast m}$ , we therefore get $u_{\ast c}=\sqrt{u_{\ast c,1}\,u_{\ast m}}=\sqrt{3.3\times 21.2}=8.4~\text{cm}~\text{s}^{-1}$ . This is identical to the shear velocity inferred from the upper current profile $u_{\ast c,2}=u_{\ast c}=8.4~\text{cm}~\text{s}^{-1}$ . As will be discussed in § 3.1.2, the experimental values of $u_{\ast c}$ and $u_{\ast m}$ can be directly determined from log-profile fitting instantaneous velocity profiles. For this test, the direct measurements of $u_{\ast c}$ and $u_{\ast m}$ are $8.8~\text{cm}~\text{s}^{-1}$ and $21.0~\text{cm}~\text{s}^{-1}$ , respectively, which are in excellent agreement with the fitted values (8.4 and $21.2~\text{cm}~\text{s}^{-1}$ ). Therefore, the internal relationship among various shear velocities suggested by the GM model is validated by this representative test. Other sinusoidal-wave–current tests yield similar quantitative comparisons.

The first-harmonic velocity amplitude profile and the lower current velocity profile shown in figure 5(a,b) give values of $k_{b}$ of 18.6 and 19.4 mm. These two values agree well with the predetermined Nikuradse equivalent sand grain roughness of the ceramic-marble bottom, $k_{N}=20$  mm, so both the lower current velocity and the wave velocity profiles are controlled by the physical bottom roughness. In most theoretical models, the no-slip boundary condition is applied at $z=z_{0}=k_{b}/30$ for both waves and currents. Our experimental results justify the validity of this choice. The upper current velocity profile in figure 5(a) gives an apparent roughness $k_{a}=30z_{0a}$ of 189.3 mm which is an order of magnitude larger than the physical bottom roughness. The improved GM model gives the following formula for $k_{a}$ :

(3.10) $$\begin{eqnarray}\frac{k_{a}}{k_{b}}=\frac{1}{{\it\alpha}}\left(\frac{5{\it\delta}_{w}}{ek_{b}}\right)^{1-{\it\alpha}},\end{eqnarray}$$

where ${\it\delta}_{w}$ is a characteristic wave boundary layer thickness. Using measured $k_{b}$ and ${\it\alpha}$ and applying (A 3)–(A 5) to give ${\it\delta}_{w}$ , (3.10) predicts $k_{a}=213$  mm, which is in good agreement with the fitted value.

These analyses suggest that the GM model (with the improvements provided by Humbyrd Reference Humbyrd2012) has very good ability in predicting collinear sinusoidal waves and currents in OWTs.

3.1.2. Bottom shear stress

The instantaneous bottom shear stress is obtained by log-profile fitting the instantaneous Reynolds-averaged velocity profiles. This method assumes that in the very-near bottom region oscillatory boundary layer flows are quasi-steady, so instantaneous Reynolds-averaged velocity profiles are logarithmic and controlled by the physical bottom roughness $k_{b}$ and the instantaneous shear velocities $u_{\ast }(t)$ . The validity of this method is supported by Jensen et al. (Reference Jensen, Sumer and Fredsøe1989) who directly measured the bottom shear stress for smooth oscillatory boundary layers using hot-film probes and found that their direct measurements were in good agreement with the estimates given by log-profile fittings. YM14 compared this method to other common methods for inferring bottom shear stress from velocity measurements, and demonstrated that this is the only valid method for rough-bottom tests in OWTs. To ensure that the selected data points correspond to near-bottom conditions, we simply use the bottom-most five data points which satisfy the previously stated data selection rules. A complete log-profile fitting analysis would take both $k_{b}$ and $u_{\ast }(t)$ as unknowns to be determined, but most theoretical models, including the one presented in § 4, assume $k_{b}$ as time-invariant. Therefore, here we adopt a modified log-profile fitting which has a predetermined $k_{b}$ given by preliminary wave-alone and current-alone tests (see § 2.2), so that it is consistent to compare the experimental results with model predictions. However, it should be noted that the difference between the modified and the complete log-profile analyses can be up to 10–15 % for the maxima of bottom shear stress, which indicates the uncertainty in the experimental results and should be considered when interpreting potential differences between predicted and measured net sediment transport rates in the WCS.

Figure 6. Instantaneous Reynolds-averaged velocity profiles of test SP400a_C40_sa and the associated modified log-profile fittings with $k_{b}=3.7$  mm (grey crosses (red online): measurements with questionable log-profile fittings; dots: measurements with good log-profile fittings; dashed lines: fitted logarithmic profiles): (a) first half-period; (b) second half-period.

Figure 6 shows the instantaneous velocity profiles and the modified log-profile fittings for every ${\rm\Delta}{\it\omega}t=22.5^{\circ }$ for a representative test, SP400a_C40_sa (the strongest sinusoidal waves plus the strongest current over the sandpaper bottom with $k_{b}=3.7$  mm). For most velocity profiles, the fitted logarithmic profiles reasonably represent the bottom-most measurements, and the coefficient of determination $R^{2}$ is ${>}0.99$ . For two short time windows when the free-stream velocity $u_{\infty }(t)$ crosses zero, i.e. the velocity profiles marked by crosses in figure 6, the fitted logarithmic profiles cannot reasonably approximate the measurements. As $u_{\infty }(t)$ decreases to zero, $\partial p/\partial x$ has the same sign as $u_{\infty }(t)$ and generally increases with ${\it\omega}t$ . Therefore, the instantaneous flow experiences an increasingly adverse pressure gradient, which eventually leads to flow separation, e.g. at ${\it\omega}t=87.8^{\circ }$ or $245.3^{\circ }$ . The same feature is also observed for experiments over the smooth and ceramic-marble bottoms, as well as many previous turbulent wave-alone boundary layers, e.g. Jensen et al. (Reference Jensen, Sumer and Fredsøe1989). Nevertheless, since $u_{\ast }(t)$ at these instants is close to zero, it is still acceptable to apply the modified log-profile fitting to the bottom-most five points and use the fitted $u_{\ast }(t)$ , which is also close to zero, as a rough estimate of the actual value. This will at least give the correct direction of the instantaneous bottom shear stress, e.g. the fitted profile at ${\it\omega}t=87.8^{\circ }$ in figure 6 has a negative slope. The duration when the log-profile fitting is not good is roughly 20–30 % of the total wave period, so this method works well for most of a wave period.

Figure 7. The bottom shear stress of test SP400a_C40_sa: modified log-profile fitting with bottom roughness $k_{b}=3.7$  mm, and smoothed time series for modified log-profile fitting.

Figure 7 shows the time series of bottom shear stress of test SP400a_C40_sa obtained from the modified log-profile fitting analysis, as well as a smoothed time series which only contains the mean and first three harmonics of the original time series. The smoothed time series nicely represents the data points, indicating that higher-order harmonics are negligibly small and have similar behaviour to experimental noise, but it does not represent the data points around zero-crossings very well, which is due to the effect of boundary layer separation. The current bottom shear stress ${\it\tau}_{cb}$ is obtained by period-averaging, and the maximum bottom shear stress ${\it\tau}_{bm}$ is obtained from the smoothed time series. The current bottom shear stress, ${\it\tau}_{cb}$ , is usually considered an input parameter for experiments and theoretical models, so we shall not discuss it here but use it for model validations in § 5. The difference between ${\it\tau}_{bm}$ and ${\it\tau}_{cb}$ is considered to be the maximum wave bottom shear stress, ${\it\tau}_{wm}={\it\tau}_{bm}-{\it\tau}_{cb}$ , which can be expressed by a wave friction factor $f_{wc}$ :

(3.11) $$\begin{eqnarray}f_{wc}=\frac{2{\it\tau}_{wm}}{{\it\rho}{U_{bm}}^{2}}.\end{eqnarray}$$

The improved GM model developed by Humbyrd (Reference Humbyrd2012) suggests that $f_{wc}$ depends on the relative bottom roughness $A_{bm}/k_{b}$ and the current condition represented by the parameter $C_{{\it\mu}}$ :

(3.12) $$\begin{eqnarray}C_{{\it\mu}}={\it\tau}_{bm}/{\it\tau}_{wm}=(1-{\it\alpha}^{2})^{-1}.\end{eqnarray}$$

By approximating the exact analytical solution, Humbyrd (Reference Humbyrd2012) obtained a simple explicit expression for $f_{wc}$ in the region $10<X_{{\it\mu}}<10^{5}$ :

(3.13) $$\begin{eqnarray}\frac{f_{wc}}{C_{{\it\mu}}}=\exp \left\{5.70(X_{{\it\mu}})^{-0.101}-7.46\right\},\end{eqnarray}$$

where

(3.14) $$\begin{eqnarray}X_{{\it\mu}}=C_{{\it\mu}}\frac{A_{bm}}{k_{b}}.\end{eqnarray}$$

The maximum wave bottom shear stress ${\it\tau}_{wm}$ leads the free-stream velocity by a phase lead ${\it\varphi}_{{\it\tau}wm}$ . For $10<X_{{\it\mu}}<10^{5}$ , the improved GM model suggests that ${\it\varphi}_{{\it\tau}wm}$ can be approximately expressed as

(3.15) $$\begin{eqnarray}{\it\varphi}_{{\it\tau}wm}\,(\text{deg.})=(0.649X_{{\it\mu}}^{-0.160}+0.118)\,180/{\rm\pi}.\end{eqnarray}$$

Figure 8. (a) Wave friction factors and (b) phase leads of the maximum wave bottom shear stress of combined sinusoidal-wave–current flows: smooth-bottom tests, sandpaper-bottom tests, ceramic-marble-bottom tests and the improved GM model).

With $C_{{\it\mu}}=1$ , (3.13) and (3.15) give the friction factor and phase lead for turbulent sinusoidal-wave boundary layers. We can plot the measured variations of $f_{wc}/C_{{\it\mu}}$ and ${\it\varphi}_{{\it\tau}wm}$ against $X_{{\it\mu}}$ and compare this with the predictions of (3.13) and (3.15), as shown in figure 8. The measurements of $f_{wc}/C_{{\it\mu}}$ have little scatter and form a consistent decreasing trend with increasing $X_{{\it\mu}}$ . The improved GM model reasonably captures the observed trend, although it slightly overestimates the friction factors by roughly 10 %. Such a small error is insignificant compared to the potential uncertainty in the determination of the bottom roughness, so the model’s performance can be considered excellent. This is not surprising, since the model’s time-invariant turbulent eddy viscosity is scaled by the maximum shear velocity $u_{\ast m}$ , which works ideally for predicting the maximum bottom shear stress. The measurements of ${\it\varphi}_{{\it\tau}wm}$ exhibit noticeable scatter, especially for the sandpaper-bottom tests (up to 5°). Nevertheless, the data points suggest that ${\it\varphi}_{{\it\tau}wm}$ decreases with increasing $X_{{\it\mu}}$ in a manner that is reasonably represented by the improved GM model, i.e. the data points appear to be distributed evenly on both sides of the curve. Given the fact that the amplitude is also well predicted, the performance of the improved GM model for prediction of the maximum wave bottom shear stress associated with combined collinear wave–current flows can be considered excellent.

3.2. The effect of wave nonlinearity on wave–current interaction

For pure nonlinear waves in OWTs, YM14 showed that a very weak mean velocity, the TA-streaming, is present. It is negative (opposing the wave direction) in the very near-bottom region, but becomes positive (following the wave direction) at higher elevations to balance the total volume flux. In this section, we briefly present some experimental results for nonlinear-wave–current flows in the WCS with focus on the current velocity profile, since the wave velocity is found to experience little influence from the superimposed current.

Figure 9. Magnitudes of measured current velocity profiles of the ST400a wave plus currents over the sandpaper bottom (dots: positive currents (in the wave direction); crosses: negative currents (against the wave direction): (a) strong (C40) currents; (b) weak (C13) currents.

Unlike sinusoidal waves, nonlinear waves have direction dependence, so we consider collinear currents with different directions. Figure 9 shows measured current velocity profiles of four typical tests over the sandpaper bottom. They have the same wave conditions, i.e. the strongest Stokes waves ST400a (first-harmonic velocity amplitude of approximately $160~\text{cm}~\text{s}^{-1}$ and 6.25 s wave period), but different current conditions. Figure 9(a) shows the current velocity profiles of the strong superimposed current C40 (roughly $46~\text{cm}~\text{s}^{-1}$ cross-section average velocity) following or opposing the wave direction, while figure 9(b) shows the same comparison for the weak superimposed current C13 (roughly $15~\text{cm}~\text{s}^{-1}$ cross-section average velocity). Here the magnitudes of measured current velocities are shown for easy comparison of currents in the two directions. The measurements suggest a pronounced difference between currents following or opposing nonlinear waves. In the very near-bottom region, the negative currents (shown by crosses) are stronger than the positive currents (shown by dots), but the discrepancy decreases with increasing elevation $z$ . For the two strong-current tests, the two current velocity profiles do not intersect within the range of observation, but there is a clear indication of crossing at a higher level. The GM-type model’s two-log-profile feature is observed for both of them, but is less pronounced for the negative current. For the weak-current tests, the two velocity profiles intersect at approximately $z=65~\text{mm}$ , and the positive current becomes increasingly stronger than the negative current above this level. The two-log-profile structure is completely obliterated for the negative current, while according to the GM model this feature should be very significant for weak currents plus strong waves, as demonstrated by the velocity profile of the positive C13 current.

The observed differences between currents following and opposing nonlinear waves agree qualitatively with the numerical results given by Holmedal et al. (Reference Holmedal, Johari and Myrhaug2013), and can be explained by TA-streaming produced by wave nonlinearity. The mean velocity profile embedded in the pure wave test ST400a_sa, which indicates the TA-streaming, is indicated by the solid line in figure 10. Such a mean velocity profile was first reported by Ribberink & Al-Salem (Reference Ribberink and Al-Salem1995). Thus, if added to the two-log-profile structure suggested by the GM model, in the very near-bottom region this mean velocity will enhance a negative superimposed current but suppress a positive superimposed current, while the situation reverses at higher elevations. As a confirmation, we subtract the magnitude of the negative current profile from the magnitude of the positive current profile shown in figure 9(a,b) and divide the difference by 2. The obtained results shown in figure 10 agree well with the measured mean velocity embedded in test ST400a_sa, i.e. the three profiles are all negative in the very near-bottom region and reach minima of comparable magnitudes ( $-3$ to $-4~\text{cm}~\text{s}^{-1}$ ) at approximately the same level $z\approx 10~\text{mm}$ , and the negative values decrease and eventually become positive at higher elevations.

Figure 10. Estimate of the embedded mean velocity of test ST400a_sa (dotted line: based on C13 currents; dashed line: based on C40 current; solid line: measured mean velocity of test ST400a_sa).

The GM model is entirely blind to the direction of currents, since it is based on sinusoidal waves, and it will always give a two-log-profile structure of the current velocity profile, which is completely absent for the test with a weak negative current, ST400a_C13r_sa. Therefore, the GM model is incapable of modelling nonlinear-wave–current flows, especially for a wave with strong nonlinearity plus a weak current. If one were to blindly apply the GM model to analyse field measurements obtained in relatively shallow waters where the waves are nonlinear, the obtained results, e.g. current shear velocity $u_{\ast c}$ , apparent roughness $k_{a}$ or physical bottom roughness $k_{b}$ , would be unreliable. In such situations, it is necessary to have a more elaborate wave–current boundary layer model which can consider the wave nonlinearity characterized by a second-harmonic free-stream velocity and a time-varying turbulent eddy viscosity.

4.1. Governing equations

Oscillatory flows in OWTs are uniform in the longitudinal direction of the test channel, so the governing momentum equation is

(4.1) $$\begin{eqnarray}\frac{\partial u}{\partial t}=-\frac{\partial (p/{\it\rho})}{\partial x}+\frac{\partial ({\it\tau}_{zx}/{\it\rho})}{\partial z},\end{eqnarray}$$

where $u$ is the Reynolds-averaged velocity in the horizontal direction, $x$ and $z$ are the horizontal and vertical coordinates, $t$ is time, ${\it\rho}$ is fluid density, $p$ is pressure and ${\it\tau}_{zx}$ is the component of the Reynolds stress which can be related to the Reynolds-averaged velocity through the concept of a turbulent eddy viscosity ${\it\nu}_{t}$ . The pressure gradient $\partial p/\partial x$ is considered depth-invariant, and its oscillatory part is related to the free-stream flow as

(4.2) $$\begin{eqnarray}-\frac{1}{{\it\rho}}\frac{\partial \tilde{p}}{\partial x}=\frac{\partial \tilde{u} _{\infty }}{\partial t},\end{eqnarray}$$

where $\tilde{u} _{\infty }$ is the oscillatory part of the free-stream velocity. The time-invariant part of the pressure gradient can be neglected if we only consider the very near-bottom part of the current boundary layer. For the currents in the WCS, we have discussed that the mean pressure gradient will not invalidate the logarithmic profile and hence can be generally neglected in the region less than 100 mm from the bottom. For mathematical convenience, we define a velocity deficit:

(4.3) $$\begin{eqnarray}u_{d}=u-\tilde{u} _{\infty },\end{eqnarray}$$

so (4.1) can be rewritten as

(4.4) $$\begin{eqnarray}\frac{\partial u_{d}}{\partial t}=\frac{\partial }{\partial z}\left({\it\nu}_{t}\frac{\partial u_{d}}{\partial z}\right).\end{eqnarray}$$

It is customary to split $u_{d}$ into a wave (time-dependent) velocity $\tilde{u} _{d}$ and a current (time-invariant) velocity $\bar{u}$ as

(4.5) $$\begin{eqnarray}u_{d}=\tilde{u} _{d}+\bar{u}\end{eqnarray}$$

where $\tilde{u} _{d}$ is expressed in the following general form:

(4.6) $$\begin{eqnarray}\tilde{u} _{d}(z,t)=\text{Re}\mathop{\sum }_{n=1}^{\infty }U_{d}^{(n)}(z)\text{e}^{\text{i}n{\it\omega}t}.\end{eqnarray}$$

The wave velocity must converge to the free-stream value asymptotically and satisfy the no-slip boundary condition which is prescribed at $z=z_{0}=k_{b}/30$ . Therefore, the following boundary conditions are prescribed for the complex amplitude of $n$ th-harmonic velocity deficit $U_{d}^{(n)}(z)$ :

(4.7) $$\begin{eqnarray}U_{d}^{(n)}=-U_{\infty }^{(n)},\quad z=z_{0},\end{eqnarray}$$

(4.8) $$\begin{eqnarray}U_{d}^{(n)}\rightarrow 0,\quad z\rightarrow \infty ,\end{eqnarray}$$

where $U_{\infty }^{(n)}$ is the complex amplitude of $n$ th-harmonic free-stream velocity. For $\bar{u}$ , the no-slip boundary condition also applies:

(4.9) $$\begin{eqnarray}\bar{u}=0,\quad z=z_{0}.\end{eqnarray}$$

The other boundary condition is specified as a reference current velocity $u_{r}$ at a reference elevation $z_{r}$ :

(4.10) $$\begin{eqnarray}\bar{u}=u_{r},\quad z=z_{r}.\end{eqnarray}$$

4.2. Turbulent eddy viscosity

For oscillatory turbulent boundary layers, Trowbridge & Madsen (Reference Trowbridge and Madsen1984a ,Reference Trowbridge and Madsen b ) suggested that the flow can be considered quasi-steady in the very near-bottom region, so the following expression for the turbulent eddy viscosity ${\it\nu}_{t}$ , obtained by generalizing the well-known results for steady turbulent boundary layer flows, is applicable in the immediate vicinity of the bottom:

(4.11) $$\begin{eqnarray}{\it\nu}_{t}={\it\kappa}\sqrt{\frac{|{\it\tau}_{b}(t)|}{{\it\rho}}}z={\it\kappa}|u_{\ast }(t)|z.\end{eqnarray}$$

Therefore, the temporal variation of ${\it\nu}_{t}$ is characterized by the temporal variation of $|u_{\ast }(t)|$ , which can be expressed as

(4.12) $$\begin{eqnarray}|u_{\ast }(t)|=\bar{u}_{\ast }\,f(t),\end{eqnarray}$$

where $\bar{u}_{\ast }$ is the period average of $|u_{\ast }(t)|$ and $f(t)$ is a dimensionless temporal variation function which is always positive and has an period-averaged value of 1. Thus, (4.11) can be rewritten as

(4.13) $$\begin{eqnarray}{\it\nu}_{t}=\bar{{\it\nu}}_{t}(z)\,f(t),\end{eqnarray}$$

where $\bar{{\it\nu}}_{t}(z)$ is the time-averaged turbulent eddy viscosity. Although at higher elevations the temporal variation of ${\it\nu}_{t}$ may not be synchronized with that in the very near-bottom region, which means that $f(t)$ also varies with $z$ , Trowbridge & Madsen (Reference Trowbridge and Madsen1984a ,Reference Trowbridge and Madsen b ) ignored this for mathematical convenience. Such a simplification will certainly deteriorate the prediction in the upper part of the boundary layer, but the flow in the very near-bottom region, which is critical for the determination of bottom shear stress, is still well-predicted. Therefore, in our theoretical model ${\it\nu}_{t}$ is treated as the product of a temporal variation function $f(t)$ and a mean turbulent eddy viscosity $\bar{{\it\nu}}_{t}(z)$ over the entire boundary layer.

We can always express $f(t)$ in terms of a Fourier series:

(4.14) $$\begin{eqnarray}f(t)=1+a_{1}\cos ({\it\omega}t+{\it\psi}_{1})+a_{2}\cos (2{\it\omega}t+{\it\psi}_{2})+\cdots \!.\end{eqnarray}$$

For sinusoidal waves, only the even harmonics of $f(t)$ exist since the two half-periods of $|u_{\ast }(t)|$ must be identical. For nonlinear waves, the two half-periods of $|u_{\ast }(t)|$ are asymmetric, e.g. the positive half-period of a Stokes wave will have larger maximum $|u_{\ast }(t)|$ than the negative half period, so odd harmonics of $f(t)$ are produced. The leading odd harmonic is the first harmonic $a_{1}\cos ({\it\omega}t+{\it\psi}_{1})$ , which interacts with the first-harmonic velocity gradient to produce a non-zero period-averaged shear stress. This is the origin of the TA-streaming (Trowbridge & Madsen Reference Trowbridge and Madsen1984b ). It should be noted that a superimposed current can also give rise to odd harmonics of $f(t)$ , so TA-streaming also exists for sinusoidal-wave–current boundary layer flows.

By only retaining the first and second harmonics of $f(t)$ , Trowbridge & Madsen (Reference Trowbridge and Madsen1984a ,Reference Trowbridge and Madsen b ) used perturbation methods to analytically study sinusoidal turbulent wave boundary layers (first-order analysis) and weakly nonlinear progressive waves (second-order analysis), but external currents were not considered in their study. Following their work, Gonzalez-Rodriguez & Madsen (Reference Gonzalez-Rodriguez and Madsen2011) studied nonlinear waves and wave–current flows in OWTs. Their predictions of TA-streaming for Stokes waves were in good agreement with measurements, but their predictions for forward-leaning waves and wave–current flows were not very successful. This is because they assumed that the first and second harmonics of $f(t)$ are very small, i.e. $a_{1},a_{2}\ll 1$ in (4.14), which is necessary for the neglect of higher-order terms. However, the predicted values of $a_{1}$ and $a_{2}$ contradict this assumption, e.g. the value of $a_{2}$ is around 0.4 which is hardly a very small number. Therefore, in § 4.3 we will develop a semi-analytical approach which does not make prior assumptions on the magnitudes of certain parameters.

The three-layer turbulent eddy viscosity of the improved GM model, i.e. (A 1), may be questionable for very weak currents, i.e. for very small ${\it\alpha}$ , the ${\it\delta}_{ct}$ given by (A 7) will be larger than ${\it\delta}_{w}$ , meaning that the wave-produced turbulence still controls the current even far outside the wave boundary layer, in violation of the basic assumption of the GM model. To remove this inconsistency, we propose the following four-layer vertical structure for $\bar{{\it\nu}}_{t}(z)$ :

(4.15) $$\begin{eqnarray}\bar{{\it\nu}}_{t}(z)=\left\{\begin{array}{@{}ll@{}}\displaystyle {\it\kappa}\bar{u}_{\ast }z, & z_{0}\leqslant z<{\it\delta}_{I},\\ \displaystyle {\it\kappa}\bar{u}_{\ast }{\it\delta}_{I}, & {\it\delta}_{I}\leqslant z<{\it\delta}_{J},\\ \displaystyle {\it\kappa}\bar{u}_{\ast }{\it\delta}_{I}\exp \left\{-{\it\gamma}\frac{z-{\it\delta}_{J}}{{\it\kappa}\bar{u}_{\ast }/{\it\omega}}\right\}\!, & {\it\delta}_{J}\leqslant z<{\it\delta}_{K},\\ \displaystyle {\it\kappa}u_{\ast c}z, & {\it\delta}_{K}\leqslant z.\end{array}\right.\end{eqnarray}$$

Figure 11. Vertical structure of the time-averaged turbulent eddy viscosity for oscillatory turbulent boundary layer flows: (a) weak current; (b) strong current.

A graphical representation of this vertical structure is shown in figure 11. Without the top-most current-dominated layer, the first three layers collectively give the $\bar{{\it\nu}}_{t}(z)$ of a pure wave boundary layer, which approximates a parabolic vertical structure scaled by the wave boundary layer thickness ${\it\delta}_{w}$ :

(4.16) $$\begin{eqnarray}\bar{{\it\nu}}_{t}(z)={\it\kappa}\bar{u}_{\ast }z\left(1-\frac{z}{{\it\delta}_{w}}\right),\quad z_{0}\leqslant z\leqslant {\it\delta}_{w}.\end{eqnarray}$$

Here ${\it\delta}_{w}$ is defined as the level where the amplitude of the first-harmonic velocity deficit reaches 1 % of the amplitude of the free-stream first-harmonic velocity. The choice of a parabolic structure as the target structure of the approximation is because this structure is commonly used for modelling steady open-channel turbulent flows. The transition levels ( ${\it\delta}_{I}$ and ${\it\delta}_{J}$ ) and the parameter ${\it\gamma}$ are determined by requiring that the total area under the parabolic curve is maintained, which gives

(4.17a−c ) $$\begin{eqnarray}{\it\delta}_{I}=0.21{\it\delta}_{w},\quad {\it\delta}_{J}=0.79{\it\delta}_{w},\quad {\it\gamma}=9.5\frac{{\it\kappa}\bar{u}_{\ast }/{\it\omega}}{{\it\delta}_{w}}.\end{eqnarray}$$

The third layer makes the wave-produced $\bar{{\it\nu}}_{t}(z)$ decay exponentially with elevation, so the singularity $\bar{{\it\nu}}_{t}({\it\delta}_{w})=0$ of the parabolic structure is avoided. The level ${\it\delta}_{K}$ is where ${\it\kappa}u_{\ast c}z$ intersects the underlying three-layer structure. With the exponential-decaying layer, ${\it\delta}_{K}$ will be generally less than ${\it\delta}_{w}$ for most current conditions, e.g. figure 11(a), so the possible inconsistency of the improved GM model is avoided. When the superimposed current is sufficiently strong, i.e. $u_{\ast c}/\bar{u}_{\ast }>{\it\delta}_{I}/{\it\delta}_{J}=0.27$ , the intersection will be within the second layer, as shown in figure 11(b), so the exponential-decaying layer does not appear and the three-layer structure of the improved GM model is reproduced.

4.3. A semi-analytical approach

In this subsection, a semi-analytical approach for solving the governing equation introduced in § 4.1 with the turbulent eddy viscosity introduced in § 4.2 is briefly presented. With (4.13) and (4.5), (4.4) can be written as

(4.18) $$\begin{eqnarray}\frac{\partial \tilde{u} _{d}}{\partial t}=\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\,f(t)\frac{\partial \tilde{u} _{d}}{\partial z}\right)+\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\,f(t)\frac{\partial \bar{u}}{\partial z}\right).\end{eqnarray}$$

Instead of using the perturbation method as Gonzalez-Rodriguez & Madsen (Reference Gonzalez-Rodriguez and Madsen2011), we divide all terms by $f(t)$ to get

(4.19) $$\begin{eqnarray}\frac{\partial \tilde{u} _{d}}{f(t)\partial t}=\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\frac{\partial \tilde{u} _{d}}{\partial z}\right)+\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\frac{\partial \bar{u}}{\partial z}\right).\end{eqnarray}$$

Following Lavelle & Mofjeld (Reference Lavelle and Mofjeld1983), we define a new time variable:

(4.20) $$\begin{eqnarray}{\it\tau}=\int f(t)\text{d}t=t+h(t),\end{eqnarray}$$

where, with $f(t)$ given by (4.14):

(4.21) $$\begin{eqnarray}h(t)=\mathop{\sum }_{n=1}^{\infty }\frac{a_{n}}{n{\it\omega}}\sin (n{\it\omega}t+{\it\psi}_{n}).\end{eqnarray}$$

We change the temporal variable of (4.19) to give

(4.22) $$\begin{eqnarray}\frac{\partial \tilde{u} _{d}}{\partial {\it\tau}}=\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\frac{\partial \tilde{u} _{d}}{\partial z}\right)+\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\frac{\partial \bar{u}}{\partial z}\right).\end{eqnarray}$$

This is the governing equation with a new temporal variable ${\it\tau}$ . For clarity, we denote the solution in ${\it\tau}$ -space by $V(z,{\it\tau})$ , and split (4.22) into a wave equation:

(4.23) $$\begin{eqnarray}\frac{\partial \tilde{V}}{\partial {\it\tau}}=\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\frac{\partial \tilde{V}}{\partial z}\right)\end{eqnarray}$$

and a current equation:

(4.24) $$\begin{eqnarray}\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\frac{\partial \bar{V}}{\partial z}\right)=0.\end{eqnarray}$$

The oscillatory velocity in ${\it\tau}$ -space can be expressed in terms of a Fourier series:

(4.25) $$\begin{eqnarray}\tilde{V}(z,{\it\tau})=\text{Re}\left\{\mathop{\sum }_{n=1}^{\infty }V^{(n)}(z)\text{e}^{\text{i}n{\it\omega}{\it\tau}}\!\right\},\end{eqnarray}$$

where $V^{(n)}$ is the complex amplitude of $n$ th-harmonic velocity in ${\it\tau}$ -space:

(4.26) $$\begin{eqnarray}V^{(n)}(z)=V_{n}(z)\exp (\text{i}{\it\varphi}_{vn}(z)),\end{eqnarray}$$

with $V_{n}(z)$ and ${\it\varphi}_{vn}(z)$ being the amplitude and phase of $V^{(n)}$ , respectively. The wave equation is solved harmonic by harmonic:

(4.27) $$\begin{eqnarray}\text{i}n{\it\omega}V^{(n)}=\frac{\partial }{\partial z}\left(\bar{{\it\nu}}_{t}(z)\frac{\partial V^{(n)}}{\partial z}\right),\end{eqnarray}$$

and, with $\bar{{\it\nu}}_{t}$ given by (4.15), both (4.24) and (4.27) are ordinary differential equations that can be solved analytically. The general analytical solutions obtained are converted back to $t$ -space using (4.20), and the conversion requires some simple numerical Fourier analyses of given time series. It should be highlighted that a mean velocity, $\bar{u}_{V}(z)$ , is produced when converting the wave solution from ${\it\tau}$ -space to $t$ -space, which essentially represents the origin of the TA-streaming. Matching boundary conditions in $t$ -space produces a set of linear equations for the unknown constants in the general solution. The reader is referred to Yuan (Reference Yuan2013) for further details on approximately solving the wave equation using change of temporal variables.

The current equation is analytically solved as introduced in appendix B. The solution suggests that the total current velocity profile is the sum of the TA-streaming $\bar{u}_{s}$ (related to $\bar{u}_{V}(z)$ ), i.e. (B 5), and another mean velocity $\bar{u}_{c}$ given by solving the current equation in ${\it\tau}$ -space, i.e. (B 8). The latter is named the ‘basic current’ hereafter, because it is shown that the current bottom shear stress is solely produced by this basic current, while waves in OWTs will not create a mean bottom shear stress, regardless of the wave shape (see appendix B for details).

Since the solution requires ${\it\nu}_{t}$ which contains parameters given by model predictions, an iterative solution procedure is developed as follows. With the free-stream oscillatory velocity, the bottom roughness and the current reference velocity at $z=100$  mm obtained from the measurements, we first apply the improved GM model (appendix A) to predict the current bottom shear stress ${\it\tau}_{cb}$ and the maximum wave bottom shear stress ${\it\tau}_{wm}$ , which together give a time series of the bottom shear stress as

(4.28) $$\begin{eqnarray}{\it\tau}_{b}(t)={\it\tau}_{wm}\cos ({\it\omega}t+{\it\varphi}_{{\it\tau}wm})+{\it\tau}_{cb}.\end{eqnarray}$$

We then analyse this time series to give initial values for $\{a_{n}\}$ , $\{{\it\psi}_{n}\}$ and $\bar{u}_{\ast }$ , which define the mean turbulent eddy viscosity $\bar{{\it\nu}}_{t}$ given by (4.15) and the temporal variation function $f(t)$ given by (4.14). The initial value of the boundary layer thickness ${\it\delta}_{w}$ is also given by the improved GM model. For simulating all tests in this study, this initiation never results in a divergence of the model results, and the solution converges quickly.

The governing equations are then solved as introduced above, so the remaining task is to obtain new estimates of the iterative parameters through a closure hypothesis. According to (4.12), we have

(4.29) $$\begin{eqnarray}\bar{u}_{\ast }f(t)=\sqrt{\frac{|{\it\tau}_{b}(t)|}{{\it\rho}}}=\sqrt{\left|{\it\kappa}\bar{u}_{\ast }z_{0}f(t)\frac{\partial u}{\partial z}\bigg|_{z=z_{0}}\right|}.\end{eqnarray}$$

Since $\bar{u}_{\ast }$ and $f(t)$ are always positive by definition, this equation can be rewritten as

(4.30) $$\begin{eqnarray}\bar{u}_{\ast }f(t)={\it\kappa}z_{0}\left|\left.\frac{\partial u}{\partial z}\right|_{z=z_{0}}\right|.\end{eqnarray}$$

The right-hand side can be easily calculated from the solution for $u$ , so the time series of $\bar{u}_{\ast }\,f(t)$ over a wave period is obtained. We numerically Fourier analyse it to give the mean and the Fourier components which define the values of $\bar{u}_{\ast }$ , $\{a_{n}\}$ and $\{{\it\psi}_{n}\}$ . In addition, the wave boundary layer thickness ${\it\delta}_{w}$ can also be determined using the solution of the first-harmonic velocity in $t$ -space. These new values will replace the previous estimates and the procedure is repeated until the following convergence criteria are met:

(4.31) $$\begin{eqnarray}\displaystyle \left|\frac{x_{new}-x_{old}}{x_{old}}\right|<1\,\%,\quad x\sim {\it\delta}_{w},u_{\ast c},\bar{u}_{\ast },a_{n}\text{e}^{\text{i}{\it\psi}_{n}}.\end{eqnarray}$$

Most of the iterative algorithms are expressed using matrices and vectors (see Yuan Reference Yuan2013 for details), so the solution procedure can be easily programmed in a manner that allows us to consider any number of terms of $\{a_{n}\}$ , $\{{\it\psi}_{n}\}$ and $\{U^{(n)}\}$ by simply changing the sizes of these matrices and vectors. This is the advantage of this theoretical model over the analytical model of Gonzalez-Rodriguez & Madsen (Reference Gonzalez-Rodriguez and Madsen2011), because we do not need to make prior assumptions on flow conditions to neglect higher-order terms, but simply keep adding terms until the predictions for the leading Fourier components of velocity and bottom shear stress converge. Thus, we essentially let the model tell us which higher-order terms can be neglected. For modelling the experimental results in this study, calculations show that it is only necessary to consider up to the fourth harmonic of ${\it\nu}_{t}$ and the fifth harmonic of velocity, since adding more harmonics will cause only minor changes to the quantities of primary interest, e.g. for Stokes waves over the range $10<A_{bm}/k_{b}<10^{5}$ the first three harmonics of bottom shear stress are changed by less than 0.5 %, 1.5 % and 5.0 % in amplitude, and less than 0.1°, 0.3° and 5.5° in phase, respectively, by adding higher-order harmonics.

5.1. Wave velocity and wave bottom shear stress

For simplicity, the validation of our theoretical model for pure wave boundary layer flows are presented based on three representative tests, SP400a_sa, ST400a_sa, and FL320_sa, corresponding to sinusoidal, Stokes and forward-leaning waves, respectively. These tests are all over the sandpaper bottom and have wave periods of 6.25 s. Their first-harmonic free-stream velocity amplitudes are the highest among tests with the same wave shape, i.e. the values of $U_{\infty ,1}$ are roughly $160~\text{cm}~\text{s}^{-1}$ , $160~\text{cm}~\text{s}^{-1}$ and $126~\text{cm}~\text{s}^{-1}$ for SP400a_sa, ST400a_sa, and FL320_sa, respectively. For sinusoidal wave boundary layers, the top-most current-dominated layer of the four-layer mean turbulent eddy viscosity is not considered since no current shear stress exists, and the model’s algorithm is accordingly simplified. Because the two half-wave-periods of a sinusoidal wave are symmetric, the even-order harmonics of the velocity vanish, and therefore the following comparisons are based on the first- and third-harmonic velocities, as shown in figure 12. The model successfully predicts the first-harmonic velocity for sinusoidal waves in the immediate vicinity of the bottom (figure 12 a,b), i.e. the logarithmic amplitude profile is predicted with a relative error less than 2 % below $z=6$  mm, and the disagreement for the first-harmonic phase is within 1°–5°. The model’s performance at higher elevations, however, is not as good as in the very near-bottom region, i.e. both the magnitudes and the locations of amplitude overshoots are not predicted very accurately. This is most likely due to oversimplification of the mean turbulent eddy viscosity $\bar{{\it\nu}}_{t}$ , i.e. (4.15), in the upper part of wave boundary layer. Nevertheless, it is the very near-bottom region that plays a crucial role in determining bottom shear stress and sediment transport, so our model’s prediction of the first-harmonic velocity of sinusoidal waves is still considered good. The third-harmonic velocity embedded in sinusoidal waves, as shown in figure 12(c,d), is a result of the temporal variation of turbulent eddy viscosity. The model to some extent predicts the overshoot structure of the third-harmonic velocity amplitude in the very near-bottom region, but the magnitudes are significantly underestimated. The third-harmonic phases are overpredicted in the very near-bottom region by 30°–50°, and the predictions in the upper part of the water column deviate even more significantly from the measurements. Therefore, it should be concluded that the model only qualitatively predicts the third-harmonic velocity. Nevertheless, the third-harmonic velocity is a very small component of the entire oscillatory velocity (3–5 %), and its contribution to bottom shear stress is negligible, so there is no necessity to predict it very accurately.

Figure 12. Prediction of the first- and third-harmonic velocities of sinusoidal wave test SP400a_sa (solid lines: predictions; dots: measurements): (a) first-harmonic amplitude, $|u^{(1)}(z)|$ ; (b) first-harmonic phase, $\text{Arg}(u^{(1)}(z))$ ; (c) third-harmonic amplitude, $|u^{(3)}(z)|$ ; (d) third-harmonic phase, $\text{Arg}(u^{(3)}(z))$ .

Figure 13. Prediction of the second-harmonic velocity of test ST400a_sa (Stokes waves) and FL320a_sa (forward-leaning waves) (solid lines: predictions; dots: measurements): (a) amplitude of test ST400a_sa, $|U^{(2)}(z)|$ ; (b) phase of test ST400a_sa, $\text{Arg}(U^{(2)}(z))$ ; (c) amplitude of test FL320_sa, $|U^{(2)}(z)|$ ; (d) phase of test FL320_sa, $\text{Arg}(U^{(2)}(z))$ .

For Stokes and forward-leaning waves, the predictions of the first- and third-harmonic velocities lead to the same conclusions as for sinusoidal waves. The prediction of the mean velocity will be specifically discussed in the § 5.2, so here we will only discuss the predictions for second-harmonic velocities, as shown in figure 13. Generally speaking, the model’s performance is very similar to that for the first-harmonic velocity of sinusoidal waves. For the Stokes-wave test (figure 13 a,b), the near-bottom logarithmic amplitude profile is excellently predicted, and second-harmonic phase is also reasonably predicted with an error less than $2^{\circ }$ , but in the upper part of the boundary layer the agreement between measurements and predictions deteriorates, i.e. the vertical elevation of the maximum amplitude overshoot is still overestimated. For the forward-leaning-wave test (figure 13 c,d), the prediction is not as good as for the Stokes-wave test. The prediction of the second-harmonic amplitude is no longer excellent in the very near-bottom region, i.e. the model overpredicts the amplitude by roughly 5 % below $z=3$  mm, and the prediction of the second-harmonic phase is also worse, i.e. the measurements are nearly depth-invariant in the region a few millimetres above the bottom, but the predictions do not pick up this feature. Given that the second-harmonic velocity is only approximately 25 % of the first-harmonic velocity in magnitude, such an inaccuracy is still acceptable.

Figure 14. Prediction of wave bottom shear stress (solid lines: predictions; dots: measurements from modified log-profile fitting analysis with $k_{b}=3.7$  mm): (a) test SP400a_sa (sinusoidal wave); (b) test ST400a_sa (Stokes wave); (c) test FL320_sa (forward-leaning wave).

Although the overall performance of the theoretical model for predicting wave velocity is not excellent, it still accurately predicts the dominant harmonics of wave velocity in the very near-bottom region, which is essential for predicting wave bottom shear stress. The predicted time series of bottom shear stress for the three representative tests are compared with measurements obtained from modified log-profile fittings in figure 14. For all tests, the model predicts the positive and negative maxima of the bottom shear stresses with an accuracy better than 10 %. The temporal variation is also well predicted, e.g. the model accurately captures the asymmetry between the two wave half-periods for the two nonlinear-wave tests. Thus, the model can accurately predict bottom shear stress in the WCS, which is critical for future studies on fluid–sediment interaction.

5.2. Mean velocity embedded in pure nonlinear-wave tests

For six representative nonlinear-wave tests, the predictions of $\bar{u}$ are compared with measurements in figure 15. For the two ceramic-marble-bottom tests (figure 15 a,b), the model accurately predicts the mean velocity in the very-near bottom region, i.e. the maximum negative velocity and the thickness of the negative-velocity region are both reasonably predicted. Above $z=100$  mm, the location where the reference value of $\bar{u}$ is specified as a boundary condition, the neglected effect of the mean pressure gradient becomes significant, and the prediction deviates from measurements above $z=100$  mm for some tests. The model prediction is also good for the two sandpaper tests (figure 15 c,d). For the two smooth-bottom tests (figure 15 e,f) the prediction is only qualitative, i.e. the overshoot is significantly overestimated (by a factor over 2). This is possibly because the viscous effect is not completely negligible as evidenced by the thickness of buffer layer estimated from (3.6) to be up to 2 mm. Therefore, it is concluded that the model can quantitatively predict the mean velocity embedded in nonlinear waves over rough bottoms, but its performance for smooth-bottom scenarios is only qualitatively acceptable.

Figure 15. Prediction of the boundary layer streaming of the nonlinear waves (black dots: measurements; grey zones: standard deviations given by spatial averaging; solid lines: prediction): (a) ST400a_ce; (b) FL320a_ce; (c) ST400a_sa; (d), FL320a_sa; (e) ST400a_sm; (f) FL320a_sm.

The total mean velocity $\bar{u}$ embedded in a nonlinear wave in OWTs is the sum of a TA-streaming, $\bar{u}_{s}$ , and a return current, $\bar{u}_{c}$ . The present model can separately predict $\bar{u}_{s}$ and $\bar{u}_{c}$ , so it is of interest to investigate them individually. For simplicity, the following discussion is based on the two ceramic-marble-bottom tests ST400a_ce and FL320a_ce. In figure 16, the predicted total mean velocities $\bar{u}$ shown in figure 15(a,b) are decomposed into $\bar{u}_{s}$ and $\bar{u}_{c}$ . For the nonlinear waves in this study, model predictions suggest that $\bar{u}_{s}$ can be approximately written as

(5.1) $$\begin{eqnarray}\bar{u}_{s}(z)\approx \bar{u}_{s,1}(z)+\bar{u}_{s,2}(z),\end{eqnarray}$$

where $\bar{u}_{s,n}$ can be approximated by (see Yuan Reference Yuan2013 for details)

(5.2) $$\begin{eqnarray}\bar{u}_{s,n}(z)\approx \text{Re}\left\{-{\textstyle \frac{1}{2}}a_{n}\text{e}^{-\text{i}{\it\psi}_{n}}u^{(n)}(z)\right\},\quad n=1,2.\end{eqnarray}$$

It can be readily seen that $\bar{u}_{s,1}$ and $\bar{u}_{s,2}$ represent the TA-streaming produced by the first- and second-harmonic interactions of the time-varying velocity gradient and the time-varying turbulent eddy viscosity, since $a_{n}$ represents the relative magnitude of the $n$ th-harmonic turbulent eddy viscosity and $u^{(n)}(z)$ is the complex amplitude of the $n$ th-harmonic velocity.

Figure 16. Decomposing of the mean velocity profile embedded in nonlinear waves (thick dashed lines: turbulence asymmetry streaming; thin dashed lines: main components of turbulence asymmetry streaming; thick dash-dot lines: basic (return) current): (a) Stokes wave (ST400a_ce); (b) forward-leaning wave (FL320a_ce).

For the Stokes wave shown in figure 16(a), $\bar{u}_{s,1}$ and $\bar{u}_{s,2}$ are both negative over the entire boundary layer; $\bar{u}_{s,1}$ reaches $-19.3~\text{cm}~\text{s}^{-1}$ at $z=100$  mm and is much stronger than $\bar{u}_{s,2}$ which is only $-2.3~\text{cm}~\text{s}^{-1}$ at the same level. Therefore, the total $\bar{u}_{s}$ is $-22.5~\text{cm}~\text{s}^{-1}$ at $z=100$  mm, which is much stronger than the total mean velocity ( $\bar{u}=1.7~\text{cm}~\text{s}^{-1}$ ) at $z=100$  mm (see figure 15 a). Thus, the superimposed return current $\bar{u}_{c}$ must be positive and have a magnitude comparable to $\bar{u}_{s}$ , i.e. $\bar{u}_{c}$ is $24.2~\text{cm}~\text{s}^{-1}$ at $z=100$  mm. The predicted current shear velocity $u_{\ast c}$ which controls $\bar{u}_{c}$ is $4.64~\text{cm}~\text{s}^{-1}$ , and is in excellent agreement with the measurement $u_{\ast c}=4.67~\text{cm}~\text{s}^{-1}$ obtained from log-profile fittings of the instantaneous velocity profiles (§ 3.1.2). Therefore, the observed mean velocity embedded in a Stokes wave is essentially a small imbalance between a sizeable TA-streaming and a sizeable return current. The reason for having such a strong $\bar{u}_{s,1}$ is that the first-harmonic turbulent eddy viscosity and the first-harmonic free-stream velocity are roughly in phase. For this test the predicted $a_{1}$ and ${\it\psi}_{1}$ are 0.28 and $-22^{\circ }$ , respectively, and the amplitude of $U_{\infty }^{(1)}$ is roughly $160~\text{cm}~\text{s}^{-1}$ , so the asymptotic value of $\bar{u}_{s,1}$ , according to (5.2), is approximately $-20~\text{cm}~\text{s}^{-1}$ , which is in agreement with the actual prediction ( $-19.3~\text{cm}~\text{s}^{-1}$ ).

The TA-streaming of forward-leaning waves is dramatically different from that of Stokes waves. As shown in figure 16(b), $\bar{u}_{s,1}$ is mostly positive over the entire boundary layer and its magnitude is slightly smaller than that of $\bar{u}_{s,2}$ which is always negative. Thus, the magnitude of $\bar{u}_{s}$ is quite small, and consequently the return current required to balance the total volume flux is weak, i.e. $\bar{u}_{c}$ is only $1.02~\text{cm}~\text{s}^{-1}$ at $z=100$  mm and the predicted $u_{\ast c}$ is only $0.95~\text{cm}~\text{s}^{-1}$ (the measured value is $1.36~\text{cm}~\text{s}^{-1}$ ). The reason for such a small and positive $\bar{u}_{s,1}$ for the forward-leaning wave is that the first-harmonic turbulent eddy viscosity is almost 90 $^{\circ }$ out of phase with the first-harmonic velocity. The predicted ${\it\psi}_{1}$ is $-96^{\circ }$ , while the predicted $a_{1}$ is still 0.28. Thus, the absolute value of the complex number given by the right-hand side of (5.2) is still large ( $18~\text{cm}~\text{s}^{-1}$ ), but its real part is only $2~\text{cm}~\text{s}^{-1}$ . Gonzalez-Rodriguez & Madsen (Reference Gonzalez-Rodriguez and Madsen2011) only considered the first-harmonic interaction $\bar{u}_{s,1}$ , which is the dominant one for Stokes waves but not for forward-leaning waves. Therefore, their model works for Stokes waves, but since $\bar{u}_{s,1}$ is dramatically different from the total $\bar{u}_{s}$ , their prediction fails for forward-leaning waves. This suggests that the TA-streaming is very sensitive to the actual flow condition.

The predicted and measured mean bottom shear stresses are compared in terms of the current shear velocity $u_{\ast c}$ in figure 17. For the Stokes-wave tests (figure 17 a), the predictions on average deviate from the measurements by a mere 0.5 %. For the forward-leaning-wave tests, the agreement is still reasonable but not as good as for the Stokes-wave tests, i.e. the fitted dashed line in figure 17(b) corresponds to a slope of 0.65, indicating that the model underestimates $u_{\ast c}$ by 35 %. This is partly because the forward-leaning-wave tests have much smaller $u_{\ast c}$ than the Stokes-wave tests, so the same absolute experimental and model inaccuracies will result in a larger relative error. At the same time, the model’s performance in predicting the mean velocity is generally worse for forward-leaning waves, as shown in figure 15, consequently the model’s inaccuracy in predicting the mean bottom shear stress for forward-leaning waves is expected to be larger.

Figure 17. Comparison of the predicted and measured current shear velocity for nonlinear waves (solid line: perfect agreement; dashed lines: least-square fit to the data): (a) Stokes waves; (b) forward-leaning waves.

The analytical solution for the mean velocity indicates that the mean bottom shear stress is solely related to the return current, so it will not be the same if the waves are produced in real coastal waters, and hence has little value for practical applications. However, this mean bottom shear stress is very important for interpreting measurements of net sediment transport under pure nonlinear waves in OWTs, such as the experimental results reported by Ribberink & Al-Salem (Reference Ribberink and Al-Salem1995), O’Donoghue & Wright (Reference O’Donoghue and Wright2004), van der A, O’Donoghue & Ribberink (Reference van der A, O’Donoghue and Ribberink2010) and Ruessink et al. (Reference Ruessink, Michallet, Abreu, Sancho, Van der A, Van der Werf and Silva2011). Generally speaking, the wave nonlinearities make the maximum onshore (positive) bottom shear stress ${\it\tau}_{bm+}$ stronger than the maximum offshore (negative) bottom shear stress ${\it\tau}_{bm-}$ , so a net onshore sediment transport, which is closely related to the difference between the two maxima, is produced. When studying this phenomenon using OWTs, the additional positive mean bottom shear stress produced by the return current will enhance this difference, and therefore increase the net onshore sediment transport rate. Based on actual measurements of bottom shear stress, the mean bottom shear stress of test ST400a_ce increases the ratio ${\it\tau}_{bm+}/{\it\tau}_{bm-}$ from 1.54 to 1.75, which means that the difference between the two maxima is increased by roughly 40 %. As a rule of thumb, the bedload sediment transport rate is proportional to the magnitude of bottom shear stress to the power of $3/2$ , so the net bedload sediment transport rate is roughly proportional to a period-averaged quantity $I_{BS}$ defined as

(5.3) $$\begin{eqnarray}I_{BS}=\overline{|{\it\tau}_{b}(t)|^{3/2}\frac{{\it\tau}_{b}(t)}{|{\it\tau}_{b}(t)|}},\end{eqnarray}$$

where the overbar indicates period average. For Stokes-wave tests in this study, $I_{BS}$ can be increased by roughly 50 % if the mean bottom shear stress is added to the total bottom shear stress. Therefore, it is possible that half of the measured net sediment transport rates under Stokes waves in OWTs is due to a facility-produced return current, which must be considered when deriving or validating empirical sediment transport formulae using OWT measurements.

5.3. Currents in the presence of sinusoidal waves

Comparisons of the measured and predicted current velocity profiles for six representative sinusoidal-wave–current tests are shown in figure 18. These tests are all with the SP400a wave (the strongest sinusoidal wave) but have different current and bottom conditions. The predictions given by the improved GM model (indicated by the dashed lines) are also provided. It should be noted that this model was derived with a three-layer turbulent eddy viscosity, but the predicted current velocity profile was further simplified by neglecting the smooth transition layer between the two logarithmic layers, so the curves have a sharp kink which we shall not consider when evaluating this model’s performance. Generally speaking, the present model very accurately predicts the upper current profiles, while the accuracy deteriorates slightly in the very near-bottom region. The test SP400a_C13_ce shown in figure 18(a) has the worst agreement between measurements and predictions. This is possibly because it suffers from the most significant effect of the mean pressure gradient, since it has the largest apparent roughness (see discussions in § 3.1.1). The predictions given by the improved GM model closely follow the predictions of the present model. The predicted and measured current shear velocities are compared in figure 19. The present model generally overestimates the current shear velocity by approximately $6.1\,\%\pm 2.7\,\%$ , while the improved GM model has a slightly better performance, i.e. the overestimate is $3.3\,\%\pm 6.1\,\%$ . Thus, we conclude that both models accurately predict the current shear velocity.

Figure 18. Prediction of current velocity profiles in the presence of sinusoidal waves (solid lines: the present model; dashed lines: the improved GM model; dots: measurements): (a) SP400a_C13_ce; (b) SP400a_C40_ce; (c) SP400a_c13_sa; (d) SP400a_c40_sa; (e) SP400a_c13_sm; (f) SP400a_c40_sm.

Figure 19. Comparison of the predicted and measured current shear velocity (solid line: perfect agreement; dashed lines: least-square fit to the data): (a) present model; (b) the improved GM model.

Despite their equally good performance, the two models differ substantially in their formulation of the effect of waves on currents. For the GM-type model, only the basic wave–current interaction is considered, i.e. the waves influence the coexisting currents by increasing the turbulent eddy viscosity in the very near-bottom region. However, for the present model, the total current velocity is the sum of a basic current velocity $\bar{u}_{c}$ and a TA-streaming $\bar{u}_{s}$ . In the very near-bottom region, $\bar{u}_{c}$ is controlled by a turbulent eddy viscosity which is scaled by the mean shear velocity $\bar{u}_{\ast }$ , which is generally larger than $u_{\ast c}$ , so the effect of an enhanced turbulent eddy viscosity still exist. However, it is weaker than that suggested by the GM model, since $\bar{u}_{\ast }$ is less than $u_{\ast m}$ , which is used in the GM model to scale the mean bottom turbulent eddy viscosity, by a factor of roughly $\sqrt{2/{\rm\pi}}=0.8$ (assuming a weak current and a strong sinusoidal wave). Therefore, it is the existence of the TA-streaming, $\bar{u}_{s}$ , which is entirely neglected in the improved GM model, that makes the GM model’s predictions agree with those obtained from the present model. To illustrate this, the total current velocity $\bar{u}$ predicted by the present model is decomposed into $\bar{u}_{c}$ and $\bar{u}_{s}$ . For simplicity, the following discussion is based on a representative test SP400a_C40_sa (figure 18 d) which has the best agreement between model predictions and measurements.

Figure 20. Decomposing of the total current velocity profile of test SP400a_C40_sa into a boundary layer streaming and a basic current velocity (solid line: total mean velocity; dashed line: boundary layer streaming; dash-dot line: basic current).

As shown in figure 20, the predicted $\bar{u}_{s}$ is opposing $\bar{u}_{c}$ , and its magnitude is roughly a quarter of $\bar{u}_{c}$ at the level $z=100$  mm, so it is a non-negligible component of the total current velocity profile. For sinusoidal-wave–current flows, $\bar{u}_{s}$ is mainly due to the interaction of the first-harmonic turbulent eddy viscosity and the first-harmonic velocity gradient. Thus, $\bar{u}_{s}$ can be roughly written as

(5.4) $$\begin{eqnarray}\bar{u}_{s}(z)\approx \text{Re}\left\{-{\textstyle \frac{1}{2}}a_{1}\text{e}^{-\text{i}{\it\psi}_{1}}u^{(1)}(z)\right\}.\end{eqnarray}$$

For a superimposed basic current $\bar{u}_{c}$ in the positive direction, model predictions suggest that ${\it\psi}_{1}$ is always close to zero. The phase of $u^{(1)}(z)$ , as suggested by measurements, changes by no more than 25 $^{\circ }$ over the entire boundary layer. Thus, $\bar{u}_{s}$ can be further approximated as

(5.5) $$\begin{eqnarray}\bar{u}_{s}(z)\approx -{\textstyle \frac{1}{2}}a_{1}|u^{(1)}(z)|.\end{eqnarray}$$

This suggests that $\bar{u}_{s}$ is in the negative direction for a positive $\bar{u}_{c}$ . For a negative $\bar{u}_{c}$ , it can be easily proved that ${\it\psi}_{1}$ is changed to roughly $-180^{\circ }$ , so $\bar{u}_{s}$ becomes positive. Therefore, the TA-streaming embedded in the sinusoidal-wave–current flows is always against the direction of the superimposed current. In the very-near bottom region, the amplitude of the first-harmonic velocity $|u^{(1)}(z)|$ follows a logarithmic profile approximately controlled by the maximum wave shear velocity $u_{\ast wm}$ , so the total mean velocity $\bar{u}$ can be written as

(5.6) $$\begin{eqnarray}\bar{u}(z)=\bar{u}_{c}+\bar{u}_{s}=\frac{1}{{\it\kappa}}\left(-\frac{1}{2}a_{1}u_{\ast wm}+\frac{{u_{\ast c}}^{2}}{\bar{u}_{\ast }}\right)\ln \left(\frac{z}{z_{0}}\right).\end{eqnarray}$$

This suggests that $\bar{u}(z)$ is still logarithmic in the very near-bottom region. The shear velocity of the basic current velocity profile, $u_{\ast c}^{2}/\bar{u}_{\ast }$ , is reduced by $a_{1}u_{\ast wm}/2$ , which corresponds to the effect of the TA-streaming. Thus, an ‘effective’ shear velocity is given by

(5.7) $$\begin{eqnarray}u_{\ast c,a}=-\frac{1}{2}a_{1}u_{\ast wm}+\frac{{u_{\ast c}}^{2}}{\bar{u}_{\ast }}.\end{eqnarray}$$

In the GM-type models, the effective shear velocity for the lower current velocity profile is $u_{\ast c,b}=u_{\ast c}^{2}/u_{\ast m}$ . For a simple yet realistic approximation, the bottom shear stress of the sinusoidal-wave–current flows can be taken as the sum of a sinusoidal wave bottom shear stress and a current bottom shear stress. With this approximation, we can evaluate the difference between $u_{\ast c,a}$ and $u_{\ast c,b}$ for various current conditions specified by $u_{\ast c}/u_{\ast wm}$ . The results suggest that the ratio $u_{\ast c,a}/u_{\ast c,b}$ is between 0.93 and 1.15 for $u_{\ast c}/u_{\ast wm}$ from 0 to 2 (very weak to very strong current). This explains why the two models give virtually identical predictions for the lower current velocity profile. Thus, by using $u_{\ast m}$ as the scaling shear velocity, the GM model implicitly includes the effect of the TA-streaming which can only be predicted with a time-varying turbulent eddy viscosity.

5.4. Currents in the presence of nonlinear waves

For currents in the presence of nonlinear waves, both the wave nonlinearity and the current can produce TA-streaming, so the total $\bar{u}_{s}$ is the sum of two components: a current-related TA-streaming $\bar{u}_{sc}$ , which is always against the current direction; and a wave-related TA-streaming $\bar{u}_{sw}$ , which always opposes the wave direction. Therefore, $\bar{u}_{sw}$ and $\bar{u}_{sc}$ are co-directional for currents following nonlinear waves, resulting in an enhanced $\bar{u}_{s}$ that opposes currents. However, $\bar{u}_{sw}$ and $\bar{u}_{sc}$ are against each other for currents opposing nonlinear waves, resulting in a reduced $\bar{u}_{s}$ in the direction of the stronger of $\bar{u}_{sw}$ and $\bar{u}_{sc}$ . For a given working frequency of the WCS pump, the total net volume discharge is maintained, so if added to a nonlinear wave, the basic current $\bar{u}_{c}$ of a positive current should be stronger (in terms of $|u_{\ast c}|$ ) than that of a negative current with the same pump frequency, so the difference in $\bar{u}_{s}$ is compensated. This explains why the vertical profiles of the total mean velocity $\bar{u}$ depend on current direction (e.g. see figure 9).

Figure 21. Predicted and measured current velocity profile of a C13 current following and opposing the ST400a wave over the sandpaper bottom: (a) predicted and measured current velocity profiles for ST400a_C13_sa, (b) decomposition of the predicted current velocity profile for test ST400a_C13_sa, (c) predicted and measured current velocity profiles for ST400a_C13r_sa, (d) decomposition of the predicted current velocity profile for test ST400a_C13r_sa.

Figure 22. Predicted and measured current velocity profile of a C40 current following and opposing the ST400a wave over the sandpaper bottom: (a) predicted and measured current velocity profiles for ST400a_C40_sa, (b) decomposition of the predicted current velocity profile for test ST400a_C40_sa, (c) predicted and measured current velocity profiles for ST400a_C40r_sa, (d) decomposition of the predicted current velocity profile for test ST400a_C40r_sa.

Figures 21 and 22 compare the model’s predictions with the measurements for the four typical tests shown in figure 9, i.e. C13 currents (roughly $15~\text{cm}~\text{s}^{-1}$ cross-section average velocity) and C40 currents (roughly $46~\text{cm}~\text{s}^{-1}$ cross-section average velocity) following and opposing the strongest Stokes waves ST400a (first-harmonic velocity amplitude of $160~\text{cm}~\text{s}^{-1}$ and 6.25 s wave period) over the sandpaper bottom. The model accurately predicts the mean velocity profiles for all four tests. For the two tests with the C40 current shown in figure 22, the predictions are in excellent agreement with the measurements. For test ST400a_C13_sa (figure 21 a), the model reasonably predicts the very steep slope in the very near-bottom region, but the transition to the upper current profile takes place above the observed level. For test ST400a_C13r_sa (figure 21 c), the current velocity profile does not have the conventional two-log-profile structure suggested by the GM model. The model successfully predicts this current velocity profile and offers a clear explanation for this phenomenon. As shown in figure 21(d), the two-log-profile structure is predicted for the basic current profile, but it is totally obliterated by the coexisting TA-streaming which is much stronger than the basic current.

The predicted current velocity profile for each test is also decomposed into $\bar{u}_{c}$ and $\bar{u}_{s}$ . The two components of $\bar{u}_{s}$ , $\bar{u}_{sc}$ and  $\bar{u}_{sw}$ , cannot be further isolated, but their mutual interaction can still be traced. Here we used the $\bar{u}_{s}$ for a pure Stokes wave ST400a_sa over the sandpaper bottom as a rough estimate of $\bar{u}_{sw}$ . Figure 16 shows that $\bar{u}_{sw}$ is negative and its magnitude can reach the order of $-20~\text{cm}~\text{s}^{-1}$ at $z=100$  mm. For a positive superimposed C13 current, $\bar{u}_{sw}$ enhances the co-directional $\bar{u}_{sc}$ , so a sizeable total $\bar{u}_{s}$ opposing the current is produced, i.e. $\bar{u}_{s}$ is approximately $-30~\text{cm}~\text{s}^{-1}$ at  $z=100$  mm, as shown in figure 21(b). However, a negative superimposed C13 current gives a positive $\bar{u}_{sc}$ opposing $\bar{u}_{sw}$ , so the total $\bar{u}_{s}$ is the difference between its two components. As shown in figure 21(d), $\bar{u}_{s}$ is only $-12~\text{cm}~\text{s}^{-1}$ at $z=100$  mm, and it follows the current, indicating that $\bar{u}_{sw}$ is stronger than $\bar{u}_{sc}$ . Since the pump maintains a constant and preset discharge, the facility has to develop an additional current to cope with the different net discharges due to TA-streaming, which is analogous to the facility-generated return current balancing the TA-streaming for pure nonlinear wave tests. Therefore, the basic currents in the two tests are significantly different. A strong basic current is predicted for the test with a positive C13 current, i.e. $\bar{u}_{c}$ at $z=100$  mm is over $50~\text{cm}~\text{s}^{-1}$ in figure 21(b), but a much weaker basic current is required for the test with a negative C13 current, i.e. $\bar{u}_{c}$ is only $-6~\text{cm}~\text{s}^{-1}$ in figure 21(d).

The two tests with C40 currents give the same conclusions. For test ST400a_C40r_sa (figure 22 d), the positive $\bar{u}_{sc}$ is slightly stronger than the negative $\bar{u}_{sw}$ , so the total $\bar{u}_{s}$ is still against the current, but its magnitude ( ${<}6~\text{cm}~\text{s}^{-1}$ at $z=100$  mm) is much smaller than that of test ST400a_C40_sa (figure 22 b) which is over $-50~\text{cm}~\text{s}^{-1}$ at $z=100$  mm. Accordingly, a strong basic current, i.e. $u_{c}=100~\text{cm}~\text{s}^{-1}$ at $z=100$  mm, is predicted for test ST400a_C40_sa, while the basic current is much weaker for test ST400a_C40r_sa, i.e. $u_{c}$ is only $-65~\text{cm}~\text{s}^{-1}$ at $z=100$  mm.

Since the current (mean) bottom shear stress is only related to the basic current, the significant difference in the basic current velocity implies that the currents following and opposing a nonlinear wave in the WCS are not comparable, despite the total mean discharge being the same. Another way to look at this implication is as follows. If we only use a target current velocity at a certain reference level or a total current discharge to specify a current, we cannot distinguish between boundary layer streaming and basic current, so such a specification of currents is inappropriate when the current bottom shear stress is the primary concern. In such situations, we suggest using the current bottom shear stress or the slope of the current velocity profile outside the wave boundary layer (see (B 8)) to specify currents, although this may make the experimental procedure more tedious.

We can validate the prediction of the basic current using the measured current bottom shear stress, since appendix B shows that the total current bottom shear stress is just the mean bottom shear stress given by the basic current. This can also indirectly validate the prediction of TA-streaming given that the total current velocity is accurately predicted. As shown in table 4, except for test ST400a_C13r_sa which has a very small $u_{\ast c}$ and consequently a large relative error, the predicted $u_{\ast c}$ is in excellent agreement with the measurements, i.e. the relative error for the remaining three tests is only approximately 5 % (overestimate), which is in agreement with the 6 % overestimate for the $u_{\ast c}$ of sinusoidal-wave–current tests (figure 19). Both measurements and predictions suggest significant differences between the current shear velocities for currents following and opposing nonlinear waves, i.e. the experimental $u_{\ast c}$ in test ST400a_C13_sa is more than four times larger than that of ST400a_C13r_sa in magnitude. Therefore, the present model not only reasonably predicts the total current velocity profile, but also accurately predicts its two components. If the improved GM model is applied to these tests, very similar values for $u_{\ast c}$ will be predicted for following and opposing currents with the same pump setting, since the model does not consider the TA-streaming associated with the wave nonlinearity. Using approximated current velocities at $z=100$  mm ( $20~\text{cm}~\text{s}^{-1}$ and $55~\text{cm}~\text{s}^{-1}$ for C13 and C40 currents, respectively) for specifying the current condition, the current shear velocities given by the improved GM model are $3.27~\text{cm}~\text{s}^{-1}$ and $6.26~\text{cm}~\text{s}^{-1}$ for tests with C13 and C40 currents, respectively. Therefore, $u_{\ast c}$ will be dramatically overpredicted for the opposing currents, e.g. a factor of 3 for test ST400a_C13r_sa, but dramatically underpredicted for the following currents, e.g. 31 % for test ST400a_C13_sa. This shows that a GM-type model is not applicable in shallow waters where waves can be highly nonlinear.

Table 4. Predicted and measured current shear velocities for nonlinear-wave–current tests.

These results also have important implications in the study of sediment transport under nonlinear-wave–current flows in OWTs. Since both the nonlinear wave and the current depend on direction, it is natural to compare the net sediment transports of currents following and opposing the nonlinear wave. If the total current discharge or a reference current velocity at a reference elevation is used to specify the current condition, the two basic currents obtained are not equivalent, i.e. the magnitudes of the current bottom shear stresses are not the same. Consequently, the two net sediment transports are not comparable. In such cases, the correct thing to do is to adjust the total current discharge to maintain the same magnitude of current bottom shear stress.