2.1. Linear invariants and waves

In the linear approximation let the system obey the equation

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{t}=\unicode[STIX]{x1D64D}\,\boldsymbol{f},\end{eqnarray}$$

where $\unicode[STIX]{x1D647}$ and $\unicode[STIX]{x1D64D}$ are some linear matrix operators, $\unicode[STIX]{x1D647}$ being independent of the time $t$ . For simplicity let the state vector $\boldsymbol{f}$ satisfy the homogeneous boundary condition at time-independent boundary ${\it\Gamma}$ ,

(2.2a ) $$\begin{eqnarray}\boldsymbol{B}\,\boldsymbol{f}|_{{\it\Gamma}}=0,\end{eqnarray}$$

and the initial condition,

(2.2b ) $$\begin{eqnarray}\boldsymbol{f}|_{t=0}=\boldsymbol{f}_{I}.\end{eqnarray}$$

The linear operator $\boldsymbol{B}$ can depend on time; here and below the subscript ‘ $I$ ’ denotes initial values.

Our consideration is based on the two following assumptions.

(i) There exists a time-independent linear operator $\unicode[STIX]{x1D648}\neq 0$ that nullifies the right-hand side of (2.1), i.e.

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D648}\unicode[STIX]{x1D64D}=0.\end{eqnarray}$$

In this case in view of the time independence of $\unicode[STIX]{x1D647}$ and $\unicode[STIX]{x1D648}$ , (2.1) possesses the linear invariant

(2.4) $$\begin{eqnarray}{\it\Lambda}={\it\Lambda}(\boldsymbol{x})=\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f},\end{eqnarray}$$

since

(2.5) $$\begin{eqnarray}\frac{\partial {\it\Lambda}}{\partial t}=\frac{\partial (\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f})}{\partial t}=0.\end{eqnarray}$$

The operator $\unicode[STIX]{x1D648}$ will be called the invariant operator. The value of the invariant ${\it\Lambda}$ is determined by the initial conditions (2.2b ):

(2.6) $$\begin{eqnarray}{\it\Lambda}=\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{I}.\end{eqnarray}$$

(ii) The equations (2.1) and (2.2a ) possess harmonic wave solutions of the form

(2.7) $$\begin{eqnarray}\boldsymbol{f}_{w}=\boldsymbol{A}_{w}(\boldsymbol{x})\text{e}^{\text{i}{\it\sigma}t},\quad {\it\sigma}=\text{const.}\neq 0.\end{eqnarray}$$

The key point is that for the wave solution (2.7) the linear invariant (2.4) is identically zero since the product $\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}$ is linear and time-independent; therefore we have

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{w}=\text{e}^{\text{i}{\it\sigma}t}\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\boldsymbol{A}_{w}=\text{const.}\quad \Rightarrow \quad \unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{w}=0.\end{eqnarray}$$

Physically, this means that the waves do not transport and affect the invariant  ${\it\Lambda}$ .

In this case it is natural to represent the solution $\boldsymbol{f}$ to the problem (2.1) and (2.2) as a sum of two solutions:

(2.9) $$\begin{eqnarray}\boldsymbol{f}=\bar{\boldsymbol{f}}(\boldsymbol{x})+\tilde{\boldsymbol{f}}(\boldsymbol{x},t).\end{eqnarray}$$

The stationary solution $\boldsymbol{f}=\bar{\boldsymbol{f}}(\boldsymbol{x})$ possesses the non-zero invariant ${\it\Lambda}$ and satisfies the equations

(2.10a−c ) $$\begin{eqnarray}\unicode[STIX]{x1D64D}\bar{\boldsymbol{f}}=0,\quad \unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\bar{\boldsymbol{f}}={\it\Lambda},\quad \boldsymbol{B}\bar{\boldsymbol{f}}|_{{\it\Gamma}}=0,\end{eqnarray}$$

and the solution $\tilde{\boldsymbol{f}}$ obeys (2.1) and has the zero invariant, i.e. 

(2.11a−d ) $$\begin{eqnarray}\unicode[STIX]{x1D647}\tilde{\boldsymbol{f}}_{t}=\unicode[STIX]{x1D64D}\tilde{\boldsymbol{f}},\quad \unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\tilde{\boldsymbol{f}}=0,\quad \boldsymbol{B}\tilde{\boldsymbol{f}}|_{{\it\Gamma}}=0,\quad \tilde{\boldsymbol{f}}_{I}=\boldsymbol{f}_{I}-\bar{\boldsymbol{f}}.\end{eqnarray}$$

If the system of waves (2.7) is complete, then the solution $\tilde{\boldsymbol{f}}$ is a linear superposition of these waves; therefore it is natural to call $\tilde{\boldsymbol{f}}$ the wave component. We now introduce the time average

(2.12) $$\begin{eqnarray}\langle a\rangle =\lim _{T_{0}\rightarrow \infty }\frac{1}{T_{0}}\int _{0}^{T_{0}}a\,\text{d}t,\end{eqnarray}$$

and assume the operators $\unicode[STIX]{x1D64D}$ and $\boldsymbol{B}$ to be such that

(2.13a,b ) $$\begin{eqnarray}\langle \unicode[STIX]{x1D64D}\tilde{\boldsymbol{f}}\rangle =\unicode[STIX]{x1D64D}\langle \tilde{\boldsymbol{f}}\rangle =0,\quad \langle \boldsymbol{B}\tilde{\boldsymbol{f}}\rangle =\boldsymbol{B}\langle \tilde{\boldsymbol{f}}\rangle =0.\end{eqnarray}$$

This takes place, for example, if $\unicode[STIX]{x1D64D}$ and $\boldsymbol{B}$ are either time-independent or include only time derivatives (as in the examples below). In this case, averaging (2.11) in time, we obtain

(2.14a−c ) $$\begin{eqnarray}\unicode[STIX]{x1D64D}\langle \tilde{\boldsymbol{f}}\rangle =0,\quad \unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\langle \tilde{\boldsymbol{f}}\rangle =0,\quad \boldsymbol{B}\langle \tilde{\boldsymbol{f}}\rangle |_{{\it\Gamma}}=0.\end{eqnarray}$$

The system (2.14) coincides with (2.10) for ${\it\Lambda}=0$ . Therefore if the problem (2.10) has a unique solution, then we have from (2.14) that

(2.15) $$\begin{eqnarray}\langle \tilde{\boldsymbol{f}}\rangle =0,\end{eqnarray}$$

i.e. the wave solution oscillates near the stationary one. If the waves (2.7) are able to radiate away from the initial perturbation domain, then the wave component $\tilde{\boldsymbol{f}}$ gradually decays at any fixed point $\boldsymbol{x}$ and the solution (2.9) tends to the stationary component $\bar{\boldsymbol{f}}$ with increasing time – so-called wave adjustment takes place.

Thus the sufficient condition for the state vector $\boldsymbol{f}$ of the linear system to split into stationary and oscillating parts is that the system possesses linear invariants and a complete set of linear waves.

An alternative way of achieving the linear fast–slow splitting exploits the decomposition of the state vector into normal modes of the system considered. If the normal modes can be divided into two sets consisting, respectively, of slow (low-frequency) modes and fast (high-frequency) modes, then the state vector can be decomposed into the slow and fast components, which are the vector’s projections onto the slow and fast sets, respectively (Leith Reference Leith1980; Salmon Reference Salmon1998). This kind of splitting is used as a first step in initialization (see e.g. Leith Reference Leith1980; Daley Reference Daley1981; Temperton Reference Temperton1988). This way is more general than the way that we use because it works in the absence of linear invariants, too. The price paid for this universality is that dynamical equations should be represented in phase space (space of the normal modes) instead of the physical one. Such a representation may be useful for numerical modelling, but for our analytical treatment the equations in physical space are much more convenient. The splitting in physical space based on linear invariants also was used for initialization (Temperton Reference Temperton1988) and for derivation of balanced models (Lynch Reference Lynch1989) as applied to the rotating shallow-water (RSW) model.

2.2. Weak nonlinearity and balanced models

We now consider the perturbed model (2.1) written in a non-dimensional form:

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{t}=\unicode[STIX]{x1D64D}\,\boldsymbol{f}+{\it\varepsilon}\unicode[STIX]{x1D649}\boldsymbol{f},\quad {\it\varepsilon}\ll 1.\end{eqnarray}$$

The term ${\it\varepsilon}\unicode[STIX]{x1D649}\boldsymbol{f}$ is a small perturbation, $\unicode[STIX]{x1D649}$ is a linear or nonlinear operator, which can depend on ${\it\varepsilon}$ , and the boundary and initial conditions (2.2a,b ) are assumed to be ${\it\varepsilon}$ -independent. Applying the invariant operator $\unicode[STIX]{x1D648}$ to (2.16) one obtains

(2.17) $$\begin{eqnarray}\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{t}={\it\varepsilon}\unicode[STIX]{x1D648}\unicode[STIX]{x1D649}\boldsymbol{f},\end{eqnarray}$$

i.e. the ‘former’ linear invariant $\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}$ ceases to be invariant but changes slowly in time. There are two time scales in the system – the slow one ${\sim}1/{\it\varepsilon}$ and the fast one ${\sim}1$ . Therefore, a solution to the problem is sought in the form of a multiple-time-scale asymptotic expansion:

(2.18) $$\begin{eqnarray}\boldsymbol{f}=\boldsymbol{f}^{(0)}(\boldsymbol{x},t,T_{1},T_{2},\ldots )+{\it\varepsilon}\boldsymbol{f}^{(1)}(\boldsymbol{x},t,T_{1},T_{2},\ldots )+\cdots \,,\quad T_{n}={\it\varepsilon}^{n}t,\end{eqnarray}$$

where $T_{n},~n=1,2,\ldots$ , are the slow times.

Substitution of (2.18) into (2.16) gives in the lowest order:

(2.19) $$\begin{eqnarray}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{t}^{(0)}=\unicode[STIX]{x1D64D}\,\boldsymbol{f}^{(0)}.\end{eqnarray}$$

Applying the invariant operator $\unicode[STIX]{x1D648}$ to (2.19), one obtains that the quantity $\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}^{(0)}$ does not depend on the fast time $t$ :

(2.20) $$\begin{eqnarray}\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}^{(0)}={\it\Lambda}(\boldsymbol{x},T_{1}).\end{eqnarray}$$

We now represent $\boldsymbol{f}^{(0)}$ as a sum of slow and fast components (cf. (2.9)),

(2.21) $$\begin{eqnarray}\boldsymbol{f}^{(0)}=\bar{\boldsymbol{f}}^{(0)}(\boldsymbol{x},T_{1})+\tilde{\boldsymbol{f}}^{(0)}(\boldsymbol{x},t,T_{1}),\end{eqnarray}$$

which obey the equations (cf. (2.10) and (2.11))

(2.22a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D64D}\bar{\boldsymbol{f}}^{(0)}=0,\quad \unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\bar{\boldsymbol{f}}^{(0)}={\it\Lambda}(\boldsymbol{x},T_{1}),\end{eqnarray}$$

(2.23a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D647}\tilde{\boldsymbol{f}}_{t}^{(0)}=\unicode[STIX]{x1D64D}\tilde{\boldsymbol{f}}^{(0)},\quad \unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\tilde{\boldsymbol{f}}^{(0)}=0.\end{eqnarray}$$

Initial conditions for $\bar{\boldsymbol{f}}^{(0)}$ and $\tilde{\boldsymbol{f}}^{(0)}$ are determined from the equations following from (2.22) and (2.23):

(2.24a−c ) $$\begin{eqnarray}\unicode[STIX]{x1D64D}\bar{\boldsymbol{f}}_{I}^{(0)}=0,\quad \unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\bar{\boldsymbol{f}}_{I}^{(0)}=\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{I},\quad \tilde{\boldsymbol{f}}_{I}^{(0)}=\boldsymbol{f}_{I}-\bar{\boldsymbol{f}}_{I}^{(0)}.\end{eqnarray}$$

The function ${\it\Lambda}(\boldsymbol{x},T_{1})$ in (2.22b ) has yet to be found; the slow evolution is determined from boundedness of the first-order correction, which obeys the equation

(2.25) $$\begin{eqnarray}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{t}^{(1)}-\unicode[STIX]{x1D64D}\,\boldsymbol{f}^{(1)}=-\unicode[STIX]{x1D647}\,\boldsymbol{f}_{T_{1}}^{(0)}+\unicode[STIX]{x1D649}^{(0)}\boldsymbol{f}^{(0)},\end{eqnarray}$$

where $\unicode[STIX]{x1D649}^{(0)}$ is the lowest-order term in the expansion of the operator $\unicode[STIX]{x1D649}$ in ${\it\varepsilon}$ . We now apply the invariant operator to (2.25), taking into account (2.22) and (2.23):

(2.26) $$\begin{eqnarray}\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\,\boldsymbol{f}_{t}^{(1)}=-\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\bar{\boldsymbol{f}}_{T_{1}}^{(0)}+\unicode[STIX]{x1D648}\unicode[STIX]{x1D649}^{(0)}\boldsymbol{f}^{(0)}.\end{eqnarray}$$

It readily follows from (2.26) that the first-order correction $\boldsymbol{f}^{(1)}$ is bounded in the fast time $t$ if

(2.27) $$\begin{eqnarray}\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\bar{\boldsymbol{f}}_{T_{1}}^{(0)}-\langle \unicode[STIX]{x1D648}\unicode[STIX]{x1D649}^{(0)}\boldsymbol{f}^{(0)}\rangle =0,\end{eqnarray}$$

where the angle brackets denote time averaging (2.12).

Equation (2.27) determines the slow evolution of the zero-order invariant $\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\bar{\boldsymbol{f}}^{(0)}$ . To close the problem for $\bar{\boldsymbol{f}}^{(0)}$ , one should express $\langle \unicode[STIX]{x1D648}\unicode[STIX]{x1D649}^{(0)}\boldsymbol{f}^{(0)}\rangle$ in terms of $\bar{\boldsymbol{f}}^{(0)}$ ; in this case $\bar{\boldsymbol{f}}^{(0)}$ could be determined from (2.22a ) and (2.27). It is unlikely that the problem can be solved in the general case, but if one assumes that the wave component $\tilde{\boldsymbol{f}}^{(0)}$ can be neglected, at least, on times ${\sim}1/{\it\varepsilon}$ , i.e. 

(2.28) $$\begin{eqnarray}\langle \unicode[STIX]{x1D648}\unicode[STIX]{x1D649}^{(0)}\boldsymbol{f}^{(0)}\rangle =\unicode[STIX]{x1D648}\unicode[STIX]{x1D649}^{(0)}\bar{\boldsymbol{f}}^{(0)},\end{eqnarray}$$

then the time evolution of the slow component is described by the equation

(2.29) $$\begin{eqnarray}\unicode[STIX]{x1D648}\unicode[STIX]{x1D647}\bar{\boldsymbol{f}}_{T_{1}}^{(0)}-\unicode[STIX]{x1D648}\unicode[STIX]{x1D649}^{(0)}\bar{\boldsymbol{f}}^{(0)}=0.\end{eqnarray}$$

In addition, $\bar{\boldsymbol{f}}^{(0)}$ should obey the corresponding boundary conditions and equations (2.24a,b ), which filter out the fast waves from the initial field $\boldsymbol{f}_{I}$ . The set of equations (2.22a ), (2.29) and (2.24a,b ) and the boundary conditions determine the balanced model (in which the fast waves are absent) related to the system (2.16).

Of course, the important question is under what conditions (2.28) is valid. Analysis of various examples from fluid dynamics shows that (2.28) is fulfilled (as could be expected) if the fast waves are absent in initial fields or their group speeds greatly exceed the velocities typical of the slow component. The corresponding scenario of wave adjustment will be referred to as ‘fast radiation’ (see § 3). More surprisingly, (2.28) is valid also if the group speeds are of the order of or smaller than the slow velocities. In this case the waves are trapped by the slow component but their self-interaction gives no contribution to the left-hand side of (2.28), at least, on times ${\sim}1/{\it\varepsilon}$ . The scenario of such adjustment will be called ‘nonlinear trapping’ (see § 4). In § 5 we consider the situation when non-stationary boundary layers arise near the boundary ${\it\Gamma}$ or inside the motion domain at large times $t$ in the linear problem (2.1) and (2.2). Each boundary layer consists of an infinite number of very short harmonic waves (2.7) with vanishingly small group speeds; therefore the corresponding scenario of adjustment is called ‘imperfect splitting’.

In the theory developed we have used the basic equations in the operator form (2.1), (2.2) and (2.16). This is convenient for general theory but for ‘practical’ calculations in the examples considered below the matrix operators $\unicode[STIX]{x1D647}$ , $\unicode[STIX]{x1D64D}$ and $\unicode[STIX]{x1D649}$ are too cumbersome and we examine the equations written in ‘standard’ form. Operator forms of the models in use are given in appendix A.

3.1. Statement of the problem

The corresponding equations of motion are written in the form:

(3.1a,b ) $$\begin{eqnarray}\boldsymbol{u}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0.\end{eqnarray}$$

Here $\boldsymbol{u}=(u,v,w)$ is the velocity, with components $u,v,w$ along the axes $x,y,z$ , respectively; $p$ is the deviation of pressure from the hydrostatic one divided by the fluid density. At the free surface $z={\it\varsigma}$ and bottom $z=-H$ the dynamic and kinematic conditions are fulfilled:

(3.2a−c ) $$\begin{eqnarray}w|_{z={\it\varsigma}}={\it\varsigma}_{t}+(u{\it\varsigma}_{x}+v{\it\varsigma}_{y})_{z={\it\varsigma}},\quad p|_{z={\it\varsigma}}=g{\it\varsigma},\quad w|_{z=-H}=0.\end{eqnarray}$$

The velocity and elevation are given at the initial moment $t=0$ :

(3.2d ) $$\begin{eqnarray}(\boldsymbol{u},{\it\varsigma})_{t=0}=(\boldsymbol{u}_{I},{\it\varsigma}_{I}).\end{eqnarray}$$

In what follows all fields are assumed to decay at infinity, i.e. as $\sqrt{x^{2}+y^{2}}\rightarrow \infty$ .

3.2. Linear invariants and waves

In the linearized problem (3.1) and (3.2),

(3.3a−d ) $$\begin{eqnarray}\boldsymbol{u}_{t}=-\boldsymbol{{\rm\nabla}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad (w,p)_{z=0}=({\it\varsigma}_{t},g{\it\varsigma}),\quad w|_{z=-H}=0,\end{eqnarray}$$

the vorticity vector $\text{curl }\boldsymbol{u}=({\it\Omega}^{x},{\it\Omega}^{y},{\it\Omega}^{z})$ is conserved in view of (3.3a ), i.e. there are three linear invariants:

(3.4a−c ) $$\begin{eqnarray}{\it\Omega}^{x}=w_{y}-v_{z}={\it\Omega}_{I}^{x},\quad {\it\Omega}^{y}=u_{z}-w_{x}={\it\Omega}_{I}^{y},\quad {\it\Omega}^{z}=v_{x}-u_{y}={\it\Omega}_{I}^{z}.\end{eqnarray}$$

Following § 2, we represent the solution in the form (2.9), i.e. 

(3.5) $$\begin{eqnarray}(\boldsymbol{u},p,{\it\varsigma})=(\bar{\boldsymbol{u}},\bar{p},\bar{{\it\varsigma}})+(\tilde{\boldsymbol{u}},\tilde{p},\tilde{{\it\varsigma}}),\end{eqnarray}$$

where the stationary solution $(\bar{\boldsymbol{u}},\bar{p},\bar{{\it\varsigma}})$ obeys the equations

(3.6a−c ) $$\begin{eqnarray}\bar{p}=\bar{{\it\varsigma}}=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\bar{\boldsymbol{u}}=0,\quad \text{curl }\bar{\boldsymbol{u}}=({\it\Omega}^{x},{\it\Omega}^{y},{\it\Omega}^{z})\end{eqnarray}$$

and the boundary conditions

(3.6d ) $$\begin{eqnarray}\bar{w}|_{z=0,-H}=0.\end{eqnarray}$$

For the wave solution $(\tilde{\boldsymbol{u}},p,{\it\varsigma})$ we have

(3.7a−e ) $$\begin{eqnarray}\tilde{\boldsymbol{u}}_{t}=-\boldsymbol{{\rm\nabla}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=0,\quad \text{curl }\tilde{\boldsymbol{u}}=0,\quad (\tilde{w},p)_{z=0}=({\it\varsigma}_{t},g{\it\varsigma}),\quad \tilde{w}|_{z=-H}=0.\end{eqnarray}$$

To find the stationary solution, one obtains the vertical velocity $\bar{w}$ from the equation ${\rm\Delta}\bar{w}={\it\Omega}_{Iy}^{x}-{\it\Omega}_{Ix}^{y}$ and boundary conditions (3.6d ), and then the horizontal velocities from the equations ${\rm\Delta}_{2}\bar{u}=-\bar{w}_{xz}-{\it\Omega}_{Iy}^{z}$ and ${\rm\Delta}_{2}\bar{v}=-\bar{w}_{yz}+{\it\Omega}_{Ix}^{z}$ , where ${\rm\Delta}=\partial _{xx}+\partial _{yy}+\partial _{zz}$ and ${\rm\Delta}_{2}=\partial _{xx}+\partial _{yy}$ . Knowing the stationary component one can determine the initial conditions for the wave component:

(3.8) $$\begin{eqnarray}(\tilde{\boldsymbol{u}},{\it\varsigma})_{t=0}=(\boldsymbol{u}_{I}-\bar{\boldsymbol{u}},{\it\varsigma}_{I}).\end{eqnarray}$$

In view of (3.7c ) the system (3.7a,b,d,e ) can be represented in terms of the potential ${\it\Phi}$ :

(3.9a−e ) $$\begin{eqnarray}\tilde{\boldsymbol{u}}=\boldsymbol{{\rm\nabla}}{\it\Phi},\quad p=-{\it\Phi}_{t},\quad {\rm\Delta}{\it\Phi}=0,\quad (g{\it\Phi}_{z}+{\it\Phi}_{tt})_{z=0}=0,\quad {\it\Phi}_{z}|_{z=-H}=0.\end{eqnarray}$$

The solution for the potential is written as a superposition of the irrotational surface gravity waves,

(3.10a−c ) $$\begin{eqnarray}{\it\Phi}=\int A({\bf\kappa})\cosh [{\it\kappa}(z+H)]\text{e}^{\text{i}({\bf\kappa}\boldsymbol{x}_{2}-{\it\sigma}t)}\,\text{d}{\bf\kappa},\quad {\bf\kappa}=(k,l),\quad {\it\sigma}^{2}=g{\it\kappa}\tanh ({\it\kappa}H),\end{eqnarray}$$

where the Fourier amplitude $A({\bf\kappa})$ is determined by the initial conditions (3.8) and $\boldsymbol{x}_{2}=(x,y)$ . If the initial fields are localized in space, then ${\it\Phi}$ decays at a fixed point $x,y$ as $t\rightarrow \infty$ because of geometrical scattering and dispersion of the waves:

(3.11a,b ) $$\begin{eqnarray}{\it\Phi}=O(1/t),\quad \tilde{\boldsymbol{u}}=O(1/t),\quad t\rightarrow \infty .\end{eqnarray}$$

Thus in the process of wave adjustment the system tends with increasing time to the vortical three-dimensional stationary state $(\bar{\boldsymbol{u}},\bar{p},\bar{{\it\varsigma}})$ specified by (3.4) and (3.6). Importantly, the stationary residual flow does not perturb the free surface in view of (3.6a,d ) (see figure 1).

3.3. Weak nonlinearity and balanced equations

We now write the system (3.1) and (3.2) in non-dimensional form using the scales of length $L=H$ , time $T=L/\sqrt{gH}$ , velocity $U$ , pressure $U\sqrt{gH}$ and elevation ${\it\varepsilon}H$ , where ${\it\varepsilon}=U/\sqrt{gH}$ :

(3.12a,b ) $$\begin{eqnarray}\boldsymbol{u}_{t}+{\it\varepsilon}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

(3.13a−c ) $$\begin{eqnarray}w|_{z={\it\varepsilon}{\it\varsigma}}={\it\varsigma}_{t}+{\it\varepsilon}(u{\it\varsigma}_{x}+v{\it\varsigma}_{y})_{z={\it\varepsilon}{\it\varsigma}},\quad p|_{z={\it\varepsilon}{\it\varsigma}}={\it\varsigma},\quad w|_{z=-1}=0.\end{eqnarray}$$

Let the nonlinearity be weak, i.e. the parameter ${\it\varepsilon}\ll 1$ . According to § 2 the solution is sought in the asymptotic form (2.18) and the lowest-order solution $(\boldsymbol{u}^{(0)},p^{(0)},{\it\varsigma}^{(0)})$ obeys equations coinciding with (3.3) in which $g$ and $H$ are replaced by 1. Therefore the formulae (3.4)–(3.11) are valid (mutatis mutandis) for the lowest-order solution if we assume that the invariants ${\it\Omega}^{x}$ , ${\it\Omega}^{y}$ and ${\it\Omega}^{z}$ in (3.6) and the amplitude $A$ in (3.10) depend on the slow time, i.e.

(3.14a,b ) $$\begin{eqnarray}\text{curl }\boldsymbol{u}^{(0)}=({\it\Omega}^{x},{\it\Omega}^{y},{\it\Omega}^{z})(\boldsymbol{x},T_{1},\ldots ),\quad A=A({\bf\kappa},T_{1},\ldots ).\end{eqnarray}$$

As in § 2, this dependence is determined from the boundedness of higher-order corrections.

For the first-order correction $(\boldsymbol{u}^{(1)},p^{(1)},{\it\varsigma}^{(1)})$ we have

(3.15) $$\begin{eqnarray}\boldsymbol{u}_{t}^{(1)}+\boldsymbol{{\rm\nabla}}p^{(1)}=-\boldsymbol{u}_{T_{1}}^{(0)}-(\boldsymbol{u}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}^{(0)}.\end{eqnarray}$$

Applying the curl operator to (3.15) one finds

(3.16a,b ) $$\begin{eqnarray}\partial _{t}\,\text{curl }\boldsymbol{u}^{(1)}=-[D_{T_{1}}^{(0)}\,\text{curl }\boldsymbol{u}^{(0)}-(\text{curl }\boldsymbol{u}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}^{(0)}],\quad D_{T_{1}}^{(0)}=\partial _{T_{1}}+\boldsymbol{u}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}.\end{eqnarray}$$

The fast time averaging (2.12) of (3.16a ) taking into account the asymptotic equation (3.11) gives

(3.17) $$\begin{eqnarray}D_{T_{1}}^{(0)}\,\text{curl }\bar{\boldsymbol{u}}^{(0)}-(\text{curl }\bar{\boldsymbol{u}}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\bar{\boldsymbol{u}}^{(0)}=0.\end{eqnarray}$$

Equation (3.17) describes the evolution of the slow vortex component and should be solved for the initial conditions coinciding with the stationary component $\bar{\boldsymbol{u}}$ of the linear problem (3.6) and (3.4) considered above. In addition, the boundary conditions (3.6d ) should be satisfied:

(3.18) $$\begin{eqnarray}\bar{w}^{(0)}|_{z=0,-1}=0.\end{eqnarray}$$

Equation (3.17) together with the no-flux boundary condition (3.18) and initial conditions constitute a closed balanced model for calculation of the slow vortex component of motion. Analysis of the fast first-order component (very similar to that in Reznik et al. (Reference Reznik, Zeitlin and Ben Jelloul2001)) shows that the amplitude (3.14b ) does not depend on the slow time $T_{1}$ , i.e. it can be considered as constant on times $t\sim 1/{\it\varepsilon}$ .

In view of (3.18) the lowest-order vortex component does not perturb the free surface on times $t\sim 1/{\it\varepsilon}$ (see figure 1). Time averaging of (3.15) gives the equation for $\bar{p}^{(1)}$ :

(3.19) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\bar{p}^{(1)}=-[\bar{\boldsymbol{u}}_{T_{1}}^{(0)}+(\bar{\boldsymbol{u}}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\bar{\boldsymbol{u}}^{(0)}].\end{eqnarray}$$

Knowing $\bar{p}^{(1)}$ one can find the slow part of the first-order elevation correction $\bar{{\it\varsigma}}^{(1)}=\bar{p}^{(1)}(x,y,0,T_{1})$ , which, generally, is non-zero, i.e. the first-order slow component of motion contributes to the elevation. In turn, the first-order fast component ceases to be irrotational, as readily follows from the vortex equation (3.16a ).

Thus the nonlinear wave adjustment in a non-rotating homogeneous fluid with free surface follows the ‘fast radiation’ scenario when the leading-order slow and wave components evolve independently of one another, at least on times ${\sim}1/{\it\varepsilon}$ . This takes place because the assumed smallness of the nonlinearity parameter ${\it\varepsilon}$ means that the typical fluid particle velocity $U$ is small as compared to the typical group speed $\sqrt{gH}$ of surface gravity waves. In this case the waves rapidly radiate away from the initial perturbation and do not interact effectively with the slow component. We have considered the case $L=H$ ; analysis of the long-wave motion $L\gg H$ (not given here) also confirms the validity of this scenario. In the short-wave case $L\ll H$ , in view of (3.10c ) the frequency ${\it\sigma}\sim \sqrt{g/L}$ and the group speed $c_{gr}\sim \sqrt{gL}$ , i.e.  ${\it\sigma}$ increases and $c_{gr}$ decreases with decreasing $L$ . However, the typical phase speed ${\it\sigma}L$ is of the order of $c_{gr}$ and the nonlinearity parameter here is equal to $U/\sqrt{gL}$ . Smallness of the parameter means again that the wave group speed greatly exceeds the typical velocity $U$ , similarly to the case $L\geqslant H$ . Therefore one can think (we did not examine this case in detail) that the ‘fast radiation’ scenario takes place in the short-wave limit too.

Note that, in the absence of rotation, the frequencies of the gravity waves are not bounded from below; the dispersion relation (3.10c ) gives for long waves ${\it\sigma}\cong \sqrt{gH}{\it\kappa}$ for ${\it\kappa}\rightarrow 0$ , i.e. there is no frequency gap between the waves and the stationary flow in this case. Nevertheless, the fast time average works well. The point is that the result of the averaging depends strongly on the structure of the initial perturbation. If the perturbation is localized, then the contribution of very long waves with very small frequencies to its Fourier integral is vanishingly small and the wave part tends to zero with increasing time at a fixed point of space. In moving coordinates the wave perturbation can be non-decaying, for example in the case of one-dimensional non-dispersive waves any localized perturbation propagates at a constant speed, conserving its shape. Obviously, the time average of such a wave perturbation at a fixed space point is zero. A similar situation takes place for the geostrophic adjustment in RSW on a half-plane with rigid lateral boundary (Reznik & Grimshaw Reference Reznik and Grimshaw2002; Reznik & Sutyrin Reference Reznik and Sutyrin2005) where the wave spectrum contains Kelvin waves whose frequency tends to zero in the long-wave limit.

The ‘fast radiation’ scenario also occurs in one- and two-layer RSW models if the dominating horizontal scale $L$ is smaller than or of the order of the corresponding Rossby scale $L_{R}$ (for details, see Reznik et al. Reference Reznik, Zeitlin and Ben Jelloul2001; Zeitlin et al. Reference Zeitlin, Reznik and Ben Jelloul2003). In the next section we will see that a different wave adjustment scenario is realized in this scale range in a three-dimensional homogeneous rotating fluid.

4.1. Statement of the problem

Let the barotropic layer considered above rotate at a constant angular speed ${\it\Omega}=f/2$ as shown in figure 2. The modified equations of motion are written as follows (see also appendix A):

(4.1a,b ) $$\begin{eqnarray}\boldsymbol{u}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}+f\boldsymbol{e}_{\boldsymbol{z}}\times \boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0.\end{eqnarray}$$

The boundary and initial conditions (3.2) remain unchanged. Here $\boldsymbol{e}_{\boldsymbol{z}}$ is the vertical unit vector.

4.2. Linear invariants and waves

Without nonlinear terms equations (4.1) and (3.2a−c ) take the form (cf. (3.3)):

(4.2a−d ) $$\begin{eqnarray}\boldsymbol{u}_{t}+f\boldsymbol{e}_{\boldsymbol{z}}\times \boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad (w,p)_{z=0}=({\it\varsigma}_{t},g{\it\varsigma}),\quad w|_{z=-H}=0.\end{eqnarray}$$

The system (4.2) conserves the non-local invariant ${\it\Pi}$ , which can be naturally treated as a PV for a rotating three-dimensional barotropic layer with free surface:

(4.3) $$\begin{eqnarray}{\it\Pi}=\int _{-H}^{0}(v_{x}-u_{y})\,\text{d}z-f{\it\varsigma}={\it\Pi}_{I}(x,y)=\int _{-H}^{0}(v_{Ix}-u_{Iy})\,\text{d}z-f{\it\varsigma}_{I}.\end{eqnarray}$$

In the representation (3.5) the stationary part now obeys the geostrophic relations

(4.4a−d ) $$\begin{eqnarray}\bar{u}=-\bar{p}_{y}/f,\quad \bar{v}=\bar{p}_{x}/f,\quad \bar{w}=0,\quad \bar{p}=g\bar{{\it\varsigma}}\end{eqnarray}$$

and is determined from the equation following from (4.3) and (4.4),

(4.5) $$\begin{eqnarray}{\rm\Delta}_{2}\bar{{\it\varsigma}}-\frac{f^{2}}{gH}\bar{{\it\varsigma}}=\frac{f}{gH}{\it\Pi}_{I}.\end{eqnarray}$$

The wave component in (3.5) is a solution to the system (4.2) with the zero invariant (4.3) and is a superposition of the wave harmonics

(4.6) $$\begin{eqnarray}(u^{(w)},v^{(w)},w^{(w)},p^{(w)},{\it\varsigma}^{(w)})=(U,V,W,P,Z)\exp [\text{i}(kx+ly-{\it\sigma}t)],\end{eqnarray}$$

where the amplitudes $U$ , $V$ , $W$ , $P$ and $Z$ depend on $z$ and ${\bf\kappa}(k,l)$ . Substitution of (4.6) into (4.2) and exclusion of $U$ , $V$ , $P$ and $Z$ gives the following eigenvalue problem for $W(z)$ and ${\it\sigma}$ :

(4.7a ) $$\begin{eqnarray}W_{zz}-\text{sgn}({\it\sigma}-f)q^{2}W=0,\end{eqnarray}$$

(4.7b,c ) $$\begin{eqnarray}[W_{z}-\text{sgn}({\it\sigma}-f)(gq^{2}/{\it\sigma}^{2})W]_{z=0}=0,\quad W|_{z=-H}=0,\end{eqnarray}$$

where $q={\it\sigma}{\it\kappa}/\sqrt{|{\it\sigma}^{2}-f^{2}|}$ . The wave spectrum consists of super-inertial gravity waves with the eigenfunctions and dispersion relation

(4.8a−c ) $$\begin{eqnarray}W=\sinh [q(z+H)],\quad qH\tanh (qH)={\it\sigma}^{2}H/g,\quad {\it\sigma}={\it\sigma}^{g}\geqslant f\end{eqnarray}$$

and sub-inertial gyroscopic waves for which

(4.9a−c ) $$\begin{eqnarray}W=\sin [q(z+H)],\quad qH\tan (qH)=-{\it\sigma}^{2}H/g,\quad {\it\sigma}={\it\sigma}^{gir}\leqslant f.\end{eqnarray}$$

In the case of low frequencies

(4.10) $$\begin{eqnarray}{\it\sigma}\ll \sqrt{g/H},\end{eqnarray}$$

one obtains the following approximate dispersion relations:

(4.11a,b ) $$\begin{eqnarray}{\it\sigma}^{g}\cong (f^{2}+gH{\it\kappa}^{2})^{1/2},\quad {\it\sigma}^{gir}\cong f/(1+{\it\kappa}^{2}H^{2}/n^{2}{\rm\pi}^{2})^{1/2},\quad n=1,2,\ldots .\end{eqnarray}$$

For typical geophysical parameters $f=10^{-4}~\text{s}^{-1}$ and $H=10^{5}~\text{cm}$ , we have $f^{2}H/g=10^{-6}\ll 1$ , and therefore for gyroscopic waves whose frequencies obey (4.9c ) the dispersion relation (4.11b ) is a good approximation throughout the range of wavenumbers. This is not the case for the gravity waves since (4.11a ) is fulfilled only for the long waves with wavenumbers

(4.12) $$\begin{eqnarray}{\it\kappa}H\ll 1.\end{eqnarray}$$

For the gravity waves with moderate and short wavelengths we have ${\it\sigma}^{g}\gg f$ and the dispersion relation (4.8b ) takes the form (3.10c ):

(4.13) $$\begin{eqnarray}{\it\sigma}^{g}=[g{\it\kappa}\tanh ({\it\kappa}H)]^{1/2},\quad {\it\kappa}H\geqslant 1.\end{eqnarray}$$

In the range of horizontal scales from $H$ to the Rossby scale $L_{R}=\sqrt{gH}/f$ ,

(4.14) $$\begin{eqnarray}H\leqslant L\leqslant L_{R},\end{eqnarray}$$

one obtains from (4.11) and (4.13) the following estimates for the group velocities of the gravity and gyroscopic waves:

(4.15a,b ) $$\begin{eqnarray}c_{gr}^{g}\sim \sqrt{gH},\quad c_{gr}^{gir}\sim {\it\delta}fH\ll c_{gr}^{g},\quad {\it\delta}=H/L.\end{eqnarray}$$

4.3. Weak nonlinearity and balanced equations

Nonlinear effects are considered in the long-wave approximation when

(4.16a,b ) $$\begin{eqnarray}H\ll L\leqslant L_{R},\quad {\it\delta}=H/L\ll 1.\end{eqnarray}$$

There exist three time scales in the problem:

(4.17a−c ) $$\begin{eqnarray}T^{g}=L/\sqrt{gH},\quad T^{r}=f^{-1},\quad T^{ad}=({\it\varepsilon}f)^{-1}.\end{eqnarray}$$

Here $T^{g}$ , $T^{r}$ and $T^{ad}$ are the typical times of gravity waves, rotation and advection, respectively; the Rossby number ${\it\varepsilon}=U/fL$ is assumed to be small,

(4.18) $$\begin{eqnarray}{\it\varepsilon}\ll 1.\end{eqnarray}$$

We now write (4.1) in non-dimensional form. In the range (4.16a ), $T^{g}\leqslant T^{r}$ , and therefore $T^{g}$ is chosen as the time scale. Also we use the scales of horizontal and vertical velocities $U$ and $W={\it\delta}U$ , pressure $U\sqrt{gH}$ and elevation ${\it\varepsilon}BH$ , where $B=L/L_{R}\leqslant 1$ . The resulting equations are written as

(4.19a−c ) $$\begin{eqnarray}\boldsymbol{u}_{t}+{\it\varepsilon}B(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}+B\boldsymbol{e}_{\boldsymbol{z}}\times \boldsymbol{u}=-\hat{\boldsymbol{{\rm\nabla}}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \hat{\boldsymbol{{\rm\nabla}}}=(\partial _{x},\partial _{y},\partial _{z}/{\it\delta}^{2}).\end{eqnarray}$$

The boundary conditions (3.2a−c ) take the form

(4.20a−c ) $$\begin{eqnarray}w|_{z={\it\varepsilon}B{\it\varsigma}}={\it\varsigma}_{t}+{\it\varepsilon}B(u{\it\varsigma}_{x}+v{\it\varsigma}_{y})_{z={\it\varepsilon}B{\it\varsigma}},\quad p|_{z={\it\varepsilon}B{\it\varsigma}}={\it\varsigma},\quad w|_{z=-1}=0.\end{eqnarray}$$

There are two small parameters ${\it\varepsilon}$ and ${\it\delta}$ in (4.19) and (4.20). We put

(4.21) $$\begin{eqnarray}{\it\delta}^{2}\leqslant {\it\varepsilon},\end{eqnarray}$$

which, in view of (4.15), is equivalent to the assumption $U\geqslant {\it\delta}fH\sim c_{gr}^{gir}$ . Since $B\leqslant 1$ and, therefore, ${\it\varepsilon}B=U/\sqrt{gH}\ll 1$ , this choice of scales means that

(4.22) $$\begin{eqnarray}c_{gr}^{g}\gg U\geqslant c_{gr}^{gir},\end{eqnarray}$$

i.e. the typical flow velocity is much smaller than the group speed of gravity waves and exceeds the group speed of gyroscopic waves. Therefore, as we will see (cf. § 3), the balanced slow component affects weakly the gravity waves and strongly the gyroscopic ones.

The solution is again sought in the asymptotic form (2.18). Here only results are given; details can be found in appendix B. The lowest-order equations are a simplified version of the set (4.2), and (4.3) takes the form

(4.23) $$\begin{eqnarray}{\it\Pi}^{(0)}=\int _{-1}^{0}(v_{x}^{(0)}-u_{y}^{(0)})\,\text{d}z-B{\it\varsigma}^{(0)}={\it\Pi}^{(0)}(x,y,T_{1},\ldots ).\end{eqnarray}$$

The slow solution $(\bar{\boldsymbol{u}}^{(0)},\bar{p}^{(0)},\bar{{\it\varsigma}}^{(0)})$ in the representation (2.21) obeys the geostrophic equations analogous to (4.4),

(4.24a−d ) $$\begin{eqnarray}\bar{v}^{(0)}=\bar{{\it\psi}}_{x}^{(0)},\quad \bar{u}^{(0)}=-\bar{{\it\psi}}_{y}^{(0)},\quad \bar{w}^{(0)}=0,\quad \bar{p}^{(0)}=\bar{{\it\varsigma}}^{(0)}=B\bar{{\it\psi}}^{(0)},\end{eqnarray}$$

and the non-dimensional version of (4.5),

(4.25) $$\begin{eqnarray}{\rm\Delta}_{2}\bar{{\it\psi}}^{(0)}-B^{2}\bar{{\it\psi}}^{(0)}={\it\Pi}^{(0)}(x,y,T_{1},\ldots ).\end{eqnarray}$$

The wave solution $(\tilde{\boldsymbol{u}}^{(0)},\tilde{p}^{(0)},\tilde{{\it\varsigma}}^{(0)})$ is characterized by the zero invariant (4.23), i.e.

(4.26) $$\begin{eqnarray}\tilde{{\it\Pi}}^{(0)}=\int _{-1}^{0}(\tilde{v}_{x}^{(0)}-\tilde{u} _{y}^{(0)})\,\text{d}z-B\tilde{{\it\varsigma}}^{(0)}=0,\end{eqnarray}$$

and can be written as the sum

(4.27) $$\begin{eqnarray}(\tilde{\boldsymbol{u}}^{(0)},\tilde{p}^{(0)},\tilde{{\it\varsigma}}^{(0)})=(\tilde{\boldsymbol{u}}^{(g)},\tilde{{\it\varsigma}}^{(0)},\tilde{{\it\varsigma}}^{(0)})+(\hat{\boldsymbol{u}},0,0).\end{eqnarray}$$

The first term on the right-hand side of (4.27) describes the long gravity waves with depth-independent horizontal velocities:

(4.28a,b ) $$\begin{eqnarray}(\tilde{u} ^{(g)},\tilde{v}^{(g)})=\int _{-1}^{0}(\tilde{u} ^{(0)},\tilde{v}^{(0)})\,\text{d}z,\quad \tilde{w}^{(g)}=\tilde{{\it\varsigma}}_{t}^{(0)}(z+1).\end{eqnarray}$$

The elevation $\tilde{{\it\varsigma}}^{(0)}$ obeys the Klein–Gordon equation:

(4.29) $$\begin{eqnarray}\tilde{{\it\varsigma}}_{tt}^{(0)}-{\rm\Delta}_{2}\tilde{{\it\varsigma}}^{(0)}+B^{2}\tilde{{\it\varsigma}}^{(0)}=0.\end{eqnarray}$$

The second term on the right-hand side of (4.27) corresponds to the inertial oscillations (i.e. long gyroscopic waves) in which the pressure, elevation and depth-averaged horizontal velocities are zero:

(4.30a,b ) $$\begin{eqnarray}\hat{u} +\text{i}\hat{v}=A(x,y,z,T_{1},\ldots )\text{e}^{-\text{i}Bt},\quad \int _{-1}^{0}A\,\text{d}z=0,\end{eqnarray}$$

(4.30c ) $$\begin{eqnarray}{\hat{w}}=-\frac{1}{2}\text{e}^{-\text{i}Bt}\int _{-1}^{z}(A_{x}-\text{i}A_{y})\,\text{d}z+\text{c.c.}\end{eqnarray}$$

For localized initial fields the gravity waves decay at a fixed point $x,y$ proportionally to $1/t$ because of the dispersion (see e.g. Reznik et al. Reference Reznik, Zeitlin and Ben Jelloul2001):

(4.31a,b ) $$\begin{eqnarray}\tilde{\boldsymbol{u}}^{(0)}=O(1/t),\quad \tilde{{\it\varsigma}}^{(0)}=O(1/t),\quad t\rightarrow \infty .\end{eqnarray}$$

At the same time the inertial oscillations do not propagate (in this approximation) and remain to be localized in the initial perturbation domain.

Knowing the PV ${\it\Pi}^{(0)}(x,y,T_{1},\ldots )$ and the amplitude of inertial oscillations $A(x,y,z,T_{1},\ldots )$ , one can calculate all other fields. Slow evolution of ${\it\Pi}^{(0)}$ and $A$ is determined from the condition of boundedness of the first-order solution (see appendix B).

The evolution equation for the slow QG component is written for the geostrophic streamfunction $\bar{{\it\psi}}^{(0)}$ :

(4.32a,b ) $$\begin{eqnarray}{\it\Pi}_{T_{1}}^{(0)}+B\text{J}(\bar{{\it\psi}}^{(0)},{\it\Pi}^{(0)})=0,\quad {\it\Pi}^{(0)}={\rm\Delta}_{2}\bar{{\it\psi}}^{(0)}-B^{2}\bar{{\it\psi}}^{(0)},\end{eqnarray}$$

where $\text{J}$ is the Jacobian. Equations (4.32) together with the initial streamfunction $\bar{{\it\psi}}_{I}^{(0)}$ (which is determined as the solution of (4.25) with the right-hand side equal to ${\it\Pi}_{I}$ in (4.3)) constitute the balanced model for the rotating layer. Obviously, (4.32) express conservation of the QG PV ${\it\Pi}^{(0)}$ and coincide exactly with the corresponding equations derived in the model of RSW (see e.g. Reznik et al. Reference Reznik, Zeitlin and Ben Jelloul2001), in which the horizontal velocities and pressure are assumed to be $z$ -independent, i.e. the gyroscopic waves (4.30) are neglected.

The equation describing the slow evolution of the inertial oscillation amplitude $A$ on times $t\sim 1/{\it\varepsilon}$ has the form:

(4.33) $$\begin{eqnarray}A_{T_{1}}+B\text{J}(\bar{{\it\psi}}^{(0)},A)+\frac{\text{i}B}{2}{\rm\Delta}_{2}\bar{{\it\psi}}^{(0)}A+\text{i}\frac{{\it\delta}^{2}B}{2{\it\varepsilon}}{\rm\Delta}_{2}\left(\int _{-1}^{z}A\,\text{d}z-z\int _{-1}^{z}A\,\text{d}z+\frac{1}{2}\int _{-1}^{0}z^{2}A\,\text{d}z\right)=0.\end{eqnarray}$$

Without the last term on the left-hand side of (4.33) we have

(4.34) $$\begin{eqnarray}|A|_{T_{1}}+B\text{J}(\bar{{\it\psi}}^{(0)},|A|)=0,\end{eqnarray}$$

i.e. the inertial oscillations are trapped by the slow QG component. The term proportional to ${\it\delta}^{2}/{\it\varepsilon}$ counteracts the trapping and arises due to slow dispersion of the long gyroscopic waves. The dispersion is effective on times $t\sim 1/{\it\varepsilon}$ if ${\it\delta}^{2}/{\it\varepsilon}\sim 1$ or, in the dimensional form,

(4.35) $$\begin{eqnarray}U\sim {\it\delta}fH\sim c_{gr}^{gir},\end{eqnarray}$$

i.e. the group speed of inertial oscillations and the typical QG velocity are of the same order.

In the barotropic, two-layer and stably continuously stratified RSW models, only gravity waves exist, the gyroscopic ones are absent, but the ‘nonlinear trapping’ scenario occurs if the dominating horizontal scale $L$ greatly exceeds the corresponding Rossby scale $L_{R}$ (Reznik et al. Reference Reznik, Zeitlin and Ben Jelloul2001; Zeitlin et al. Reference Zeitlin, Reznik and Ben Jelloul2003). Physically, this is due to the fact that the surface and internal gravity waves with $L\gg L_{R}$ are close to inertial oscillations with small group speeds (see e.g. (4.11a )).

Surprisingly, the fast inertial oscillations, being trapped by the slow component, do not contribute to the latter on times ${\sim}1/{\it\varepsilon}$ as readily seen from (4.32). Self-interaction of the oscillations gives a contribution to the first-order PV equation (B 7) only on a double inertial frequency; the corresponding terms with zero frequency cancel each other. The same is valid in the RSW models (Reznik et al. Reference Reznik, Zeitlin and Ben Jelloul2001; Zeitlin et al. Reference Zeitlin, Reznik and Ben Jelloul2003), in a three-dimensional rotating barotropic fluid confined between two rigid lids (Reznik Reference Reznik2014a ) and in a rotating stably neutrally stratified (SNS) fluid (Reznik Reference Reznik2014b ). In the RSW models, the trapped oscillations do not affect the slow motion even on longer times ${\sim}1/{\it\varepsilon}^{2}$ . At the same time, in a three-dimensional rotating fluid, the long gyroscopic waves affect the slow component on times ${\sim}1/{\it\varepsilon}^{2}$ (Reznik Reference Reznik2014a ).

5.1. Statement of the problem

We now consider a layer of continuously stratified fluid of constant depth $H$ confined between two rigid lids $z=0,-H$ (figure 3). In the Boussinesq approximation, the system is described by the equations (see also appendix A)

(5.1a−c ) $$\begin{eqnarray}\boldsymbol{u}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}+\boldsymbol{e}_{z}b=-\boldsymbol{{\rm\nabla}}p,\quad b_{t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}b-N^{2}w=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

with the no-flux boundary conditions and initial conditions

(5.2a,b ) $$\begin{eqnarray}w|_{z=0,-H}=0,\quad (u,v,b)_{t=0}=(u_{I},v_{I},b_{I}).\end{eqnarray}$$

Here $b=g{\it\rho}/{\it\rho}_{0}$ is the buoyancy, ${\it\rho}/{\it\rho}_{0}$ and $p$ are the variations of density and pressure from their hydrostatic profiles divided by the reference density ${\it\rho}_{0}$ , and $N$ is the buoyancy frequency.

5.2. Linear invariants and linear waves

The linearized system (5.1),

(5.3a−c ) $$\begin{eqnarray}\boldsymbol{u}_{t}+\boldsymbol{e}_{z}b=-\boldsymbol{{\rm\nabla}}p,\quad b_{t}-N^{2}w=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

conserves the vertical component of vorticity,

(5.4) $$\begin{eqnarray}{\it\Omega}^{z}=v_{x}-u_{y}={\it\Omega}_{I}^{z}=v_{Ix}-u_{Iy},\end{eqnarray}$$

and, in view of (5.3b ) and (5.2), the buoyancy at the boundaries $z=0,-H$ ,

(5.5) $$\begin{eqnarray}b|_{z=0,-H}=b_{I}|_{z=0,-H}.\end{eqnarray}$$

In the corresponding representation (2.9), the stationary solution obeys the equations

(5.6a−c ) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\bar{p}+\boldsymbol{e}_{z}\bar{b}=0,\quad \bar{w}=0,\quad \bar{u}_{x}+\bar{v}_{y}=0.\end{eqnarray}$$

Assuming the fields to decay at infinity in any horizontal plane, one finds that

(5.7a−c ) $$\begin{eqnarray}\bar{p}=\bar{b}=0,\quad \bar{u}=-\bar{{\it\psi}}_{y},\quad \bar{v}=\bar{{\it\psi}}_{x},\end{eqnarray}$$

where the streamfunction $\bar{{\it\psi}}$ is found from the equation

(5.8) $$\begin{eqnarray}{\rm\Delta}_{2}\bar{{\it\psi}}={\it\Omega}_{I}^{z},\end{eqnarray}$$

following from (5.6c ) and (5.4).

Lighthill (Reference Lighthill1996) examined the linear wave adjustment of localized perturbation in an unbounded (vertically and horizontally) stratified ocean and showed that the motion is split into a steady vortical flow obeying (5.6)–(5.8) and radiated internal waves propagating away from the initial perturbation. In our case the fluid is bounded by two rigid lids, which results in the existence of two additional ‘boundary’ invariants (5.5) absent in the vertically unbounded fluid. These invariants are not used when determining the stationary solution (5.7), which does not satisfy (5.5) if the initial buoyancy $b_{I}$ is non-zero at the boundaries. In this case the time-independent non-zero boundary conditions (5.5) should be satisfied by the wave part of the solution $\tilde{\boldsymbol{u}}=(\tilde{u} ,\tilde{v},w),p,b$ , which obeys the equations

(5.9a−d ) $$\begin{eqnarray}\tilde{\boldsymbol{u}}_{t}+\boldsymbol{e}_{z}b=-\boldsymbol{{\rm\nabla}}p,\quad b_{t}-N^{2}w=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=0,\quad \tilde{{\it\Omega}}^{z}=\tilde{v}_{x}-\tilde{u} _{y}=0,\end{eqnarray}$$

(5.10a,b ) $$\begin{eqnarray}w|_{z=0,-H}=0,\quad (\tilde{u} ,\tilde{v},b)_{t=0}=(u_{I}-\bar{u},v_{I}-\bar{v},b_{I}).\end{eqnarray}$$

As a rule, the equations of internal waves (5.9) and (5.10) are reduced to a problem for the vertical velocity $w$ (see e.g. Miropol’sky Reference Miropol’sky2001). From (5.9a–c ) one derives the equation

(5.11) $$\begin{eqnarray}{\rm\Delta}w_{tt}+N^{2}{\rm\Delta}_{2}w=0,\end{eqnarray}$$

which should be solved under the no-flux conditions (5.10a ) and initial conditions

(5.12a,b ) $$\begin{eqnarray}(w,w_{t})_{t=0}=(w_{I},{\dot{w}}_{I}),\quad w_{I}=-\int _{-H}^{z}(\partial _{x}u_{I}+\partial _{y}v_{I})\,\text{d}z.\end{eqnarray}$$

The initial field ${\dot{w}}_{I}$ in (5.12a ) can be expressed in terms of $u_{I}$ , $v_{I}$ and $b_{I}$ (see e.g. Reznik Reference Reznik2013).

The solution to the problem (5.11), (5.10a ) and (5.12a ) is well known (see e.g. Miropol’sky Reference Miropol’sky2001). The velocity $w$ is written as the Fourier integral in the horizontal coordinates,

(5.13) $$\begin{eqnarray}w=\int {\hat{w}}({\bf\kappa},z,t)\text{e}^{\text{i}(kx+ly)}\,\text{d}{\bf\kappa},\quad {\bf\kappa}=(k,l),\end{eqnarray}$$

and for the Fourier amplitude ${\hat{w}}$ we have

(5.14a−c ) $$\begin{eqnarray}({\hat{w}}_{zz}-{\it\kappa}^{2}{\hat{w}})_{tt}-{\it\kappa}^{2}N^{2}{\hat{w}}=0,\quad {\hat{w}}|_{z=0,-H}=0,\quad ({\hat{w}},{\hat{w}}_{t})_{t=0}=({\hat{w}}_{I},\hat{{\dot{w}}}_{I}).\end{eqnarray}$$

Here and below the hat denotes the corresponding Fourier amplitude. The solution for ${\hat{w}}$ can be written as the series

(5.15a,b ) $$\begin{eqnarray}{\hat{w}}=\mathop{\sum }_{n=1}^{\infty }[{\hat{w}}_{In}\cos ({\it\sigma}_{n}t)+(\hat{{\dot{w}}}_{In}/{\it\sigma}_{n})\sin ({\it\sigma}_{n}t)]G_{n}(z,{\it\kappa}),\quad {\it\sigma}_{n}=1/\sqrt{{\it\lambda}_{n}},\end{eqnarray}$$

where $G_{n}$ and ${\it\lambda}_{n}$ are the eigenfunctions and eigenvalues of the Sturm–Liouville problem:

(5.16a,b ) $$\begin{eqnarray}G_{zz}-{\it\kappa}^{2}G=-{\it\lambda}{\it\kappa}^{2}N^{2}G,\quad G|_{z=0,-H}=0.\end{eqnarray}$$

In the case $N=\text{const.}$ we have

(5.17a,b ) $$\begin{eqnarray}{\it\lambda}_{n}=(1+n^{2}{\rm\pi}^{2}/{\it\kappa}^{2}H^{2})/N^{2},\quad G_{n}=\sqrt{2/H}\sin (n{\rm\pi}z/H).\end{eqnarray}$$

In view of the conditions

(5.18) $$\begin{eqnarray}(w_{I},{\dot{w}}_{I})_{z=0,-H}=0,\end{eqnarray}$$

the coefficients ${\hat{w}}_{In}$ and $\hat{{\dot{w}}}_{In}$ decay rapidly with increasing $n$ :

(5.19a,b ) $$\begin{eqnarray}{\hat{w}}_{In}=O(n^{-3}),\quad \hat{{\dot{w}}}_{In}=O(n^{-3}),\quad n\rightarrow \infty .\end{eqnarray}$$

Since ${\it\sigma}_{n}=O(n^{-1}),~n\rightarrow \infty$ , the series (5.15a ) converges absolutely and uniformly in the domain $[-H\leqslant z\leqslant 0]\times [0\leqslant t\leqslant \infty ]$ . Solutions for the horizontal velocities $\tilde{u} ,\tilde{v}$ and pressure $p$ are given in appendix C.

The problem for the buoyancy is conveniently reduced to the form

(5.20a−c ) $$\begin{eqnarray}{\rm\Delta}c_{tt}+N^{2}{\rm\Delta}_{2}c=0,\quad c|_{z=0,-H}=(b_{I}/N^{2})_{z=0,-H},\quad (c,c_{t})_{t=0}=(c_{I},w_{I}),\end{eqnarray}$$

where $c=b/N^{2}$ is the normalized buoyancy and $c_{I}=b_{I}/N^{2}$ . If the initial buoyancy is zero at the boundaries, i.e. $b_{I}|_{z=0,-H}=0$ , then (5.20) is analogous to the problem (5.11) and (5.12) and the buoyancy Fourier amplitude can be represented as an absolutely and uniformly converging series of the form (5.15). However, for

(5.21) $$\begin{eqnarray}b_{I}|_{z=0,-H}\neq 0,\end{eqnarray}$$

the situation becomes somewhat more complicated. First, in view of (5.20a ), the time average $\langle c\rangle =0$ inside the interval $[-H,0]$ , but $\langle c\rangle \neq 0$ at the boundaries by virtue of (5.20b ) and (5.21). Second, the vertical gradients of horizontal velocities grow linearly in time near the boundaries since it readily follows from (5.9a,b ) and (5.10) that

(5.22) $$\begin{eqnarray}(\tilde{u} _{z},\tilde{v}_{z})_{z=0,-H}=t{\rm\nabla}_{2}(b_{I})_{z=0,-H}+(\tilde{u} _{Iz},\tilde{v}_{Iz})_{z=0,-H}.\end{eqnarray}$$

These peculiarities indicate that in close vicinities of the boundaries the solutions for $\tilde{u} ,\tilde{v},c$ have a boundary-layer structure when, for example, near $z=0$ the normalized buoyancy $c$ at large times is represented as

(5.23) $$\begin{eqnarray}c=c(x,y,zt),\end{eqnarray}$$

where the right-hand side tends to zero at any $z<0$ , $t\rightarrow \infty$ but $c(x,y,0)\neq 0$ .

To understand better the physics of the boundary layer, we consider the solution to (5.20) in terms of the Fourier amplitude ${\hat{c}}$ , which can be written in the forms (see appendix C)

(5.24) $$\begin{eqnarray}{\hat{c}}={\hat{c}}_{I}+\mathop{\sum }_{n=1}^{\infty }\{{\hat{c}}_{In}[\cos ({\it\sigma}_{n}t)-1]+({\hat{w}}_{In}/{\it\sigma}_{n})\sin ({\it\sigma}_{n}t)\}G_{n}(z,{\it\kappa})\end{eqnarray}$$

or

(5.25) $$\begin{eqnarray}{\hat{c}}=\mathop{\sum }_{n=1}^{\infty }[{\hat{c}}_{In}\cos ({\it\sigma}_{n}t)+({\hat{w}}_{In}/{\it\sigma}_{n})\sin ({\it\sigma}_{n}t)]G_{n}(z,{\it\kappa}).\end{eqnarray}$$

In the case (5.21) we have (appendix C)

(5.26) $$\begin{eqnarray}{\hat{c}}_{In}=O(n^{-1}),\quad n\rightarrow \infty ,\end{eqnarray}$$

and therefore the terms of the series in (5.24) behave as $O(n^{-3})$ , $n\rightarrow \infty$ , at a fixed $t$ , i.e. the series converges absolutely and uniformly in the domain $[-H,0]\times [0,T_{0}]$ , where $T_{0}$ is an arbitrary finite time. However, the convergence gets worse with increasing time $t$ since in (5.24) the terms ${\hat{c}}_{In}[\cos ({\it\sigma}_{n}t)-1]\sim t^{2}/n^{3}$ and $({\hat{w}}_{In}/{\it\sigma}_{n})\sin ({\it\sigma}_{n}t)\sim t/n^{3}$ as $n\rightarrow \infty$ . The series (5.25) has the same form as (5.15) but it converges much slower because of (5.26). Any partial sum in the series has zero time average (2.12); the sum does not have ‘boundary layer’ form (5.23) and its contribution to the corresponding Fourier integral of the form (5.13) for $c$ tends to zero with increasing time because of dispersion. This means that the non-stationary boundary layers developing at large times near the surface and bottom arise as a result of joint impact of the low-frequency wave harmonics with numbers $1\ll n\leqslant \infty$ in the series (5.24) and (5.25), characterized by a very weak dispersion.

5.3. Non-stationary boundary layer

To describe the buoyancy at large times, one uses the new variable (cf. Il’in Reference Il’in1970, Reference Il’in1972; Kamenkovich & Kamenkovich Reference Kamenkovich and Kamenkovich1993)

(5.27a,b ) $$\begin{eqnarray}C=\frac{1}{t}\int _{0}^{t}c\,\text{d}t,\quad c=(tC)_{t}.\end{eqnarray}$$

The meaning of the variable $C$ is that the contribution of the rapidly oscillating part of the field $c$ to $C$ becomes negligible at large times $t\gg 1$ , as seen from (5.25). Equations (5.20) are rewritten in non-dimensional form using the length scale $H$ and the time scale $1/N$ :

(5.28a−c ) $$\begin{eqnarray}{\rm\Delta}C_{tt}+\frac{2}{t}{\rm\Delta}C_{t}+{\rm\Delta}_{2}C=\frac{1}{t}{\rm\Delta}w_{I},\quad C|_{z=0,-1}=c_{I}|_{z=0,-1},\quad (C,C_{t})_{t=0}=(c_{I},w_{I}/2).\end{eqnarray}$$

For simplicity, $N$ is assumed to be constant.

Let us consider the problem (5.28) at large times. Outside the boundary layers the vertical scale is assumed to be of the order of the layer depth and the solution here is sought in the following asymptotic form:

(5.29) $$\begin{eqnarray}C=\frac{1}{t}C_{1}(x,y,z)+\frac{1}{t^{2}}C_{2}(x,y,z)+\cdots \,.\end{eqnarray}$$

Substitution of (5.29) into (5.28a ) gives

(5.30) $$\begin{eqnarray}{\rm\Delta}_{2}C_{1}={\rm\Delta}w_{I}.\end{eqnarray}$$

Generally, $C_{1}/t$ does not satisfy the boundary conditions (5.28b ); this discrepancy is corrected by boundary layers confined to narrow vicinities near the boundaries. The solution near $z=0$ is sought in the form of the expansion

(5.31) $$\begin{eqnarray}C=D_{0}(x,y,{\it\xi})+\frac{1}{t}D_{1}(x,y,{\it\xi})+\cdots \,,\quad {\it\xi}=-zt,\end{eqnarray}$$

where ${\it\xi}$ is the boundary layer stretched variable.

In the leading order we have the equation

(5.32) $$\begin{eqnarray}{\it\xi}^{2}D_{0}^{(4)}+6{\it\xi}D_{0}^{(3)}+6D_{0}^{\prime \prime }+{\rm\Delta}_{2}D_{0}=0,\end{eqnarray}$$

where the superscripts and primes mean differentiation with respect to ${\it\xi}$ . The function $D_{0}$ should satisfy the boundary conditions

(5.33a,b ) $$\begin{eqnarray}D_{0}|_{{\it\xi}=0}=c_{I}|_{z=0},\quad D_{0}\rightarrow 0,\quad {\it\xi}\rightarrow \infty .\end{eqnarray}$$

The corresponding equations for the Fourier amplitude $\hat{D}_{0}$ are written as

(5.34a−c ) $$\begin{eqnarray}{\it\xi}^{2}\hat{D}_{0}^{(4)}+6{\it\xi}\hat{D}_{0}^{(3)}+6\hat{D^{\prime \prime }}_{0}-{\it\kappa}^{2}\hat{D}_{0}=0,\quad \hat{D}_{0}|_{{\it\xi}=0}={\hat{c}}_{I}|_{z=0},\quad \hat{D}_{0}\rightarrow 0,\quad {\it\xi}\rightarrow \infty .\end{eqnarray}$$

The solution to (5.34) has the form (Kamke Reference Kamke1976)

(5.35) $$\begin{eqnarray}\hat{D}_{0}=\frac{{\hat{c}}_{I}|_{z=0}}{\sqrt{{\it\kappa}{\it\xi}}}\text{J}_{1}(2\sqrt{{\it\kappa}{\it\xi}}),\end{eqnarray}$$

where here and below $\text{J}_{n}$ is the Bessel function of $n$ th order. The lowest-order normalized buoyancy and vertical velocity are found from (5.27b ), (5.35) and (5.9b ) as

(5.36a,b ) $$\begin{eqnarray}{\hat{c}}_{0}=(t\hat{D}_{0})_{t}={\hat{c}}_{I}|_{z=0}\text{J}_{0}(2\sqrt{{\it\kappa}{\it\xi}}),\quad {\hat{w}}_{0}={\hat{c}}_{0t}=-\hat{b}_{I}|_{z=0}\frac{\sqrt{{\it\kappa}{\it\xi}}}{t}\text{J}_{1}(2\sqrt{{\it\kappa}{\it\xi}}).\end{eqnarray}$$

Using (5.36b ) and (5.9c,d ) we find the Fourier amplitudes of the lowest-order boundary layer horizontal velocities:

(5.37) $$\begin{eqnarray}(\hat{\tilde{u} }_{0},\hat{\tilde{v}}_{0})=\frac{\text{i}}{{\it\kappa}}(k,l)\hat{b}_{I}|_{z=0}\text{J}_{0}(2\sqrt{{\it\kappa}{\it\xi}}).\end{eqnarray}$$

One can readily check that (5.37) agrees with (5.22) and the horizontal boundary layer velocity is time-independent at $z=0$ . The boundary layer near the bottom $z=-H$ is analysed in the same way.

5.4. Weak nonlinearity and balanced equations

We now write (5.1) in non-dimensional form using the scales of length $L=H$ , time $T=1/N_{0}$ , velocity $U$ , pressure $UHN_{0}$ and buoyancy $UN_{0}$ :

(5.38a−c ) $$\begin{eqnarray}\boldsymbol{u}_{t}+{\it\varepsilon}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}+\boldsymbol{e}_{z}b=-\boldsymbol{{\rm\nabla}}p,\quad b_{t}+{\it\varepsilon}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}b-N^{2}w=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0.\end{eqnarray}$$

Here $N_{0}$ is the typical buoyancy frequency and the nonlinearity parameter ${\it\varepsilon}$ is assumed to be small:

(5.39) $$\begin{eqnarray}{\it\varepsilon}=U/HN_{0}\ll 1.\end{eqnarray}$$

The boundary and initial conditions (5.2) keep their form with $H$ replaced by 1.

Again, the solution is sought in the asymptotic form (2.18); details are given in appendix C. The lowest-order solution $(u^{(0)},v^{(0)},w^{(0)},b^{(0)},p^{(0)})$ obeys the linear equations (5.3) and conditions (5.2), and therefore the slow component of the solution $\bar{u}^{(0)}(\boldsymbol{x},T_{1},\ldots ),\ldots$ satisfies (5.6) and (5.7),

(5.40a−c ) $$\begin{eqnarray}\bar{p}^{(0)}=\bar{b}^{(0)}=\bar{w}^{(0)}=0,\quad \bar{u}^{(0)}=-\bar{{\it\psi}}_{y},\quad \bar{v}^{(0)}=\bar{{\it\psi}}_{x},\end{eqnarray}$$

and (5.8) takes the form

(5.41) $$\begin{eqnarray}{\rm\Delta}_{2}\bar{{\it\psi}}={\it\Omega}^{z}(\boldsymbol{x},T_{1}).\end{eqnarray}$$

The initial conditions for the slow component follow from (5.41):

(5.42a−c ) $$\begin{eqnarray}{\rm\Delta}_{2}\bar{{\it\psi}}_{I}={\it\Omega}_{I}^{z},\quad \bar{u}_{I}^{(0)}=-\bar{{\it\psi}}_{Iy},\quad \bar{v}_{I}^{(0)}=\bar{{\it\psi}}_{Ix}.\end{eqnarray}$$

The buoyancy at $z=0,-1$ also depends on the slow time $T_{1}$ , i.e. instead of (5.5) we have

(5.43a,b ) $$\begin{eqnarray}b^{(0)}|_{z=0}=b_{0}(x,y,T_{1},\ldots ),\quad b^{(0)}|_{z=-1}=b_{-1}(x,y,T_{1},\ldots ).\end{eqnarray}$$

The functions $b_{0,-1}$ are known only at the initial moment:

(5.44) $$\begin{eqnarray}b_{0,-1}|_{t=0}=b_{I}|_{z=0,-1}.\end{eqnarray}$$

The dependence of ${\it\Omega}^{z}$ and $b_{0,-1}$ on $T_{1}$ is determined from the condition of boundedness of first-order fields (see appendix C).

The lowest-order wave component is described by (5.9) and (5.10); in this case the solution (5.15) for the vertical velocity remains unchanged, at least up to times $t=O(1/{\it\varepsilon})$ (see Reznik et al. Reference Reznik, Zeitlin and Ben Jelloul2001). As for the buoyancy and horizontal velocity, the representations (5.24), (5.25) and (C 3b ) remain relevant for $t\ll 1/{\it\varepsilon}$ ; for $t=O(1/{\it\varepsilon})$ the solution for the boundary layer becomes relevant, with the quantities $b_{I}|_{z=0,-1}$ replaced by $b_{0,-1}(x,y,T_{1},\ldots )$ in (5.22), (5.36b ) and (5.37), and $c_{I}|_{z=0,-1}$ replaced by $b_{0,-1}/N^{2}|_{z=0,-1}$ in (5.28b ), (5.33a ), (5.34b ), (5.35) and (5.36a ).

The streamfunction $\bar{{\it\psi}}$ obeys the equation for two-dimensional hydrodynamics,

(5.45) $$\begin{eqnarray}{\rm\Delta}_{2}\bar{{\it\psi}}_{T_{1}}+\text{J}(\bar{{\it\psi}},{\rm\Delta}_{2}\bar{{\it\psi}})=0,\end{eqnarray}$$

while the surface buoyancy is described by the equations

(5.46a ) $$\begin{eqnarray}z=0,-1:\quad b_{T_{1}}^{(0)}+\text{J}(\bar{{\it\psi}},b^{(0)})+\tilde{u} _{0}b_{x}^{(0)}+\tilde{v}_{0}b_{y}^{(0)}=0.\end{eqnarray}$$

The ‘wave’ velocities $\tilde{u} _{0}$ and $\tilde{v}_{0}$ do not depend on the fast time at the boundaries in view of (5.37) and in terms of the Fourier amplitudes are given by

(5.46b ) $$\begin{eqnarray}(\hat{\tilde{u} }_{0},\hat{\tilde{v}}_{0})_{z=0,-1}=\left.\frac{\text{i}}{{\it\kappa}}(k,l)\hat{b}^{(0)}\right|_{z=0,-1}.\end{eqnarray}$$

Equations (5.45) and (5.46) determine the evolution of the slow component and should be solved under the initial conditions following from (5.42a ) for the streamfunction and from (5.43) and (5.44) for the boundary buoyancy. The slow vertical velocity $\bar{w}^{(0)}$ is zero in view of (5.40a ), and therefore the slow motion is quasi-two-dimensional, i.e. the motion at any plane $z=z_{0}=\text{const.}$ does not depend on the motion at other planes, whereas the motions at different horizons differ from one another if the initial vorticity ${\it\Omega}_{I}^{z}$ varies with  $z$ .

During the wave adjustment (see figure 3) the wave parts of the vertical velocity and pressure decay with increasing time at any fixed point $x$ , $y$ , $z$ . As for the buoyancy and wave part of the horizontal velocity, they turn out to be confined to the non-stationary boundary layers near the surface and bottom at large times. Thus any localized initial perturbation tends with increasing time to the sum of a quasi-two-dimensional state determined by the vortex equation (5.45) and the boundary layer fields. The boundary layers are characterized by sharp vertical gradients of the buoyancy and horizontal velocity growing proportionally to the time $t$ ; the layer thicknesses decrease proportionally to $t$ . The time dependence of the buoyancy and wave horizontal velocities slows down when approaching the boundaries; at the boundaries the quantities do not depend on the fast time $t$ . It is impossible to split the slow and wave components even in the linear approximation; therefore we call this the ‘imperfect splitting’ scenario.

The same scenario takes place for the wave adjustment in an SNS fluid consisting of a stratified upper layer and a homogeneous lower layer, the density and other fields being continuous at the interface between the layers (Reznik Reference Reznik2013). In this model the vertical component of vorticity is invariant in the upper layer and the full three-dimensional vorticity components are invariants in the lower one. Generally, the stationary solution related to these invariants is discontinuous at the interface between the layers. The result of this is that, in parallel with the boundary layers near rigid surface, an internal boundary layer arises at large times near the interface between the stratified and homogeneous fluids. The layer is characterized by sharp vertical gradients of the horizontal velocity.

Dissipation can be important in these domains with sharp gradients. As noted by an anonymous reviewer, ‘In the presence of dissipation, mean flow can be generated due to the dissipative operator’. Probably, in our case the dissipation will transform the non-stationary boundary layers into stationary ones, and, in doing so, contributes to the slow flow. This interesting and important effect will be considered elsewhere.

5.5. Effects of rotation

Linearized equations of a rotating stratified fluid have the form (cf. (4.2) and (5.3)):

(5.47a−c ) $$\begin{eqnarray}\boldsymbol{u}_{t}+f\boldsymbol{e}_{\boldsymbol{z}}\times \boldsymbol{u}+\boldsymbol{e}_{\boldsymbol{z}}b=-\boldsymbol{{\rm\nabla}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad b_{t}-N^{2}w=0.\end{eqnarray}$$

The boundary and initial conditions (5.2) remain unchanged. The conserved quantity here is the PV,

(5.48) $$\begin{eqnarray}{\it\Pi}=v_{x}-u_{y}-f(b/N^{2})_{z}={\it\Pi}_{I}.\end{eqnarray}$$

The stationary geostrophic solution is expressed in terms of the geostrophic streamfunction $\bar{{\it\psi}}$ as

(5.49a−d ) $$\begin{eqnarray}\bar{u}=-\bar{{\it\psi}}_{y},\quad \bar{v}=-\bar{{\it\psi}}_{y},\quad \bar{b}=-f\bar{{\it\psi}}_{z},\quad \bar{{\it\psi}}=\bar{p}/f,\end{eqnarray}$$

determined from the equation

(5.50) $$\begin{eqnarray}{\rm\Delta}_{2}\bar{{\it\psi}}+f^{2}(\bar{{\it\psi}}_{z}/N^{2})_{z}={\it\Pi}_{I},\end{eqnarray}$$

following from (5.48) and (5.49). The conservation of boundary buoyancy (5.5) gives here the boundary conditions for (5.50):

(5.51) $$\begin{eqnarray}\bar{{\it\psi}}_{z}|_{z=0,-H}=-(1/f)b_{I}|_{z=0,-H}.\end{eqnarray}$$

That is, in the rotating case, the stationary solution and conditions (5.5) do not ‘contradict’ each other unlike in the non-rotating fluid.

Nevertheless, non-stationary boundary layers arise too in the rotating fluid. It follows from (5.47) and (5.5) that, at the boundaries $z=0,-H$ (M. V. Kalashnik, private communication),

(5.52a,b ) $$\begin{eqnarray}U_{zt}+\text{i}fU_{z}=b_{Ix}+\text{i}b_{Iy},\quad U=u+\text{i}v,\end{eqnarray}$$

whence one obtains that

(5.53) $$\begin{eqnarray}U_{z}=[U_{Iz}+(\text{i}/f)(b_{Ix}+\text{i}b_{Iy})]\text{e}^{-\text{i}ft}-(\text{i}/f)(b_{Ix}+\text{i}b_{Iy}),\quad z=0,-H.\end{eqnarray}$$

Thus for $U_{Iz}|_{z=0,-H}\neq 0$ and/or ${\rm\nabla}_{2}b_{I}|_{z=0,-H}\neq 0$ inertial oscillations proportional to $\text{e}^{-\text{i}ft}$ are generated in the vertical gradient of the horizontal velocity at the boundaries. Let us multiply (5.9) and (5.10a ) by the harmonic $\text{e}^{\pm \text{i}ft}$ and apply the time averaging (2.12) to the resulting equations. A simple investigation of the obtained system shows (see e.g. Reznik Reference Reznik2014b ) that inertial oscillations localized in the horizontal plane are prohibited inside a stably stratified fluid, i.e.

(5.54) $$\begin{eqnarray}\lim _{T_{0}\rightarrow \infty }\frac{1}{T_{0}}\int _{0}^{T_{0}}(\boldsymbol{u},b,p)\text{e}^{\pm \text{i}ft}\,\text{d}t=0.\end{eqnarray}$$

Therefore near the boundaries $z=0,-H$ at large times boundary layers should arise that provide a transition from the ‘boundary’ motion (5.53) to the motion in the inner domain. A description of the boundary layers is given in Reznik (Reference Reznik2014b ).

In a rotating SNS fluid the boundary layer near the interface is more intense than that near the rigid boundaries; the layer prevents penetration into the stably stratified upper fluid of ‘prohibited’ inertial oscillations (long sub-inertial gyroscopic waves) which occur in the homogeneous fluid and cannot exist in the stratified one (Reznik Reference Reznik2014b ).