2. Model equations

The model considered in this work applies to incompressible flows where, whenever a free surface is present, surface tension is ignored. The problem is readily converted to non-dimensional form using the depth or wall-to-wall distance $h$ and surface current peak velocity $U_0$, as follows:

(2.1)\begin{align} (x,y,z) & \mapsto (x,y,z)h,\quad {\boldsymbol{k}}\mapsto {\boldsymbol{k}}/h,\quad t \mapsto th/U_0, \nonumber\\ \hat{\boldsymbol{u}}_\textrm{t} & \mapsto \hat{\boldsymbol{u}}_{t}U_0,\quad \hat p_\textrm{t} \mapsto \hat p_\textrm{t} \rho U_0^2. \end{align}

Here ${\boldsymbol {k}}=(k_x,k_y)$ is the wavenumber in the surface plane, $\hat {\boldsymbol {u}}_\textrm {t}$ the fluid velocity and $\hat p_\textrm {t}$ the pressure. Hatted quantities pertain to real space with the Fourier (wave vector) space counterparts written without hats. Flows over sinusoidally shaped boundaries of low steepness are in focus in the present study, as sketched in figure 1, yet any bathymetry shape can be analysed by superposition according to Fourier's theorem. A slip velocity is assumed at the two boundary reference planes $z=0$ and $1$, here achieved by stretching commonly used velocity profiles a displacement length $\delta$ such that the stagnation points of the bulk current $U(z)$ fall outside of the interior domain. Small undulations observed in streamlines near an actual corrugated wall can then approximately be modelled as the displaced wall sketched in figure 1. Crude as this approximation may appear in terms of boundary layers, it suffices for our purposes because the nature of the vortex mechanism to be studied is kinematic, not dynamic, and does not rely on wall friction (other than that needed for generating the principal shear). The boundaries, located at $z=\hat \eta _\textrm {b}(\boldsymbol {r})$ and $z=1+\hat \eta _\textrm {s}(\boldsymbol {r})$, perturb and redirect current energy to undulate the velocity field within. Here, $\boldsymbol {r}=(x,y)$ is the horizontal coordinate.

Figure 1. Sketch of the problem set-up. (Magnitudes of $\eta _\textrm {s}$ and $\eta _\textrm {b}$ are exaggerated.) (a) Wall bounded and (b) free surface.

The problem to be solved consists of the Navier–Stokes equations, along with a kinematic condition at the boundaries:

(2.2a)\begin{gather} \left.\begin{aligned} \left[ \frac{\partial \hat{\boldsymbol{u}}_\textrm{t}}{\partial t} \right] + ( \hat{\boldsymbol{u}}_\textrm{t} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}})\hat{\boldsymbol{u}}_\textrm{t} + \hat{\boldsymbol{\nabla}} \hat p_\textrm{t} & = \left[ Re^{-1} \hat\nabla^2\hat{\boldsymbol{u}} \right] - Fr^{-2} \boldsymbol{e}_z\\ \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} \hat{\boldsymbol{u}}_\textrm{t} & =0 \end{aligned}\right\};\quad \hat\eta_\textrm{b}(\boldsymbol{r}) \leq z\leq~1+\hat\eta_\textrm{s}(\boldsymbol{r}), \end{gather}

(2.2b)\begin{gather}\hat{\boldsymbol{u}}_\textrm{t} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} \hat\eta_\textrm{b} = \hat w;\quad z = \hat\eta_\textrm{b}(\boldsymbol{r}), \end{gather}

(2.2c)\begin{gather}\hat{\boldsymbol{u}}_\textrm{t} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} \hat\eta_\textrm{s} = \hat w;\quad z = 1+\hat\eta_\textrm{s}(\boldsymbol{r}), \end{gather}

with $\hat {\boldsymbol {\nabla }}=(\partial _x,\partial _y,\partial _z)$. The Froude number is $Fr=U_0/\sqrt {g h}$ and $Re= U_0 h/\nu$ is the Reynolds number ($g$ being the gravitational acceleration and $\nu$ the kinetic eddy viscosity). The viscous and transient terms placed in square brackets are considered only in relation to the resonant second-order mode interaction—first-order solutions are assumed inviscid and at a steady state. When the upper surface is free, $\hat \eta$ is a function of the local flow and the inviscid dynamic boundary condition that pressure be constant at the surface is imposed.

3. Linearised solution

Assume that the motions attributed to the perturbed boundaries are small compared with the current velocity $U(z)\,$. Separating these, we have

(3.1a)\begin{gather} \hat{\boldsymbol{u}}_\textrm{t}(\boldsymbol{r},z)\, = \boldsymbol{U}(z) + \hat{\boldsymbol{u}}(\boldsymbol{r},z,t)\,; \quad \boldsymbol{U} = (U,0,0), \end{gather}

(3.1b)\begin{gather}\hat p_\textrm{t}(\boldsymbol{r},z)\, = Fr^{-2}(1-z) + \hat p(\boldsymbol{r},z,t)\,. \end{gather}

A constant reference pressure term has been set to zero. A steady state of the first-order motion can be reached without the influence of viscosity. Assuming such a state has been reached and the Reynolds number is large, we neglect the transient and viscous terms in the first-order motion. (Note that viscous effects in the steady state are indirectly modelled via the shape of $U(z)$.) Let the boundary topography $\hat \eta$ be described as a superposition of modes $\eta ({\boldsymbol {k}}){\textrm {e}}^{{\textrm {i}} {\boldsymbol {k}}\boldsymbol {\cdot }\boldsymbol {r}}$ (real when summed), $\eta \in \{\eta _\textrm {b},\eta _\textrm {s}\}$. The flow field $\{\hat {\boldsymbol {u}},\hat p\}(\boldsymbol {r},z)$ of the linear solution becomes a superposition of these wave modes $\{\boldsymbol {u},p\}({\boldsymbol {k}},z)$. After linearisation and insertion of (3.1), the linearised Euler equations read as

(3.2a)\begin{gather} {\textrm{i}} k_x U \boldsymbol{u} + \boldsymbol{U}'(z) w +\boldsymbol{\nabla} p = 0, \end{gather}

(3.2b)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

where $\boldsymbol {\nabla }=({\textrm {i}} k_x,{\textrm {i}} k_y,\partial _z)$. Eliminating $u$, $v$ and $p$ then yields the Rayleigh equation

(3.3)\begin{equation} w''-{ \left( k^2+U''/U \right) }w = 0, \end{equation}

$k=(k_x^2+k_y^2)^{1/2}$.

The kinematic boundary conditions (2.2b) and (3.1b), linearised about the reference planes, read as

(3.4a)\begin{gather} w(0) = {\textrm{i}} k_{x}U(0) \eta_\textrm{b}, \end{gather}

(3.4b)\begin{gather}w(1) = {\textrm{i}} k_{x}U(1) \eta_\textrm{s}. \end{gather}

(We describe in Akselsen & Ellingsen (Reference Akselsen and Ellingsen2019a) a procedure for extending the lower boundary condition to bathymetries of finite amplitude, but we will here consider only the linear conditions.) If free surface flow is considered then the additional linear dynamic condition

(3.5)\begin{equation} \eta_\textrm{s} = Fr^{2} p(1)\, \end{equation}

governs the upper boundary, determining $\eta _\textrm {s}$.

3.1. Numerical solutions for arbitrary current profiles

An arsenal of methods are at our disposal for evaluating the boundary value problem (3.3)–(3.5). For example, Phillips et al. (Reference Phillips, Wu and Lumley1996) adopted a Galerkin method at this stage. The problem (3.3)–(3.5) is, however, rather simple and so we rather opt for a simple shooting technique.

Writing the Rayleigh equation (3.3) in terms of the pressure and integrating once yields

(3.6)\begin{equation} p''-2\frac{U'}{U}p'-k^2 p = 0. \end{equation}

Next, substituting the vertical velocity component by the pressure component using (3.2a) with (3.4), the boundary conditions on $p$ read as

(3.7)\begin{equation} p'(0) = [k_x U(0)]^2 \eta_\textrm{b},\quad \begin{cases} p'(1) = [k_x U(1)]^2 \eta_\textrm{s}; & \text{[fixed boundaries],}\\ p'(1)-[Fr\, k_x U(1)]^2 p(1) =0; & \text{[free surface].} \end{cases} \end{equation}

Equations (3.6) and (3.7) reveal that $p$ is independent of the sign of $k_x$ and $k_y$ if $\eta$ is. The velocity components in terms of $p$ are

(3.8a–c)\begin{equation} u =-\frac{p}{U} - \frac{U' p'}{k_x^2 U^2},\quad v = - \frac{k_y p}{k_x U},\quad w = \frac{{\textrm{i}} p'}{k_x U}. \end{equation}

Integration of (3.6) is performed numerically using a standard ordinary differential equation (ODE) solver and the guess for $p(0)$ adjusted according to the error in (3.7). Computation time has never been found to exceed a second on a regular laptop computer, and the procedure yields $p(z)$ and its derivative, from which the velocity components are directly retrieved using (3.2a).

3.2. Analytical solution for a power-law current profile

Following Phillips & Shen (Reference Phillips and Shen1996) and Phillips et al. (Reference Phillips, Wu and Lumley1996), we now consider the particular family of power-law current profiles where $U$ is proportional to $z$ raised to a power $q$. The power law allows for analytical solutions while representing a wide variety of primary flow profiles. We assume the free surface geometry of figure 1(b). Several forms of the commonly studied current profiles are obtained with different values of $q$; the uniform current most commonly considered is then recovered by setting $q$ to zero. A linear current profile is recovered when $q=1$. The linear profile has been investigated frequently in the literature since it is—in two dimensions—the only rotational flow for which potential theory is applicable. In the intermediate range $0<q<1$, concave-up profiles of the kind sketched in figure 1(b), resembling that observed in turbulent bottom boundary layer flows reside. A range of intermediate exponent values have over the years been suggested for turbulent flows (Chen Reference Chen1991; Cheng Reference Cheng2007). Flows with a surface shear layer may be modelled with $q > 1$, resulting in a class of concave-down profiles.

Akselsen & Ellingsen (Reference Akselsen and Ellingsen2019a) used a profile $U=z^q$ with the $z$-axis defined so that the lower boundary reference plane was at $z=\delta$. This plane is in the current formalism located at $z=0$ and we instead stretch and shift the vertical axis:

(3.9a,b)\begin{equation} U(z) = [\zeta(z)]^q,\quad \zeta(z) = \frac{z+\delta}{1+\delta};\quad [0\leq z\leq 1]. \end{equation}

Inserting (3.9a,b) and adopting the substitution $w=\sqrt \zeta W(Z)$; $Z=k(z+\delta )$, the Rayleigh equation (3.3) is remoulded into the modified Bessel equation

(3.10)\begin{equation} W''(Z)+\frac{W'(Z)}{Z}-{ \left( 1+\frac{(1/2-q)^2}{Z^2} \right) }W(Z)=0 \end{equation}

whose two linearly independent homogeneous solutions are known. Written in terms of a generic flow variable $\phi$, we have

(3.11)\begin{equation} \phi = \sum_\pm b^\pm \phi^\pm \end{equation}

with

(3.12a)\begin{gather} w^\pm(z) = {\textrm{i}}\sqrt \zeta I_{\pm(q-1/2)}(\tilde k \zeta)\,, \end{gather}

(3.12b)\begin{gather}p^\pm(z) = \frac{k_x}{k} \zeta^{q+1/2} I_{\pm(q+1/2)}(\tilde k \zeta)\,, \end{gather}

(3.12c)\begin{gather}\boldsymbol{u}_{\textrm{h}}^\pm(z) = \frac{{\textrm{i}} q \zeta^{q-1} w^\pm \boldsymbol{e}_x-\tilde{\boldsymbol{k}} p^\pm}{\tilde k_x \zeta^q}, \end{gather}

where $I_\alpha$ is the modified Bessel function of the first kind of order $\alpha$ and $\boldsymbol {u}_{\textrm {h}} = (u,v)$ is the horizontal velocity vector. A stretched wave vector $\tilde {\boldsymbol {k}} = {\boldsymbol {k}}(1+\delta )$ has also been introduced. The remaining constants $a^+$ and $a^-$ are determined by the two boundary conditions; (3.4) and (3.5) yield the relationships

(3.13)\begin{equation} b^\pm =\begin{cases} {\textrm{i}} k_{x} \dfrac{\eta_\textrm{s} w^\mp(0)-U(0) \eta_\textrm{b} w^\mp(1) }{ w^\pm(1)w^\mp(0) -w^\pm(0)w^\mp(1)}, & \text{ (fixed boundaries)},\\ {\textrm{i}} k_{x} U(0) \eta_\textrm{b} \left[ w^\pm(0) -\dfrac{w^\pm(1)\,-{\textrm{i}} k_x Fr^2 p^\pm(1)\,}{w^\mp(1)\, - {\textrm{i}} k_x Fr^2 p^\mp(1)\,}w^\mp(0) \right]^{-1}, & \text{ (free surface)}. \end{cases}\end{equation}

Further details on the analytical free surface solution is given in Akselsen & Ellingsen (Reference Akselsen and Ellingsen2019a).

4. The resonant second-order wave interaction

Similar to Craik (Reference Craik1970), we consider boundary undulations composed of a pair of sinusoidal waves directed symmetrically about the streamwise direction $x$, i.e.

(4.1)\begin{equation} \hat\eta = \frac{a}{4}\left[{\textrm{e}}^{{\textrm{i}}( k_{1x} x + k_{1y} y)}+{\textrm{e}}^{{\textrm{i}}( k_{1x} x - k_{1y} y)} + \mathrm{c.c.} \right] =a\cos(k_{1x} x)\cos(k_{1y} y),\end{equation}

($\eta \in \{\eta _\textrm {b},\eta _\textrm {s}\}$, $a\in \{a_\textrm {b},a_\textrm {s}\}$). An illustration is shown in figure 2. A free surface will to linear order also have the same functional form as the bathymetry, thus behaving kinematically similar to an undulated upper wall with the important difference that the mean current and secondary flow see this as a full-slip boundary. The first-order wave modes involved each have amplitudes $a/4$ and the four wave vectors $(\pm k_{1x},\pm k_{1y})$ (the sign of each component is varied independently). Second-order harmonics, in turn, have wave vectors $\boldsymbol {\kappa }$ which are sums of pairs of these. The harmonics thus come in four different types with wave vectors $\boldsymbol {\kappa }=\pm 2(k_{1x},k_{1y})$, $\boldsymbol {\kappa } = 0$, $\boldsymbol {\kappa } =(\pm 2k_{1x},0)$ and $\boldsymbol {\kappa } = (0,\pm 2k_{1y})$. The last type of harmonic is resonant—see Craik (Reference Craik1970) and Akselsen & Ellingsen (Reference Akselsen and Ellingsen2019b) for further details. The resonance will manifest in the formation of vortex pairs aligned in the current flow direction, as illustrated in figure 2.

Figure 2. Boundary topography and the orientation of vortices relative to it: (a) bottom topography with wave vectors $(k_{1x},\pm k_{1y})$ indicated; (b) schematic illustration of swirling streamlines in the cross-flow plane.

From here we adopt $\boldsymbol {u}$ and $p$ as the second-order components unless otherwise stated. The equations of motion for the second-order components of the Stokes expansion read as

(4.2a)\begin{gather} (\partial_t-Re^{-1}\nabla^2)\boldsymbol{u} + \boldsymbol{U}' w + \boldsymbol{\nabla} p = -[(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}]_{1\times}, \end{gather}

(4.2b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}=0 \end{gather}

for this particular interaction. Left-hand variables are second-order components while the right-hand cross terms are the second-order interactions of first-order components. Eliminating the horizontal velocity components and the pressure, the Orr–Sommerfeld equation for the second-order vertical velocity $w(z,t;{\boldsymbol {k}}_1)$ becomes

(4.3)\begin{gather} (\partial_t-Re^{-1}\nabla^2)\nabla^2 w = \mathcal{R}(z)\,; \end{gather}

(4.4)\begin{gather}\mathcal{R} ={\textrm{i}} \boldsymbol{\kappa} \boldsymbol{\cdot} \partial_z [(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}_{\textrm{h}}]_{1\times} + \kappa^2 [(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) w]_{1\times} \end{gather}

with $\nabla ^2=\partial _z^2-\kappa ^2$, $\kappa =2|k_{1y}|$. Due to symmetry, we have $\boldsymbol {\nabla }=(0,{\textrm {i}}\kappa _y,\partial _z)$ and $[(\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla })u]_{1\times } = 0$ which simplifies the right-hand expression (4.4) significantly; using (3.8a,c) to express the first-order velocity component in terms of the first-order pressure, we find the remarkably simple expression

(4.5)\begin{equation} \mathcal{R}(z)= 8 \frac{k_{1y}^2}{k_{1x}^2}\frac{U'}{U^3} \left[ { \left( k_{1x}^2-k_{1y}^2 \right) }p_1^2 + (p_1')^2 \right]. \end{equation}

The first-order pressure mode $p_1$ is the solution of (3.6)–(3.7) with ${\boldsymbol {k}}={\boldsymbol {k}}_1$. Note that the curvature term $(U''/U)w$ is not present in (4.3) because $\kappa _x = 0$, i.e. the resonant wave is two-dimensional and orthogonal to the shear current.

For arbitrary profiles, considered in § 3.1, the right-hand side of (4.3) is evaluated directly from the output of the ODE integrator. If the power-law profile (3.9a,b) is adopted (§ 3.2) then the right-hand side can be written explicitly as

(4.6)\begin{equation} \mathcal{R}(z) = 8 \frac{q k_{1y}^2 }{1+\delta} \sum_{s=\pm1}\frac{k_{1x}^2-sk_{1y}^2}{k_1^2} \left\{\sum_{\pm} b_1^{\pm} I_{\pm(q+s(1/2))} \big[k_1(z+\delta)\big]\right\}^2,\end{equation}

$b_1^\pm$ and its velocity components are as given in (3.11)–(3.13) with ${\boldsymbol {k}}={\boldsymbol {k}}_1$.

The resonant second-order harmonic is uniform in $x$ which allows for a streamfunction

(4.7)\begin{equation} \psi={\textrm{i}} w/\kappa_y \end{equation}

to be expressed in the $yz$-plane. An expression for the pressure component is given in the appendix.

The resonant wave interactions will come with both signs $\boldsymbol {\kappa }=(0,\kappa _y)=(0,\pm 2k_{1y})$, with the corresponding modes being complex conjugates of each other. It is therefore sufficient to consider only the positive spanwise wavenumber $\kappa _y=\kappa =2|k_{1y}|$ and write, for any flow variable $\hat \phi$,

(4.8)\begin{equation} \hat\phi = 2 \operatorname{Re}(\phi)\cos(\kappa y) - 2\operatorname{Im}(\phi) \sin(\kappa y).\end{equation}

The streamfunction $\psi$ is purely imaginary when $\eta _\textrm {b}$ and $\eta _\textrm {s}$ are real. Contour plots of $-2\,\operatorname {Im}\,\psi$ are therefore presented among the results of § 6. When $-\operatorname {Im}\,\psi >0$, the vertical flow is directed in towards lines aligned with peaks and troughs in the wall topography, with the horizontal velocity perturbation $\hat u$ being strongest along these lines. This serves to increasingly tilt vortex particle trajectories in the $xy$-plane with time. Vortex rotation is in the opposite direction when $-\operatorname {Im}\,\psi <0$.

In what follows, two limits of (4.3) shall be investigated, namely the inviscid transient problem $Re^{-1} = 0$ and the viscous steady-state problem $\partial _t\boldsymbol {\cdot } \rightarrow 0$.

4.1. Inviscid transient problem

We now set $Re^{-1} = 0$. Because the harmonic in question is resonant it will not reach a time-independent state without the intervention of viscosity or nonlinear dynamics. The solution of (4.3) can be obtained using the method of variation of parameters which results in

(4.9a,b)\begin{equation} w = \sum_\pm d^\pm {\textrm{e}}^{\pm \kappa z} + w_{\times}(z)\,;\quad w_{\times}(z)\, = \frac{t}{\kappa} \int_0^z \textrm{d} \xi\mathcal{R}(\xi)\, \sinh \kappa (z-\xi). \end{equation}

We have here imposed circulation quiescence at $t=0$. The function $w_{\times }(z)$ is evaluated at all $z$ with simple numerical integration.

The second-order boundary conditions for fixed wall boundaries, Taylor expanded about $z=0$ and $z=1$, reduce to

(4.10)\begin{equation} w(0,t) = w(1,t) = 0. \end{equation}

A little more caution is required for a free surface boundary; kinematic and dynamic conditions then reduce to

\[ \partial_t \eta_\textrm{s}-w =0;\quad Fr^2 p -\eta_\textrm{b} = Fr^4 p_1 p_1',\quad [z=1]. \]

Combining these with the $y$-component of the second-order momentum equation (4.2a) and the continuity equation (4.2b) yield

(4.11)\begin{equation} \partial_t^2\partial_z w +\kappa_y^2 Fr^{-2}w = 0,\quad [z=1]. \end{equation}

Vertical velocity $w$ is only linear in time so that the boundary conditions become (4.10) also at a free surface. Fixing the integration constants of (4.9a,b), our inviscid solution reads as

(4.12)\begin{equation} w = w_{\times}(z)\, - w_{\times}(1)\frac{\sinh\kappa z}{\sinh \kappa}. \end{equation}

The streamwise velocity component $u$ is retrieved directly from inserting (4.12) into the momentum equation (4.2a). The interaction term disappears due to symmetry and we obtain

(4.13)\begin{equation} u(z,t)=-\frac{t}{2}U'(z) w(z,t), \end{equation}

which is proportional to $t^2$.

4.2. Viscous steady-state problem

We now turn to the viscous problem assuming a steady state has been reached. It should be noted here that the second-order solution constructed here is driven by the interaction kinematics of first-order waves which are themselves assumed to be inviscid; we assume that the periodic dynamics of the first-order waves are affected little by viscosity at moderate to large Reynolds numbers, and that the role of viscosity at that order is in all essentials accounted for through the prescribed shear of $U(z)$. The resonant mode, on the other hand, continues to grow in time until it becomes strong enough for viscosity to become significant, and the flow reaches a steady state. The solution to (4.3) when denying any time dependency can be found by applying the solution formula of the previous problem twice. It is

(4.14)\begin{gather} w = \sum_\pm (d_0^\pm + z d_1^\pm) {\textrm{e}}^{\pm \kappa z} + w_{\times}(z)\,; \end{gather}

(4.15a,b)\begin{gather}w_{\times}(z)\, = \frac{Re}{2\kappa^3 } \int_0^z \textrm{d} \xi \mathcal{R}(\xi)\, G[ \kappa (z-\xi)];\quad G(Z)=\sinh(Z)-Z\cosh(Z). \end{gather}

A full-slip condition is assumed at the upper boundary if considering a free surface, as is often reasonable for an atmospheric interface between water and air. We find that all first-order cross terms cancel also for these stress conditions, even after Taylor expansion about $z=1$ (see also Craik Reference Craik1970). Appropriate full-slip stress conditions are $v'+ {\textrm {i}} \kappa _y w = 0$ about the reference plane. The other boundary conditions discussed in § 4.1 still hold; employing the continuity equation (3.2b) yields

(4.16)\begin{equation} w=w''=0 \quad[\text{full-slip}] \end{equation}

at $z=0$ or $1$ if considering the upper or lower boundary, respectively. We will also consider no-slip boundaries, which with (4.2b) implies that

(4.17)\begin{equation} w=w'= 0\quad[\text{no-slip}] \end{equation}

at $z=0$ or $1$. With these conditions we have derived the integration constants $d_0^+$, $d_0^-$, $d_1^+$ and $d_1^-$ in (4.14) for the three cases, full-slip/full-slip, no-slip/no-slip and full-slip/no-slip, at the upper/lower boundary. Explicit solutions for $w(z)$ are presented in appendix A.

The horizontal velocity component of the vortex motion is again obtained by integrating the $x$-component of the linearised momentum equation (4.2a). One finds that

(4.18)\begin{equation} u(z) = \sum_\pm d_u^\pm {\textrm{e}}^{\pm\kappa z}+ \frac{Re}{\kappa} \int_0^z\textrm{d}\xi U'(\xi) w(\xi)\sinh\kappa (z- \xi). \end{equation}

Appropriate no-slip and full-slip boundary conditions are $u=0$ and $u'=0$, respectively, evaluated at the appurtenant reference planes. (Cross terms merely represent the second-order correction of a first-order motion which itself does not support stress conditions.) Coefficients $d_u^\pm$ are given in appendix A.

5. Comparison with strong shear CL2 instability and the theory of the generalized Lagrangian mean

We will start this section by relating our findings to GLM theory before considering features of strong shear CL1 instability compared with strong shear CL2 instability.

Andrews & McIntyre's (Reference Andrews and McIntyre1978) theory of the GLM has in recent years dominated in the study of CL2-type Langmuir circulation (Leibovich Reference Leibovich1980; Craik Reference Craik1982a,Reference Craikb; Phillips Reference Phillips1998, Reference Phillips2005). This theory is commonly used to derive the eigenvalue problem that governs the stability of unidirectional waves to flows of longitudinal vortex form. In order to connect our work with the GLM-based literature, we shall briefly discuss how key observations from the GLM theory relates to the CL1-type circulation considered in the present paper.

The GLM Navier–Stokes equations, averaged in $x$, read as

(5.1)\begin{equation} \overline{D}{}^\textrm{L} (\overline{\hat u}{}^\textrm{L}_i - \hat\wp_i) + \overline{\hat u}{}^\textrm{L}_{j,i}(\overline{\hat u}{}^\textrm{L}_j- \hat\wp_j)+\hat\pi_{,i} = Re^{-1} \left[ \overline{( \hat\nabla^2\hat u_i)}{}^\textrm{L} + \overline{\hat\xi_{j,i}(\hat\nabla^2\hat u_j)^l} \right]. \end{equation}

Except for the symbol for the pseudomomentum

(5.2)\begin{equation} \hat\wp_i\equiv -\overline{\hat\xi_{j,i} \hat u_{j}^l} \end{equation}

and the use of hats to denote physical space variables, we have here adopted the common notation. The operator $\overline {(\cdot )}{}^\textrm {L}$ is the average along a Lagrangian trajectory described by the displacement $\boldsymbol {\xi }$ of fluid particles. The fluctuating part of the Lagrangian velocity is denoted by $\hat {\boldsymbol {u}}^l=\hat {\boldsymbol {u}} (\hat {\boldsymbol {x}} + \hat {\boldsymbol {\xi }},t)-\overline {\hat {\boldsymbol {u}}}{}^\textrm {L}$. Gradient terms are lumped into $\hat \pi$ which may be thought of as an effective pressure. Overlines denote the streamwise average, indices $\{1,2,3\}$ refer respectively to the spatial dimensions $\{x,y,z\}$, repeated indices are summed and comma denotes differentiation. We refer the reader to the above references, in particular Andrews & McIntyre (Reference Andrews and McIntyre1978) and Craik (Reference Craik1982a,Reference Craikb), for a more complete description.

Envisage a wave field of $O(\epsilon )$ generated by the wall topography. Only $O(\epsilon ^2)$ terms can survive the streamwise averaging. To $O(\epsilon ^2)$, (5.1) reduces to

(5.3)\begin{equation} \partial_t (\overline{\hat{\boldsymbol{u}}}{}^\textrm{L} - \hat{\boldsymbol{\wp}})+\hat{\boldsymbol{\nabla}}[U \overline{\hat u}{}^\textrm{L}_1 + \hat\pi-\tfrac12 U^2] -\hat\wp_1 \hat{\boldsymbol{\nabla}} U+ \overline{\hat u}{}^\textrm{L}_3 \boldsymbol{U}_{,3}- Re^{-1}\hat\nabla^2\overline{\hat{\boldsymbol{u}}}=0. \end{equation}

The streamwise averaged Eulerian velocity, pseudomomentum and differential drift (Stokes drift) $\overline {\hat {\boldsymbol {u}}}{}^\textrm {S}=\overline {\hat {\boldsymbol {u}}}{}^\textrm {L}-\overline {\hat {\boldsymbol {u}}}$ will at $O(\epsilon ^2)$ have modes

\[ [\boldsymbol{u}(z,t),\boldsymbol{\wp}(z,t),\boldsymbol{u}^\textrm{S}(z,t)]{\textrm{e}}^{{\textrm{i}} \kappa_y y}. \]

The streamwise component of $\boldsymbol {u}$ obeys, by (5.3),

(5.4)\begin{equation} (\partial_t - Re^{-1} \nabla^2 )u_{1}= \partial_t ( \wp_1 - u_1^\textrm{S}) - (u_3+u^\textrm{S}_3) U_{,3}. \end{equation}

The other components of (5.3), together with continuity $\boldsymbol {\nabla }\boldsymbol {\cdot }\overline {\boldsymbol {u}}=0$, can at $O(\epsilon ^4)$ be reduced into an Orr–Sommerfeld type equation

(5.5)\begin{equation} (\partial_t-Re^{-1}\nabla ^2)\nabla ^2 u_3 = - \kappa^2 \wp_1 U_{,3} - \partial_{t}[{\textrm{i}} \kappa_y (\wp_2-u_2^\textrm{S})_{,3} + \kappa^2 (\wp_3-u_3^\textrm{S})], \end{equation}

$\nabla ^2=\partial _{zz}-\kappa ^2$. The second right-hand term vanishes when we next consider a time-independent primary flow.

For evaluating (5.5), one must first find the linear displacement $\boldsymbol {\xi }$. Following Craik (Reference Craik1982a), we have derived the particle trajectories for time independent $U$, adopting Craik's method for matching trajectories to the appropriate streamlines. Writing $\hat {\boldsymbol {\xi }}$ as a superposition of modes ${\boldsymbol {\xi }}{\textrm {e}}^{{\textrm {i}}(k_x x + k_y y)}$ and using (3.8a–c) to express the $O(\epsilon )$ wave field, we eventually find the simple relation

(5.6)\begin{equation} \boldsymbol{\xi} = \frac{\boldsymbol{\nabla} p}{k_x^2 U^2}+O(\epsilon^2), \end{equation}

$\boldsymbol {\nabla } = ({\textrm {i}} k_x,{\textrm {i}} k_y,\partial _z)$. Expression (5.6) implicitly assumes the absence of critical layers where $U=0$, as is common. Adopting this to compute the pseudomomentum (5.2), we get

(5.7a)\begin{gather} \hat\wp_1 = -\frac{4}{k_x^2 U^3} \left[ { \left( k_x^2+k_y^2 \right) } p^2 + ( p')^2 \right] -\frac{4\cos(2k_y y)}{k_x^2 U^3} \left[ { \left( k_x^2-k_y^2 \right) } p^2 + ( p')^2 \right] + O(\epsilon^3), \end{gather}

(5.7b)\begin{gather}\hat\wp_2 = \hat\wp_3 = 0. \end{gather}

Inserting (5.7) into (5.5) and matching it to (4.3) successfully recovers (4.5). Generalised Lagrangian mean theory thus accords with the more primitive mode coupling theory utilized in § 4, although the benefits of GML theory are not as patent as in studies of CL2 instability. The horizontally uniform part of the pseudomomentum can, for $k_y\to \infty$, be compared with Craik's (Reference Craik1982a) equation (3.4a) with his $\phi =-4p'/\epsilon k_x^2 U$, whence we find that the spanwise-periodic topography generates, for the same amplitude $a$, half the pseudomomentum of the spanwise uniform topography, in addition to the spanwise periodic part.

When the shear intensity is of $O(1)$, the study of CL2-type instabilities turns more complicated than its weak shear counterpart due to the influence of wave distortion (Craik Reference Craik1982b; Phillips Reference Phillips2005). Wave distortion comes about as the pseudomomentum generated by the spanwise-periodic wave perturbation (here called the Langmuir wave) reaches a magnitude great enough to re-enter the eigenvalue problem governing the stability of the spanwise-periodic wave. The pseudomomentum generated by the spanwise-periodic wave itself must then be estimated in a separate analysis, first performed by Craik (Reference Craik1982b), yielding a coupled pair of differential equations governing the stability to longitudinal vortex form. This equation set has since been termed the ‘generalized Craik–Leibovich equations’ (CLg), with the range of validity extending beyond instability from wind-driven surface waves and into boundary layer flows.

Wave distortion enters the CL1 problem less directly due to three prominent differences between the two mechanisms. First, where the CL2 mechanism primarily consists of a distortion of the streamwise flow generating a spanwise-periodic wave, the CL1 mechanism works by directly imposing a spanwise-periodic wave of prescribed wavelength. Second, CL2 wave distortion is present immediately from the scaling of the initial perturbation of the primary flow, while the instability of CL1 is algebraic, growing linearly from out of a completely unperturbed state. If not sufficiently suppressed by viscosity, the wave may eventually grow until wave distortion takes effect or other longitudinal vortex instabilities occur. Third, where the spanwise-periodic wave in the CL2 mechanism generates pseudomomentum at the same wavelength, the CL1 equivalent is an interaction between two spanwise-periodic wave fields, generating higher and lower spanwise harmonics.

Wave distortion will eventually arise in CL1 if the viscosity is insufficient to restrain the circulation to a moderate intensity. This then occurs via pseudomomentum generated by the Langmuir wave. Only streamwise-periodic particle displacements and velocity perturbations can combine in (5.2) to produce non-zero pseudomomentum (by construction, $\overline {\hat {\boldsymbol {\xi }}}=0$). The next order of pseudomomentum relevant in terms of wave distortion is then generated by two primary $O(\epsilon )$ harmonics interacting with the Langmuir wave. The Langmuir wave must therefore be of $O(1)$ (and not $O(\epsilon )$) for this distortion to scale with the principal pseudomomentum (5.7). Such wave distortions will then be spanwise periodic with wavenumbers $k_{1y}$ and $3k_{1y}$. The streamwise component of the Langmuir wave is, similar to CL2, the first to become significantly large as this is proportional to $Re^2$. Flow simulations are presented later in § 7. Examining velocity fields from these simulations (examining $\overline {\hat u}$ minus its spanwise mean) reveals that the streamwise component of the Langmuir wave is about $5\,\%$ the magnitude of the principal flow at $Re_\tau =30$, $a=0.0635$, and significantly less for smaller Reynolds numbers and amplitudes. It seems quite possible that wave distortion generated by the Langmuir wave plays a role in the unstable behaviour and emergence of higher wavenumbers observed at larger Reynolds numbers. The theoretical model accords with the observed magnitude of the streamwise Langmuir wave component and indicates that the no-slip boundary condition for the Langmuir wave is active in reducing its magnitude below $O(\epsilon ^2 Re^2)$.

7. Simulated laminar flow

The perturbation theory presented so far sheds light on the subject of structurally generated Langmuir circulation, yet a vital question still remains to be answered: What happens to the boundary layer? The presence of Langmuir circulation beneath free surfaces is of course well established, but these boundaries allow for a slip velocity whereas flow stagnation takes place at a rigid boundary. Will Langmuir circulation even appear near a no-slip boundary? We have so far circumvented the issue by assuming that the bulk current velocity $U(z)$ does have some slip velocity at the reference planes $z=0,1$. This assumption is perhaps less alien when considering a turbulent boundary layer – assuming a slip governed by the ‘law of the wall’ is a common approach within turbulence modelling – but is not easily justified in laminar flows. We must then consider wall undulation amplitudes that penetrate some depth into the flow field, and these amplitudes must in turn somehow be related to the undulation amplitude $a$ and the displacement height $\delta$ introduced in the perturbation theory. A sketch is presented in figure 17 where we regard a streamsurface within the flow as the ‘seen’ surface. The closer this surface is to the actual wall (the smaller the $\delta$) the more the boundary undulations will resemble the actual known undulation of the wall. Further away the ‘seen’ amplitude will be smaller and less sinusoidal as also the opposing boundary contributes to compress the streamlines. On the other hand, the perturbation theory rests on the assumption that the variation in current velocity across the perturbation is small relative to the velocity at the reference plane, which improves in validity at larger $\delta$.

Figure 17. Conceptual representation near the boundary with no-slip at the walls with a slip velocity reference boundary located some distance $\delta$ within the flow. Streamlines resemble the actual boundary curvature more closely when $\delta$ is smaller than the amplitude $a$ (a), although this makes the perturbation amplitudes large relative to the assumptions of the linear theory. The situation is the opposite if $\delta >a$ (b).

We have at our disposal a code for numerical flow simulations using the lattice–Boltzmann method, with which we have generated simulation results within the laminar flow regime for the purpose of demonstrating that the artefacts of our perturbation theory do in fact manifest in real, nonlinear flows. Although the lattice–Boltzmann method is usually not first choice in terms of accuracy and computational efficiency, it proved easily adaptable to the present problem since the undulated boundaries could be represented simply by marking nodes outside the flow domain as ‘solid’ nodes. The code has been validated against two- and three-dimensional lid-driven cavity flow cases. These benchmarks are presented in the supplementary material available at https://doi.org/10.1017/jfm.2020.490.

In order to compare the simulated flow field with our theoretical results, we compute an approximation of $\psi$ from the streamwise averaged velocity field based on the principle of volume flux. Provided the vortices are not skewed ($\vartheta$ is a right or straight angle),

(7.1)\begin{equation} \tilde\psi(z) = \int_{0}^z \textrm{d} \xi \overline{\hat v}(y_0,\xi),\end{equation}

approximates $\psi$, $\overline {\hat v}$ being the streamwise mean of the spanwise velocity component. Spanwise position $y_0$ is the spanwise coordinate of the vortex centre. We here assume $y_0=\pi /4k_{1y}$. Simulation results for a range of wall Reynolds numbers $Re_\tau = (h/2)u_* /\nu$, where $u_*=\sqrt {\tau _{\textrm {w}}/\rho }$, are shown in figures 18 and 19. The wall Reynolds number is here the appropriately fixed parameter (a measure of the flow forcing) and relates to $Re=h U_0/\nu$ as $Re=Re_\tau ^2$ in the case of flat-walled, laminar Poiseuille flow. The figures show, when compared with the viscous, stationary results of figures 14–16, qualitative agreement between simulation and theory in terms of vortex orientation and rotational direction. Normalised steady-state quantities $\hat \psi /(Re\, a^2)$ and $\hat u/(Re\, a)^2$ are in the theory invariant to the Reynolds number, although Reynolds number dependency is evident in the simulation results, notably in the $\vartheta =\pi /2$ case where the value of the Reynolds number affects the degree to which vortices are merged to cover the cross-sections or exist as separate, co-rotating vortex pairs. It should also be noted that, contrary to the $\vartheta =0$ topography, the $\vartheta =\pi /2$ topography is not antisymmetric about the horizontal plain $z=1/2$. This generates some minor vertical asymmetry in each vortex due to dynamic effects. The asymmetry differs between the four vortices in each spanwise period; we here display the vortex centred at $y=\pi /4k_{1y}$.

Figure 18. A double wall-bounded flow. Normalised streamfunctions $\tilde \psi /(Re\,a^2)$ estimated from the steady states of lattice–Boltzmann simulations using (7.1): (a) $\vartheta =0$, $a= 0.0391$, (b) $\vartheta =0$, $a= 0.0625$, (c) $\vartheta =0$, $a= 0.0938$, (d) $\vartheta =\pi /2$, $a= 0.0391$, (e) $\vartheta =\pi /2$, $a= 0.0625$, (f) $\vartheta =\pi /2$, $a= 0.0938$. Here $k_1 \approx 2\pi$, $\theta \approx \pi /8$.

Figure 19. Numerically computed streamlines of the two-dimensional velocity field obtained when averaging along the streamwise dimension: (a) $\vartheta =0$, $a = 0.0391$, (b) $\vartheta =0$, $a = 0.0625$, (c) $\vartheta =0$, $a = 0.0938$, (d) $\vartheta =\pi /2$, $a = 0.0391$, (e) $\vartheta =\pi /2$, $a = 0.0625$, (f) $\vartheta =\pi /2$, $a = 0.0938$. Here $Re_\tau = 20$, $k_1 \approx 2\pi$, $\theta \approx \pi /8$.

Figure 20 shows slices of the velocity field in the $yz$-plane at a prescribed $x$-location. These are not averaged in the streamwise direction and so give an impression of the circulation strength relative to the first-order undulating motion. In the cases studied, the two types of motion are found to be of the same order of magnitude, with first-order motion dominating near peaks and troughs, and circulation dominating in the intervening regions. The slices shown in the plots of figure 20 are at the somewhat arbitrary location $x=\pi /6k_{1x}$. Circulation intensity relative to the undulating motion increases with increasing Reynolds numbers.

Figure 20. Streamline slices of a full three-dimensional flow field: (a) $Re_\tau =15$, $\vartheta =0$, (b) $Re_\tau =30$, $\vartheta =0$, (c) $Re_\tau =15$, $\vartheta =\pi /2$, (d) $Re_\tau =30$, $\vartheta =\pi /2$. Here $a\approx 0.0625$, $k_1 \approx 2\pi$, $\theta \approx \pi /8$. Slice at $x=\pi /6k_{1x}$.

Quantitatively, the simulations most closely resemble the no-slip boundary condition cases where the value for $\delta$ is small. This dependency escalates with increased wall undulation amplitude. For the sake of briefness and clarity, we have not attempted to adjust the simulation parameters to better accommodate our theoretical construct. Indeed, the undulation amplitudes used in the normalisation are the actual amplitudes and not an estimate of the streamsurface amplitudes seen by the flow above some displacement thickness. (Recall that the circulation is theoretically proportional to these amplitudes squared.) Likewise, the duct height $h$ is taken as the average distance between the actual walls.

Earlier, in § 6.1, we made the observation that large values of $\theta$ (i.e. boundary undulations that are stretched out in the streamwise direction rather than the spanwise) can generate vortices which rotate in the opposite direction relative to the boundary topography (pushing fluid away from peaks and troughs as opposed to in towards them). Simulation results for such a case, where $\theta = 3\pi /8$ and $k_1=\pi$, are presented in figure 21 along with some corresponding theoretical streamfunctions. We observe that the rotation is indeed reversed in the bulk interior of the cross-section (compared with, for example, figure 18), although thin additional vortices appear close to the walls which still rotate in an unreversed direction. This may be an artefact of the dynamic stresses caused by the undulating boundaries which has not been incorporated into our kinematic theory. (The only means with which other harmonics are filtered out field is through averaging in the streamwise direction.) Viscosity serves to modify the first-order wave field and give a spanwise perturbation to the primary flow. Dynamic effects may generate circulation with its own rotational preference, as observed in turbulent flows by Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015), Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2018), Yang & Anderson (Reference Yang and Anderson2018), Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015) and Kevin et al. (Reference Kevin, Monty, Bai, Pathikonda, Nugroho, Barros, Christensen and Hutchins2017). Presumably, a dynamic rotational effect is always present, either assisting or impeding the Langmuir circulation.

Figure 21. ‘Rotating the boundary topography’; $\theta = 3\pi /8$. This will, according to theory, cause rotation in the opposite direction relative to the boundary undulations. Here $k_1=\pi$, $\vartheta = 0$. (a) Theoretical, viscous, no-slip and (b) simulated, $a = 0.0625$.

Secondary vortex motion, aligned with the flow, can be generated via spanwise intermittent roughness patches (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015) or streamwise-aligned obstacles (Sirovich & Karlsson Reference Sirovich and Karlsson1997; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Kevin et al. Reference Kevin, Monty, Bai, Pathikonda, Nugroho, Barros, Christensen and Hutchins2017; Yang & Anderson Reference Yang and Anderson2018). Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015) demonstrated that these structures are related to Prandtl's secondary flow of the second kind, driven and sustained by spatial gradients in the Reynolds-stress components.

Our simulations start showing signs of vortex instability at wall Reynolds numbers above 30. An example is shown in figure 22. The degree of instability also increases with increased boundary undulation amplitude and seems more unstable for $\vartheta =\pi /2$ than $\vartheta =0$ as the former generate vortices which are more stretched out over the cross-section. Vortex instability manifests as the spawning of higher wavenumber vortices. These are naturally less steady as the Reynolds number increases. Longitudinal vortex instability is a well-studied phenomenon since their discovery by Craik (Reference Craik1977), and their effect on strongly sheared flows such as those considered here thoroughly investigated by Phillips and co-workers (Phillips & Wu Reference Phillips and Wu1994; Phillips & Shen Reference Phillips and Shen1996; Phillips et al. Reference Phillips, Wu and Lumley1996; Phillips& Dai Reference Phillips and Dai2014). Vortex instability is likely to play a decisive role in the evolution of flows over undulating topographies as we move towards the turbulent regime.

Figure 22. Unstable vortex field: $Re_\tau = 40$, $\vartheta =\pi /2$, $\theta = \pi /8$, $k_1 = 2\pi$.

8. Summary and closing remarks

We have studied the creation of streamwise vortices in shear flows over a wall with a particular topography. The roll-like flow structures are closely analogous to Langmuir circulation created through a strong shear ‘CL1’-type mechanism for the algebraic instability of longitudinal vortices. Instead of a shear flow beneath the wavy free surface normally associated with Langmuir circulation, the same mechanism can be imposed on a flow by furnishing walls with a criss-crossing wavy pattern. Under sufficient influence of viscosity, the flow structures due to the CL1 instability will stabilise and form coherent longitudinal vortices whose transverse size matches the wavelength of the chosen wall topography. This mode of circulation requires no external forcing other than maintaining the ambient shear flow, but functions by redirecting the vorticity inherently present in the current itself.

Assuming a shear profile of arbitrary form $U(z)$ when unperturbed, a resonant wave–current interaction is computed in two limiting cases: transient inviscid circulation and steady-state viscous circulation. The vortical motion in the cross-flow plane grows proportional to $a^2 t$ in the early, essentially inviscid stage ($a$, wave-pattern amplitude; $t$: time), and is proportional to $Re\,a^2$ in the latter, viscous steady state. The instability is algebraic in nature, growing linearly with time in the former case and to the Reynolds number in the latter. The streamwise velocity perturbations are respectively proportional to these quantities squared.

After the onset of the instability, at first purely kinematically driven and essentially inviscid, viscosity serves to push the vortex centres out away form the undulating boundaries. Vortex intensity is observed to be greatest in the wavelength range where vortex height and width are comparable. Circulation is also typically strongest when the half-angle $\theta$ between the crossing plane waves making up the sinusoidal wall corrugations lies in the range $10^\circ <\theta <25^\circ$. A reversal of the direction of circulation is possible within the range $\theta >45^\circ$, albeit with weaker rotation strength. Several weak vortices can coexist vertically in the region where this switch of rotational direction takes place. For a channel flow between two similarly corrugated walls, the circulation generated near each wall can merge into single vortices spanning the cross-section, provided the two boundaries are suitably shifted relative to each other. Such a merger is more pronounced the closer the boundaries are situated relative to vortex wavelength.

The investigation has been supplemented with laminar flow simulations using a lattice–Boltzmann method. These simulations show that the forced circulation phenomenon does indeed take place also in laminar flows over no-slip boundaries, and that the resulting circulation is qualitatively similar to the results of the perturbation theory even though the latter conceptually assumes a current velocity at the walls. The laminar simulations furthermore show that circulation is quickly established relative to the time it takes for the primary flow to stabilise. The circulation intensity is comparable to the intensity of the first-order motion. Dynamic effects, excluded from the theoretical model, were here also observed in near-wall regions. Quantitatively, a degree of uncertainty is related to the closure assumption of a displacement height $\delta$ adopted in the theory.

The practical potential of the mechanism here studied is as yet unexplored and subject to conjecture. Literature, as presented in § 1, suggests that longitudinal vortices can be beneficial in terms of turbulent drag reduction, although the divers mechanisms hitherto employed to create such vortices differ fundamentally from the Langmuir mechanism suggested herein. Therefore, studies of the continued development of CL1-type Langmuir circulation in the turbulent regime are desirable, particularly regarding the stability of the generated structures for increasing Reynolds numbers. The optimal scale of induced Langmuir vortices is an open question as vortices could be made to exist either as a tiny structure deep within the boundary layer, or as large structures reaching into the bulk of the flow, depending on the scale of the wall undulations. As indicated by the literature, tiny vortex structures could conceivably act to stabilise the transition to turbulence while large vortices could act to reduce drag in flows which are already fully turbulent.