2.1. Model formulation

We consider the evolution of 3-D localized finite-amplitude perturbations of a steady parallel unidirectional boundary-layer shear flow $\boldsymbol {U}^{\star }(z^{\star })$ adjacent to an infinite horizontal boundary. Throughout the paper we denote dimensional quantities by superscript stars. The boundary layer has a weak density stratification $N^{\star 2}(z^{\star })=g^{\star }\,{\rm d}\rho ^{\star }_{0}(z^{\star })/{\rm d} z^{\star }/\rho ^{\star }_{0}(z^{\star })$, where $N^{\star 2}$ is the buoyancy frequency, $g^{\star }$ is gravitational acceleration, $\rho ^{\star }_{0}$ is the reference density that depends only on the vertical coordinate $z^{\star }$. The full density $\rho ^{\star }$ is sum of the reference density $\rho _{0}^{\star }(z)$ and density perturbation $\bar {\rho }^{\star }$

(2.1)\begin{equation} \rho^{{\star}}=\rho_{0}^{{\star}}(z^{{\star}})+ \bar{\rho}^{{\star}}(x^{{\star}},y^{{\star}},z^{{\star}}, t^{{\star}});\quad |\bar{\rho}^{{\star}}|\ll\rho_{0}^{{\star}}, \end{equation}

the perturbation density $\bar {\rho }^{\star }$ depends on horizontal coordinates $x^{\star }$, $y^{\star }$, vertical coordinate $z^{\star }$ and time $t^{\star }$. Both the shear and stratification are confined to the boundary layer as sketched in figure 1. There are no other assumptions regarding the profiles of the shear and stratification.

Figure 1. Sketch of geometry of a generic free-surface boundary-layer profile with the shear and stratification localized in a layer of thickness $d$. The vertical scales of shear and stratification are assumed be of the same order as the boundary-layer thickness $d$. No other assumptions regarding the profiles of shear and stratification are required. A typical free-surface boundary layer has the maximum of velocity at the surface $U_{max}=U(0)$.

In the Cartesian frame with the fluid occupying the half-space $z^{\star }>0$ and with $x^{\star }$ and $y^{\star }$ directed streamwise and spanwise, respectively, the Navier–Stokes equations complemented by the mass conservation and incompressibility equations take the form

(2.2a)\begin{gather} \rho_{0}^{{\star}}[D_{t^{{\star}}}u^{{\star}}+w^{{\star}}U^{{\star}\prime}]+ p_{x^{{\star}}}^{{\star}} ={-}\bar{\rho}^{{\star}}[D_{t^{{\star}}}u^{{\star}}+ w^{{\star}}U^{{\star}\prime}] -\rho_{0}^{{\star}}(\boldsymbol{u}^{{\star}}\boldsymbol{\cdot}\boldsymbol{\nabla})u^{{\star}}- \bar{\rho}^{{\star}}(\boldsymbol{u}^{{\star}}\boldsymbol{\cdot}\boldsymbol{\nabla})u^{{\star}} \nonumber\\ +\,\mu^{{\star}}\nabla^{2}u^{{\star}}, \end{gather}

(2.2b)$$\begin{gather} \rho_{0}^{{\star}} D_{t^{{\star}}}v^{{\star}}+p_{y^{{\star}}}^{{\star}}={-}\bar{\rho}^{{\star}} D_{t^{{\star}}}v^{{\star}}-\rho_{0}^{{\star}} (\boldsymbol{u}^{{\star}}\boldsymbol{\cdot}\boldsymbol{\nabla})v^{{\star}}-\bar{\rho}^{{\star}} (\boldsymbol{u}^{{\star}}\boldsymbol{\cdot}\boldsymbol{\nabla})v^{{\star}}+ \mu^{{\star}}\nabla^{2}v^{{\star}}, \end{gather}$$

(2.2c)$$\begin{gather}\rho_{0}^{{\star}} D_{t^{{\star}}}w^{{\star}}+p_{z^{{\star}}}^{{\star}}+ \bar{\rho}^{{\star}}g^{{\star}}={-}\bar{\rho}^{{\star}} D_{t^{{\star}}}w^{{\star}}- \rho_{0}^{{\star}}(\boldsymbol{u}^{{\star}}\boldsymbol{\cdot}\boldsymbol{\nabla})w^{{\star}}- \bar{\rho}^{{\star}}(\boldsymbol{u}^{{\star}}\boldsymbol{\cdot}\boldsymbol{\nabla})w^{{\star}}+ \mu^{{\star}}\nabla^{2}w^{{\star}}, \end{gather}$$

(2.2d)$$\begin{gather}D_{t^{{\star}}}\bar{\rho}^{{\star}}+w^{{\star}} {\rho_{0}^{{\star}}}^{\prime} =K^{{\star}}\nabla^{2}\bar{\rho}^{{\star}}- (\boldsymbol{u}^{{\star}}\boldsymbol{\nabla}) \bar{\rho}^{{\star}}, \end{gather}$$

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

where $\boldsymbol {U}^{\star }=(U^{\star }(z^{\star }),0,0)$ is the basic flow localized in the boundary layer, $\boldsymbol {u}^{\star }=(\boldsymbol {q}^{\star },w^{\star })=(u^{\star },v^{\star },w^{\star })$ and $p^{\star }$ are, respectively, the velocity and pressure perturbations, $D_{t^{\star }}=\partial _{t^{\star }}+U^{\star }\partial _{x^{\star }}$ is the material derivative, $\mu ^{\star }$ is the dynamic viscosity, while $K^{\star }$ is the mass diffusivity coefficient, $\nabla ^{2}$ stands for the 3-D Laplacian. Here, the prime denotes the derivatives with respect to $z^{\star }$. When the model is applied to turbulent boundary layers, viscosity and diffusivity are understood as eddy viscosity and diffusivity and, for simplicity, assumed to be constant. We impose no restrictions on $\boldsymbol {U}^{\star }(z^{\star })$, apart from the assumption that the flow is plane parallel (non-parallel effects and 3-D boundary layers will be considered elsewhere). In contrast to the original derivations in Shrira (Reference Shrira1989) and Voronovich et al. (Reference Voronovich, Shrira and Stepanyants1998), we do not require ${U}^{\star }(z^{\star })$ not to have inflection points.

The boundary conditions for the perturbations $\boldsymbol {u}^{\star }$ at the undisturbed fluid surface $z^{\star }=0$ are: the ‘no-flux’ condition

(2.3)\begin{equation} w^{{\star}}(z^{{\star}}=0)=0, \end{equation}

complemented by the no-stress condition

(2.4)\begin{equation} \left.\frac{\partial u^{{\star}}}{\partial z^{{\star}}}\right|_{z^{{\star}}=0}= \left.\frac{\partial v^{{\star}}}{\partial z^{{\star}}}\right|_{z^{{\star}}=0}=0. \end{equation}

The boundary condition at infinity is that of vanishing perturbation velocity

(2.5)\begin{equation} \boldsymbol{u^{{\star}}} \to 0, \quad\mbox{as}\ z^{{\star}} \to \infty. \end{equation}

We complete the formulation of our initial problem by specifying the perturbation velocity field at the initial moment, $\boldsymbol {u}^{\star }(\boldsymbol {x}^{\star }, 0)$. We are primarily interested in localized initial perturbations

(2.6)\begin{equation} \boldsymbol{u^{{\star}}}(\boldsymbol{x^{{\star}}}, 0)\to 0, \quad\mbox{as}\ x^{{\star}},y^{{\star}} \to \infty. \end{equation}

The Navier–Stokes equations (2.2) with the initial and boundary conditions (2.5), (2.3) and (2.4), constitute the mathematical formulation of the problem.

2.2. Scaling

We begin by specifying the scaling. The basic boundary-layer shear flow $U^{\star }(z^{\star })$ has a maximal velocity, usually at the surface, which we denote as $V_{0}^{\star }$. As a characteristic value of the buoyancy frequency $N^{\star }(z^{\star })$ we take its value $N_{0}^{\star }$ at the surface. If $N^{\star }(z^{\star }=0)$ vanishes, any other reference point can be chosen. The characteristic streamwise and spanwise scales of perturbations we denote as $L^{\star }$, while for the vertical scale we choose the boundary-layer thickness $d^{\star }$. We are considering long perturbations with $L^{\star }\gg d^{\star }$. We non-dimensionalize the variables as follows:

(2.7)\begin{equation} \left.\begin{gathered} \tilde{u}=\frac{u^{{\star}}}{V_{0}^{{\star}}},\quad \tilde{v}=\frac{v^{{\star}}}{V_{0}^{{\star}}},\quad \tilde{w}=\frac{w^{{\star}}}{V_{0}^{{\star}}},\quad \tilde{x}=\frac{x^{{\star}}}{L^{{\star}}},\\ \tilde{y}=\frac{y^{{\star}}}{L^{{\star}}},\quad \tilde{z}=\frac{z^{{\star}}}{d^{{\star}}}, \quad \tilde{U}=\frac{U^{{\star}}}{V_{0}^{{\star}}},\quad \tilde{t}=\frac{V_{0}^{{\star}}}{L^{{\star}}}t^{{\star}},\\ \tilde{N}=\frac{N^{{\star}}}{N_{0}^{{\star}}},\quad \tilde{p}= \frac{p^{{\star}}}{\rho_{0}^{{\star}}{V_{0}^{{\star}}}^{2}},\quad \tilde{\rho}=\frac{\rho^{{\star}}}{\rho_{0}^{{\star}}}, \end{gathered}\right\} \end{equation}

where tildes denote non-dimensional quantities. To proceed, we first estimate the magnitudes of the perturbations from the governing equations (2.2). With tildes omitted the magnitude of the vertical velocity perturbation $[\,w\,]$ expressed in terms of the magnitudes of the horizontal components $[\,q\,]$,

(2.8)\begin{equation} [\,w\,]=\frac{d^{{\star}}}{L^{{\star}}}[\,q\,]. \end{equation}

The momentum equation (2.2a) yields the characteristic time scale,

(2.9)\begin{equation} [\,t\,]=\frac{L^{{\star}}}{V_{0}^{{\star}}}. \end{equation}

On multiplying mass conservation equation (2.2d) by $g^{\star }$ and dividing it by the reference density $\rho _{0}^{\star }$, we substitute buoyancy frequency $N^{\star 2}$ for $({g^{\star }}/{\rho _{0}^{\star }}) ({\partial \rho _{0}^{\star }}/{\partial z^{\star }})$ using the Boussinesq approximation, which yields

(2.10)\begin{equation} g^{{\star}}D_{t^{{\star}}}\left(\frac{\bar{\rho}^{{\star}}}{\rho_{0}^{{\star}}}\right) + w^{{\star}}N^{\star2}(z^{{\star}})=g^{{\star}}\frac{\nu^{{\star}}}{Sc\rho_{0}^{{\star}}} \nabla^{2}\bar{\rho}^{{\star}}-\frac{g^{{\star}}}{\rho_{0}^{{\star}}} (\boldsymbol{u}^{{\star}}\boldsymbol{\nabla})\bar{\rho}^{{\star}}, \end{equation}

where $Sc=\nu ^{\star }/{K^{\star }}$ is the Schmidt number, the ratio of the kinematic viscosity $\nu ^{\star }={\mu ^{\star }}/ {\rho _{0}^{\star }}$ and mass diffusivity ${K^{\star }}$.

To proceed we express (2.10) in non-dimensional form and substitute equations (2.8) and (2.9), to get an estimate characteristic magnitude of density perturbations (tildes and star dropped)

(2.11)\begin{equation} \left[\frac{\bar{\rho}}{\rho_{0}}\right]=\frac{N_{0}^{2}d}{g}\frac{[\,q\,]}{V_{0}}. \end{equation}

In the derivation, we employ the Boussinesq approximation, whereby, if not multiplied by $g$, we neglect the terms in the equations due to the equilibrium density $\rho _{0}$ variations. The Boussinesq approximation holds under the assumption

(2.12)\begin{equation} \frac{N_{0}^{2}}{g/d}\ll \frac{[\,q\,]}{V_{0}}, \end{equation}

which is valid in most oceanic and atmospheric contexts where the equilibrium density variations are small. This inequality allows us to drop small nonlinear buoyancy terms. In nature the parameters of the surface layer of the ocean vary very widely. Still, we provide a few characteristic values: $V_0^{\star }\sim 10^{-1}$ m s$^{-1}$, $d^{\star } \sim 50\,\textrm {m},\ N_0 \sim 10^{-4}\,\textrm {s}^{-1}$.

The regimes we study are characterized by four independent non-dimensional small/large parameters specifying, respectively,

1. (i) the smallness of nonlinearity: ${[\,q^{\star }\,]}/{V_{0}^{\star }}=\varepsilon \ll 1$;

2. (ii) the smallness of characteristic wavenumbers, which we also refer to as the dispersion parameter: ${d^{\star }}/{L^{\star }}=\varepsilon _{D}\ll 1$;

3. (iii) the weakness of the boundary-layer stratification characterized by the Richardson number: $Ri=({N_{0}^{\star }d^{\star }}/{V_{0}^{\star }})^{2}\ll 1$;

4. (iv) the weakness of dissipative effects: $Re^{-1}= {\nu ^{\star }}/{V_{0}^{\star }d^{\star }}\ll 1$.

From now on we will be using only the non-dimensional tilde variables (2.7) omitting the tildes. Under the Boussinesq approximation upon non-dimensionalization the Navier–Stokes equations (2.2) take the form

(2.13a)$$\begin{gather} D_{t}u+wU^{\prime}+p_{x}={-}\varepsilon(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})u+ \frac{1}{Re}\left[\frac{1}{\varepsilon_{D}}\partial_{zz}+ \varepsilon_{D}\nabla_{{\perp}}^{2}\right]u, \end{gather}$$

(2.13b)$$\begin{gather}D_{t}v+p_{y}={-}\varepsilon(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})v+ \frac{1}{Re}\left[\frac{1}{\varepsilon_{D}}\partial_{zz}+ \varepsilon_{D}\nabla_{{\perp}}^{2}\right]v, \end{gather}$$

(2.13c)$$\begin{gather}D_{t}w+\frac{1}{\varepsilon_{D}}p_{z}+\varepsilon\varepsilon_{D}Ri \bar{\rho}={-}\varepsilon(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})w+ \frac{1}{Re}\left[\frac{1}{\varepsilon_{D}}\partial_{zz}+ \varepsilon_{D}\nabla_{{\perp}}^{2}\right]w, \end{gather}$$

(2.13d)$$\begin{gather}D_{t}\bar{\rho}+ \frac{1}{\varepsilon_{D}}N^{2}(z)w=\frac{\varepsilon_{D}}{Re\, Sc}\left[\frac{1}{\varepsilon_{D}}\partial_{zz}+ \varepsilon_{D}\nabla_{{\perp}}^{2}\right]\bar{\rho}- \varepsilon(\boldsymbol{u}\boldsymbol{\nabla})\bar{\rho}, \end{gather}$$

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

where $\nabla _{\perp }^{2}=\partial _{xx}^{2}+\partial _{yy}^{2}$.

After some algebra, upon eliminating pressure from (2.13), we get a single equation for the vertical component of velocity ${w}$

(2.14)\begin{equation} D_{t}^{2}\partial_{zz}^2w-U^{\prime\prime}D_{t}\partial_{x}w={\mathcal{N}} -Ri\,N^{2}\nabla_{{\perp}}^{2}w+\left(\varepsilon_{D}Re\right)^{{-}1}D_{t} \partial_{zzzz}^{4}w+{\mathcal{M}}, \end{equation}

where

(2.15)\begin{equation} \left.\begin{gathered} {\mathcal{N}}=\varepsilon D_{t}\partial_{z} \boldsymbol{\nabla}_{{\perp}}[(\boldsymbol{u}\boldsymbol{\nabla})\boldsymbol{q}]- \varepsilon_{D}^{2}D_{t}^{2}\nabla_{{\perp}}^{2}w-\varepsilon_{D}^{2} \varepsilon D_{t}\nabla_{{\perp}}^{2}(\boldsymbol{u}\boldsymbol{\nabla})w,\\ {\mathcal{M}}=Ri \varepsilon\nabla_{{\perp}}^{2}(\boldsymbol{u}\boldsymbol{\nabla})\bar{\rho}+ \frac{1}{Sc\,Re}Ri\varepsilon\nabla_{{\perp}}^{2}\bar{\rho}^{\prime\prime}. \end{gathered}\right\} \end{equation}

Note that the equation is closed only in the linear approximation and just the leading-order viscous and diffusion terms are retained. The contribution of diffusion of mass is characterized by the inverse Schmidt number. Under assumption $Sc\sim O(1)$ diffusion of mass proves to be negligible in our further analysis.

Aiming to describe the dynamics of 3-D perturbations in the boundary layer taking into account nonlinearity, long-wave dispersion, stratification and viscous effects in the distinguished limit we set

(2.16a–d)\begin{equation} \varepsilon=\frac{[\,q\,]}{V_{0}}\ll 1,\quad \varepsilon_{D}=\frac{d}{L}= O(\varepsilon),\quad Ri=\left(\frac{N_{0}d}{V_{0}}\right)^{2}=O(\varepsilon),\quad Re=O(\varepsilon^{{-}2}). \end{equation}

The adopted scaling can be also expressed in terms of the Reynolds number

(2.17a–c)\begin{equation} \varepsilon \sim Re^{{-}1/2}, \quad \varepsilon_{D} \sim Re^{{-}1/2}, \quad Ri \sim Re^{{-}1/2}. \end{equation}

2.3. Asymptotic expansion

To rationalize our choice of asymptotic expansion in powers of $\varepsilon$, it is helpful to consider first the inviscid non-stratified reduction of the Navier–Stokes equations (2.13) linearized around the steady profile $U(z)$ with pressure $p$, the streamwise and transverse perturbation velocities, $u$ and $v$, expressed in terms of vertical velocity $w$ (see e.g. Miropol'sky (Reference Miropol'Sky2001), (2.58)). For long-wave perturbations such a reduction reads

(2.18a)$$\begin{gather} \partial_{z}D_{t}u={-}[\partial_{z}(wU^{\prime}) -\partial_{x}(D_{t}w)], \end{gather}$$

(2.18b)$$\begin{gather}\partial_{z}D_{t}v=\partial_{y}(D_{t}w), \end{gather}$$

(2.18c)$$\begin{gather}\partial_{z}p=D_{t}w . \end{gather}$$

Noting that, under adopted assumptions specifying the wavelength scale of the perturbations in terms of the nonlinearity parameter, $\varepsilon _{D}= O(\varepsilon )$, we have

(2.19a,b)\begin{equation} \partial_{x}\sim \partial_{y}\sim O(\varepsilon), \quad\partial_{z}\sim O(1). \end{equation}

Since in our scaling $U,U^{\prime }\sim O(1)$, it is easy to see that to leading order the material derivative $D_{t}=(U-c)\partial _{x}\sim O(\varepsilon )$, where $c$ is an unspecified yet phase velocity of long perturbations. By virtue of our definition of $\varepsilon (\varepsilon =[u]/V_0)$, $u$ is $O(\varepsilon )$. Therefore, upon omitting the higher-order terms, the relations (2.18) reduce to

(2.20a)$$\begin{gather} (U-c)\partial_{x} u={-}\underbrace{wU^{\prime}}_{O(\varepsilon^{2})},\quad \implies w\sim O(\varepsilon^{2}), \end{gather}$$

(2.20b)$$\begin{gather}\partial_{z}[(U-c)v]=\underbrace{\partial_{y}(U-c)w}_{O(\varepsilon^{3})}, \quad\implies v\sim O(\varepsilon^{3}), \end{gather}$$

(2.20c)$$\begin{gather}\partial_{z}p=\underbrace{(U-c) \partial_{x}w}_{O(\varepsilon^{3})},\quad \implies p \sim O(\varepsilon^{3}). \end{gather}$$

In addition to the above perturbation velocities scaling, we derive below the scaling of density from the linearized inviscid reduction of the Navier–Stokes equation. First, it is straightforward to express perturbation density $\bar {\rho }$ in terms of vertical velocity $w$ from the linearized equation (2.10) according to

(2.21a,b)\begin{equation} gD_{t}\left(\frac{\bar{\rho}}{\rho_{0}}\right)=N^{2}w,\quad \implies g(U-c)\partial_{x}\bar{\rho}=\underbrace{\rho_{0}}_{O(1)} \underbrace{N^{2}w}_{O(\varepsilon^{3})},\quad \bar{\rho}\sim O(\varepsilon^{2}). \end{equation}

Thus, assuming the streamwise perturbation velocity $u$ to be $O(\varepsilon )$ and setting $\varepsilon _{D}={d}/{L} = O(\varepsilon )$ uniquely dictates the scaling (2.20), (2.10) of all other dependent variables inside the boundary layer. Taking into account the nonlinear and viscous terms neglected in this analysis does not affect the found scaling. We note that the above scaling has no self-consistent alternatives. Therefore, for the full Navier–Stokes equations (2.13) in the boundary layer we adopt the following asymptotic expansion:

(2.22a)$$\begin{gather} u=U(z)+\varepsilon u_{1}+\varepsilon^{2}u_{2}+\varepsilon^{3}u_{3}+\cdots, \end{gather}$$

(2.22b)$$\begin{gather}w=\varepsilon^{2}w_{2}+\varepsilon^{3}w_{3}+\varepsilon^{4}w_{4}+\cdots, \end{gather}$$

(2.22c)$$\begin{gather}v=\varepsilon^{3}v_{3}+\varepsilon^{4}v_{4}+\cdots, \end{gather}$$

(2.22d)$$\begin{gather}p=\varepsilon^{3}p_{3}+\varepsilon^{4}p_{4}+\cdots, \end{gather}$$

(2.22e)$$\begin{gather}\bar{\rho}=\varepsilon^{2}\bar{\rho}_{2}+\varepsilon^{3}\bar{\rho}_{3}+\cdots, \end{gather}$$

where $u_i, v_i, w_i, p_i$ and $\bar {\rho }_{i}$ are $O(1)$ functions of $x,y,z, t$. The scaling (2.22) will be employed only inside the boundary layer, while outside the boundary layer and in the immediate vicinity of the boundary, in the viscous sub-layer, the scaling is different and will be specified in the next section.

3.1. Preliminary consideration. Layout of the triple-deck scheme

In this section we derive the nonlinear evolution equation (1.1) for long-wave 3-D perturbations employing a version of the ‘triple deck’ asymptotic approach (e.g. Neiland Reference Neiland1969; Stewartson Reference Stewartson1969; Messiter Reference Messiter1970; Van Dyke Reference Van Dyke1975). As is common for this approach, we distinguish three domains or ‘decks’ in $z$, with different balance between the terms in (2.14). The domains are sketched in figure 2. Following the convention the bulk of the boundary layer is referred to as the ‘main deck’ or ‘middle deck’. There we balance nonlinearity, dispersion and viscous effects and employ the asymptotic expansion (2.22). Immediately adjacent to the boundary lies much thinner viscous sub-layer or the ‘first deck’, where viscous terms are dominant, while the nonlinearity is small, but not negligible. The semi-infinite domain outside the boundary layer, where both viscous and nonlinear terms are negligible and the flow is potential, is referred to as the ‘outer flow’ or the ‘third deck’. We adopt this terminology introduced by Stewartson (Reference Stewartson1969). In contrast to the widely followed convention we do not a priori scale our variables in terms of powers of inverse Reynolds number, in our context we prefer to use the scaling in powers of $\varepsilon$. The derivation largely follows that in Voronovich et al. (Reference Voronovich, Shrira and Stepanyants1998), with three key differences:

1. (i) In contrast to Voronovich et al. (Reference Voronovich, Shrira and Stepanyants1998), the presence of weak stratification in the boundary layer (‘main deck’) is taken into account; we consider stratification to be weak and localized inside the boundary layer.

2. (ii) The Reynolds number here is allowed to be smaller: $Re\sim \varepsilon ^{-2}$, not $\sim \varepsilon ^{-3}$.

3. (iii) No assumption regarding the absence of inflection points is required.

   

   Figure 2. Sketch of triple-deck structure illustrating three different regions. The first region adjacent to the boundary is the viscous sub-layer where viscosity dominates, the middle region is the main deck, where nonlinearity, dispersion, stratification and viscosity are all balanced. The third semi-infinite region is the outer flow region where the motion is irrotational.

   The rationale for the assumption (iii) was to exclude a priori the possibility of linear inflectional instability, which, in principle, might interfere with the nonlinear dynamics we are studying. However, a closer look at the now available results on linear instabilities in boundary layers (e.g. Healey Reference Healey2017) enabled us to lift the restriction. Although the profiles with inflection points can be indeed unstable, these inflectional instabilities are long-wave instabilities with the maximal growth rates $O(\varepsilon ^{3})$ and, thus, do not affect the processes with the $O(\varepsilon ^{-2})$ characteristic periods which we focus upon.

3.2. Inside the boundary layer. The main deck

We begin with analysis of the motion in the main deck. The scaling (2.22) based upon the distinguished limit which balances nonlinearity, weak dispersion, stratification and viscosity, provides the basis of our asymptotic analysis inside the main deck. On substituting the already adopted relations between the small parameters: $\varepsilon _{D} \sim Ri\sim \varepsilon$, and introducing re-scaled Richardson and Reynolds numbers denoted as $Ri_{\ast }$ and $Re_{\ast }$, such that $Ri=\varepsilon Ri_{\ast }$ and $Re_{\ast }=\varepsilon ^{2} Re$, into equation (2.14), we make the scaling of each term more explicit

(3.1)\begin{equation} D_{t}^{2}\partial_{zz}^2 w-U^{\prime\prime}D_{t}\partial_{x}w= \varepsilon D_{t}\partial_{z}{\mathcal{N}}_{L}-\varepsilon Ri_{{\ast}} N^{2}\nabla_{{\perp}}^{2}w +\frac{\varepsilon}{Re_{{\ast}}}D_{t}w^{iv}-\varepsilon^{2}D_{t}^{2}\nabla_{{\perp}}^{2}w+{\mathcal{M}}_{1}, \end{equation}

where

(3.2)\begin{equation} {\mathcal{M}}_{1}=\varepsilon^{2}Ri_{{\ast}}\nabla_{{\perp}}^{2} (\boldsymbol{u}\boldsymbol{\nabla})\bar{\rho}-\varepsilon^{3} D_{t}\nabla_{{\perp}}^{2}(\boldsymbol{u}\boldsymbol{\nabla})w+ \varepsilon^{4}\frac{1}{Sc\,Re_{{\ast}}}Ri_{{\ast}}\nabla_{{\perp}}^{2}\bar{\rho}^{\prime\prime}, \end{equation}

and

(3.3)\begin{align} {\mathcal{N}}_{L} &= [\partial_{x}(u\partial_{x}u+v\partial_{y}u+w\partial_{z}u)+ \partial_{y}(u\partial_{x}v+v\partial_{y}v+w\partial_{z}v)], \nonumber\\ &=[\varepsilon\partial_{x}(\varepsilon^{3}u\partial_{x}u+ \varepsilon^{5}v\partial_{y}u+\varepsilon^{3}w\partial_{z}u)+ \varepsilon\partial_{y}(\varepsilon^{5}u\partial_{x}v+ \varepsilon^{7}v\partial_{y}v+\varepsilon^{5}w\partial_{z}v)] \nonumber\\ &=\varepsilon^{4}[\partial_{x}(u\partial_{x}u+w\partial_{z}u)] + \varepsilon^{6}[\partial_{x}(v\partial_{y}u) +\partial_{y} (u\partial_{x}v +w\partial_{z}v +\varepsilon^{2} v\partial_{y}v)]. \end{align}

Although the streamwise and spanwise scales are assumed to be comparable, according to (2.22) the spanwise velocity is two orders of magnitude smaller, which enables us to split the nonlinear term ${\mathcal {N}}_{L}$ into three parts and neglect the $O(\varepsilon ^{6})$, $O(\varepsilon ^{7})$ and $O(\varepsilon ^{8})$ contributions. The terms $\varepsilon ^{2}Ri_{\ast }\nabla _{\perp }^{2}(\boldsymbol {u}\boldsymbol {\nabla })\bar {\rho }$, $\varepsilon ^{4}({1}/{Sc\,Re_{\ast }})Ri_{\ast }\nabla _{\perp }^{2}\bar {\rho }^{\prime \prime }$ and $\varepsilon ^{3} D_{t}\nabla _{\perp }^{2}((\boldsymbol {u}\boldsymbol {\nabla })w)$, on the right-hand side of (3.1) are $O(\varepsilon ^{7})$ and above, these terms will also be neglected in our further analysis. Upon these simplifications, recalling that $w=O(\varepsilon ^{2})$, we re-write (3.1) explicitly pulling out $\varepsilon$

(3.4)\begin{align} \underbrace{D_{t}^{2}\partial_{zz}^2w-U^{\prime\prime}D_{t} \partial_{x}w}_{O(\varepsilon^{4}{terms})} &= \varepsilon^{5}D_{t}\partial_{z}\partial_{x} [u\partial_{x}u+w\partial_{z}u]-\varepsilon^{5}Ri_{{\ast}}N^{2}\nabla_{{\perp}}^{2}w+\cdots \nonumber\\ &\quad +\frac{\varepsilon^{5}}{Re_{{\ast}}}D_{t}w^{iv}- \varepsilon^{6}D_{t}^{2}\nabla_{{\perp}}^{2}w. \end{align}

Since $v= O(\varepsilon ^{3})\ \mbox {and}\ u=O(\varepsilon ^{1})$, it is easy to see that, by virtue of the continuity equation, $\partial _{x}u= -\partial _{z}w + O(\varepsilon ^{4})$. To solve (3.4) we adopt a moving coordinate frame $x\to x-ct$ and use standard multiple-scale method by introducing fast and slow non-dimensional independent variables

(3.5a,b)\begin{equation} Z_{1}=\varepsilon z,\quad \tau=\varepsilon t, \end{equation}

where we recall that $c$ is the speed of perturbations in the long-wave limit which will be specified later. Then the squared material derivative $D_{t}^{2}$ and $\partial _{zz}^{2}$ take the form

(3.6a,b)\begin{equation} D_{t}^{2}=\frac{V_{0}^{2}}{L^{2}}[(U-c)\partial_{x}+\varepsilon\partial_{\tau}]^{2}= \frac{V_{0}^{2}}{L^{2}}D_{t}^{2}\sim\varepsilon^{2}D_{t}^{2},\quad \partial_{zz}^{2}=\frac{1}{d^{2}}(\partial_{z}+\varepsilon\partial_{Z_{1}})^{2}\approx \frac{1}{d^{2}}\partial_{zz}^{2}\sim\partial_{zz}^{2}. \end{equation}

Equation (3.4) for the vertical component of velocity, $w$, reduces to

(3.7)\begin{equation} D_{t}^{2}\partial_{zz}^2w-U^{\prime\prime}D_{t}\partial_{x}w=\varepsilon D_{t}\partial_{z}\partial_{x} [u\partial_{x}u+w\partial_{z}u]-\varepsilon Ri_{{\ast}} N^{2}\nabla_{{\perp}}^{2}w+\frac{\varepsilon}{Re_{{\ast}}} D_{t}w^{\prime\prime\prime\prime}-\varepsilon^{2}D_{t}^{2}\nabla_{{\perp}}^{2}w, \end{equation}

where

(3.8a,b)\begin{equation} D_{t}=(U-c)\partial_{x}+\varepsilon\partial_{\tau},\quad \partial_{z}=\partial_{z}+\varepsilon\partial_{Z_{1}}. \end{equation}

We recall that the mean flow $U(z)$ and stratification $N(z)$ are assumed to be localized in the boundary layer and, correspondingly, to depend only on the fast scale $z$.

First, we impose the no-flux condition at $z=0$. We will deal with the complete boundary conditions at the boundary in § 3.4. We also require vanishing of the velocity as $Z_{1}\to \infty$. In addition, we introduce ‘inner boundary conditions’, the condition of matching at the outer boundary of the boundary layer. Thus,

(3.9a–c)\begin{equation} w(z=0)=0;\quad w(z\to \infty)=const=w(Z_{1}\to 0),\quad w({Z_{1}\to \infty})\to 0. \end{equation}

We will seek an asymptotic solution to the boundary value problem (3.7), (3.8a,b) and (3.8a,b) employing power series in $\varepsilon$ (2.22). On finding the solution to (3.7) for $w$ at a certain order in $\varepsilon$, we find $u,v$ and $p$ with the corresponding accuracy from the basic equations

(3.10a)$$\begin{gather} D_{t}u+wU^{\prime}+p_{x}={-}\varepsilon(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})u+ \frac{\varepsilon}{Re_{{\ast}}}u^{\prime\prime}, \end{gather}$$

(3.10b)$$\begin{gather}D_{t}v+p_{y}={-}\varepsilon(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})v+ \frac{\varepsilon}{Re_{{\ast}}}v^{\prime\prime}, \end{gather}$$

(3.10c)$$\begin{gather}\partial_{x}u+\partial_{y}v+\partial_{z}w=0. \end{gather}$$

The solutions for $u$, $v$ and $p$, accurate to a corresponding order in $\varepsilon$, are further used for the derivation of the next order terms for $w$, the cycle is repeated as many times as necessary.

At the first step, on substituting (2.22) into (3.7) and setting $\varepsilon =0$, we see that in the leading-order nonlinearity, stratification and viscous dissipation drop out. Taking into account (3.8a,b) we obtain for $w_{2}$ the long-wave limit of the Rayleigh equation

(3.11)\begin{equation} (U-c)\partial_{xx}[(U-c)w_{2}^{\prime\prime}-U^{\prime\prime}w_{2}]=0, \end{equation}

where derivatives with respect to the fast variable $z$ are denoted by primes. It can be easily seen from (3.11) that the $x,y$ and $z$ dependencies can be separated. Assuming the disturbance to be localized or periodic in the streamwise direction, the general solution to equation (3.11) is convenient to present in the form

(3.12)\begin{equation} w_{2}=[f(x,y,Z_{1})\ast\partial_{x}A(x,y,\tau)]\cdot(U(z)-c), \end{equation}

where $\ast$ designates the convolution of two functions

(3.13)\begin{equation} \varphi\ast\psi=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\varphi({x}^{\prime},{y}^{\prime}) \psi(x-{x}^{\prime},y-{y}^{\prime})\,\mathrm{d}{x}^{\prime}\,\mathrm{d}{y}^{\prime}. \end{equation}

Here, as will be made explicit a few lines below, $A(x,y,\tau )$ is the amplitude of the $x$-component of velocity perturbation, while the arbitrary function $f(x,y,Z_{1})$ in (3.12) is the general representation of a function of $(x,y,Z_{1},\tau )$ localized or periodic in $(x,y)$. The fundamental properties of the presentation (3.12) become more transparent on performing the Fourier transform of ansatz (3.12). The standard way of representing a function of $(x,y,Z_{1},\tau )$ is to decompose it into a set of spatial orthogonal functions with time-dependent amplitudes. We chose Fourier in $(x,y)$ and a particular function $f$ specifying each Fourier mode that depends on the slow vertical spatial variable $Z_{1}$. The specific dependence on the slow spatial scale $f(Z_{1})$ will be found a few lines below. The boundary condition $w_{2}(z=0)=0$ specifies the eigenvalue $c$

(3.14)\begin{equation} c=U|_{z=0}\equiv U_0. \end{equation}

The mode we are considering is to leading order a vorticity wave, modified in the next orders by stratification and viscosity. To leading order its speed is the mean flow velocity at $z=0$, which usually is its maximal value for the typical ‘no-stress’ flows. In any case, it plays no role in our further analysis since its only significance is in specifying the reference frame, the surface velocity $U_0$ will be removed by the Galilean transformation at the next step.

To proceed further, we first find the other components of perturbation wave field from (3.10) using the asymptotic expansion (2.22). On substituting the leading-order solution for $w_{2}$ into (3.10) we find the other components of the perturbation field

(3.15a)$$\begin{gather} u_{1}={-}U^{\prime}(f\ast A), \end{gather}$$

(3.15b)$$\begin{gather}v_{2}=0, \end{gather}$$

(3.15c)$$\begin{gather}p_{2}=0. \end{gather}$$

Note that (3.15a) clarifies the physical sense of amplitude $A$, it is indeed, up to a factor $-U^{\prime }$, the amplitude of the $x$-component of perturbation velocity. The above relations show that to leading order the motion is extremely simple: the particles of the vorticity wave motion just oscillate in the streamwise direction, while the vertical and, especially, spanwise velocities and pressure perturbations are much smaller in the adopted long-wave approximation. This peculiar feature is specific for long vorticity waves (see Voronovich et al. Reference Voronovich, Shrira and Stepanyants1998). Such a simplicity of the motion of interest in the leading order is the key element enabling for a remarkably simple description of the 2-D nonlinear dynamics of vorticity waves which we will discuss below.

Substituting expressions (3.12) for $w_{2}$ and (3.15a) for $u_1$ into (3.7) we get an equation for $w_{3}$

(3.16)\begin{align} w_{3}^{\prime\prime}-\frac{U^{\prime\prime}}{U-c}w_{3} &={-}M\frac{U^{\prime\prime}}{U-c}+\frac{P}{Re_{{\ast}}} \frac{U^{\prime\prime\prime\prime}}{U-c}-\hat{G}_{4}[Q]\frac{N^{2}(z)}{U-c} \nonumber\\ &\quad -T\frac{[(U-c)^{2}]^{\prime}}{U-c}+R \frac{[(U^{\prime})^{2}-(U-c)U^{\prime\prime}]^{\prime}}{U-c}, \end{align}

where

(3.17a–e)\begin{gather} \begin{gathered} M=(f\ast A_{\tau}),\quad P=(f\ast A),\quad Q=(f\ast A_{x}),\quad T=(f_{Z_{1}}\ast A_{x}),\nonumber\\ R=(f\ast A)(f\ast A_{x}), \end{gathered} \end{gather}

and $\hat {G}_{4}$ is an integral operator which in the Fourier space is equivalent to multiplication by ${(k_{x}^{2}+k_{y}^{2})}/{k_{x}^{2}}$. The subscripts $x$,$\tau$ and $Z_{1}$ stand for the corresponding partial derivatives. Recall that the buoyancy frequency $N(z)$ is confined to the boundary layer, i.e.

(3.18)\begin{equation} N(z)\rightarrow 0, \quad\mbox{as}\ z\rightarrow\infty . \end{equation}

At the boundary $N(z=0)=N_{0}$. Usually, $N_{0}$ is the maximum of $N(z)$, but this point is immaterial for our consideration. The general solution to the inhomogeneous Rayleigh equation (3.16) can be written as

(3.19)\begin{align} w_{3} &= M-T(U-c)\int^{\infty}_{z}\mathrm{d}\xi+B(U-c)\int^{\infty}_{z} \frac{\mathrm{d}\xi}{(U-c)^{2}} \nonumber\\ &\quad +\frac{P}{Re_{{\ast}}}(U-c)\int^{\infty}_{z} \left[\frac{U^{\prime\prime\prime}}{(U-c)^{2}}\right]\mathrm{d}\xi- \hat{G}_{4}[Q](U-c)\int_{z}^{\infty}\left[\frac{S(\xi)}{(U-c)^{2}}\right]\mathrm{d}\xi-RU^{\prime}. \end{align}

The contribution of stratification is given by the term with $\hat {G}_{4}[Q]$, where function $S(z)$ in the integrand is proportional to the difference of the local reference density and that of the homogeneous fluid; it is given by the integral

(3.20)\begin{equation} S(z)=\int_{z}^{\infty}N^{2}(\xi)\,\mathrm{d}\xi. \end{equation}

In (3.19) $B$ is an arbitrary constant specifying the ‘amplitude’ of the diverging fundamental solution to the homogeneous Rayleigh equation in the long-wave limit. To eliminate singularity in the integrals in (3.19) we chose $B$ in such a way that the equation for $w_{3}$ takes the form

(3.21)\begin{align} w_{3} &= M-T(U-c)\int_{z}^{\infty}[1-Y(U-c)^{{-}2}]\,\mathrm{d}\xi+ \frac{P}{Re_{{\ast}}}(U-c)\int_{z}^{\infty} \left[\frac{U^{\prime\prime\prime}}{(U-c)^{2}}\right]\mathrm{d}\xi \nonumber\\ &\quad -\hat{G}_{4}[Q](U-c)\int_{z}^{\infty} \left[\frac{S(z)}{(U-c)^{2}}\right]\mathrm{d}\xi-RU^{\prime} , \end{align}

where the integration constant $Y$ has been chosen to be

(3.22)\begin{equation} Y=\lim_{z\rightarrow\infty}(U-c)^{2}=c^{2}=U_{0}^{2}. \end{equation}

To evaluate the singular integrals above in (3.21) we assume that $U^{\prime \prime \prime }(0)\ne 0$ and $U^{\prime }(0)\ne 0$ and expand $U-c$ near the boundary as $U^{\prime }(0)z=U_0^{\prime }z$.

The derived main-deck solution (3.21) ought to be matched with the lower-deck solution, the matching is carried out in Appendix A. Here, we just note that, upon formally applying to solution (3.21), the no-flux condition at the boundary, i.e. requiring $w|_{z=0,Z_{1}=0}=0$ we get, after some algebra, a closed equation for the amplitude $A$ containing so far unspecified function $f(Z_{1})$ of the slow vertical variable $Z_{1}$

(3.23)\begin{align} & (f(0)\ast A_{\tau})-U^{\prime}(0)\left(\left(f(0)\ast A\right) \left(f(0)\ast A_{x}\right)\right)+\left(\frac{U_{0}^{2}}{U^{\prime}(0)}\right) \left(f_{Z_{1}}(0)\ast A_{x}\right) \nonumber\\ &\quad -\frac{1}{2}U_{0}Ri\hat{G}_{2}[f(0)\ast A_{x}]+ \frac{1}{Re_{{\ast}}}\left(\frac{U^{\prime\prime\prime}(0)}{U^{\prime}(0)}\right) \left(f(0)\ast A\right)=0 , \end{align}

where $f(0)\equiv f(x,y,Z_{1}=0)$. The a posteriori justification of our use of the no-flux boundary condition is provided in Appendix A. Of course, the small $z$ expansion becomes invalid in the immediate vicinity of the boundary, a composite uniformly valid expansion is derived in Appendix A. Operator $\hat {G}_{2}$, which appears in (3.23), is related to $\hat {G}_{4}$ as follows: in the Fourier space $\hat {G}_{2}$ and $\hat {G}_{4}$ correspond, respectively, to $k_y^2/k_x^2$ and $(k_y^2 +k_x^2)/k_x^2$. Going from (3.21) to (3.23) we have subtracted 1 from the kernel in $\hat {G}_{4}$ to get $\hat {G}_{2}$, which means that the small correction to the long-wave velocity caused by stratification $c_{1}=-\frac {1}{2}U_{0}Ri$ was taken into account by adjusting the speed of the moving coordinate frame, the adjustment eliminates the term $c_{1}A_{x}$ appearing in the old frame due to correction to the long-wave velocity.

3.3. Matching the outer flow

In this subsection in order to obtain our final evolution equation for the 3-D perturbation amplitude $A$ we proceed to find the unknown function $f(Z_1)$ for the outer flow slow vertical motion and then substitute it into (3.23). To find the as yet unspecified dependence of the found solution (3.12) and (3.21) on the slow vertical variable $Z_{1}$, we first proceed to the next order in $\varepsilon$ in equation (3.7) for $w_{2}$. Following the same asymptotic procedure and using (3.12) and (3.21) we express $u, v, p$ in terms of amplitude $A$

(3.24)\begin{align} u_{2} &={-}Y(f_{Z_{1}}\ast \nabla_{{\perp}}^{{-}2}A_{xx}) (U-c)^{{-}1}+(f_{Z_{1}}\ast A)U^{\prime}\int_{z}^{\infty} \left[1-Y(U-c)^{{-}2}\right]\mathrm{d}\xi \nonumber\\ &\quad +\frac{1}{2}(f\ast A)^{2}U^{\prime\prime}+\hat{G}_{4}[P] \frac{S(z)}{U-c}+\hat{G}_{4}[P]U^{\prime}\int_{z}^{\infty} \left(\frac{S(z)}{(U-c)^{2}}\right)\mathrm{d}\xi \nonumber\\ &\quad -\frac{\hat{P}U^{\prime}}{Re_{{\ast}}}\int_{z}^{\infty} \left(\frac{U^{\prime\prime\prime}}{(U-c)^{2}}\right)\mathrm{d}\xi -\frac{\hat{P}}{Re_{{\ast}}}\frac{U^{\prime\prime\prime}}{(U-c)}, \end{align}

(3.25)$$\begin{align} v_{3}&={-}Y(f_{Z_{1}}\ast \nabla_{{\perp}}^{{-}2}A_{xy})(U-c)^{{-}1}+\hat{G}_{3}[(f\ast A)]S(z)(U-c)^{{-}1}, \end{align}$$

(3.26)$$\begin{align}p_{3}&=Y(f_{Z_{1}}\ast \nabla_{{\perp}}^{{-}2}A_{xx})-\hat{G}_{4}[(f\ast A)]S(z), \end{align}$$

where $\hat {G}_{3}$ is a pseudo-differential operator specified below, while $\hat {P}=(f\ast \partial _{x}^{-1}A)$. In the Fourier space the $\textbf{r}$-space operators $\hat {G}_{3}$ correspond to $({k_{y}^{2}}/{k_{x}^{2}}) ({k^{2}}/{k_{x}^{2}})$. Further on we will use symbol $\Leftrightarrow$ to denote such correspondence. Expressions for $P$, $Y$ and $S(z)$ are given by (3.17a–e), (3.22) and (3.20).

Now consider the next-order term for the vertical velocity, $w_{4}$. After some algebra it can be brought to the form

(3.27)\begin{align} \partial_{x}[(U-c)w_{4}^{\prime\prime}-U^{\prime\prime}w_{4}] &={-}(f_{Z_{1}Z_{1}}\ast A_{xx})(U-c)^{2} \nonumber\\ &\quad -(f\ast\nabla_{{\perp}}^{2}A_{xx})(U-c)^{2}+F(x,y,\tau,z,Z_{1}), \end{align}

where $F(x,y,\tau,z,Z_{1})$, being specified by a very bulky expression, is not given here. Crucially, it tends to zero as $z\to \infty$ faster than $|z|^{-1}$. In contrast, for an arbitrary function of $f(x,y,Z_{1})$ the first two terms on the right-hand side of equation (3.27) do not tend to zero as $z\to \infty$. As a result, the integration of (3.27) yields secular growth of the correction $w_{4}$ as $z\to \infty$, which does not allow the matching condition $w(z\to \infty,Z_{1})\to \mbox {const}$ to be satisfied. We put these secular terms to zero, which gives us an equation determining the dependence of function $f$ on the slow variable $Z_{1}$

(3.28)\begin{equation} (f_{Z_{1}Z_{1}}\ast A_{xx})+(f\ast\nabla_{{\perp}}^{2}A_{xx})=0. \end{equation}

The equation is complemented by the boundary condition at infinity

(3.29)\begin{equation} w(Z_{1}\to\infty) =0. \end{equation}

To find $f(Z_{1})$ we, exploiting the homogeneity of the problem with respect to horizontal coordinates, perform the Fourier transform with respect to $x,y$ in the boundary problem (3.28) and (3.29). Making the Fourier transform of the convolution and pulling out the terms with amplitude $A$, we get the boundary value problem

(3.30)\begin{equation} \partial_{Z_{1}Z_{1}}^{2}\hat{f}({\boldsymbol k})-k^{2}\hat{f}({\boldsymbol k})=0, \end{equation}

with the boundary conditions

(3.31a,b)\begin{equation} \hat{f}({\boldsymbol k})|_{Z_{1} \to \infty}=0,\quad \hat{f}({\boldsymbol k})|_{Z_{1}=0}=1. \end{equation}

Here, $\hat {f}({\boldsymbol k},Z_{1})$ is the Fourier transform of $f(x,y,Z_{1})$

(3.32)\begin{equation} f(x,y,Z_{1})=\frac{1}{4{\rm \pi}^{2}}\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty}\hat{f}(k_{x},k_{y},Z_{1}) \exp({{\rm i}(k_{x}x+k_{y}y)})\,\mathrm{d}k_{x}\,\mathrm{d}k_{y}, \end{equation}

and $k=|\boldsymbol k|$, $k^{2}=k_{x}^{2}+k_{y}^{2}$. The simplest normalization, $\hat {f}({\boldsymbol k})|_{Z_{1}=0}=1$, has been introduced for convenience to normalize the motion in the outer deck $\hat {f}(k_{x},k_{y},Z_{1})$ near the boundary. As one could easily anticipate, the motion in the outer layer is potential and satisfies the Laplace equation (3.30). Hence, it is easy to find the solution of the boundary value problem (3.30), (3.31a,b) satisfying the boundary condition (3.31a,b)

(3.33)\begin{equation} \hat{f}(\boldsymbol{k}, Z_{1})=\exp({-kZ_{1}})\quad Z_{1}\to \infty. \end{equation}

Next, we designate $\partial _{Z_{1}}\hat {f}(\boldsymbol {k},0)\equiv \check {F}(k_{x},k_{y})$ and take into account that, at the boundary, ${f}(0)=\delta (x)\delta (y)$. The solution to (3.30) readily yields the kernel of the integral operator $\check {F}(k)=-k$. Substituting these findings into (3.23), we obtain the nonlinear evolution equation for the amplitude of 3-D long-wave perturbations in the distinguished limit

(3.34)\begin{equation} A_{\tau}-\alpha_{n} AA_{x}-\beta_{1}\hat{G}_{1}[A_{x}]- \beta_{2}\hat{G}_{2}[A_{x}]+\gamma A=0, \end{equation}

where the non-local operators $\hat {G}_{1}$ and $\hat {G}_{2}$ are

(3.35)$$\begin{gather} \hat{G}_{1}[\varphi(\boldsymbol{r})]=\frac{1}{4{\rm \pi}^{2}}\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty}|\boldsymbol{k}|\varphi(\boldsymbol{r}_{1}) \exp({{\rm i}\boldsymbol{k}(\boldsymbol{r}-\boldsymbol{r}_{1})})\,\mathrm{d}\boldsymbol{k}\,\mathrm{d}\boldsymbol{r}_{1}, \end{gather}$$

(3.36)$$\begin{gather}\hat{G}_{2}[\psi(\boldsymbol{r})]=\frac{1}{4{\rm \pi}^{2}} \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \frac{k_y^{2}}{k_{x}^{2}}\psi(\boldsymbol{r}_{1})\exp({{\rm i}\boldsymbol{k} (\boldsymbol{r}-\boldsymbol{r}_{1})})\,\mathrm{d}\boldsymbol{k}\,\mathrm{d}\boldsymbol{r}_{1}. \end{gather}$$

Recall that $k=|\boldsymbol {k}|=\sqrt {k_{x}^{2}+k_{y}^{2}}$. The coefficients are specified by the parameters of the flow at the surface: $\alpha _{n}=U^{\prime }(0), \beta _{1}={U_{0}^{2}}/{U^{\prime }(0)}, \gamma =({1}/{Re_{\ast }}){U^{\prime \prime \prime }(0)}/{U^{\prime }(0)}, \beta _{2}=U_0Ri/2$. The explicit account of viscous dissipation and weak stratification yield, respectively, the Rayleigh friction-type term $\gamma A$ and a very specific dispersion term $\beta _{2}\hat {G}_{2}[A_{x}]$ due to buoyancy. Coefficients $\alpha _{n}$ and $\beta _{1}$ of (3.34) can be removed by re-scaling to obtain the ‘universal’ form of the evolution equation. The equation is universal in the sense that the specific profiles of the boundary layer and stratification are immaterial. The residual specificity of the boundary layer retained in the coefficients $\alpha _{n}, \beta _{1}, \gamma, \beta _{2}$ in (3.34) can be further reduced by re-scaling the variables

(3.37a–d)\begin{equation} \tau_{1}=\tau U^{\prime}(0),\quad d=\frac{U_{0}}{U^{\prime}(0)},\quad x_{1}=x/d,\quad A_{1}={-}\frac{1}{d}A, \end{equation}

which upon setting $d=1$ and dropping the subscripts yields the ultimate form of the equation

(3.38)\begin{equation} A_{\tau}+AA_{x}-\hat{G}_{1}[A_{x}]-\tilde{\beta}_{2}\hat{G}_{2}[A_{x}]+\tilde{\gamma}A=0. \end{equation}

The non-local operators $\hat {G}_{1}$ and $\hat {G}_{2}$ remain the same and are given by (3.35) and (3.36). The only remaining two coefficients

(3.39a–c)\begin{equation} \tilde{\beta}_{2}=\frac{N_{0}^{2}}{2\varepsilon^2(U^{\prime}(0))^{2}},\quad \tilde{\gamma}= \frac{1}{Re_{{\ast}}} \frac{U^{\prime\prime\prime}(0)}{(U^{\prime}(0))^{2}},\quad Re_{{\ast}}=\varepsilon^{2}Re, \end{equation}

characterize the relative importance of the novel dispersion due to weak stratification and the Rayleigh friction-type term compared with the nonlinear one in the evolution equation. In this context the very weak ${Ri}=O(\varepsilon )$ stratification proves to be essential. Although the flows with ${Ri}<1/4$ might be linearly unstable, in our context linear instabilities play no role since we focus on linearly stable situations; linear instabilities will be considered elsewhere.

3.4. Viscous sub-layer

3.4.1. Divergence of the asymptotic expansion and Tollmien's rescaling

Consider more closely behaviour of the main-deck solutions $u_{2}$ and $v_{3}$ given by (3.24) and (3.25) as $z\to 0$. Near the boundary as $z\to 0$, we apply Taylor series to expand $(U-c)\approx U^{\prime }(0)z$ to obtain in physical space both streamwise and spanwise velocities

(3.40)$$\begin{gather} u_{2}={-}\frac{U_{max}^{2}}{U^{\prime}(0)}(f_{Z_1}\ast \nabla_{{\perp}}^{{-}2}A_{xx})\frac{1}{z}+\frac{U_{max}^{2}}{U^{\prime}(0)}(f_{Z_1}\ast A)\frac{1}{z}+\frac{1}{2}(f\ast A)^{2}U^{\prime\prime}(0), \end{gather}$$

(3.41)$$\begin{gather}v_{3}={-}\frac{U_{max}^{2}}{U^{\prime}(0)}(f_{Z_{1}}\ast \nabla_{{\perp}}^{{-}2}A_{xy})\frac{1}{z}+\frac{1}{2}U_{0}Ri\hat{G}_{3}(f\ast A) \frac{1}{z}. \end{gather}$$

First, it is important to note that the last four terms due to stratification and viscosity of (3.24) all cancel each other. In contrast to the streamwise perturbation velocity $u_{2}$, the spanwise component yields a new extra singular term due to non-zero stratification.

Let us take a closer look at the divergence of the streamwise velocity component $u_{2}$ given by (3.40), by carrying out some algebra we simplify the first two terms on the right-hand side with singularities at $z=0$. To this end we transform them into the Fourier space

(3.42) \begin{align} & -\frac{U_{max}^{2}}{U^{\prime}(0)}\frac{k_{x}^{2}}{k_{x}^{2}+k_{y}^{2}} (\,\hat{f}_{Z_{1}}\cdot\hat{A})\frac{1}{z}+\frac{U_{max}^{2}}{U^{\prime}(0)} (\,\hat{f}_{Z_{1}}\cdot\hat{A})\frac{1}{z}=\frac{U_{max}^{2}}{U^{\prime}(0)} \left[\frac{k_{x}^{2}+k_{y}^{2}-k_{x}^{2}}{k_{x}^{2}+k_{y}^{2}}\right] (\,\hat{f}_{Z_{1}}\cdot\hat{A})\frac{1}{z}, \nonumber\\ &\quad =\frac{U_{max}^{2}}{U^{\prime}(0)}\left[\frac{k_{y}^{2}}{k_{x}^{2}+k_{y}^{2}}\right] (\,\hat{f}_{Z_{1}}\cdot\hat{A})\frac{1}{z}. \end{align}

Substituting this simplified form of the singular terms and taking into account their inverse Fourier transform, we obtain a more compact form of the expression for $u_{2}$

(3.43)\begin{equation} u_{2}=\frac{U_{max}^{2}}{U^{\prime}(0)}(f_{Z_1}(0)\ast \nabla_{{\perp}}^{{-}2}A_{yy})\frac{1}{z}+ \frac{1}{2}(f\ast A)^{2}U^{\prime\prime}(0). \end{equation}

It is easy to see that uniformity of the main-deck asymptotic expansion breaks down as $z\to 0$. Indeed, according to (3.43) and (3.41) unless the spanwise component dependence vanishes, velocity components $u_{2}$ and $v_{3}$ diverge

(3.44a,b)\begin{equation} u_{2}\sim 1/z,\quad v_{3}\sim 1/z. \end{equation}

To get rid of the singularities, in Appendix A we re-scale the variables appropriately and solve the Navier–Stokes equations in the region immediately adjacent to the boundary (the lower deck) and then match these solutions with the solutions already obtained for the main deck. The main results of the detailed analysis carried out in Appendix A can be summarized as follows: (i) under the chosen scaling a uniformly valid explicit asymptotic solution has been found for the no-stress boundary layers, although in a quite complicated form; (ii) in the context of derivation of evolution equation the specific form of the solution in the viscous sub-layer is immaterial, since under adopted scaling the inner solution anyway does not affect the motion in the main deck and, hence, does not contribute to the evolution equation.