3.1.1 Vertical structure functions

Evidently, the vertical structure functions $\hat{\unicode[STIX]{x1D702}}_{0}$ , $\hat{\unicode[STIX]{x1D702}}_{1}$ , $\hat{\unicode[STIX]{x1D713}}_{0}$ and $\hat{\unicode[STIX]{x1D713}}_{1}$ cannot be chosen independently. In fact, by considering the linearized (in $\unicode[STIX]{x1D6FC}$ ) governing equations (2.1) at zeroth order in $\unicode[STIX]{x1D700}$ , it can be shown that $\hat{\unicode[STIX]{x1D702}}_{0}=\hat{\unicode[STIX]{x1D713}}_{0}$ , and both vertical structure functions satisfy the same modal equation:

(3.4) $$\begin{eqnarray}{\displaystyle \frac{\text{d}^{2}\,\hat{\unicode[STIX]{x1D702}}_{0}}{\text{d}z^{2}}}+\left({\displaystyle \frac{N^{2}}{\unicode[STIX]{x1D714}^{2}}}-1\right)k^{2}\hat{\unicode[STIX]{x1D702}}_{0}=0.\end{eqnarray}$$

For stratification prescribed by $N^{2}(z)$ , this differential eigenvalue problem can be solved for $\hat{\unicode[STIX]{x1D702}}_{0}$ and the corresponding dispersion relation $\unicode[STIX]{x1D714}(k)$ , both of which depend upon the vertical mode number. Crucially, equation (3.4) must be solved subject to no-flow boundary conditions at the two confining horizontal walls to ensure the solution propagates only horizontally. At the next order in $\unicode[STIX]{x1D700}$ , we obtain the differential equation describing the vertical structure function $\hat{\unicode[STIX]{x1D702}}_{1}$ , as well as $\hat{\unicode[STIX]{x1D713}}_{1}$

(3.5) $$\begin{eqnarray}\left({\displaystyle \frac{1}{k^{2}}}{\displaystyle \frac{\text{d}^{2}}{\text{d}z^{2}}}+{\displaystyle \frac{N^{2}}{\unicode[STIX]{x1D714}^{2}}}-1\right)\hat{\unicode[STIX]{x1D702}}_{1}=2\left((1-\unicode[STIX]{x1D712}){\displaystyle \frac{N^{2}}{\unicode[STIX]{x1D714}^{2}}}-1\right)\hat{\unicode[STIX]{x1D702}}_{0},\end{eqnarray}$$

where we have defined $\unicode[STIX]{x1D712}\equiv c_{g}/c_{p}$ . We further find that $\hat{\unicode[STIX]{x1D702}}_{1}=\hat{\unicode[STIX]{x1D713}}_{1}-(1-\unicode[STIX]{x1D712})\hat{\unicode[STIX]{x1D702}}_{0}$ . Alternatively, the compatibility condition (3.5) can be regarded as a definition of the group velocity (see appendix B for details). Below, we will show that only the vertical structure function $\hat{\unicode[STIX]{x1D702}}_{0}$ is required to compute the Stokes drift and the induced Eulerian flow.

3.1.2 Stokes drift

One can follow the methodology first performed by Stokes (Reference Stokes1847) to derive the order amplitude-squared displacement of fluid parcels by a horizontally periodic internal mode. Generally, the Stokes drift is given by

(3.6) $$\begin{eqnarray}\boldsymbol{u}_{S}\equiv \overline{\boldsymbol{u}_{x}^{(1)}\unicode[STIX]{x1D709}^{(1)}+\boldsymbol{u}_{z}^{(1)}\unicode[STIX]{x1D702}^{(1)}},\end{eqnarray}$$

in which the $x$ and $z$ subscripts denote partial derivatives, $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D702}$ are the horizontal and vertical displacements, respectively, and the overline denotes averaging over the fast time scales of the waves. From the polarization relations for internal modes in table 1 and using (3.4), we find the following leading order

(3.7a,b ) $$\begin{eqnarray}u_{S}={\displaystyle \frac{1}{4}}c_{p}|A_{0}|^{2}{\displaystyle \frac{\text{d}^{2}\hat{\unicode[STIX]{x1D702}}_{0}^{2}}{\text{d}z^{2}}},\quad w_{S}=-{\displaystyle \frac{1}{4}}c_{p}(1+\unicode[STIX]{x1D712})\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{X}|A_{0}|^{2}{\displaystyle \frac{\text{d}\hat{\unicode[STIX]{x1D702}}_{0}^{2}}{\text{d}z}}.\end{eqnarray}$$

Equation (3.7a ) corresponds to that derived by Thorpe (Reference Thorpe1968) (his appendix 6).

We note in passing that the Stokes drift (3.7) is divergent ( $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{S}\neq 0$ ), which is generally true for Lagrangian velocities even in incompressible fluids. Indeed, we can readily confirm that (3.7) satisfies the identity for volume conservation $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{S}=(1/2)\unicode[STIX]{x2202}_{tzz}\overline{(\unicode[STIX]{x1D702}^{(1)})^{2}}$ from generalized Lagrangian-mean theory (equation (9.4) of Andrews & McIntyre (Reference Andrews and McIntyre1978)), shown here correct to leading order for our case (see also McIntyre (Reference McIntyre1988) for a discussion of the analogous effect for surface waves).

4.1.1 Linear solutions: $O(\unicode[STIX]{x1D6FC})$

Although the polarization relationships in table 1 hold in each layer, we must still identify explicit solutions for the vertical structure functions, including the matching condition between the layers. To do so, it is convenient to express the linear solution in terms of a velocity potential, which has the form

(4.3) $$\begin{eqnarray}\unicode[STIX]{x1D719}^{(1)}=\{C_{0}(X,T)\hat{\unicode[STIX]{x1D719}}_{0}(z)+\unicode[STIX]{x1D700}C_{1}(X,T)\hat{\unicode[STIX]{x1D719}}_{1}(z)\}\exp [\imath (kx-\unicode[STIX]{x1D714}t)]+O(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D700}^{2}),\end{eqnarray}$$

where we can show that $C_{0}=-\imath c_{p}A_{0}$ , $\hat{\unicode[STIX]{x1D719}}_{0}(z)=\hat{\unicode[STIX]{x1D702}}_{0}^{\prime }/k$ , $C_{1}=(c_{p}/k)A_{0,X}$ and $\hat{\unicode[STIX]{x1D719}}_{1}(z)=\left(\hat{\unicode[STIX]{x1D702}}_{1}^{\prime }+(2-\unicode[STIX]{x1D712})\hat{\unicode[STIX]{x1D702}}_{0}^{\prime }\right)/k$ (cf. table 1). We obtain for the vertical structure function

(4.4) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D702}}_{0}(z)=\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{\sinh (k(D-z))}{\sinh (kD_{+})}} & \unicode[STIX]{x1D6FF}\leqslant z\leqslant D,\\ {\displaystyle \frac{\sinh \left(k(D+z)\right)}{\sinh (kD_{-})}} & \unicode[STIX]{x1D6FF}>z\geqslant -D,\\ \end{array}\right.\end{eqnarray}$$

which is continuous in $z$ in order to satisfy the kinematic interfacial boundary condition (4.1). From the linearized dynamic interfacial boundary condition (4.2), we obtain

(4.5) $$\begin{eqnarray}g^{\prime }\unicode[STIX]{x1D702}_{I}^{(1)}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{+}^{(1)}}{\unicode[STIX]{x2202}t}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{-}^{(1)}}{\unicode[STIX]{x2202}t}}\quad \text{for }z=\unicode[STIX]{x1D6FF}.\end{eqnarray}$$

Hence $\hat{\unicode[STIX]{x1D719}}_{0}(z)$ is not continuous across the interface. We thus recover the well-known linear dispersion relation at $O(\unicode[STIX]{x1D6FC}^{1}\unicode[STIX]{x1D700}^{0})$ :

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=g^{\prime }k{\displaystyle \frac{1}{CT_{+}+CT_{-}}},\end{eqnarray}$$

in which we have defined $CT_{\pm }=\coth (kD_{\pm })$ for convenience.

Furthermore, it is readily evident that the pressure jump across the interface $\unicode[STIX]{x0394}p^{(1)}\equiv p_{+}^{(1)}(z=\unicode[STIX]{x1D6FF})-p_{-}^{(1)}(z=\unicode[STIX]{x1D6FF})=-\unicode[STIX]{x1D70C}_{ref}g^{\prime }A_{0}$ (from table 1 and (4.4)) and the linearized total pressures infinitesimally above and below the interface are equal ( $p_{tot,+}^{(1)}|_{z=\unicode[STIX]{x1D6FF}+\unicode[STIX]{x1D702}_{I}}=p_{+}^{(1)}|_{z=\unicode[STIX]{x1D6FF}}-\unicode[STIX]{x1D70C}_{+}g\unicode[STIX]{x1D702}_{I}^{(1)}=p_{tot,-}^{(1)}|_{z=\unicode[STIX]{x1D6FF}+\unicode[STIX]{x1D702}_{I}}=p_{-}^{(1)}|_{z=\unicode[STIX]{x1D6FF}}-\unicode[STIX]{x1D70C}_{-}g\unicode[STIX]{x1D702}_{I}^{(1)}$ ). The explicit $O(\unicode[STIX]{x1D6FC}^{1}\unicode[STIX]{x1D700}^{1})$ solution is given in appendix B.1 for completeness.

4.1.2 Stokes drift and transport

By substituting in the vertical structure function (4.4) into the general expression for the Stokes drift (3.7), we obtain

(4.7a,b ) $$\begin{eqnarray}\displaystyle & u_{S,\pm }={\displaystyle \frac{1}{2}}c_{p}k^{2}|A_{0}|^{2}{\displaystyle \frac{\cosh (2k(D\mp z))}{\sinh ^{2}(kD_{\pm })}},\quad w_{S,\pm }=\pm {\displaystyle \frac{1}{4}}c_{p}k(1+\unicode[STIX]{x1D712})\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{X}|A_{0}|^{2}{\displaystyle \frac{\sinh (2k(D\mp z))}{\sinh ^{2}(kD_{\pm })}}. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

We define the Stokes transport in each layer $Q_{ST,\pm }$ as a second-order Eulerian quantity given by vertically integrating the linear Eulerian velocity from the boundary to the time-varying linear interface $\unicode[STIX]{x1D702}_{I}$ . Explicitly,

(4.8) $$\begin{eqnarray}Q_{ST,\pm }\equiv \pm \overline{\int _{\unicode[STIX]{x1D6FF}+\unicode[STIX]{x1D702}_{I}^{(1)}}^{\pm D}u_{\pm }^{(1)}\,\text{d}z}=\mp \overline{\left.u_{\pm }^{(1)}\right|_{z=\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D702}_{I}^{(1)}}={\displaystyle \frac{1}{2}}\unicode[STIX]{x1D714}\coth (kD_{\pm })|A_{0}|^{2},\end{eqnarray}$$

in which the overline denotes averaging over the fast scales of the waves and we have only retained $O(\unicode[STIX]{x1D6FC}^{2})$ terms. To derive (4.8) we have used the polarization relationships in table 1 with the vertical structure function (4.4) and its derivative evaluated at $z=\unicode[STIX]{x1D6FF}$ , as given explicitly in table 2. Consequently, we note that the Stokes transport is equal to the vertical integral of the horizontal Stokes drift (4.7a ) in each layer.

Table 2. For interfacial waves in a two-layer Boussinesq fluid, expressions for the amplitude of the linear, $O(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D700}^{0})$ , fields in the top (plus sign) and bottom (minus sign) layers evaluated at $z=\unicode[STIX]{x1D6FF}$ and the relevant leading-order wave-averaged products between pairs of fields. Values are given in terms of the amplitude envelope of the interfacial displacement $A_{0}$ . The actual fields at the interface are found by multiplying by $\exp [\imath (kx-\unicode[STIX]{x1D714}t)]$ and taking the real part of the result.

4.1.3 Wave-induced Eulerian flow

The mean-flow forcing equation (2.2) simply reduces to the Laplace equation in both layers (cf. (3.8)). The forcing of an Eulerian mean flow instead comes from the two nonlinear interfacial boundary conditions. Following a Stokes expansion of the kinematic interfacial boundary condition $D\unicode[STIX]{x1D702}_{I}/Dt=w$ at $z=\unicode[STIX]{x1D6FF}+\unicode[STIX]{x1D702}_{I}$ , we have at second order in amplitude:

(4.9) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D719}}_{\pm }^{(2)}}{\unicode[STIX]{x2202}z}}-{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D702}}_{I}^{(2)}}{\unicode[STIX]{x2202}t}}=\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{X}\overline{u_{\pm }^{(1)}\unicode[STIX]{x1D702}_{I}^{(1)}}=\mp {\displaystyle \frac{1}{2}}\unicode[STIX]{x1D714}CT_{\pm }\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{X}|A_{0}|^{2}\quad \text{for }z=\unicode[STIX]{x1D6FF}.\end{eqnarray}$$

The forcing on the right-hand side corresponds to the divergence of the Stokes transport (4.8). Similarly, from a Stokes expansion of Bernoulli’s equation and equality of total pressures on either side of the interface, we have

(4.10) $$\begin{eqnarray}\displaystyle g^{\prime }\overline{\unicode[STIX]{x1D702}}_{I}^{(2)}+{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D719}}_{-}^{(2)}}{\unicode[STIX]{x2202}t}}-{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D719}}_{+}^{(2)}}{\unicode[STIX]{x2202}t}} & = & \displaystyle \overline{{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{+}^{(2)}}{\unicode[STIX]{x2202}z\unicode[STIX]{x2202}t}}\unicode[STIX]{x1D702}_{I}^{(1)}}+\overline{{\displaystyle \frac{1}{2}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{+}^{(1)}|^{2}}-\overline{{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{-}^{(2)}}{\unicode[STIX]{x2202}z\unicode[STIX]{x2202}t}}\unicode[STIX]{x1D702}_{I}^{(1)}}-\overline{{\displaystyle \frac{1}{2}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{-}^{(1)}|^{2}}\quad \text{for }z=\unicode[STIX]{x1D6FF},\nonumber\\ \displaystyle & = & \displaystyle {\displaystyle \frac{1}{4}}\unicode[STIX]{x1D714}^{2}(CT_{+}^{2}-CT_{-}^{2})|A_{0}|^{2}\quad \text{for }z=\unicode[STIX]{x1D6FF},\end{eqnarray}$$

in which the correlations in table 2 have been used to evaluate the forcing at the interface. We note that our (4.9), (4.10) are equivalent to (2.9) in Keady (Reference Keady1971) and (3.2)–(3.3) in Grimshaw & Pullin (Reference Grimshaw and Pullin1985). Because we assume a Boussinesq fluid, equations (4.9), (4.10) imply in the symmetric case that the set-down or set-up of the wave-averaged interface is zero ( $\overline{\unicode[STIX]{x1D702}}_{I}^{(2)}=0$ ). It is therefore not possible to recover from (4.9), (4.10) the case of surface gravity waves, in which the overlying density is zero (cf. non-Boussinesq) and the wave-averaged free surface sets down.

Although a general solution for the mean flow in both layers $\overline{\unicode[STIX]{x1D719}}_{+}^{(2)}$ and $\overline{\unicode[STIX]{x1D719}}_{-}^{(2)}$ and the wave-averaged interface $\overline{\unicode[STIX]{x1D702}}_{I}^{(2)}$ can be found, we additionally assume that the wavepacket extent is broad compared to the layer depths ( $\unicode[STIX]{x1D70E}\gg D_{\pm }$ ). In this limit, the Eulerian return flow is shallow and horizontally local to the overlying packet, unlike the more general non-local return flow that arises when this assumption is not made. This non-local return flow can be much wider than the overlying packet (cf. McIntyre (Reference McIntyre1981) for such non-local return flow for surface gravity waves). In the limit $\unicode[STIX]{x1D70E}\gg D_{\pm }$ , it is more convenient to find solutions in terms of the streamfunction. Explicitly, we assume solutions of the form $\overline{\unicode[STIX]{x1D713}}_{\pm }^{(2)}(\tilde{x},z)=(B_{\pm }\unicode[STIX]{x1D714}/H)(D\mp z)|A_{0}|^{2}$ and $\overline{\unicode[STIX]{x1D702}}_{I}^{(2)}(\tilde{x})=(B_{0}/H)|A_{0}|^{2}$ , in which $B_{\pm }$ and $B_{0}$ are constants to be determined and we have introduced the translating coordinate $\tilde{x}\equiv x-c_{g}t$ , which is related to the slow variable $X$ through $X=\unicode[STIX]{x1D700}\tilde{x}$ . Substituting these into (4.9), in which $\unicode[STIX]{x2202}_{z}\overline{\unicode[STIX]{x1D719}}_{\pm }^{(2)}=\unicode[STIX]{x2202}_{x}\overline{\unicode[STIX]{x1D713}}_{\pm }^{(2)}$ , gives two equations relating the coefficients $B_{0}$ and $B_{\pm }$ . Taking the $x$ -derivative of (4.10), using $\unicode[STIX]{x2202}_{x}\overline{\unicode[STIX]{x1D719}}_{\pm }^{(2)}=-\unicode[STIX]{x2202}_{z}\overline{\unicode[STIX]{x1D713}}_{\pm }^{(2)}$ , and substituting for $\overline{\unicode[STIX]{x1D713}}_{\pm }^{(2)}$ and $\overline{\unicode[STIX]{x1D702}}_{I}^{(2)}$ yields a third equation in the three coefficients. Solving the resulting system of simultaneous equations gives an explicit expression for $\overline{\unicode[STIX]{x1D713}}_{\pm }^{(2)}$ from which the horizontal induced Eulerian flow is found to be

(4.11) $$\begin{eqnarray}\overline{u}_{\pm }^{(2)}=\unicode[STIX]{x1D714}kf_{\pm }(kD_{+},kD_{-})|A_{0}|^{2},\end{eqnarray}$$

where

(4.12) $$\begin{eqnarray}f_{\pm }(kD_{+},kD_{-})\equiv {\displaystyle \frac{1}{4kH{\mathcal{D}}}}(CT_{+}+CT_{-})\left(2\unicode[STIX]{x1D712}^{2}+\unicode[STIX]{x1D712}kD_{\mp }(CT_{\mp }-CT_{\pm })-2kD_{\mp }CT_{\pm }\right),\end{eqnarray}$$

in which $\unicode[STIX]{x1D712}\equiv c_{g}/c_{p}$ ,

(4.13) $$\begin{eqnarray}{\mathcal{D}}=k\overline{D}(CT_{+}+CT_{-})-\unicode[STIX]{x1D712}^{2},\end{eqnarray}$$

and $\overline{D}=D_{+}D_{-}/H$ . We can also obtain the wave-averaged interface:

(4.14) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D702}}_{I}^{(2)}={\displaystyle \frac{1}{4H{\mathcal{D}}}}\left(2\unicode[STIX]{x1D712}k(D_{-}\,CT_{+}-D_{+}\,CT_{-})+k^{2}\overline{D}H(CT_{+}^{2}-CT_{-}^{2})\right)|A_{0}|^{2}.\end{eqnarray}$$

In the symmetric case ( $D\equiv D_{+}=D_{-}$ ), it is clear from (4.14) that $\overline{\unicode[STIX]{x1D702}}_{I}^{(2)}=0$ . In this case $f_{\pm }(kD,kD)=-1/(2kD\tanh (kD))$ and hence the total volume flux $Q_{S}+Q_{E}$ is zero in each layer.

In the asymmetric case, the induced Eulerian flow is changed significantly, particularly in the long-wave limit $kD_{\pm }\ll 1$ , in which case $CT_{\pm }\simeq 1/(kD_{\pm })+(kD_{\pm })/3$ , $\unicode[STIX]{x1D712}\simeq 1-k^{2}\overline{D}H/3$ and, from (4.13), ${\mathcal{D}}\simeq k^{2}\overline{D}H$ , in which $\overline{D}=D_{+}D_{-}/H$ and $H=D_{+}+D_{-}$ is the total domain depth. Thus the induced Eulerian flow becomes singular in this limit and likewise the Eulerian transport in each layer is singular:

(4.15) $$\begin{eqnarray}{\displaystyle \frac{Q_{E,\pm }}{|A_{0}|^{2}\sqrt{g^{\prime }/H}}}=\mp {\displaystyle \frac{12}{(kH)^{2}}}{\displaystyle \frac{\unicode[STIX]{x1D6FF}/H}{\left(1-4(\unicode[STIX]{x1D6FF}/H)^{2}\right)^{3/2}}}+O(\hat{k}^{0}).\end{eqnarray}$$

It is evident from (4.15) that the long-wave singularity only occurs for unequal layer depths ( $\hat{\unicode[STIX]{x1D6FF}}\neq 0$ ). As a consequence, the induced Eulerian transport becomes very large for shallow interfacial wavepackets in an asymmetric two-layer fluid.

This is illustrated in figure 2, which plots the transport in the top layer as a function of $kH$ and $\unicode[STIX]{x1D6FF}/H$ , with $\unicode[STIX]{x1D6FF}\equiv (D_{-}-D_{+})/2$ denoting the interfacial position away from the centre (see figure 1 a). While the Stokes transport (figure 2 a) shows a modest increase with decreasing $kH$ and increasing $\unicode[STIX]{x1D6FF}/H$ , the induced Eulerian flow becomes significantly negative if $\unicode[STIX]{x1D6FF}/H>0$ (shallower upper-layer depth) and significantly positive if $\unicode[STIX]{x1D6FF}/H<0$ (deeper upper-layer depth). Correspondingly, the wave-averaged interface displaces significantly into the shallowest layer. The set-up or set-down of the wave-averaged interface thus enhances the (negative) return flow in the shallowest layer.

Figure 2. For interfacial waves in a two-layer Boussinesq fluid, regime diagram for the total volume flux $(Q_{S}+Q_{E})/(\sqrt{g^{\prime }/H}|A_{0}|^{2})$ in the top (a) and bottom (b) layers, their sum (c), as well as the magnitude of the wave-averaged interface (d). Note that the colour scale is chosen to saturate for large negative values.

4.1.4 Lagrangian displacements

As a wavepacket passes horizontally, Lagrangian fluid parcels are displaced according to the sum of the Stokes drift $u_{S}$ and the Eulerian flow $u_{E}$ . For the particular case of a wavepacket with a horizontally Gaussian amplitude envelope $A_{0}(\tilde{x})=a_{0}\exp [-\tilde{x}^{2}/(2\unicode[STIX]{x1D70E}^{2})]$ , the two associated net displacements can be obtained from integrating (4.7a ) and (4.11) with respect to time between $t\rightarrow \pm \infty$

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}x_{S}=u_{S}\left(\tilde{x}=0\right)\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}/c_{g}={\displaystyle \frac{\unicode[STIX]{x1D70E}}{2}}\sqrt{\unicode[STIX]{x03C0}}(k|a_{0}|)^{2}{\displaystyle \frac{\cosh 2k(D_{\pm }\mp z)}{\unicode[STIX]{x1D712}\sinh ^{2}kD_{\pm }}}, & \displaystyle\end{eqnarray}$$

(4.17) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}x_{E}=u_{E}\left(\tilde{x}=0\right)\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}/c_{g}={\displaystyle \frac{\unicode[STIX]{x1D70E}}{2}}\sqrt{\unicode[STIX]{x03C0}}(k|a_{0}|)^{2}{\displaystyle \frac{f_{\pm }(kD_{+},kD_{-})}{\unicode[STIX]{x1D712}}}. & \displaystyle\end{eqnarray}$$

4.2.1 Linear solutions: $O(\unicode[STIX]{x1D6FC})$

As for the two-layer case, the polarization relationships in table 1 hold in each of the two unstratified layers as well as the stratified layer. We proceed to identify explicit solutions for the vertical structure functions, including the matching condition between the layers. Using (3.4), the vertical structure function $\hat{\unicode[STIX]{x1D702}}_{0}$ satisfies

(4.19a,b ) $$\begin{eqnarray}\left.\begin{array}{@{}ll@{}}\hat{\unicode[STIX]{x1D702}}_{0}^{\prime \prime }-k^{2}\hat{\unicode[STIX]{x1D702}}_{0}=0\quad & \text{for }d<|z|<D,\\ \hat{\unicode[STIX]{x1D702}}_{0}^{\prime \prime }+\unicode[STIX]{x1D6FE}^{2}\hat{\unicode[STIX]{x1D702}}_{0}=0\quad & \text{for }|z|\leqslant d,\end{array}\right\}\end{eqnarray}$$

in which $\unicode[STIX]{x1D6FE}=k\,(N_{0}^{2}/\unicode[STIX]{x1D714}^{2}-1)^{1/2}$ . Explicit solutions to (4.19a,b ) are found by applying matching conditions at the ‘interface’ between the stratified and unstratified fluid, namely at $z=\pm d$ in a linear approximation. Explicitly, the kinematic and dynamic interface conditions require that both $\hat{\unicode[STIX]{x1D702}}_{0}$ and its vertical derivative $\hat{\unicode[STIX]{x1D702}}_{0}^{\prime }$ are continuous at $z=\pm d$ , as is evident from the entries for the vertical displacement and pressure in table 1, respectively. These two pairs of matching conditions at the interfaces together with the requirement that $w=0$ at $z=\pm D$ give the following vertical structure function for even modes:

(4.20) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D702}}_{0}=\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{\cos \unicode[STIX]{x1D6FE}d}{\sinh k\unicode[STIX]{x0394}}}\sinh k(D-z)\quad & d<z\leqslant D,\\ \cos \unicode[STIX]{x1D6FE}z\quad & |z|\leqslant d,\\ {\displaystyle \frac{\cos \unicode[STIX]{x1D6FE}d}{\sinh k\unicode[STIX]{x0394}}}\sinh k(D+z)\quad & -d>z\geqslant -D,\end{array}\right.\end{eqnarray}$$

in which $\unicode[STIX]{x0394}=D-d$ . Furthermore, they give the dispersion relation, given implicitly by the condition

(4.21) $$\begin{eqnarray}k\coth k\unicode[STIX]{x0394}=\unicode[STIX]{x1D6FE}\tan \unicode[STIX]{x1D6FE}d.\end{eqnarray}$$

In the limit of uniform stratification ( $d\rightarrow D$ , $\unicode[STIX]{x0394}\rightarrow 0$ ), we have $\unicode[STIX]{x1D6FE}D=(2j-1)\unicode[STIX]{x03C0}/2$ for positive integers $j$ . Thus for the lowest even mode ( $j=1$ ) we have $\unicode[STIX]{x1D6FE}\rightarrow m=\unicode[STIX]{x03C0}/H$ . In the limit of a two-layer fluid ( $d\rightarrow 0$ , $\unicode[STIX]{x0394}\rightarrow D$ ), we have $\tan \unicode[STIX]{x1D6FE}d\simeq \unicode[STIX]{x1D6FE}d$ . Thus (4.21) becomes

(4.22) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}\sim (k\coth kD)^{1/2}d^{-1/2}\end{eqnarray}$$

and the dispersion relation for the two-layer fluid ((4.6) with $D=D_{+}=D_{-}$ ) can be recovered. Although we will not examine the case of odd modes in detail, it is relevant for the discussion below to note that the dispersion relation for odd modes (whose streamfunction varies as $\sin \unicode[STIX]{x1D6FE}z$ for $|z|\leqslant d$ ) satisfies

(4.23) $$\begin{eqnarray}k\coth k\unicode[STIX]{x0394}=-\unicode[STIX]{x1D6FE}\cot \unicode[STIX]{x1D6FE}d.\end{eqnarray}$$

For completeness, the explicit $O(\unicode[STIX]{x1D6FC}^{1}\unicode[STIX]{x1D700}^{1})$ solution for the vertical structure function $\hat{\unicode[STIX]{x1D702}}_{1}$ is given in appendix B.2.

4.2.2 Stokes drift and transport

By substituting the vertical structure function (4.20) into the general expression for the Stokes drift (3.7), we can obtain the Stokes drift for even modes in top-hat stratification

(4.24) $$\begin{eqnarray}u_{S}={\displaystyle \frac{1}{2}}c_{p}k^{2}|A_{0}|^{2}\left\{\begin{array}{@{}ll@{}}\left[{\displaystyle \frac{\cos \unicode[STIX]{x1D6FE}d}{\sinh k\unicode[STIX]{x0394}}}\right]^{2}\cosh 2k(D-z)\quad & d<z\leqslant D,\\ -{\displaystyle \frac{\unicode[STIX]{x1D6FE}^{2}}{k^{2}}}\cos 2\unicode[STIX]{x1D6FE}z\quad & |z|\leqslant d,\\ \left[{\displaystyle \frac{\cos \unicode[STIX]{x1D6FE}d}{\sinh k\unicode[STIX]{x0394}}}\right]^{2}\cosh 2k(D+z)\quad & -d>z\geqslant -D.\end{array}\right.\end{eqnarray}$$

Even though the vertical structure function $\hat{\unicode[STIX]{x1D702}}_{0}$ in (4.20) is continuously differentiable at $|z|=d$ , the second derivative is discontinuous. Thus the Stokes drift has jump discontinuities where the background density gradient changes discontinuously.

The (Eulerian) Stokes transport (as defined in (4.8)) in each layer is found in general to be

(4.25) $$\begin{eqnarray}Q_{ST}=\left\{\begin{array}{@{}ll@{}}{\textstyle \frac{1}{4}}\unicode[STIX]{x1D6FE}c_{p}\sin \left(2\unicode[STIX]{x1D6FE}d\right)|A_{0}|^{2}\quad & d<|z|\leqslant D,\\ -{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FE}c_{p}\sin \left(2\unicode[STIX]{x1D6FE}d\right)|A_{0}|^{2}\quad & |z|\leqslant d,\end{array}\right.\end{eqnarray}$$

which, using (4.21), equals the vertical integral of the Stokes drift (4.24).

In the two-layer limit for which $d\ll D$ with $N_{0}^{2}(2d)=g^{\prime }$ kept constant, we find $\unicode[STIX]{x1D6FE}\sim \sqrt{k\coth (kD)/d}$ . Thus where $d<|z|<D$ , equation (4.24) recovers the prediction (4.7) for the Stokes drift in a symmetric two-layer fluid. However, unlike the case of a two-layer fluid, the total Stokes transport is zero. This indicates that the Stokes transport in the stratified layer is retrograde to the flows above and below and that it becomes infinitely large as the layer becomes infinitesimally thin. Indeed, evaluating (4.24) at $z=0$ assuming $d/D\ll 1$ gives $u_{S}(z=0)\simeq -[g^{\prime }k/(4\unicode[STIX]{x1D714}d)]\,|A_{0}|^{2}$ , although the total Stokes transport in this layer $Q_{ST}=-[\unicode[STIX]{x1D714}/\!\tanh kD]|A_{0}|^{2}$ is finite. This is an indication of singular behaviour associated with transport properties in the limit of an infinitesimally thin interface. Specifically, the mode-2-like structure evident in (4.24) is not permitted in a two-layer fluid.

Finally, we note that in the limit of uniform stratification ( $d\rightarrow D$ , $\unicode[STIX]{x1D6FE}\rightarrow m$ ), equation (4.24) gives

(4.26) $$\begin{eqnarray}u_{S}=-{\displaystyle \frac{1}{2}}{\displaystyle \frac{m^{2}}{|\boldsymbol{k}|}}\,N_{0}|A_{0}|^{2}\cos (2mz),\end{eqnarray}$$

in which $\unicode[STIX]{x1D6FE}\rightarrow m=\unicode[STIX]{x03C0}/H$ is the vertical mode number of the lowest even mode. This is the result previously found by Thorpe (Reference Thorpe1968).

4.2.3 Wave-induced Eulerian flow

We focus our attention on the case in which the wavepacket is long relative to the depth of the unstratified layer ( $\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D6E5}\gg 1$ ). In that case, we can seek a separable solution for the Eulerian mean flow of the form $\overline{\unicode[STIX]{x1D713}}^{(2)}=N_{0}|A_{0}(X)|^{2}\overline{\unicode[STIX]{x1D6F9}}(z)$ , in which $N_{0}$ is the characteristic buoyancy frequency and the non-dimensional vertical structure $\overline{\unicode[STIX]{x1D6F9}}(z)$ is given by the solution of (3.13). Solving separately in each layer specifically for even-mode wavepackets gives the general solution

(4.27) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D6F9}}=\left\{\begin{array}{@{}ll@{}}b_{+}\,\unicode[STIX]{x1D707}\,(D-z) & d<z\leqslant D,\\ b_{0}\sin \unicode[STIX]{x1D707}z+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}\sin 2\unicode[STIX]{x1D6FE}z & |z|\leqslant d,\\ b_{-}\,\unicode[STIX]{x1D707}\,(D+z) & -d>z\geqslant -D,\\ \end{array}\right.\end{eqnarray}$$

in which

(4.28) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}={\displaystyle \frac{1}{4}}{\displaystyle \frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D707}}}{\displaystyle \frac{1}{\unicode[STIX]{x1D712}}}\left(1+2\unicode[STIX]{x1D712}\right)(1-\unicode[STIX]{x1D712})\left[\left({\displaystyle \frac{2\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D707}}}\right)^{2}-1\right]^{-1},\end{eqnarray}$$

$\unicode[STIX]{x1D707}\equiv N_{0}/c_{g}$ and, as before, $\unicode[STIX]{x1D712}\equiv c_{g}/c_{p}$ . The constants $b_{0}$ and $b_{\pm }$ in (4.27) are found by kinematic and dynamic matching conditions at $z=\pm d$ evaluated at order amplitude squared.

Explicitly, the kinematic condition requires the vertical displacement to be continuous across each interface. Using (3.10) at $|z|=d$ we have

(4.29) $$\begin{eqnarray}\left.\overline{\unicode[STIX]{x1D702}}^{(2)}\right|_{z=\pm d}=\left.\left[-\unicode[STIX]{x1D707}\unicode[STIX]{x1D6F9}+{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{1}{\unicode[STIX]{x1D712}}}-1\right){\displaystyle \frac{\text{d}\hat{\unicode[STIX]{x1D702}}_{0}^{2}}{\text{d}z}}\right]\right|_{z=\pm d}|A_{0}|^{2},\end{eqnarray}$$

where we can ignore the singularity in $\text{d}N^{2}/\text{d}z$ at $|z|=d$ , as $\text{d}N^{2}/\text{d}z=0$ in both layers on either side of the interface. Given that $\hat{\unicode[STIX]{x1D702}}_{0}$ and $\hat{\unicode[STIX]{x1D702}}_{0}^{\prime }$ vary continuously, we must have from (4.29) that $\overline{\unicode[STIX]{x1D6F9}}$ is continuous at $z=\pm d$ .

The dynamic matching condition requires that total pressure at order amplitude squared is continuous. From the horizontal momentum equation, expressed here in terms of the total pressure for convenience, we determine the order amplitude-squared horizontal pressure gradient at $z=\pm d+\unicode[STIX]{x1D702}$ from a Stokes expansion around $z=\pm d$

(4.30) $$\begin{eqnarray}\displaystyle & {\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{ref}}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left.\left[p_{tot}^{(2)}+{\displaystyle \frac{\unicode[STIX]{x2202}p_{tot}^{(1)}}{\unicode[STIX]{x2202}z}}\unicode[STIX]{x1D702}^{(1)}\right]\right|_{z=\pm d}=\left.\left[-{\displaystyle \frac{\unicode[STIX]{x2202}u^{(2)}}{\unicode[STIX]{x2202}t}}-{\displaystyle \frac{\unicode[STIX]{x2202}^{2}u^{(1)}}{\unicode[STIX]{x2202}z\,\unicode[STIX]{x2202}t}}\unicode[STIX]{x1D702}^{(1)}-{\displaystyle \frac{\unicode[STIX]{x2202}u^{(1)}u^{(1)}}{\unicode[STIX]{x2202}x}}-{\displaystyle \frac{\unicode[STIX]{x2202}u^{(1)}w^{(1)}}{\unicode[STIX]{x2202}z}}\right]\right|_{z=\pm d}. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Using the polarization relations in table 1, averaging over the fast scales and integrating over one slow $X$ derivative gives after some manipulation

(4.31) $$\begin{eqnarray}{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{ref}}}\left.\left[\overline{p}_{tot}^{(2)}+\overline{{\displaystyle \frac{\unicode[STIX]{x2202}p_{tot}^{(1)}}{\unicode[STIX]{x2202}z}}\unicode[STIX]{x1D702}^{(1)}}\right]\right|_{z=\pm d}=\left.\left[-c_{g}N_{0}\overline{\unicode[STIX]{x1D6F9}}^{\prime }+{\displaystyle \frac{1}{4}}c_{p}^{2}\left(2\unicode[STIX]{x1D712}\hat{\unicode[STIX]{x1D702}}_{0}\hat{\unicode[STIX]{x1D702}}_{0}^{\prime \prime }-(\hat{\unicode[STIX]{x1D702}}_{0}^{\prime })^{2}\right)\right]\right|_{z=\pm d^{\mp }}|A_{0}|^{2},\end{eqnarray}$$

where we note that only the zeroth-order vertical structure function $\hat{\unicode[STIX]{x1D702}}_{0}$ plays a role, as for the mean-flow forcing equation (3.8). Now, insisting that the left-hand side of (4.31) is continuous at the interfaces and noting that $\hat{\unicode[STIX]{x1D702}}_{0}^{\prime \prime }$ is discontinuous at $z=\pm d$ , we find that $\overline{\unicode[STIX]{x1D6F9}}^{\prime }$ is discontinuous at the interfaces such that the jump in $\overline{\unicode[STIX]{x1D6F9}}^{\prime }$ from the unstratified to stratified side of each interface is

(4.32) $$\begin{eqnarray}J\equiv \overline{\unicode[STIX]{x1D6F9}}^{\prime }|_{\text{d}^{+}}-\overline{\unicode[STIX]{x1D6F9}}^{\prime }|_{\text{d}^{-}}=\overline{\unicode[STIX]{x1D6F9}}^{\prime }|_{-d^{-}}-\overline{\unicode[STIX]{x1D6F9}}^{\prime }|_{-d^{+}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{N_{0}}{c_{p}}}\cos ^{2}\unicode[STIX]{x1D6FE}d.\end{eqnarray}$$

The discontinuity of $\overline{\unicode[STIX]{x1D6F9}}^{\prime }$ indicates that, like the Stokes drift, the Eulerian induced flow jumps discontinuously at the interfaces. From (4.24), the jump in the Stokes drift is of equal but opposite sign, so that the Lagrangian velocity $u_{L}=u_{S}+u_{E}$ is continuous across the interface.

Applying these two matching conditions to the general solution for $\overline{\unicode[STIX]{x1D6F9}}$ in (4.27), we find

(4.33) $$\begin{eqnarray}\displaystyle & \displaystyle b_{0}=-{\displaystyle \frac{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}(\sin 2\unicode[STIX]{x1D6FE}d+2\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394}\cos 2\unicode[STIX]{x1D6FE}d)+J\unicode[STIX]{x0394}}{\sin \unicode[STIX]{x1D707}d+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\cos \unicode[STIX]{x1D707}d}}, & \displaystyle\end{eqnarray}$$

(4.34) $$\begin{eqnarray}\displaystyle & \displaystyle b_{\pm }=\pm {\displaystyle \frac{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}(\sin 2\unicode[STIX]{x1D6FE}d\cos \unicode[STIX]{x1D707}d-(2\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D707})\cos 2\unicode[STIX]{x1D6FE}d\sin \unicode[STIX]{x1D707}d)-(J/\unicode[STIX]{x1D707})\sin \unicode[STIX]{x1D707}d}{\sin \unicode[STIX]{x1D707}d+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\cos \unicode[STIX]{x1D707}d}}. & \displaystyle\end{eqnarray}$$

With these expressions the predicted induced Eulerian flow is

(4.35) $$\begin{eqnarray}u_{E}=c_{g}\,\unicode[STIX]{x1D707}^{2}|A_{0}|^{2}\left\{\begin{array}{@{}ll@{}}\pm b_{\pm } & |d|<z\leqslant D,\\ -\!\left(b_{0}\cos \unicode[STIX]{x1D707}z+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}(2\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D707})\cos 2\unicode[STIX]{x1D6FE}z\right) & |z|\leqslant d.\end{array}\right.\end{eqnarray}$$

Profiles of the induced Eulerian flow and Stokes drift together with their sum (the Lagrangian flow), are plotted in three cases with symmetric top-hat stratification in the first row of figure 3. In all three cases, $u_{E}$ and $u_{S}$ have opposite sign at the middle, top and bottom of the domain, though it is not always the case that these flows act destructively to produce the Lagrangian flow. Generally, it is true that $u_{S}|_{z=0}<0$ and $u_{S}|_{z=\pm D}>0$ and that $u_{S}|_{z=0}$ becomes increasingly negative as $d/D$ becomes smaller. For long waves ( $kH\ll 1$ ), $u_{E}$ is opposite-signed to $u_{S}$ at the middle, top and bottom of the domain. However, $u_{E}$ can change sign and indeed become very large for moderate $kH$ . This dynamics is examined in detail below.

Figure 3. For top-hat stratification, vertical profiles of the Stokes drift (red), the induced Eulerian flow (blue) and the total Lagrangian flow (black) with indicated values of non-dimensional wavenumber $kH$ , half-depth of the stratified layer $d$ and upward displacement of the mid-point of the stratified layer from $z=0$ , $\unicode[STIX]{x1D6FF}$ . All velocities are normalized by $|A_{0}|^{2}N_{0}/H$ , as denoted by the tilde. Panels (a–c) correspond to symmetric cases, whereas (d–f) show asymmetric cases.

In the uniform stratification limit, $d\rightarrow D$ , equation (4.35) recovers the result previously found by McIntyre (Reference McIntyre1973) and Grimshaw (Reference Grimshaw1977):

(4.36) $$\begin{eqnarray}u_{E}={\displaystyle \frac{1}{2}}{\displaystyle \frac{m^{2}}{|\boldsymbol{k}|}}{\displaystyle \frac{(1+2m^{2}/|\boldsymbol{k}|^{2})(1-m^{2}/|\boldsymbol{k}|^{2})}{1-4m^{6}/|\boldsymbol{k}|^{6}}}\,N_{0}|A_{0}|^{2}\cos (2mz),\end{eqnarray}$$

in which $m=\unicode[STIX]{x03C0}/H$ is the vertical wavenumber of the mode-1 wavepacket. Different from the Stokes drift, $u_{E}$ exhibits a singularity at a critical horizontal wavenumber $k_{c}\equiv (4^{1/3}-1)^{1/2}m\simeq 0.766m$ . At this value, corresponding to $\unicode[STIX]{x1D707}\equiv N_{0}/c_{g}=2\unicode[STIX]{x1D6FE}\rightarrow 2m$ , the group velocity of the wavepacket is equal to the horizontal phase and group speed of the long mode-2 induced flow, $N_{0}/(2m)$ . Thus the singularity is an indication of the induced flow being resonantly excited (McIntyre Reference McIntyre1973). For $k>k_{c}$ ( $k<k_{c}$ ) the mode-2 disturbance flow has faster (slower) speed than the wavepacket. Crucially, (4.36) predicts a change in sign of the induced flow as $k$ changes from being greater than $k_{c}$ . For $k>k_{c}$ , $u_{S}$ and $u_{E}$ are opposite signed. If $k<k_{c}$ , on the other hand, the flows add constructively to make up the Lagrangian flow.

In the two-layer fluid limit, it follows that $b_{\pm }\rightarrow \mp N_{0}d/(c_{p}\unicode[STIX]{x1D707}D)$ and $u_{E}\rightarrow -\unicode[STIX]{x1D714}|A_{0}|^{2}/(2D\tanh kD)$ . The induced Eulerian flow for $|z|>d$ in this limit is thus equal to that found for the Eulerian flow in a two-layer fluid (4.11) (provided $D\ll \unicode[STIX]{x1D70E}$ , in both cases). On the other hand, the Eulerian velocity for $|z|<d$ becomes infinitely large, while the volume flux remains finite, as follows because the vertically integrated Eulerian flow in any top-hat stratification is always zero. Like the Stokes drift, taking $d\rightarrow 0$ is a singular limit owing to the mode-2-like vertical structure of (4.27) and (4.35), which is not permitted in a two-layer fluid.

Interestingly, while $\unicode[STIX]{x1D707}=2m$ gave rise to a singularity in the case of uniformly stratified fluid, the corresponding apparent singularity occurring for $\unicode[STIX]{x1D707}=2\unicode[STIX]{x1D6FE}$ disappears if $d<D$ . Nevertheless, new singularities appear corresponding to zeros of the denominator, $\sin \unicode[STIX]{x1D707}d+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\cos \unicode[STIX]{x1D707}d$ , appearing in the coefficients $b_{0}$ and $b_{\pm }$ in (4.33) and (4.34). By examination of the odd-mode dispersion relation (4.23) in the long-wave limit (as appropriate for the horizontally long, odd-mode structures induced by the wavepacket, cf. (4.27)), it is clear that the singularities occur when the induced flow corresponds to an odd mode with $\unicode[STIX]{x1D6FE}$ in (4.23) replaced by $\unicode[STIX]{x1D707}=N_{0}/c_{g}$ . These singularities occur if $kH$ is larger than $1$ , in which case $c_{g}$ is small and $\unicode[STIX]{x1D707}$ is correspondingly large.

In general, the structure of the induced Eulerian flow as it depends upon $kH$ can be quite complex, as illustrated figure 4, which plots versus $kH$ values of $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D707}$ , the predicted flow at the centre and top of the domain, as well as the velocity jump at $|z|=d$ in two circumstances for which the relative depth of the stratified region is $d=0.5D$ and $0.1D$ . The induced flows at $z=0$ and $z=D$ show the resonance peaks occurring for $kH\simeq 3$ and $5$ for $d=0.5D$ . Multiple resonance peaks occur for still larger $kH$ in both cases with $d=0.5D$ and $0.1D$ .

Figure 4. For symmetric top-hat stratification, (a,b) dispersion relation coefficients showing $\unicode[STIX]{x1D707}$ (solid lines) and $\unicode[STIX]{x1D6FE}$ (long-dashed lines). Below are plotted normalized values associated with the predicted induced Eulerian flow (long-dashed lines, blue) and the Stokes drift (short-dashed lines, red) and their sum (solid black lines) (c,d) at the vertical centre of the domain, (e,f) at the top of the domain, and (g,h) also showing the velocity jump at the interfaces at $z=\pm d$ where $N^{2}$ changes from $N_{0}^{2}$ to $0$ . Flow values are normalized by $|A_{0}|^{2}N_{0}/H$ , as denoted by the tilde. In (a,c,e,g) $d=0.5D=0.25H$ and in (b,d,f,h) $d=0.1D=0.05H$ . The upward arrows on the bottom axes indicate particular values of $kH$ examined in theory and numerical simulations (see § 5).

As in the case of the Stokes transport for wavepackets in this three-layer stratification, the transport due to the Eulerian flow in the middle layer ( $|z|<d$ ), $Q_{E}=-2c_{g}\unicode[STIX]{x1D707}^{2}b_{+}\unicode[STIX]{x1D6E5}|A_{0}|^{2}$ , is equal and opposite to the sum of the Eulerian transport in the upper and lower layers. However, the depth-integrated flows $Q_{S}$ and $Q_{E}$ in each layer are opposite but not equal and opposite, reflecting displacements of the wave-averaged interfaces. Hence there is differential Lagrangian transport over the depth of the passing wavepacket. Such effects will be considered in more detail in the following examination of flows induced in asymmetric top-hat stratification.

4.3.1 Linear solutions: $O(\unicode[STIX]{x1D6FC})$

The vertical structure of the waves is given generally by

(4.37) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D702}}_{0}=\left\{\begin{array}{@{}ll@{}}B_{+}\sinh k(D-z) & d+\unicode[STIX]{x1D6FF}<z\leqslant D,\\ B_{0}\cos \unicode[STIX]{x1D6FE}z+B_{1}\sin \unicode[STIX]{x1D6FE}z & -d+\unicode[STIX]{x1D6FF}\leqslant z\leqslant d+\unicode[STIX]{x1D6FF},\\ B_{-}\sinh k(D+z) & -D\leqslant z<-d+\unicode[STIX]{x1D6FF}.\end{array}\right.\end{eqnarray}$$

Seeking solutions for even modes we take $B_{0}=1$ . As will be shown, $B_{1}$ is small if $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}$ is small.

For notational convenience, we define the following:

(4.38) $$\begin{eqnarray}\displaystyle \hspace{-3.0pt}\left.\begin{array}{@{}llll@{}}C_{\pm }=\cosh [k\unicode[STIX]{x1D6E5}(1\pm \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}})], & C_{0}=\cosh (k\unicode[STIX]{x1D6E5}), & S_{\pm }=\sinh [k\unicode[STIX]{x1D6E5}(1\pm \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}})], & S_{0}=\sinh (k\unicode[STIX]{x1D6E5}),\\ c_{\pm }=\cos [\unicode[STIX]{x1D6FE}d(1\pm \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6E5}/d)], & c_{0}=\cos (\unicode[STIX]{x1D6FE}d), & s_{\pm }=\sin [\unicode[STIX]{x1D6FE}d(1\pm \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6E5}/d)], & s_{0}=\sin (\unicode[STIX]{x1D6FE}d).\end{array}\right\}\hspace{-15.00002pt} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Continuity of $\hat{\unicode[STIX]{x1D702}}_{0}$ and its derivative gives the implicit dispersion relation

(4.39) $$\begin{eqnarray}\displaystyle 0 & = & \displaystyle (kC_{-}c_{+}-\unicode[STIX]{x1D6FE}S_{-}s_{+})(kC_{+}s_{-}+\unicode[STIX]{x1D6FE}S_{+}c_{-})+(kC_{-}s_{+}+\unicode[STIX]{x1D6FE}S_{-}c_{+})(kC_{+}c_{-}-\unicode[STIX]{x1D6FE}S_{+}s_{-})\nonumber\\ \displaystyle & \simeq & \displaystyle 2(kC_{0}c_{0}-\unicode[STIX]{x1D6FE}S_{0}s_{0})(kC_{0}s_{0}+\unicode[STIX]{x1D6FE}S_{0}c_{0})+O(\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}^{2}).\end{eqnarray}$$

In the last expression the first term in parenthesis gives the dispersion relation for even modes, while the second term gives the dispersion relation for odd modes in the case $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}=0$ . That the correction for finite $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}$ enters at order $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}^{2}$ demonstrates that the dispersion relation is relatively insensitive to breaking symmetry.

Imposing continuity of $\hat{\unicode[STIX]{x1D702}}_{0}$ and its derivatives gives the following expressions for the coefficients in terms of $B_{0}$ :

(4.40) $$\begin{eqnarray}\displaystyle & \displaystyle B_{1}=B_{0}{\displaystyle \frac{kC_{+}c_{-}-\unicode[STIX]{x1D6FE}S_{+}s_{-}}{kC_{+}s_{-}+\unicode[STIX]{x1D6FE}S_{+}c_{-}}}\simeq \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}B_{0}{\displaystyle \frac{\unicode[STIX]{x1D6E5}(k^{2}+\unicode[STIX]{x1D6FE}^{2})S_{0}c_{0}}{kC_{0}s_{0}+\unicode[STIX]{x1D6FE}S_{0}c_{0}}}+O(\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}^{2}), & \displaystyle\end{eqnarray}$$

(4.41) $$\begin{eqnarray}\displaystyle & \displaystyle B_{\pm }=B_{0}{\displaystyle \frac{\unicode[STIX]{x1D6FE}}{kC_{\mp }s_{\pm }+\unicode[STIX]{x1D6FE}S_{\mp }c_{\pm }}}\simeq B_{0}{\displaystyle \frac{\unicode[STIX]{x1D6FE}}{k_{0}C_{0}s_{0}+\unicode[STIX]{x1D6FE}S_{0}c_{0}}}\left[\!1\mp \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}{\displaystyle \frac{\unicode[STIX]{x1D6E5}(k^{2}+\unicode[STIX]{x1D6FE}^{2})S_{0}s_{0}}{kC_{0}s_{0}+\unicode[STIX]{x1D6FE}S_{0}c_{0}}}\!\right]+O(\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}^{2}),\qquad & \displaystyle\end{eqnarray}$$

in which we note that the denominator in these expressions is non-zero for even modes in the long-wave limit. Hence $B_{1}$ is small if $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}$ is small, as expected. Consequently, from (3.7a ), the Stokes drift and transport should be similar to that for the symmetric case given by (4.24).

4.3.2 Wave-induced Eulerian flow

The vertical structure of the Eulerian induced flow is given generally by

(4.42) $$\begin{eqnarray}\displaystyle \overline{\unicode[STIX]{x1D6F9}}=\left\{\begin{array}{@{}ll@{}}b_{+}\unicode[STIX]{x1D707}(D-z) & d+\unicode[STIX]{x1D6FF}<z\leqslant D,\\ b_{0}\sin \unicode[STIX]{x1D707}z+b_{1}\cos \unicode[STIX]{x1D707}z+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}[(B_{0}^{2}-B_{1}^{2})\sin 2\unicode[STIX]{x1D6FE}z-2B_{0}B_{1}\cos 2\unicode[STIX]{x1D6FE}z] & -d+\unicode[STIX]{x1D6FF}\leqslant z\leqslant d+\unicode[STIX]{x1D6FF},\\ b_{-}\unicode[STIX]{x1D707}(D+z) & -D\leqslant z<-d+\unicode[STIX]{x1D6FF}.\\ \end{array}\right.\hspace{-15.00002pt} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

We require continuity of $\overline{\unicode[STIX]{x1D6F9}}$ and jumps in $\overline{\unicode[STIX]{x1D6F9}}^{\prime }$ at the interfaces by analogy with (4.32) in which $J_{\pm }=N_{0}/(2c_{p})B_{\pm }^{2}S_{\mp }$ are the jumps at the upper (upper signs) and lower (lower signs) interface. For $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}$ small, these are given approximately by $J_{\pm }\simeq B_{0}^{2}(J_{0}\pm \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}J_{1})$ , in which

(4.43a,b ) $$\begin{eqnarray}\displaystyle J_{0}={\displaystyle \frac{1}{2}}N_{0}c_{p}S_{0}^{2}(kC_{0}s_{0}+\unicode[STIX]{x1D6FE}S_{0}c_{0})^{-2}\quad \text{and}\quad J_{1}=2J_{0}\left({\displaystyle \frac{(k^{2}+\unicode[STIX]{x1D6FE}^{2})\unicode[STIX]{x0394}S_{0}s_{0}}{kC_{0}s_{0}+\unicode[STIX]{x1D6FE}S_{0}c_{0}}}-k\unicode[STIX]{x0394}{\displaystyle \frac{C_{0}}{S_{0}}}\right). & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

It is straightforward, though algebraically cumbersome, to solve for the coefficients $b_{0}$ , $b_{1}$ and $b_{\pm }$ . In particular, in the limit of small $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}$ , we arrive at the following approximate expressions for the coefficients $b_{0}$ and $b_{1}$ in (4.42):

(4.44) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}rcl@{}}b_{0}\,\, & \simeq \,\, & -{\displaystyle \frac{B_{0}^{2}}{\unicode[STIX]{x1D707}\unicode[STIX]{x0394}{\mathcal{C}}_{0}+{\mathcal{S}}_{0}}}[\unicode[STIX]{x0394}J_{0}+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}[2\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394}(c_{0}^{2}-s_{0}^{2})+2s_{0}c_{0}]]+O(\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}^{2}),\\ b_{1}\,\, & \simeq \,\, & -\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}{\displaystyle \frac{B_{0}^{2}}{\unicode[STIX]{x1D707}\unicode[STIX]{x0394}{\mathcal{S}}_{0}-{\mathcal{C}}_{0}}}[\unicode[STIX]{x1D6E5}(J_{0}-J_{1})-(\unicode[STIX]{x1D707}\unicode[STIX]{x0394}^{2}){\mathcal{S}}_{0}(\unicode[STIX]{x0394}J_{0}+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}[2\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394}(c_{0}^{2}-s_{0}^{2})+2s_{0}c_{0}])]\\ \,\, & \,\, & +\,\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D6E5}}\left[(2\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394})^{2}(2s_{0}c_{0})-2{\displaystyle \frac{\unicode[STIX]{x1D6E5}(k^{2}+\unicode[STIX]{x1D6FE}^{2})S_{0}c_{0}}{kC_{0}s_{0}+\unicode[STIX]{x1D6FE}S_{0}c_{0}}}(2\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394}(2s_{0}c_{0})+(c_{0}^{2}-s_{0}^{2}))\right]\\ \,\, & \,\, & +\,O(\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}^{2}),\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

in which ${\mathcal{C}}_{0}=\cos (\unicode[STIX]{x1D707}d)$ and ${\mathcal{S}}_{0}=\sin (\unicode[STIX]{x1D707}d)$ . Of course, if $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}=0$ , then $b_{1}=0$ and we recover the fact that the induced flow has odd structure in $z$ . However, the denominator in $b_{1}$ turns out to be small for long waves. Explicitly, we find $\unicode[STIX]{x1D707}\simeq \unicode[STIX]{x1D6FE}(1+(3/2)k^{2}/\unicode[STIX]{x1D6FE}^{2})$ , ${\mathcal{C}}_{0}\simeq c_{0}-(3/2)(k^{2}/\unicode[STIX]{x1D6FE}^{2})(\unicode[STIX]{x1D6FE}d)s_{0}$ and ${\mathcal{S}}_{0}\simeq s_{0}+(3/2)(k^{2}/\unicode[STIX]{x1D6FE}^{2})(\unicode[STIX]{x1D6FE}d)c_{0}$ . Together with the dispersion relation (4.39), which suggests $\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394}c_{0}-s_{0}+O(\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}^{2})=0$ , we find that the denominator varies as $k^{2}$ . Also in the limit $kH\rightarrow 0$ , $\unicode[STIX]{x1D6FD}\rightarrow 0$ as a consequence of $\unicode[STIX]{x1D712}\rightarrow 1$ . Hence for long waves in moderately asymmetric top-hat stratification we have

(4.45) $$\begin{eqnarray}b_{1}\simeq -{\displaystyle \frac{\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}}{(kH)^{2}}}{\displaystyle \frac{2\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394}H^{2}(J_{0}-J_{1})}{3d(\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394}c_{0}+s_{0})}}.\end{eqnarray}$$

Likewise, by continuity of $\overline{\unicode[STIX]{x1D6F9}}$ , this implies that $b_{\pm }$ is large if $kH\ll \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}^{1/2}$ . Hence, modulated wave packets containing long waves are expected to induce large Eulerian flows in the unstratified regions at the top and bottom of the domain.

This is illustrated in figure 3(d–f) in all of which $kH=0.3$ and $d=0.5D$ , and $\unicode[STIX]{x1D6FF}/H$ equals $0.001$ ( $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}=0.004$ ), $0.01$ ( $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}=0.04$ ) and $0.1$ ( $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}=0.4$ ). A slight breaking of symmetry is observed in the first case, but this becomes pronounced in the second case for which $(kH)^{2}=2.25\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}$ , in which case $u_{E}$ adopts a dominant mode-1 structure. In the final case with $(kH)^{2}\ll \unicode[STIX]{x1D700}_{\unicode[STIX]{x1D6FF}}$ the induced Eulerian flow is an order of magnitude larger than the Stokes drift and so dominates the Lagrangian flow.

Figure 5 shows the induced Eulerian flows at the top and bottom of the domain as a function of $kH$ in cases with $d=0.5D$ and with $\unicode[STIX]{x1D6FF}=0.01H$ (a) and $\unicode[STIX]{x1D6FF}=0.1H$ (b). Both flows increase in magnitude as a power law with decreasing $kH$ . The flow in the shallower top layer is moderately larger (for small $kH$ ).

Figure 5. As in figure 4(e), but for asymmetric top-hat stratification and showing on a log–log plot the magnitude of the Eulerian induced flows at both the top and bottom of the domain. The middle of the stratified region of half-depth $d=0.5D$ is shifted vertically from $z=0$ by (a) $\unicode[STIX]{x1D6FF}=0.01H$ and (b) $\unicode[STIX]{x1D6FF}=0.1H$ . Flow values are normalized by $|A_{0}|^{2}N_{0}/H$ , as denoted by the tilde. The dotted line indicates a power law dependence of $(kH)^{-2}$ .

Figure 6 illustrates the induced Eulerian and total Lagrangian transport in the top and middle layer as it depends upon $\unicode[STIX]{x1D6FF}/H$ and $kH$ . The transport in the lower layer is given by that in the top layer when $\unicode[STIX]{x1D6FF}\rightarrow -\unicode[STIX]{x1D6FF}$ . As well as the strong induced flows occurring near the resonance about $kH\sim 3$ , there is a clear near resonance occurring for wavepackets composed of long waves in which the flow is negative in the shallower layer and negative throughout the stratified layer. Such large transport disappears when the flows are of strictly equal depth.

Figure 6. For asymmetric top-hat stratification, transport resulting from Eulerian induced flows (a,c) and the total volume flux (b,d). The middle stratified layer extends over $-D/2+\unicode[STIX]{x1D6FF}<z<D/2+\unicode[STIX]{x1D6FF}$ . The volume fluxes in the unstratified top layer ( $D/2+\unicode[STIX]{x1D6FF}<z<D$ ) are shown in (a) and (b). The volume fluxes in the stratified layer are shown in (c) and (d).

4.4.1 Linear solutions: $O(\unicode[STIX]{x1D6FC})$

For an arbitrary stratification prescribed by $N^{2}$ , the vertical structure of internal modes can be found through a Galerkin analysis. Shifting vertical coordinates so that $z$ lies in the range $0\leqslant z\leqslant H$ , $N^{2}$ is decomposed into a Fourier cosine series and $\hat{\unicode[STIX]{x1D702}}$ is decomposed into a Fourier sine series: $\hat{\unicode[STIX]{x1D702}}=\sum _{j}A_{j}\sin (jm_{0}z)$ , with $m_{0}=\unicode[STIX]{x03C0}/H$ . Thus (3.4) is reduced to a matrix eigenvalue problem for the eigenvector of coefficients $A_{m}$ with corresponding eigenvalue being the squared frequency $\unicode[STIX]{x1D714}^{2}$ . The lowest mode corresponds to that with the highest frequency. This procedure is applied to exponential stratification given by

(4.46) $$\begin{eqnarray}N^{2}=N_{0}^{2}\exp (z/\unicode[STIX]{x1D70E}_{e}),\end{eqnarray}$$

for $-H\leqslant z\leqslant 0$ . The resulting predicted vertical structure of $\hat{\unicode[STIX]{x1D702}}_{0}$ is shown in figure 7(a) for the case with $\unicode[STIX]{x1D70E}_{e}=0.1H$ and $kH=0.3$ . The predicted group velocity is found by finding the frequencies of the lowest mode with wavenumbers $1.01k$ and $0.99k$ and so estimating $c_{g}\simeq [\unicode[STIX]{x1D714}(1.01k)-\unicode[STIX]{x1D714}(0.99k)]/(0.02k)$ .

4.4.2 Wave-induced Eulerian flow

Following a similar Galerkin procedure to solve (3.13), the resulting matrix equation is used to find the vertical structure of the Eulerian induced flow $\overline{\unicode[STIX]{x1D6F9}}(z)$ . From this, we construct the streamfunction of the induced flow $\overline{\unicode[STIX]{x1D713}}^{(2)}=N_{0}|A_{0}|^{2}\overline{\unicode[STIX]{x1D6F9}}$ and the corresponding Eulerian induced flow $\overline{u}^{(2)}=-\unicode[STIX]{x2202}_{z}\overline{\unicode[STIX]{x1D713}}^{(2)}$ , which is shown in figure 7(a) for the case with $\unicode[STIX]{x1D70E}_{e}=0.1H$ and $kH=0.3$ . In this figure $\overline{\unicode[STIX]{x1D6F9}}$ has been normalized by is maximum absolute value. It is evident that the induced flow has a very similar vertical structure to the linear vertical structure function, consistent with long-wave resonant excitation. In fact, the magnitude of $|\overline{\unicode[STIX]{x1D6F9}}|$ is very large. The corresponding profiles of the Stokes drift, the induced Eulerian flow and the Lagrangian flow are shown in figure 7(b), illustrating that the Eulerian flow is indeed dominant and drives the strong negative flow at the surface.

Figure 7. For exponential stratification, (a) profiles of the background stratification with $h=0.1H$ (thin solid line), vertical structure of mode-1 waves having $kH=0.3$ (dashed line) and the vertical structure of the Eulerian induced flow normalized by its maximum absolute value (thick blue solid line), and (b) corresponding profiles of the Stokes drift (dashed red line), induced Eulerian flow (dashed blue line) and Lagrangian flow (black line), all normalized by $N_{0}|A_{0}|^{2}/H$ , as denoted by the tilde.

As anticipated by the discussion of induced flows in asymmetric two-layer and top-hat stratification, the induced Eulerian flow increases as $kH$ decreases for exponential stratification, as shown in figure 8. The induced Eulerian flows at the top and bottom vary as $(kH)^{-2}$ , for small $kH$ , as expected. For $kH=0.1$ the speeds of the induced Eulerian flow at the surface and bottom are each three orders of magnitude larger than the corresponding Stokes drift. The values are so large that the speed at the surface can be comparable with the horizontal flow due to the waves themselves if $a_{0}/H\sim (kH)^{-2}\ll 1$ . Evidently, for large $kH$ we also observe resonances.

Figure 8. As in figure 5, but showing flows induced by waves in exponential stratification with e-folding depth $\unicode[STIX]{x1D70E}_{e}=0.1H$ .