2 Internal gravity wave–mean flow interactions with zonal symmetry

We consider a fluid within a two-dimensional domain, periodic in the zonal $x$ -direction with period $L$ and semi-infinite in the vertical $z$ -direction. The bottom boundary is a vertically undulating line located on average at $z=0$ . The fluid is considered incompressible, Boussinesq, viscous with viscosity $\unicode[STIX]{x1D708}$ and linearly stratified with buoyancy frequency $N$ . For the sake of simplicity, we ignore any buoyancy diffusion process. This approximation is relevant for experimental configurations where the stratification agent is salt, given the low diffusivity $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D708}/1000$ , but it does not apply to the atmosphere and the ocean, where turbulent viscosity and diffusivity have the same order of magnitude.

Throughout this work, we solely consider monochromatic waves. Let us introduce the typical zonal wavenumber $k=2\unicode[STIX]{x03C0}/L$ , angular frequency $\unicode[STIX]{x1D714}$ and amplitude of the bottom undulation $h_{b}$ . There are three independent dimensionless numbers in the problem. The Froude number $Fr=\unicode[STIX]{x1D714}/N$ controls the angle of propagation of the wave. The wave Reynolds number $Re=\unicode[STIX]{x1D714}/(k^{2}\unicode[STIX]{x1D708})$ controls the viscous damping and the viscous boundary layer thickness of the wave field. When considering the lee wave generation case, this wave Reynolds number scales as $UL/\unicode[STIX]{x1D708}$ , where $U$ is the typical mean zonal velocity. The third parameter is the dimensionless amplitude of the wave $\unicode[STIX]{x1D716}=h_{b}k$ , which corresponds to the typical slope of the bottom boundary and controls the linearity of the wave. In numerical simulations, an additional aspect ratio $r=kH$ and a wave Péclet number $\unicode[STIX]{x1D714}/(k^{2}\unicode[STIX]{x1D705})$ have to be taken into account, because the domain has a finite height $H$ , and because it includes a buoyancy diffusivity $\unicode[STIX]{x1D705}$ . Both parameters will be much larger than one in the numerical simulations presented in this paper, and we will assume that they do not play a significant role in this limit. We use $k^{-1}$ and $\unicode[STIX]{x1D714}^{-1}$ as reference length and time for the space–time coordinates, $c=\unicode[STIX]{x1D714}/k$ as a reference velocity, $N^{2}/k$ as a reference buoyancy, and write the dynamical equation in a dimensionless form

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+Fr^{-2}b\boldsymbol{e}_{z}+Re^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u},\\ \displaystyle \unicode[STIX]{x2202}_{t}b+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}b+w=0,\\ \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{array}\right\}\end{eqnarray}$$

where $\boldsymbol{u}=(u,w)$ is the two-dimensional velocity, $p$ the renormalised pressure, $b$ the buoyancy anomaly, $\boldsymbol{e}_{z}$ the unit vector of the vertical direction pointing upwards and $\unicode[STIX]{x1D6FB}^{2}=\unicode[STIX]{x2202}_{xx}+\unicode[STIX]{x2202}_{zz}$ the standard Laplacian operator.

Previous studies in the context of acoustic streaming have investigated the effect of changing boundary conditions on mean flow properties (Xie & Vanneste Reference Xie and Vanneste2014). In this paper devoted to internal wave streaming, we discuss two different bottom boundary conditions on $z=\unicode[STIX]{x1D716}h(x,t)$ , i.e.

(2.2a-c ) $$\begin{eqnarray}\text{free-slip:}\quad w=\unicode[STIX]{x1D716}(\unicode[STIX]{x2202}_{t}h+u\unicode[STIX]{x2202}_{x}h),\quad \unicode[STIX]{x1D642}[\boldsymbol{n}_{h}]\boldsymbol{\cdot }\boldsymbol{n}_{h}^{\bot }=0;\qquad \text{no-slip:}\quad \boldsymbol{u}=\unicode[STIX]{x1D716}\unicode[STIX]{x2202}_{t}h\boldsymbol{e}_{z},\end{eqnarray}$$

where $\boldsymbol{n}_{h}=\unicode[STIX]{x1D735}(z-\unicode[STIX]{x1D716}h(x,t))$ is a local normal vector of the bottom boundary, $\boldsymbol{n}_{h}^{\bot }$ a local tangent vector and $\unicode[STIX]{x1D642}$ the velocity gradient tensor ( $G_{ij}=\unicode[STIX]{x2202}_{j}u_{i}$ ). This free-slip condition is the one implemented in the numerical model considered in this paper (see the MITgcm’s User Manual for the MIT global circulation model (GCM) (Adcroft et al. Reference Adcroft, Campin, Doddridge, Dutkiewicz, Evangelinos, Ferreira, Follows, Forget, Fox-Kemper, Heimbach, Hill, Hill, Hill, Jahn, Klymak, Losch, Marshall, Maze, Mazloff, Menemenlis, Molod and Scott2018)). It is equivalent to the stress-free condition when boundary curvature can be neglected. In the stress-free case, $\unicode[STIX]{x1D642}$ is replaced by its symmetric part only. Regarding the boundary streaming, we confirmed that the discrepancies between stress-free and our free-slip condition arise only in non-hydrostatic regimes of internal waves. Therefore, in most practical cases, the results obtained by considering the free-slip condition (2.2) will also be relevant for numerical simulation using the stress-free condition. Finally, we require all gradients with respect to $z$ to vanish as $z\rightarrow \infty$ .

When considering a progressive pattern ( $h(x,t)=h(x-t)$ ) in (2.2), a Galilean change of reference yields the case of lee-wave generation by a depth-independent mean flow passing over bottom topography. Then, the free-slip bottom boundary condition for the generation of lee waves obviates the need to treat the near-bottom critical layer induced by a more realistic no-slip condition (Passaggia et al. Reference Passaggia, Meunier and Le Dizès2014). Regarding the free-slip condition, the predictions will be compared against direct numerical simulations of monochromatic lee-wave generation using the MIT GCM (Adcroft, Hill & Marshall Reference Adcroft, Hill and Marshall1997), which specifically uses our definition for the free-slip condition. The no-slip boundary condition in (2.2) is relevant to modelling the generation of internal gravity waves in laboratory experiments using vertically oscillating bottom membranes (Plumb & McEwan Reference Plumb and McEwan1978; Semin et al. Reference Semin, Facchini, Pétrélis and Fauve2016) or a system of plates and camshafts (Gostiaux et al. Reference Gostiaux, Didelle, Mercier and Dauxois2006). We will, however, consider limiting cases where the viscous boundary layer is larger than the boundary height variations, which is not always the case in actual experiments.

We decompose any field $\unicode[STIX]{x1D719}$ into a mean flow part $\overline{\unicode[STIX]{x1D719}}$ and a wave part $\unicode[STIX]{x1D719}^{\prime }$ using the zonal averaging procedure (see Bühler Reference Bühler2014):

(2.3) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D719}}(z,t)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}x\,\unicode[STIX]{x1D719}(x,z,t),\quad \unicode[STIX]{x1D719}^{\prime }=\unicode[STIX]{x1D719}-\overline{\unicode[STIX]{x1D719}}.\end{eqnarray}$$

The averaging of the zonal momentum equation in (2.1) leads to the mean flow evolution equation

(2.4) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\overline{u}=-\unicode[STIX]{x2202}_{z}\overline{u^{\prime }w^{\prime }}+Re^{-1}\unicode[STIX]{x2202}_{zz}\overline{u}.\end{eqnarray}$$

The source of streaming is the divergence of the Reynolds stress $-\unicode[STIX]{x2202}_{z}\overline{u^{\prime }w^{\prime }}$ . To compute this term, we subtract the averaged equations from (2.1) and linearise the result assuming $(u^{\prime },w^{\prime },b^{\prime },p^{\prime })=O(\unicode[STIX]{x1D716})$ with $\unicode[STIX]{x1D716}\ll 1$ . In the limit $\unicode[STIX]{x1D716}\ll 1$ , in the absence of dissipation in the buoyancy equation, it turns out that $\overline{b}$ scales with $\unicode[STIX]{x1D716}^{2}$ at all time (see Bühler Reference Bühler2014). In the following, we will ignore $\overline{b}$ and its evolution.

At this stage, we assume that $|\overline{u}|\ll 1$ . Starting from a state of rest, at early times of its evolution, the mean flow is weak, which justifies this assumption. At later times, the feedback of the mean flow on the wave can no longer be ignored (Kataoka & Akylas Reference Kataoka and Akylas2015; Fan, Kataoka & Akylas Reference Fan, Kataoka and Akylas2018), as will be discussed in more detail in § 4 (see also (A 1) in appendix A). This case without feedback from the mean flow leads to homogeneous wave equations, which provides a simple framework to describe essential features of boundary streaming:

(2.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x2202}_{t}u^{\prime }+\unicode[STIX]{x2202}_{x}p^{\prime }-Re^{-1}\unicode[STIX]{x1D6FB}^{2}u^{\prime }=0,\\ \displaystyle \unicode[STIX]{x2202}_{t}w^{\prime }+\unicode[STIX]{x2202}_{z}p^{\prime }-Fr^{-2}b^{\prime }-Re^{-1}\unicode[STIX]{x1D6FB}^{2}w^{\prime }=0,\\ \displaystyle \unicode[STIX]{x2202}_{t}b^{\prime }+w^{\prime }=0,\\ \displaystyle \unicode[STIX]{x2202}_{x}u^{\prime }+\unicode[STIX]{x2202}_{z}w^{\prime }=0.\end{array}\right\}\end{eqnarray}$$

The coupled equations (2.4) and (2.5) form a quasi-linear model for the interaction between boundary generated viscous waves and the zonal mean flow. The Reynolds stress divergence, $-\unicode[STIX]{x2202}_{z}\overline{u^{\prime }w^{\prime }}$ , at the origin of streaming is the only nonlinear term remaining in the problem. It acts as a forcing term and is computed from the wave field.

We perform the wave–mean decomposition on the boundary conditions (2.2) and we linearise the result assuming as above a wave amplitude of order $\unicode[STIX]{x1D716}$ on an asymptotically flat boundary at $z=0$ , i.e.

(2.6a,b ) $$\begin{eqnarray}\text{free-slip:}\quad \left\{\begin{array}{@{}l@{}}\displaystyle \unicode[STIX]{x2202}_{z}\overline{u}=0,\\ \displaystyle w^{\prime }-\unicode[STIX]{x2202}_{t}h=0,\\ \displaystyle \unicode[STIX]{x2202}_{z}u^{\prime }=0;\end{array}\right.\qquad \text{no-slip:}\quad \left\{\begin{array}{@{}l@{}}\displaystyle \overline{u}=0,\\ \displaystyle w^{\prime }-\unicode[STIX]{x2202}_{t}h=0,\\ \displaystyle u^{\prime }=0.\end{array}\right.\end{eqnarray}$$

In the free-slip case, the Reynolds stress divergence vanishes at the bottom ( $\unicode[STIX]{x2202}_{z}\overline{u^{\prime }w^{\prime }}|_{z=0}=0$ ), while, in the no-slip case, the Reynolds stress itself vanishes at the bottom ( $\overline{u^{\prime }w^{\prime }}|_{z=0}=0$ ). Given that $\overline{u^{\prime }w^{\prime }}|_{z\rightarrow \infty }=0$ for damped waves, the integrated streaming in the no-slip case has to be zero: $\int _{0}^{\infty }\unicode[STIX]{x2202}_{z}\overline{u^{\prime }w^{\prime }}\,\text{d}z=0$ . Consequently, all the streaming far from the bottom boundary has to be compensated for by an opposite boundary streaming.

3 From viscous waves to boundary streaming

3.1 Viscous internal gravity waves

We describe in this section the detailed structure of the Reynolds stress divergences for both the free-slip and the no-slip boundary conditions, when the mean flow can be neglected. Inserting the ansatz $(u^{\prime },w^{\prime },b^{\prime },p^{\prime })=\text{Re}[(\tilde{u} ,\tilde{w},\tilde{b},\tilde{p})\text{e}^{\text{i}(x+mz-t)}]$ into (2.5) leads to the dispersion relation for viscous internal gravity waves, expressed here as

(3.1) $$\begin{eqnarray}m^{2}=\frac{\text{i}Re}{2}\left(1\pm \sqrt{1+\frac{4\text{i}}{Fr^{2}Re}}\right)-1.\end{eqnarray}$$

Among the four possible solutions for $m$ , we retain only the two upward-propagating ones, by discarding the solutions with a negative imaginary part. To simplify the discussion, it will be useful to express these solutions in the asymptotic regime $Fr^{2}Re\gg 1$ , followed by $Fr\ll 1$ :

(3.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle m_{w}=-1/Fr+\text{i}/(2L_{Re})+O((ReFr^{3})^{-1}),\\ \displaystyle m_{bl}=(1+\text{i})/\unicode[STIX]{x1D6FF}_{Re}+O(Re^{1/2}),\end{array}\right\}\end{eqnarray}$$

with

(3.3a,b ) $$\begin{eqnarray}L_{Re}=ReFr^{3}\quad \text{and}\quad \unicode[STIX]{x1D6FF}_{Re}=\sqrt{2/Re}.\end{eqnarray}$$

The solution $m_{w}$ corresponds to the propagating solution converging towards the inviscid solution in the limit $Re\rightarrow \infty$ . Here $L_{Re}$ is the damping length scale of the wave beam, scaling linearly with the wave Reynolds number. The solution $m_{bl}$ corresponds to the wave boundary layer. The boundary layer thickness, given by $\unicode[STIX]{x1D6FF}_{Re}$ , scales as $Re^{-1/2}$ as in the classical case of a horizontally oscillating flat boundary. This last solution is needed to match the propagating solution with the viscous boundary conditions and is analogous to the one discussed in acoustic boundary streaming (Nyborg Reference Nyborg1958). Importantly the ratio $L_{Re}/\unicode[STIX]{x1D6FF}_{Re}$ diverges in the limit $Fr^{2}Re\rightarrow +\infty$ . This limit, therefore, allows for a clear separation between bulk and boundary effects.

The viscous internal wave dispersion relation has already been extensively studied. Chini & Leibovich (Reference Chini and Leibovich2003) considered a finite Prandtl number, which provides an additional branch of boundary layer solutions associated with the diffusion operator in the buoyancy equation. They also gave asymptotic expansions for large Reynolds number. Grisouard & Thomas (Reference Grisouard and Thomas2016) considered the effect of a Coriolis force, which also gives rise to an additional branch of boundary layer solutions. Although rotation, buoyancy diffusion and their associated boundary layer solutions undoubtedly impact boundary streaming, we do not consider these additional effects, to simplify the presentation.

In the case of a progressive sine-shaped bottom undulation, $h(x,t)=\text{Re}[\text{e}^{\text{i}(x-t)}]$ , the general expression of the wave field is given by the linear combination of a propagating ( $w$ ) and a boundary layer ( $bl$ ) part:

(3.4) $$\begin{eqnarray}[u^{\prime },w^{\prime },b^{\prime },p^{\prime }]=\text{Re}\{(\unicode[STIX]{x1D719}_{w}\boldsymbol{P}[m_{w}]\text{e}^{\text{i}m_{w}z}+\unicode[STIX]{x1D719}_{bl}\boldsymbol{P}[m_{bl}]\text{e}^{\text{i}m_{bl}z})\text{e}^{\text{i}(x-t)}\},\end{eqnarray}$$

with

(3.5) $$\begin{eqnarray}\boldsymbol{P}[m]=[1,-m^{-1},\text{i}Fr^{-2}m^{-1},Fr^{-2}(1+m^{2})^{-1}].\end{eqnarray}$$

Here $\boldsymbol{P}[m]$ is the polarisation of the wave obtained from (2.5), $(m_{w},m_{bl})$ are given in (3.2) and $(\unicode[STIX]{x1D719}_{w},\unicode[STIX]{x1D719}_{bl})$ are scalars determined by the boundary conditions (2.6), i.e.

(3.6a,b ) $$\begin{eqnarray}\text{free-slip:}\quad \left\{\begin{array}{@{}l@{}}\displaystyle \unicode[STIX]{x1D719}_{w}=\text{i}\unicode[STIX]{x1D716}\frac{m_{w}m_{bl}^{2}}{m_{bl}^{2}-m_{w}^{2}},\\ \displaystyle \unicode[STIX]{x1D719}_{bl}=\text{i}\unicode[STIX]{x1D716}\frac{m_{bl}m_{w}^{2}}{m_{w}^{2}-m_{bl}^{2}};\end{array}\right.\qquad \text{no-slip:}\quad \left\{\begin{array}{@{}l@{}}\displaystyle \unicode[STIX]{x1D719}_{w}=\text{i}\unicode[STIX]{x1D716}\frac{m_{w}m_{bl}}{m_{bl}-m_{w}},\\ \displaystyle \unicode[STIX]{x1D719}_{bl}=\text{i}\unicode[STIX]{x1D716}\frac{m_{w}m_{bl}}{m_{w}-m_{bl}}.\end{array}\right.\end{eqnarray}$$

The generic vertical profiles of the wave field $u^{\prime }$ are drawn in figure 1 for both boundary conditions. Most of the differences between the two profiles are located in the boundary layer close to the bottom. We will see that these different profiles lead to very different boundary streaming behaviours, by computing the Reynolds stress divergence of the corresponding wave fields.

Figure 1. (a) Example of a linear computation of the vertical profile of the fully established wave field, $u^{\prime }$ , in the absence of a mean flow, with the free-slip (blue) and no-slip (red) boundary conditions. (b) Zoom on the boundary layer of the wave. The wave damping length, $L_{Re}$ , and the boundary layer thickness, $\unicode[STIX]{x1D6FF}_{Re}$ , are represented on the graph along with the inviscid vertical wavelength, $\unicode[STIX]{x1D706}_{z}=2\unicode[STIX]{x03C0}Fr$ .

3.2 Reynolds stress divergence

The Reynolds stress $\overline{u^{\prime }w^{\prime }}$ is composed of cross-terms involving both the propagative and the boundary layer contributions. In the limit of small viscosity, the ‘self-interaction’ of the propagating contribution decreases exponentially over a scale $L_{Re}$ . This effect is responsible for bulk streaming. All the other terms involve a pairing with the boundary layer contribution that decays exponentially over the scale $\unicode[STIX]{x1D6FF}_{Re}$ . The sum of these terms induces the boundary streaming. We thus decompose the Reynolds stress into a bulk and a boundary term, i.e.

(3.7) $$\begin{eqnarray}\overline{u^{\prime }w^{\prime }}(z)=F_{w}(z)+F_{bl}(z).\end{eqnarray}$$

In the remainder of this section, the quasi-linear computations will be performed by using the exact solutions of (3.1). In order to get insights on the basic differences between the free-slip and the no-slip cases, it is useful, however, to estimate the Reynolds stress by using the asymptotic expression (3.2) for both boundary conditions in (2.6), i.e.

(3.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\text{free-slip:}\quad \left\{\begin{array}{@{}l@{}}\displaystyle F_{w}(z)=\frac{\unicode[STIX]{x1D716}^{2}}{2Fr}\exp \left\{-\frac{z}{L_{Re}}\right\},\\ \displaystyle F_{bl}(z)=\frac{\unicode[STIX]{x1D716}^{2}}{Fr^{2}2\sqrt{2Re}}\exp \left\{-\frac{z}{\unicode[STIX]{x1D6FF}_{Re}}\right\}\left(\sin \frac{z}{\unicode[STIX]{x1D6FF}_{Re}}+\cos \frac{z}{\unicode[STIX]{x1D6FF}_{Re}}\right);\end{array}\right.\\ \text{no-slip:}\quad \left\{\begin{array}{@{}l@{}}\displaystyle F_{w}(z)=\frac{\unicode[STIX]{x1D716}^{2}}{2Fr}\exp \left\{-\frac{z}{L_{Re}}\right\},\\ \displaystyle F_{bl}(z)=-\frac{\unicode[STIX]{x1D716}^{2}}{2Fr}\exp \left\{-\frac{z}{\unicode[STIX]{x1D6FF}_{Re}}\right\}\cos \frac{z}{\unicode[STIX]{x1D6FF}_{Re}}.\end{array}\right.\end{array}\!\!\right\}\end{eqnarray}$$

The bulk Reynolds stress $F_{w}$ is the same at leading order for both the free-slip and the no-slip cases. The difference lies in the boundary-driven Reynolds stress $F_{bl}$ . The corresponding asymptotic expressions for the streaming body forces are

(3.9) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\text{free-slip:}\quad \left\{\begin{array}{@{}l@{}}\displaystyle -\unicode[STIX]{x2202}_{z}F_{w}(z)=\frac{\unicode[STIX]{x1D716}^{2}}{2Fr^{4}Re}\exp \left\{-\frac{z}{L_{Re}}\right\},\\ \displaystyle -\unicode[STIX]{x2202}_{z}F_{bl}(z)=\frac{\unicode[STIX]{x1D716}^{2}}{2Fr^{2}}\exp \left\{-\frac{z}{\unicode[STIX]{x1D6FF}_{Re}}\right\}\sin \frac{z}{\unicode[STIX]{x1D6FF}_{Re}};\end{array}\right.\\ \displaystyle \text{no-slip:}\quad \left\{\begin{array}{@{}l@{}}\displaystyle -\unicode[STIX]{x2202}_{z}F_{w}(z)=\frac{\unicode[STIX]{x1D716}^{2}}{2Fr^{4}Re}\exp \left\{-\frac{z}{L_{Re}}\right\},\\ \displaystyle -\unicode[STIX]{x2202}_{z}F_{bl}(z)=-\frac{\unicode[STIX]{x1D716}^{2}\sqrt{Re}}{2Fr\sqrt{2}}\exp \left\{-\frac{z}{\unicode[STIX]{x1D6FF}_{Re}}\right\}\left(\cos \frac{z}{\unicode[STIX]{x1D6FF}_{Re}}+\sin \frac{z}{\unicode[STIX]{x1D6FF}_{Re}}\right).\end{array}\right.\end{array}\!\!\right\}\end{eqnarray}$$

In the free-slip case, the boundary forcing amplitude does not depend on the wave Reynolds number at leading order, only its e-folding height does. This amplitude decreases with the Froude number. This effect can be seen in figure 2, where the free-slip Reynolds stress divergence $\unicode[STIX]{x2202}_{z}\overline{u^{\prime }w^{\prime }}$ is plotted for three different values of Reynolds and Froude numbers. These quasi-linear calculations are successfully compared to high-resolution direct numerical simulations (DNS) of the established wave pattern generated by a depth-independent flow above a sine-shaped topography in a linearly stratified fluid.

Figure 2. Plot of the vertical profile of the Reynolds stress divergence in the absence of mean flow ( $\overline{u}=0$ ) considering the free-slip boundary condition for different pairs $(Re,Fr)$ . The symbols come from high-resolution DNS while the dashed lines come from the full linear theory without mean flow. The other dimensionless parameters for the simulation are $\unicode[STIX]{x1D716}=0.01$ (wave amplitude) and $r=6L_{Re}$ (domain aspect ratio); the resolution is $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}z=\unicode[STIX]{x1D6FF}_{Re}/50$ ; the grid is stretched above $z=6L_{Re}$ to avoid wave reflection; the simulated data have been smoothed over 10 time steps of the simulation to get rid of the fast motion coming from surface waves present in the numerical model.

Figure 3. (a) Plot of the vertical profile of the Reynolds stress divergence for the no-slip boundary condition computed using the full linear theory without mean flow. (b) Plot of the vertical profile of the mean flow at $t=10$ computed using the quasi-linear model for the no-slip boundary condition. (c) Hovmöller diagram of the mean flow, $\overline{u}(z,t)$ , computed using the quasi-linear theory for the scenario in which the lower boundary condition is no-slip. The parameters are $Re=200$ , $Fr=0.3$ and $\unicode[STIX]{x1D716}=0.005$ .

In the no-slip case, boundary forcing is opposite to (and much stronger than) the bulk forcing, as shown in figure 3(a). The underlying reason is the vanishing of the integral of the Reynolds stress divergence over the whole domain, as discussed at the end of § 2. According to (3.9), the amplitude of the boundary forcing evaluated at the bottom scales as $\unicode[STIX]{x1D716}^{2}Re^{1/2}/Fr$ . In the limit $ReFr^{2}\gg 1$ , this amplitude is much larger than in the free-slip case. In addition, it increases with the Reynolds number. However, we will see in § 3.3 that the amplitude does not blow up in a distinguished limit that is consistent with the linearisation of the equations.

3.3 Boundary flows

We now look for the mean flow response to the Reynolds stress divergences, by inserting the linear predictions for wave fields into (2.4). When ignoring the influence of the mean flow on the wave fields, (2.4) becomes a linear diffusion equation with a steady forcing, which can be decomposed into a bulk and a boundary contribution, as in (3.7).

The typical time scales $\unicode[STIX]{x1D70F}_{w}$ and $\unicode[STIX]{x1D70F}_{bl}$ for the mean flow to reach a given velocity $U$ in the presence of either bulk or boundary streaming forcing terms are obtained by balancing $\unicode[STIX]{x2202}_{t}\overline{u}$ with $\unicode[STIX]{x2202}_{z}F_{w}$ and $\unicode[STIX]{x2202}_{z}F_{bl}$ , respectively. Using the large-Reynolds-number asymptotic estimates given in (3.9) leads then to $\unicode[STIX]{x1D70F}_{bl}/\unicode[STIX]{x1D70F}_{w}\sim 1/(Fr^{2}Re)$ in the free-slip case and $\unicode[STIX]{x1D70F}_{bl}/\unicode[STIX]{x1D70F}_{w}\sim 1/(Fr^{2}Re)^{3/2}$ in the no-slip case. We thus expect the boundary streaming to dominate over the bulk streaming at the early stage of the mean flow evolution in both cases.

At a quasi-linear level, the early stage of the mean flow evolution is obtained for both the free-slip and the no-slip conditions by solving (2.4) numerically, assuming that the wave field is described by (3.1), (3.4), (3.5) and (3.6). A finite size domain is considered in the simulations with an aspect ratio $r=6L_{Re}$ . The waves are computed as if the domain were semi-infinite and a free-slip upper boundary condition is considered for the mean flow.

Figure 4. Hovmöller diagrams of the mean flow, $\overline{u}(z,t)$ , for the scenario in which the bottom boundary condition is free-slip: (a) DNS, (b) quasi-linear model and (c) quasi-linear model without the boundary streaming terms in the Reynolds stress divergence. The parameters are $Re=200$ , $Fr=0.3$ and $\unicode[STIX]{x1D716}=0.01$ , with $\text{d}x=\text{d}z=\unicode[STIX]{x1D6FF}_{Re}/15$ . The grid is exponentially stretched on the vertical axis above $z=6L_{Re}$ in the DNS. At larger time, around $t\sim 300$ , the mean flow induced by bulk streaming becomes larger than the mean flow induced by boundary streaming.

In figure 4, we compare the quasi-linear predictions for the free-slip boundary condition against direct numerical simulations. The parameters are $Re=200$ , $Fr=0.3$ and $\unicode[STIX]{x1D716}=0.01$ . For those parameters, the wave boundary layer thickness is $\unicode[STIX]{x1D6FF}_{Re}=0.1$ and the viscous damping length is $L_{Re}=5.15$ . The Hovmöller diagrams focus on an area close to the bottom boundary. We use a vertical resolution of $\text{d}z=0.0067$ , which resolves the wave boundary layer properly. In the DNS, a stretched grid has been implemented on the vertical to avoid any downward reflection. The quasi-linear model captures well the boundary streaming effect. To emphasise the crucial role of the boundary streaming term, we added a diagram in figure 4 of a quasi-linear computation where the boundary forcing has been removed in (2.4) ( $F_{bl}=0$ in (3.7)). We clearly see that the presence of boundary streaming is important to predict accurately the early evolution of the mean flow in this case.

In figure 3(c), we show a Hovmöller diagram of the mean flow computed using the quasi-linear model in the case of no-slip boundary condition. The parameters are $Re=200$ , $Fr=0.3$ and $\unicode[STIX]{x1D716}=0.005$ . As expected from the discussion following (3.9), the boundary forcing generates a strong boundary mean flow going in a direction opposite to the direction of the bulk mean flow. Consistently with our previous estimates of typical time scales for the mean flow evolution, the establishment of the bulk flow occurs at a time scale larger than the establishment of the quasi-stationary boundary flow.

In the no-slip case, the mean flow eventually reaches a stationary state given by

(3.10) $$\begin{eqnarray}\overline{u}_{\infty }(z)=Re\int _{0}^{z}\overline{u^{\prime }w^{\prime }}(z^{\prime })\,\text{d}z^{\prime }.\end{eqnarray}$$

Then, the contribution from the boundary streaming is negligible with respect to the contribution from the bulk streaming. This can be quantified by computing the order of magnitude of typical mean flow amplitudes $U_{w}$ and $U_{bl}$ obtained by splitting the Reynolds stresses in (3.10) into a bulk and a boundary contribution, respectively. Using the large-Reynolds-number asymptotic expressions obtained in (3.8) assuming $ReFr^{2}\rightarrow +\infty$ and $Fr\rightarrow 0$ , we get $U_{w}\sim (\unicode[STIX]{x1D716}ReFr)^{2}$ and $U_{bl}\sim \unicode[STIX]{x1D716}^{2}Re^{1/2}Fr^{-1}$ . Their ratio scales as $(ReFr^{2})^{3/2}$ , and thus tends to diverge: the bulk flow is dominant in the long-time limit.

In the free-slip case, no stationary regime is reached and the mean flow amplitude keeps increasing in time. It can be assessed by considering the $z$ -integrated momentum, $P(t)=\int _{0}^{\infty }\overline{u}(z,t)\,\text{d}z$ . Using the free-slip boundary condition and integrating (2.4), we get $P(t)=(\overline{u^{\prime }w^{\prime }}|_{z=0})t$ . At sufficiently large times, the mean flow varies over the characteristic length scale $\sqrt{t/Re}$ . Consequently, the mean flow amplitude $P/\sqrt{t/Re}$ increases as $t^{1/2}$ : eventually, the feedback of the mean flow on the wave will no longer be negligible. We can however use this mean flow amplitude estimate, together with the large-Reynolds-number asymptotic expressions in (3.8), to infer that $U_{bl}/U_{w}\sim 1/(FrRe^{1/2})$ . This scaling has been obtained under the assumption $Fr^{2}Re\rightarrow +\infty$ . This means that the bulk flow is dominant in the long-time limit, just as in the no-slip case. It is also possible to estimate the time scale $\unicode[STIX]{x1D70F}$ for which the mean flow induced by the bulk streaming becomes of the same order as the mean flow induced by the boundary streaming. When this occurs, the long-time limit is relevant for the estimate of the mean flow induced by boundary streaming, as above: $U_{bl}\sim F_{bl}(0)\sqrt{\unicode[STIX]{x1D70F}/Re}$ . By contrast, assuming $L_{Re}\gg \unicode[STIX]{x1D6FF}_{Re}$ , the flow induced by the bulk streaming must be estimated using an early-time limit: $U_{w}\sim \unicode[STIX]{x1D70F}\unicode[STIX]{x2202}_{z}F_{w}|_{z=0}$ . Then, using $U_{bl}\sim U_{w}$ yields $\unicode[STIX]{x1D70F}\sim Fr^{4}Re^{2}$ . Using the parameters corresponding to figure 4 yields $\unicode[STIX]{x1D70F}\sim 300$ .

3.4 Limitation of the quasi-linear model

To derive the quasi-linear model around a state of rest presented above, the only necessary assumption is $\unicode[STIX]{x1D716}\rightarrow 0$ , with all other parameters fixed. The quasi-linear numerical calculations have been made using the actual solution of the dispersion relation (3.1), but we obtained scalings by assuming simplified expressions for the wave field in the inviscid limit $Fr^{2}Re\rightarrow +\infty$ , together with the hydrostatic limit $Fr\rightarrow 0$ . These two conditions imply that $\unicode[STIX]{x1D6FF}_{Re}/L_{Re}\sim (ReFr^{2})^{-3/2}\rightarrow 0$ , and therefore make possible a clear distinction between a bulk and a boundary contribution to streaming. To establish a self-consistent distinguished limit, we write

(3.11) $$\begin{eqnarray}(\unicode[STIX]{x1D716},Fr,Re)=(\unicode[STIX]{x1D716},\unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D716}^{-\unicode[STIX]{x1D6FD}}).\end{eqnarray}$$

The two simplifying assumptions above correspond to $\unicode[STIX]{x1D6FD}>2\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FC}>0$ . With these assumptions, the wave amplitudes scale as $(\unicode[STIX]{x1D719}_{w},\unicode[STIX]{x1D719}_{bl})\sim (\unicode[STIX]{x1D716}^{1-\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D716}^{1-2\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD}/2})$ in the free-slip case and $(\unicode[STIX]{x1D719}_{w},\unicode[STIX]{x1D719}_{bl})\sim (\unicode[STIX]{x1D716}^{1-\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D716}^{1-\unicode[STIX]{x1D6FC}})$ in the no-slip case. We now list the conditions required for the validity of the linearisation procedure.

In the bulk, neglecting the nonlinear (advection) terms with respect to the viscous terms and the time derivative terms yields the conditions $\unicode[STIX]{x1D6FD}<1+\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FC}<1$ , respectively. Within the boundary layer, the full solution is composed of a bulk and a boundary term. Neglecting nonlinear terms involving both bulk and boundary components leads to an additional constraint $\unicode[STIX]{x1D6FD}<2$ . Neglecting nonlinear terms involving only boundary components leads to an additional constraint $\unicode[STIX]{x1D6FD}>4\unicode[STIX]{x1D6FC}-2$ only in the free-slip case. Neglecting nonlinear terms in the bottom boundary conditions (2.2) does not add any new constraint. There are therefore six inequalities to be satisfied for $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ in the free-slip case, five inequalities for the no-slip case. For the latter case, regimes of parameters for which these conditions are all fulfilled are presented in figure 5. The black area corresponds to regimes fulfilling all the constraints. The additional condition required for the free-slip case is also fulfilled within this black area.

Figure 5. Distinguished limits for the validity of the linear dynamics around a state of rest in the no-slip case: $(Fr,Re)=(\unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D716}^{-\unicode[STIX]{x1D6FD}})$ . Each line delimits a half-plane where one of the constraints is satisfied. The grey scale shows the number of constraints that are satisfied. The black region corresponds to the range of exponents $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ for which the asymptotic approach is self-consistent: all the constraints are satisfied for those scalings. In the free-slip case, there is an additional constraint $\unicode[STIX]{x1D6FD}<4\unicode[STIX]{x1D6FC}-2$ , which is not represented here, but which is fulfilled within the black area. The dashed red lines correspond to limit cases above which the mean flows induced by the bulk streaming and boundary streaming impact the wave field. The red dot corresponds to the regime $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=(1,2)$ where the distinguished limit is marginally satisfied, and where two-way coupling between waves and mean flow can no longer be neglected.

In all the above analysis, we have neglected the feedback of the mean flow on the wave field. This is always valid at sufficiently short times. However, we saw that this can never be satisfied at a long time in the free-slip case since the mean flow keeps increasing in time. In the no-slip case, we found that both the bulk and the boundary mean flow are indeed negligible with respect to the horizontal phase speed ( $|\unicode[STIX]{x2202}_{t}|\gg |\overline{u}\unicode[STIX]{x2202}_{x}|$ ), as $U_{w}\sim (\unicode[STIX]{x1D716}ReFr)^{2}\rightarrow 0$ and $U_{bl}\sim \unicode[STIX]{x1D716}^{2}Re^{1/2}Fr^{-1}\rightarrow 0$ in the distinguished limit.

In the no-slip case, we expect a two-way coupling between waves and mean flow when the induced flows are of order one – with $\unicode[STIX]{x1D716}ReFr\sim 1$ for the bulk flow and $\unicode[STIX]{x1D716}Re^{1/4}Fr^{-1}\sim 1$ for the boundary-driven flow – since the terms involving the mean flow can no longer be ignored to compute the wave field in that case. These additional conditions are represented by the dashed red lines in figure 5. The red dot corresponds to the regime $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=(1,2)$ satisfying marginally the distinguished limit while allowing for order-one mean flows induced by both the bulk and the boundary forcing. Within this regime the viscous terms are of the order of nonlinear terms in the bulk wave equation, thus invalidating the quasi-linear approach. This limitation can be bypassed by introducing additional dissipative terms, such as a linear friction term in the zonal flow equation or the buoyancy equation. Such terms would allow us to control the typical vertical length scale for wave attenuation, related to the intensity of the bulk internal wave streaming, without varying the Reynolds number that constrains the mean flow vertical gradients. To avoid the introduction of such additional parameters, we choose in the following to consider the quasi-linear equations as an ad hoc model for wave–mean flow interactions. This simplified model will illustrate how boundary streaming can affect mean flow properties in the bulk when the feedback of the mean flow on the wave field is taken into account.

4 Application to an idealised analogue of the quasi-biennial oscillation

We consider a standing wave pattern imposed at the bottom boundary: $h(x,t)=\cos (x)\cos (t)$ with a no-slip boundary condition. This idealised configuration is thought to capture the essential mechanism at the origin of the equatorial stratospheric quasi-biennial oscillation (Plumb Reference Plumb1977), and has been experimentally studied by Plumb & McEwan (Reference Plumb and McEwan1978). Two linear waves with equal amplitude and opposite zonal phase speeds are emitted by such a bottom excitation. The resulting Reynolds stress is simply the sum of the Reynolds stresses computed from each individual wave plus a rapidly oscillating term that can be smoothed out by averaging over this fast oscillation. The Reynolds stress divergences induced by the two waves are opposite and annihilate each other in the absence of mean flow. Above a certain value of the amplitude of the waves, a Hopf bifurcation occurs: a vacillating mean flow is generated and approaches a limit cycle (Plumb Reference Plumb1977). Plumb & McEwan (Reference Plumb and McEwan1978) reported the spontaneous generation of an oscillating mean flow in laboratory experiments when the wave amplitude exceeds a threshold, and compared their measurements against quasi-linear computations. They considered a no-slip bottom boundary condition for the mean flow but inviscid impermeability condition for the wave field, allowing them to ignore any boundary layer effect. Here, we investigate the effect of the viscous boundary layers and the associated boundary streaming on the oscillation arising with the standing wave excitation, assuming a no-slip condition for both the mean flow and the waves. We show that the inclusion of boundary streaming induces important alterations on the mean flow in this idealised model of wave–mean flow interactions.

In § 3, we ignored the effect of the mean flow on the wave field. We need here to take this feedback into account, as the initial instability arises from a perturbation of the mean flow itself. The effect of the mean flow on the wave is included by performing a WKB expansion of the wave field following the method of Muraschko et al. (Reference Muraschko, Fruman, Achatz, Hickel and Toledo2015), but including dissipative effects. The full calculation is detailed in appendix A. The Reynolds stress divergence is then computed and inserted into the mean flow equation (2.4) in order to compute the long-time evolution of $\overline{u}$ . This task is done numerically using the results of appendix A and for the no-slip boundary condition in (3.6). While Plumb & McEwan (Reference Plumb and McEwan1978) considered an asymptotic expression for the bulk solution of the dispersion relation (A 5), our numerical calculations use the actual solutions. As discussed at the end of appendix A, this solution captures important corrections close to the critical layers, where the mean flow is of the order of the wave zonal phase speed.

The resulting Hovmöller diagrams of mean flow time series are shown in figure 6 for different values of the Reynolds number. The time series used for the upper plots have been computed using the full quasi-linear model while those used for the bottom plots have been computed without the boundary layer contributions. All simulations start with the same initial perturbation. In figure 6(a), we see that the inclusion of boundary streaming has altered the critical parameter values at which the bifurcation to mean flow reversals occurs. In figure 6(b), the Reynolds number is increased and the oscillation is present in both cases. However, the oscillation period is decreased by 20 % when the boundary streaming is included. By further increasing the Reynolds number we see in figure 6(c) that the inclusion of the boundary streaming significantly changes the mean flow oscillation. This new regime exhibiting an additional frequency in the signal can actually be reached without the boundary streaming but at a larger Reynolds number. A similar regime has been reported by Kim & MacGregor (Reference Kim and MacGregor2001), and we will study these bifurcations more thoroughly in a companion paper. Our aim here is to show that the presence of the wave boundary layer has an impact on such bifurcation diagrams, in addition to significantly altering the period of oscillations in the periodic case.

Figure 6. The mean flow, $\overline{u}(z,t)$ , is generated by the streaming coming from two counter-propagating waves with the same amplitude and opposite horizontal phase speed, generated by a vertically oscillating bottom boundary with no-slip condition using the quasi-linear model. Hovmöller diagrams of the mean flow time series are shown for three different Reynolds numbers. In each panel (a–c) the upper plot corresponds to a case where the boundary streaming has been included in the computation while the lower plot corresponds to a computation with the same parameters but without the contribution of the boundary streaming terms. In all cases, $Fr=0.15$ and $\unicode[STIX]{x1D716}=0.3$ . In panel (b), the mean flow oscillates with an oscillation period of approximately $40$ and $50$ time units for the case with (upper plot) and without (lower plot) boundary streaming, respectively.

Figure 7. (a) Mean flow snapshots extracted from the two time series plotted in figure 6(b) computed with (BL) and without (NoBL) boundary streaming and taken at $t=10$ and $t=12.5$ , respectively. (b) Plot of the associated Reynolds stress divergences obtained for the wave propagating in the direction of positive mean flow (‘ $+$ ’) and for the counter-propagating one (‘ $-$ ’) considering the mean flow profile obtained by either including (BL) or ignoring (NoBL) the boundary streaming. (c) Plot of the total Reynolds stress divergences, the sum of the two counter-propagating waves contributions for the case with boundary streaming (BL) and without boundary streaming (NoBL).

For the range of parameters corresponding to figure 6(b,c), the mean flow reaches an amplitude close to the phase speed of the waves. To investigate the effect of the boundary layers in these cases, we consider two mean flow snapshots, plotted in figure 7, taken from the two time series plotted in figure 6(b). The snapshots are taken at the same stage of the oscillation cycle.

The Reynolds stress divergences computed using the two mean flow snapshots shown in figure 7(a) and considering the counter-propagating waves separately are plotted in figure 7(b). The total Reynolds stress divergences are plotted in figure 7(c). As expected, the boundary layers significantly modify the streaming vertical profiles. Interestingly, while bulk streaming is dominated by the wave travelling in the same direction as the mean flow, the main discrepancy between the case with and without boundary streaming comes from the boundary forcing associated with the wave going in a direction opposite to the mean flow near the bottom boundary.

In figure 6(b), we see that the bottom profile of the mean flow is approximately steady before a reversal. Let us call $\unicode[STIX]{x1D706}_{Re}$ the typical length over which the mean flow reaches its extremum value. This velocity is of order one as it is close to the gravity wave zonal phase speed. Using (2.4), we infer the typical scaling of $\unicode[STIX]{x1D706}_{Re}$ by balancing the viscous term $Re^{-1}\unicode[STIX]{x2202}_{zz}\overline{u}$ and the streaming term $\unicode[STIX]{x2202}_{z}\overline{u^{\prime }w^{\prime }}$ . This last term is dominated by the bulk contribution $F_{w}$ associated with the wave that propagates in the same direction as the mean flow. This yields $\unicode[STIX]{x1D706}_{Re}\sim (Re\,F_{w}(0))^{-1}$ . The order of magnitude of $F_{w}(0)$ can be obtained by using the asymptotic expression in (3.8), under the assumption $Fr\rightarrow 0$ and $Re^{2}Fr\rightarrow \infty$ . This yields $\unicode[STIX]{x1D706}_{Re}\sim Fr/(\unicode[STIX]{x1D716}^{2}Re)$ . Using the parameters of figure 7, we find $\unicode[STIX]{x1D706}_{Re}\sim 0.1$ .

The presence of the mean flow thus leads to a new ‘bulk scale’ $\unicode[STIX]{x1D706}_{Re}$ that differs from the bulk dissipative scale $L_{Re}$ for the flow at rest; in particular, in our distinguished limit (i.e. an inviscid limit), the former tends to zero, while the latter tends to infinity. In the bulk, the wave propagating in the same direction as the mean flow is controlled by $\unicode[STIX]{x1D706}_{Re}$ ; the wave propagating in the other direction is controlled by $L_{Re}$ . The typical scale of bulk flow reversals is thus given by $\unicode[STIX]{x1D706}_{Re}$ . We expect boundary streaming to have a significant impact on these reversals when they occur within the boundary layer, i.e. when $\unicode[STIX]{x1D706}_{Re}\sim \unicode[STIX]{x1D6FF}_{Re}$ or $\unicode[STIX]{x1D706}_{Re}\ll \unicode[STIX]{x1D6FF}_{Re}$ . This is indeed the case in figure 6(c), where $\unicode[STIX]{x1D6FF}_{Re}\approx 0.1$ . It is instructive to establish the range of parameters for which this condition is satisfied. Using $\unicode[STIX]{x1D6FF}_{Re}\sim 1/Re^{1/2}$ , we find that this length scale is larger than or of the order of $\unicode[STIX]{x1D706}_{Re}$ when $Fr/(\unicode[STIX]{x1D716}^{2}Re^{1/2})\rightarrow 0$ . The above analysis suggests the possibility for active control of the boundary layer on the bulk flow when $(Fr,Re)=(\unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FC}},\unicode[STIX]{x1D716}^{-\unicode[STIX]{x1D6FD}})$ with $\unicode[STIX]{x1D6FD}\geqslant 4-2\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FC}>0$ and $\unicode[STIX]{x1D6FD}>2\unicode[STIX]{x1D6FC}$ .

As seen in figure 5, the distinguished limits consistent with an active control of the boundary layer on bulk flow reversal in the ad hoc quasi-linear model do not overlap with the distinguished limits ensuring the validity of the quasi-linear dynamics around a state of rest. In fact, both sets of constraints can only be marginally satisfied at the point $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=(1,2)$ . This is the main caveat of the analysis presented in this section and illustrated in figure 7: because $\unicode[STIX]{x1D716}ReFr\sim 1$ , the nonlinear terms involving bulk waves are not negligible with respect to the bulk viscous term. Furthermore, since $\unicode[STIX]{x1D716}Re^{1/2}\sim 1$ , the boundary layer can hardly be considered as linear. Whether active control of the boundary on the bulk flow persists when nonlinear effects are reincorporated into the problem needs to be addressed in a future work.

Another caveat of the quasi-linear model presented here comes from the assumptions underlying the WKB approach used to compute the wave field. We assumed that there is a vertical scale separation between the wave and the mean flow and that the wave field reaches its steady state in a time much shorter than the typical time of evolution of the mean flow. The wave field is thus computed using a static WKB approximation with a frozen-in-time mean flow. Since the mean flows shown in figure 6 exhibit sharp shear at the bottom, and since they reach values of the order of the zonal phase speed of the waves, those hypotheses are not valid. Nevertheless, this WKB approximation is the simplest way of accounting for the mean flow effect on the wave field, and a useful first step to understand their interactions.

We should finally stress that the no-slip bottom boundary condition is most certainly irrelevant to the actual quasi-biennial oscillation occurring in the upper atmosphere and that our model has been derived in a distinguished limit for which the viscous boundary layer is much larger than the boundary height variations, which is not satisfied in laboratory experiments. However, despite these limitations, our results show that the boundary conditions and the associated wave boundary layers should not be overlooked, since boundary streaming may have a quantitative impact on mean flow reversals in the domain bulk.