2.1. Framework

The OCFs can be described by RANS equations flanked by suitable boundary conditions, such as the kinematic and dynamic conditions at the free surface, and no slip and impermeability at the bed. We will adopt the bulk longitudinal velocity $\hat {U}_0$, the channel depth $\hat {h}_0$, the corresponding hydrodynamic time scale $\hat {h}_0/\hat {U}_0$ and $\rho {\hat {U}}_0^2$ to non-dimensionalize velocities, lengths, times and pressure, respectively ($\rho$ is the water density). Body forces are non-dimensionalized by the module of the gravity acceleration. Henceforth, the hat symbol identifies dimensional quantities.

As reported in the sketch of figure 1, we consider the (dimensionless) reference frame ${x}_i=(x,y,z)$, having axes parallel to the streamwise, bed-normal and spanwise directions, respectively. The corresponding Reynolds-averaged velocity term is $\boldsymbol {u}:=(u,v,w)$. Assuming high Reynolds numbers and hence negligible viscous stresses, the non-dimensional RANS equations read as

(2.1)$$\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial{t}}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}+\boldsymbol{\nabla}{p} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}-\frac{\boldsymbol{k}}{F^2}=0, \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

where $F=\hat {U}_0/({g \hat {h}_0})^{1/2}$ is the Froude number, g is the gravity acceleration, p is pressure, $\boldsymbol {\tau }$ is the dimensionless Reynolds stress tensor, $\boldsymbol {k}=(\sin {\theta },-\cos {\theta },0)$ is the projection of the gravity unit-vector along the chosen coordinate system, with $\theta$ being the slope of the bed, with respect to a horizontal plane (see figure 1).

2.2. Secondary currents of Prandtl's second kind in OCF

Secondary currents are cellular structures with longitudinal vorticity (see figure 1$a$) that appear in the time-averaged flow. In the present paper, we consider the case of a straight open channel in uniform flow conditions, whereby secondary flows are of Prandtl's second kind as they are the result of a combined effect between the three-dimensionality of the flow imposed by lateral boundaries and the anisotropy of turbulence (Nezu & Nakagawa Reference Nezu and Nakagawa1993). Secondary currents are clearly triggered in proximity to lateral boundaries where they occur in the form of ‘corner vortices’ and are characterized by intense lateral and vertical mean velocities. Their cellular structure then replicates along the lateral direction as a result of mass and momentum conservation, but their intensity fades out while moving towards the centre of the channel, so much so that they become difficult to detect experimentally (Zampiron, Cameron & Nikora Reference Zampiron, Cameron and Nikora2020). Nevertheless, Adrian & Marusic (Reference Adrian and Marusic2012) pointed out that, even for the case of very wide channels (i.e. channels whose width is much greater than their depth), secondary currents should be present, however weak, even in the mid-cross-section.

From a hydrodynamic instability point of view, a classical result indicates that, for a sufficiently shallow rectangular channel, the footprint of secondary currents on the spanwise profile of the streamwise (time-averaged) velocity can be interpreted as a turbulence-induced instability. The structure of such an instability in the $(y,z)$ cross-section was obtained by Ikeda (Reference Ikeda1981) and Nezu & Nakagawa (Reference Nezu and Nakagawa1993) after a term-by-term order of magnitude analysis of the time-averaged vorticity budget equation. By imposing a balance between the production of vorticity due to turbulence anisotropy and the suppression of vorticity due to Reynolds shear stress, assuming that the difference in the Reynolds normal stresses is linearly $y$-distributed, and the bed shear stress is laterally periodically perturbed, Ikeda (Reference Ikeda1981) came to a solution for the velocity field imposed by secondary currents which, by employing the dimensionless parameters introduced in the present work, reads as

(2.3a,b)\begin{equation} v_0={-}\omega\beta{{\cal{G}}(y)}\cos{(\beta z)},\quad w_0=\omega\frac{\textrm{d}{\cal{G}}(y)}{\textrm{d}y}\sin{(\beta z)}. \end{equation}

An approximated derivation of (2.3a,b) from (2.1), and a general solution for the distribution of the vertical velocity ${\cal {G}}$ is detailed in Appendix A. The parameter $\omega$ is physically dependent on the Reynolds number and the channel aspect ratio (i.e. the ratio between channel width and depth), that accounts for circulation of the cell. The transverse wavenumber $\beta =2{\rm \pi} /L_z$ accounts for the lateral spacing ($L_z$) of secondary currents. According to Townsend (Reference Townsend1976), in the absence of longitudinal bed forms such as ridges, the dominant cell refers to the case $\beta ={\rm \pi}$, whereby the cell width and the depth are identical ($L_z/2=1$). A cross-stream view of the secondary current structure is shown in figure 1$(a)$.

In the results section, it will be shown that $\omega$ is of paramount importance since it accounts for the intensity of secondary cells. In fact, by depth averaging (2.3a,b) and computing the circulation along a curve ${\cal {L}}$ lying on the plane $y-z$ and embracing a whole cell, we readily obtain that

(2.4a,b)\begin{equation} \omega=\frac{1}{2{\rm \pi}}\left|\oint_{\cal{L}} (v_0,w_0)\boldsymbol{\cdot}\,\textrm{d}\kern0.06em {\boldsymbol x}\right|=\frac{\varOmega_x}{2{\rm \pi}},\quad\text{and}\quad\omega=\left|\frac{{\rm \pi}\langle{v_0}\rangle}{4 - {\rm \pi}^2}\right|_{z=0,\beta={\rm \pi}}, \end{equation}

where brackets indicate the depth-averaging operator. Thus $\omega$ is proportional to the circulation of the velocity field around a curve embracing a whole cell. From the Stokes theorem, this can be interpreted as the value of the streamwise vorticity averaged over the surface covered by a cell ($\varOmega _x$) and hence a proxy for the intensity of secondary currents. It is noteworthy that the above result is independent of the shape of the primary velocity distribution $\cal {F}(y)$. Equation (2.4a,b) also shows that this value can be easily obtained by depth averaging the (non-dimensional) vertical velocity component $v_0$ from experimental data.

2.3. Depth-averaged equations

In order to derive the equations describing the dynamics of the problem at hand, we adopt the well-known shallow-water approximation (Tan Reference Tan1992), whereby $\hat {h}_0$ is much smaller than the dominant horizontal length scales. Moreover, we assume that the longitudinal velocity components $u$ can be described by the factorization $u=\cal {F}(y)U(x,z,t)$, whereas the lateral velocity component is decomposed into a component related to the secondary current, that is with zero depth averaging, plus a vertically uniform spanwise component, namely $w=w_0(y,z)+W(x,z,t)$. Therefore, $U$ and $W$ are the depth-averaged longitudinal and lateral velocities, respectively. Under these hypotheses and upon depth averaging, (2.1) reduces to

(2.5)$$\begin{gather} \frac{\partial{\boldsymbol{U}}}{\partial{t}}+{\boldsymbol{U}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{ {\perp}}{{\boldsymbol{U}}}+\frac{\boldsymbol{\nabla}_{ {\perp}}{\zeta}}{F^2}-\frac{\boldsymbol{\nabla}_{ {\perp}} \boldsymbol{\cdot} ({h}{\boldsymbol{T}})}{h}-\frac{\boldsymbol{\tau}_b}{h}=0, \end{gather}$$

(2.6)$$\begin{gather}\frac{\partial h}{\partial t}+\boldsymbol{\nabla}_{ {\perp}} \boldsymbol{\cdot} (h\boldsymbol{U})=0, \end{gather}$$

where ${\boldsymbol {U}}(x,z,t)=\{U,W\}$ is the depth-averaged Reynolds-averaged flow field, $\boldsymbol {\nabla }_{ \perp }=\{\partial _x,\partial _z\}$, $\zeta =h - x\sin {\theta }$ is the free surface elevation, $h(x,z,t)$ is the water depth, $\boldsymbol {\tau }_b$ is the bottom shear stress and ${\boldsymbol {T}}$ is a $(2\times 2)$ tensor resulting from the sum of the depth-averaged Reynolds stress tensor and the dispersive tensor arising from depth averaging; $\boldsymbol {\tau }_b$ and ${\boldsymbol {T}}$ can be modelled by means of the Chezy and Boussinesq formulations as

(2.7)$$\begin{gather} \boldsymbol{\tau}_b={-}c_f{\boldsymbol{U}}|\boldsymbol{U}|, \end{gather}$$

(2.8)$$\begin{gather}\boldsymbol{T}=\langle\nu_t(\boldsymbol{\nabla}_{{\perp}}{\boldsymbol{u}}_{{\perp}}+\boldsymbol{\nabla}_{{\perp}}{\boldsymbol{u}}_{{\perp}}^T)\rangle-\langle(\boldsymbol{u}_{{\perp}}-\boldsymbol{U})\otimes(\boldsymbol{u}_{{\perp}}-\boldsymbol{U})\rangle, \end{gather}$$

where $c_f=\sin {\theta }/F^2$ is the friction factor, which depends on the relative bed roughness $d\boldsymbol{\cdot }{h}^{-1}$ and the Reynolds number $Re$, ${\boldsymbol {u}}_{\perp }=\{u,w\}$, $\nu _t=\sqrt {c_{f}}hU{{\cal {N}}(y)}$ is the eddy viscosity, brackets refer to depth averaging and $\otimes$ refers to the dyadic product. The shape of function ${\cal {N}}(y)$ that defines the eddy viscosity $\nu _t$ will be derived in § 4. The first term on the right-hand side of (2.8) accounts for the depth averaging of the Reynolds stresses, while the second one produces the so-called dispersive terms that arise from departure of plug flow conditions.

Two remarks are in order. Firstly, the vector equation (2.5) states the momentum conservation along the longitudinal and transverse directions and it was derived assuming that the vertical momentum conservation reduces to the hydro-static balance $p \sim {F}^{-2}(h - y)$. This condition is the hallmark of the shallow-water equations and it can be shown from a dimensional analysis of the momentum conservation along the vertical direction, that it is constrained to the condition that the channel aspect ratio is much larger than the Froude number. Equation (2.6) states mass continuity after the use of the kinematic condition at the free surface (namely $\partial _{t}{h}=w - u_{\perp } \boldsymbol {\cdot } \boldsymbol {\nabla }_{ \perp }h$). Secondly, depth averaging of the longitudinal velocity provides $\langle {\cal {F}}\rangle =1$ by construction while $\langle {w}_0\rangle \sim \langle {\textrm {d}\cal {G}}/{\textrm {d}y}\rangle \sim 0$, for $\beta \sim {\rm \pi}$. From the above remarks, we obtain that (2.8) reads as

(2.9) \begin{equation} \boldsymbol{T}=\sqrt{c_{f}}\,U\left( \begin{array}{@{}cc@{}} \dfrac{1-\cal{I}_4}{\sqrt{c_{f}}}U+2h{\cal{I}}_0\dfrac{\partial U}{\partial x} & h{\cal{I}}_0\dfrac{\partial U}{\partial y}+h{\cal{I}}_1\dfrac{\partial W}{\partial x}-\dfrac{{\cal{I}}_2 \omega\sin (\beta z)}{\sqrt{c_{f}}} \\ h{\cal{I}}_0\dfrac{\partial U}{\partial y}+h{\cal{I}}_1\dfrac{\partial W}{\partial x}-\dfrac{{\cal{I}}_2 \omega\sin (\beta z)}{\sqrt{c_{f}}} & 2h{\cal{I}}_1\dfrac{\partial W}{\partial z}+2h{\cal{I}}_3 \beta \omega \cos (\beta z)-\dfrac{{\cal{I}}_5\omega^2}{U\sqrt{c_{f}}}\left[\sin{\left(\beta{z}\right)}\right]^2 \end{array} \right), \end{equation}

with

(2.10a–f) \begin{align} {\cal{I}}_0=\langle{\cal{F}}{\cal{N}}\rangle,\enspace {\cal{I}}_1=\langle{\cal{N}}\rangle,\enspace {\cal{I}}_2=\langle{\cal{F}}{\cal{G}'}\rangle,\enspace {\cal{I}}_3 = \langle{\cal{N}}{\cal{G}'}\rangle \sim 0,\enspace {\cal{I}}_4=\langle{\cal{F}}^2\rangle,\enspace {\cal{I}}_5=\langle{\cal{G}'}^2\rangle, \end{align}

whose exact formulation was obtained upon making a specific choice for $\cal {F}$ (see § 4) and is reported in Appendix B.

At the base state ($\partial _t=\partial _x=0$), after combining (2.3a,b) with (2.5)–(2.9), and recalling that $\omega \ll {1}$, it appears that the effect of the secondary currents on the streamwise and spanwise momentum conservation is a $z$-dependent modulation to the longitudinal velocity and depth, respectively, $u_0=\cal {F}(y)U_0(z)$, and $h_0(z)$ with

(2.11)$$\begin{gather} U_0=1-\omega\varPhi_{1}\cos{(\beta z)}-\omega^{2}\left[\tfrac{1}{4}{\varPhi_{1}^2}-\varPhi_{2}\cos{(2\beta{z})}\right]+{O}(\omega^4), \end{gather}$$

(2.12)$$\begin{gather}h_0=1+\tfrac{1}{2}\omega^2F^2{\cal{I}}_5\cos{(2\beta{z})}+{O}(\omega^4), \end{gather}$$

where

(2.13a,b) \begin{align} \varPhi_{1} = \frac{\beta{\cal{I}}_2}{\sqrt{c_{f}} \left(2 \sqrt{c_{f}}+\beta ^2 {\cal{I}}_{0}\right)},\quad \varPhi_{2} = \frac{F^{2}{\cal{I}}_{5} (c_{f} - c_{f,h}) \left(2 \sqrt{c_{f}} + \beta ^2 {\cal{I}}_{0}\right)^2 + 3\beta ^2 {\cal{I}}_{2}^2}{4\sqrt{c_{f}}\left(2 \sqrt{c_{f}} - \beta ^2 {\cal{I}}_{0}\right)^2 \left(\sqrt{c_{f}} + 2 \beta ^2 {\cal{I}}_{0}\right)}, \end{align}

with $c_{f,h}={\partial {c_f}}/{\partial {h}}|_{h=1}$. Equation (2.11) is crucial for the following analysis since it will be shown that such a transverse modulation, combined with the free surface response, induces an instability of the inflectional type to the secondary currents. It is straightforward to show that the terms of order ${O}({\omega ^3})$ are zero in the above expansions. It is also worth noticing that (2.11) is reminiscent of the spanwise modulation derived by Waleffe (Reference Waleffe1995), for the case of streamwise rolls superimposed on a laminar Couette flow. In the following analysis we will neglect terms of order $\omega ^2$ or higher, since they complicate the algebraic treatment without any numerically relevant effect (see Appendix E).

3.1. Perturbative solution

In order to account for large-scale longitudinal meandering of the secondary cell axis, let us consider a time-dependent normal-mode perturbation, by adopting the following putative solution to (2.5)–(2.6):

(3.1)\begin{equation} \{U,W,h\}=\{U_0,0,h_0\}+\epsilon\{U_1(z),W_1(z),h_1(z)\}\exp({{\textrm{i}}\alpha{x}+\lambda t})+\textrm{c.c.}, \end{equation}

where the base state $U_0$ and $h_0$ are provided by (2.11) and (2.12) truncated at order $\omega$, $\epsilon$ is a small perturbative parameter, $\alpha$ is the longitudinal wavenumber of the perturbation, $\lambda$ is its growth rate and $\textrm {c.c.}$ refers to complex conjugation.

It is important to remark that under non-uniform and unsteady conditions the friction factor depends on the local value of the depth, through the ratio $d \boldsymbol{\cdot} {h}^{-1}(x,z,t)$. This aspect is accounted for through a Taylor expansion as follows:

(3.2)\begin{equation} c_f(h)\sim\left.{c}_{f}\right|_{h=1}+\left.\frac{\partial{c_f}}{\partial{h}}\right|_{h=1}(h-1)=c_{f0}+\epsilon\left.\frac{\partial{c_f}}{\partial{h}}\right|_{h=1}h_1 \textrm{e}^{\textrm{i}\alpha{x}+\lambda{t}}+\textrm{c.c.}, \end{equation}

where $c_{f0}={c}_{f}|_{h=1}$ is the friction factor computed for the unperturbed flow. After substituting the above ansatz in (2.5)–(2.6) and taking into account (2.7)–(2.13a,b), at the order $\epsilon$ we obtain the eigenvalue problem

(3.3)\begin{equation} u_*(\boldsymbol{A}+\omega\boldsymbol{B})\boldsymbol{q}={\lambda}\boldsymbol{q}, \end{equation}

where $\boldsymbol {q}=\{U_1,W_1,h_1\}^T$, while $u_*$:=$\sqrt {c_{f0}}$ is the dimensionless friction velocity, and

(3.4)\begin{gather} \boldsymbol{A}=\left( \begin{array}{@{}ccc@{}} {\cal{I}}_0{\cal{D}}^2 - \sigma_0 & \textrm{i}\,{\cal{I}}_1\alpha{\cal{D}} & \sigma_1 \\ \textrm{i}\cal{I}_0\alpha\cal{D} & 2\cal{I}_1\cal{D}^2-\sigma_4 & -\dfrac{\cal{D}}{u_*F^2}\\ -\dfrac{\textrm{i}\alpha}{u_*} & -\dfrac{\cal{D}}{u_*} & -\dfrac{\textrm{i}\alpha}{u_*} \end{array} \right), \end{gather}

(3.5)\begin{gather} \boldsymbol{B}=\left( \begin{array}{@{}ccc@{}} \sigma_2\cal{C}+\sigma_9 \cal{S}\cal{D} -\dfrac{\textrm{i}\sigma_6}{\alpha }\cal{C}\cal{D}^2 & -\dfrac{\textrm{i}\alpha\beta\sigma_8}{2}\cal{S} -\dfrac{ \sigma_7\beta}{u_*}\cal{S} + \dfrac{\textrm{i}\alpha\sigma_8}{2} \cal{C}\cal{D} & \sigma_3\cal{C} + \sigma_9 \cal{S}\cal{D} \\ \sigma_5 \cal{S} + \sigma_6\cal{C}\cal{D} & \sigma_7 \sigma_4 \cal{C} + \sigma_8(\cal{C}\cal{D} - \beta\cal{S})\cal{D} & \textrm{i}\,\alpha \sigma _9 \cal{S}\\ 0 & 0 & \dfrac{\textrm{i}\,\alpha\sigma_7}{u_*}\cal{C} \end{array} \right). \end{gather}

In these matrices $\cal {D}^n=\text {d}^n/\text {d}z^n$, $\cal {C}=\cal {C}(z)=\cos {(\beta {z})}$, $\cal {S}=\cal {S}(z)=\sin {(\beta {z})}$, while the coefficients $\sigma _i$ are reported in Appendix C. The presence of the functions $\cal {C}$ and $\cal {S}$ in $\boldsymbol {B}$ precludes an analytical solution of (3.3), unless $\omega =0$.

3.2. Spectral solution of the eigenvalue problem

We firstly consider the mapping $z \rightarrow {\xi }=\beta {z}$, so that the operators $\boldsymbol {A}$ and $\boldsymbol {B}$ become periodic in the interval $\xi =[0,2{\rm \pi} ]$. The differential eigenvalue problem (3.3) may be solved by a Fourier Galerkin spectral method (see e.g. Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2006). Accordingly, a modal representation of the approximate solution with trigonometric polynomials is adopted

(3.6)\begin{equation} \boldsymbol{q}=\sum_{k={-}N/2}^{k=N/2}\hat{\boldsymbol{q}}_k \textrm{e}^{{\text{i}}k\xi}, \end{equation}

where the cutoff parameter $N$ is an even arbitrarily large integer.

A weak formulation of the problem is obtained by selecting a set of test functions $\boldsymbol {\psi }_{l,n}$ that determine the weights of the residual

(3.7)\begin{equation} \int_0^{2{\rm \pi}}[u_*(\boldsymbol{A}+\omega\boldsymbol{B})-\lambda\boldsymbol{I}]\boldsymbol{q}{\boldsymbol{\psi}}_{l,n}\,\textrm{d}\xi=0,\quad (n={-}N/2,\ldots,N/2;\quad l=1,\ldots,6). \end{equation}

A natural choice of the test functions is $\boldsymbol {\psi }_{l,n}=(2{\rm \pi} )^{-1}\boldsymbol {v}_l{e}^{-{\text {i}}n\xi }$, where $\boldsymbol {v}_l$ is any vector of the canonical basis in $\mathbb {C}^3$ (i.e. one component is equal to 1 or $\text {i}$, all the others $0$).

In order to reduce (3.7) to an algebraic eigenvalue problem, let us first observe that, by expressing $\cos {\xi }$ and $\sin {\xi }$ in terms of exponentials $\textrm {e}^{\pm {\text {i}}\xi }$, one can decompose the matrix $\boldsymbol {B}$ in (3.7) as

(3.8)\begin{equation} \boldsymbol{B}=\textrm{e}^{-\text{i}\xi}\boldsymbol{B}_{{-}1}+\textrm{e}^{\text{i}\xi}\boldsymbol{B}_{1}=\sum_{j={\pm}{1}}\textrm{e}^{\text{i}j\xi}\boldsymbol{B}_j, \end{equation}

where $\boldsymbol {B}_j$ are the matrices with constant coefficients with respect to $\xi$, similarly to the matrix $\boldsymbol {A}$ in (3.4). Hence,

(3.9)\begin{align} [u_*(\boldsymbol{A}+\omega\boldsymbol{B}) - \lambda\boldsymbol{I}]\boldsymbol{q}&=\sum_{k={-}N/2}^{N/2}\left[u_*\left(\boldsymbol{A}+ \omega\sum_{j={\pm}{1}} \textrm{e}^{\text{i}j\xi}\boldsymbol{B}_j\right)-\lambda\boldsymbol{I}\right] \textrm{e}^{\text{i}k\xi}\hat{\boldsymbol{q}}_k\nonumber\\ &=\sum_{k={-}N/2}^{N/2}\left[u_*\left(\hat{\boldsymbol{A}}_k+ \omega\sum_{j={\pm}{1}} \textrm{e}^{\text{i}j\xi}\hat{\boldsymbol{B}}_{j,k}\right)-\lambda\boldsymbol{I}\right] \textrm{e}^{\text{i}k\xi}\hat{\boldsymbol{q}}_k\nonumber\\ &=\sum_{k={-}N/2}^{N/2}\left[u_*\left(\hat{\boldsymbol{A}}_k\hat{\boldsymbol{q}}_k \textrm{e}^{\text{i}k\xi}+ \omega\sum_{j={\pm}{1}}\hat{\boldsymbol{B}}_{j,k}\hat{\boldsymbol{q}}_k \textrm{e}^{\text{i}(k+j)\xi} \right)-\lambda\hat{\boldsymbol{q}}_k \textrm{e}^{\text{i}k\xi}\right], \end{align}

where the matrices $\hat {\boldsymbol {A}}_k$, $\hat {\boldsymbol {B}}_{j,k}$ are obtained from ${\boldsymbol {A}}$, ${\boldsymbol {B}}_{j}$ by replacing the derivative operator $\cal {D}=\text {d}/\text {d}\xi$ by $\text {i}k$ (Fourier differentiation). The entries of these matrices are reported in Appendix D.

After substituting (3.9) in (3.7), one easily obtains the algebraic eigenvalue problem of order $3(N+1)$

(3.10)\begin{equation} u_*\left(\hat{\boldsymbol{A}}_{n}\hat{\boldsymbol{q}}_n+\omega\sum_{j={\pm}{1}}\hat{\boldsymbol{B}}_{j,n-j}\hat{\boldsymbol{q}}_{n-j}\right)={\lambda}\hat{\boldsymbol{q}}_n, \quad (n={-}N/2,\ldots,N/2). \end{equation}

Note that the matrix of order $3(N + 1)$ corresponding to the left-hand side is block tridiagonal, with diagonal blocks given by $u_*\hat {\boldsymbol {A}}_{n}$ and off-diagonal blocks given by $u_*\omega \hat {\boldsymbol {B}}_{1,n-1}$ (left), and $u_*\omega \hat {\boldsymbol {B}}_{-1,n+1}$ (right). A ready-to-use set of Matlab scripts for the solution of (3.10) can be found in the supplementary material available at https://doi.org/10.1017/jfm.2021.769.