2. Basic equations

The flow is treated as incompressible, inviscid and two-dimensional, and the background rotation is neglected. The flow is then governed by

(2.1)\begin{gather} \left( \rho_0 + \rho \right) \frac{{\rm D} u}{{\rm D} t} ={-} \frac{\partial p}{\partial x}, \end{gather}

(2.2)\begin{gather}\left( \rho_0 + \rho \right) \frac{{\rm D} w}{{\rm D} t} ={-} \frac{\partial p}{\partial z} - \rho g, \end{gather}

(2.3)\begin{gather}{\frac{\partial u}{\partial x}} + {\frac{\partial w}{\partial z}} = 0, \end{gather}

(2.4)\begin{gather}{\frac{{\rm D} \rho}{{\rm D} t}} - {N^2}\,\frac{\rho_0}{g}\,w = 0, \end{gather}

where

(2.5)\begin{equation} {\frac{{\rm D}}{{\rm D} t}} = {\frac{\partial}{\partial t}} + u\,{\frac{\partial}{\partial x}} + w\, {\frac{\partial}{\partial z}}, \end{equation}

the velocity is $(u, w)$, the dynamic pressure is $p$, $\rho _0 (z)$ is the base-state density, $\rho$ is the disturbance density, and $N$ is defined by

(2.6)\begin{equation} N^2 ={-} \frac{g {\rho_0}_z}{\rho_0}. \end{equation}

It will be useful to employ the field $\zeta (x,z,t)$, which is vertical displacement of material lines. The kinematic relationship between velocity and $\zeta$ is

(2.7)\begin{equation} \frac{{\rm D} \zeta}{{\rm D} t} = w. \end{equation}

The interfacial conditions between any two layers in a stratified fluid are developed in Grimshaw & McHugh (Reference Grimshaw and McHugh2013). In general, the interface may have a jump in density or buoyancy frequency. The kinematic condition on the interface is continuity of normal velocity. This condition is treated in the usual manner, with

(2.8)\begin{equation} \eta_t + u \eta_x = w \end{equation}

on $z = \eta$, using velocity components for each layer. Since the interface is also a material line, the kinematic interfacial condition is merely

(2.9)\begin{equation} \zeta = \eta \end{equation}

on $z = \eta$. Using a Taylor series in the usual way, (2.9) becomes

(2.10)\begin{equation} \zeta + \zeta_z \eta + \tfrac{1}{2} \zeta_{zz} \eta^2 + \cdots = \eta \end{equation}

on $z=0$. This expression may be used recursively (Grimshaw & McHugh Reference Grimshaw and McHugh2013) to get

(2.11)\begin{equation} \eta = \zeta + \zeta \zeta_z + \zeta \zeta_z^2+ \tfrac{1}{2} \zeta^2 \zeta_{zz} + \cdots \end{equation}

on $z=0$. The kinematic condition may now be written as

(2.12)\begin{equation} [\zeta + \zeta \zeta_z + \zeta \zeta_z^2 + \tfrac{1}{2} \zeta^2 \zeta_{zz} + \cdots]_-^+ = 0 \end{equation}

on $z = 0$.

The dynamic condition at the interface is continuity of total pressure:

(2.13)\begin{equation} [ p_0 + p ]_-^+ = 0 \end{equation}

on $z = \eta$, where $p_0 (z)$ is the base-state hydrostatic pressure. Using a Taylor series in (2.13) produces

(2.14) \begin{equation} [( p + p_0) + ( p_z + {p_0}_z ) \eta + \tfrac{1}{2} ( p_{zz} + {p_0}_{zz} ) \eta^2 + \tfrac{1}{6} ( p_{zzz} + {p_0}_{zzz} ) + \cdots]_-^+ = 0 \end{equation}

on $z = 0$. The base-state pressure is continuous at all positions, thus

(2.15)\begin{equation} p_0|_-^+ = 0. \end{equation}

Equation (2.14) is used recursively to evaluate derivatives of disturbance pressure. Furthermore, the variable $\eta$ is eliminated from the dynamic condition using (2.11). After some manipulation, the dynamic interfacial condition becomes

(2.16)\begin{equation} \left[p - \rho_0 g \zeta - \frac{1}{2}\,\rho_0 N^2 \zeta^2 + \frac{1}{6}\,\rho_0 \left(\frac{{\rm d}N^2}{{\rm d}z} - \frac{N^4}{g} \right) \zeta^3\right]_-^+ = 0, \end{equation}

on $z = 0$, where higher-order terms have been dropped. Equation (2.16) is valid for any interface, including discontinuous $\rho _0$ and/or $N$.

With continuous density $\rho _0$ and discontinuous $N$, but $N$ constant in each layer, (2.16) reduces to

(2.17)\begin{equation} \left[p - \frac{1}{2}\,\rho_0 N^2 \zeta^2- \frac{1}{6}\,\rho_0\,\frac{N^4}{g}\,\zeta^3\right]_-^+ = 0 \end{equation}

on $z=0$.

3. The amplitude equations

3.1. Away from the interface

A packet of monochromatic waves is assumed to impinge on the interface from below. The wave amplitude $\alpha$ is assumed small, and all dependent variables are expanded in a power series in $\alpha$:

(3.1)\begin{equation} \zeta = \alpha \zeta^{(1)} + \alpha^2 \zeta^{(2)} + \alpha^2 \bar{\zeta} + \cdots, \end{equation}

and similarly for $u,w,\rho,p$, where the superscript indicates the order of the approximation. Here, $\bar {\zeta }$ is the wave-induced mean. Previous theory (Shrira Reference Shrira1981) has shown that the wave-induced mean is non-zero at second order. This mean part is separated from the oscillatory second-order solution merely for convenience.

An upwardly propagating internal wave has the leading-order solution

(3.2)\begin{equation} \zeta^{(1)} = I E \exp({-{\rm i} n_1 z}) + I^* E^* \exp({{\rm i} n_1 z}), \end{equation}

with

(3.3)\begin{equation} E = \exp({{\rm i}(kx - \sigma t)}), \end{equation}

where $I$ is the incident wave amplitude, $k$ is the horizontal wavenumber, $n_1$ is the vertical wavenumber, and $\sigma$ is the wave frequency, all assumed positive and real. However, this solution alone cannot meet the interfacial conditions, which require a reflected and transmitted wave at this order. Following Grimshaw & McHugh (Reference Grimshaw and McHugh2013) and McHugh (Reference McHugh2015), the leading-order solution including reflected and transmitted waves is

(3.4)\begin{gather} \zeta^{(1)} = I E \exp({-{\rm i} n_1 z}) + I^* E^* \exp({{\rm i} n_1 z}) + R E \exp({{\rm i}n_1 z}) + R^* E^* \exp({-{\rm i} n_1 z}),\quad z<0, \end{gather}

(3.5)\begin{gather}\zeta^{(1)} = T E \exp({-{\rm i} n_2 z}) + T^* E^* \exp({{\rm i} n_2 z}),\quad z>0, \end{gather}

where $n_1,n_2$ are the vertical wavenumbers in the lower and upper layers, respectively, and $R,T$ are the reflected and transmitted wave amplitudes, respectively. All terms proportional to $E$ will be called the primary harmonic, while the secondary harmonic has $E^2$, and so on. Thus the above expression for $\zeta ^{(1)}$ contains only the primary harmonic. The secondary harmonic will appear in $\zeta ^{(2)}$ as a result of nonlinear effects in the interfacial conditions.

The wave amplitude is chosen to be vertically modulated but horizontally uniform, as in Sutherland (Reference Sutherland2006), Grimshaw & McHugh (Reference Grimshaw and McHugh2013) and McHugh (Reference McHugh2015). The parameter $\epsilon$ is defined to be the inverse of the packet length, normalized by the horizontal wavenumber $k$. This makes $\epsilon$ proportional to the ratio of horizontal wavelength to vertical packet length. A solution is obtained by assuming that the packet is long, making $\epsilon$ small. Slow variables are introduced to achieve this solution:

(3.6)\begin{equation} \left.\begin{gathered} t_j = \epsilon^j t,\\ z_j = \epsilon^j z, \end{gathered}\right\} \end{equation}

with $j = 1,2$. The amplitude functions $I,R,T$ all depend on these slow variables,

The waves create a wave-induced mean flow, as discussed in Grimshaw & McHugh (Reference Grimshaw and McHugh2013). The mean is defined here as a horizontal average ($x$-average), denoted with an overbar. The components of the mean solution that are needed for the interfacial conditions are $\bar {\zeta }$ and $\bar {u}$, which may be determined to second order using

(3.7)\begin{equation} \left.\begin{gathered} \alpha^2 \bar{\zeta} ={-} \frac{1}{2}\,\frac{\partial}{\partial z}\langle \zeta^2 \rangle,\\ \alpha^2 \bar{u} = \frac{k}{\sigma}\,N^2 \langle \zeta^2\rangle, \end{gathered}\right\} \end{equation}

where the angle brackets mean average. Using (3.4) and (3.5) in these equations gives (Grimshaw & McHugh Reference Grimshaw and McHugh2013)

(3.8)\begin{gather} \bar{\zeta} ={-} {\rm i} 2 n_1 ( I^* R \exp({{\rm i} 2 n_1 z}) - I R^* \exp({-{\rm i} 2 n_1 z})),\quad z < 0, \end{gather}

(3.9)\begin{gather}\bar{u} = 2 N_1^2\,\frac{k}{\sigma} [I I^* + R R^* + I^* R \exp({{\rm i} 2 n_1 z}) + I R^* \exp({-{\rm i} 2 n_1 z})],\quad z < 0, \end{gather}

(3.10)\begin{gather}\bar{u} = 2 N_2^2\,\frac{k}{\sigma}\,T T^*,\quad z > 0, \end{gather}

and $\bar {\zeta } = 0$ for $z>0$.

The governing equations (2.1)–(2.7) may be written as

(3.11)\begin{gather} \rho_0 \left( u_t + Q_1 \right) ={-} \frac{\partial p}{\partial x}, \end{gather}

(3.12)\begin{gather}\rho_0 \left( w_t + Q_3 \right) ={-} \frac{\partial p}{\partial z} - \rho g, \end{gather}

(3.13)\begin{gather}u_x + w_z = 0, \end{gather}

(3.14)\begin{gather}\rho_t + Q_4 - {N^2}\,\frac{\rho_0}{g}\,w = 0, \end{gather}

(3.15)\begin{gather}\zeta_t + Q_5 = w, \end{gather}

where the $Q_j$ are the sums of the nonlinear terms in each equation. Eliminate all variables in the linear terms except $\zeta$ to achieve

(3.16)\begin{align} &\frac{\partial^3}{\partial t^3}\left( \zeta_{xx} + ( \rho_0 \zeta_z )_z \right)+ \rho_0 N^2 \zeta_{xxt} \nonumber\\ &\quad = (\rho_0 {Q_1}_{xt} )_z- \rho_0 {Q_3}_{xxt} + {Q_4}_{xx} \nonumber\\ &\qquad - \frac{\partial^2}{\partial t^2} (\rho_0 {Q_5}_{xx} + ( \rho_0 {Q_5}_z )_z )- \rho_0 N^2 {Q_5}_{xx}. \end{align}

Now assume Boussinesq flow. It may be shown that $Q_4 \approx \rho _0 N^2 Q_5$, allowing some cancellation, which in turn permits an integration with respect to time:

(3.17)\begin{equation} \frac{\partial^2}{\partial t^2} \nabla^2 \zeta + N^2 \zeta_{xx} \approx {Q_1}_{xz} - {Q_3}_{xx} - \nabla^2 {Q_5}_t. \end{equation}

Away from the interface, the nonlinear quantities may be approximated using

(3.18)\begin{equation} \left.\begin{gathered} Q_1 \approx \alpha^3 (\bar{u} u_x^{(1)} + \bar{u}_z w^{(1)} ),\\ Q_3 \approx \alpha^3 ( \bar{u} w_x^{(1)} ),\\ Q_5 \approx \alpha^3 (\bar{u} \zeta_x^{(1)} + \bar{\zeta}_z w^{(1)} ). \end{gathered}\right\} \end{equation}

Only terms that contribute to the primary harmonic are included here, and higher-order contributions, even to the primary harmonic, are neglected.

Further details of the derivation of the three amplitude equations are provided in McHugh (Reference McHugh2015). The resulting equations are

(3.19) \begin{align} & I_t + c_g I_z - {\rm i}\,\tfrac{1}{2}c_g^\prime I_{zz}\nonumber\\ &\quad + {\rm i} \alpha^2 2 \sigma [( k^2 + n_1^2 ) ( | I |^2 + | R |^2 ) + ( k^2 - \tfrac{3}{2}n_1^2 ) | R |^2 ] I = 0,\quad z<0,\end{align}

(3.20) \begin{align} & R_t - c_g R_z - {\rm i}\,\tfrac{1}{2}c_g^\prime R_{zz} \nonumber\\ &\quad + {\rm i} \alpha^2 2 \sigma [( k^2 + n_1^2 ) ( | I |^2 + | R |^2 ) + ( k^2 - \tfrac{3}{2} n_1^2 ) | I |^2 ] R = 0,\quad z<0, \end{align}

(3.21)\begin{equation} T_t + c_g T_z - {\rm i}\,\tfrac{1}{2} c_g^\prime T_{zz}+ {\rm i} \alpha^2 2 \sigma ( k^2 + n_2^2 ) | T |^2 T = 0,\quad z>0. \end{equation}

The definitions of $I,R,T$ in McHugh (Reference McHugh2015) are slightly different than here. McHugh (Reference McHugh2015) defined $I,R,T$ in terms of vertical velocity rather than vertical displacement, as here and in Grimshaw & McHugh (Reference Grimshaw and McHugh2013). Furthermore, McHugh (Reference McHugh2015) neglected the contribution made by the mean buoyancy $\bar {b}$, reflected in the mean displacement $\bar {\zeta }$. This effect is included here, and results in the $3/2$ that appears in (3.19) and (3.20). All other aspects of the evolution equations match McHugh (Reference McHugh2015).

3.2. At the interface

3.2.1. The second harmonic

The interfacial conditions (2.12) and (2.17) must be satisfied at $z = 0$, and they must be expressed in terms of the wave amplitudes in the two layers. The interface generates second harmonics, and the interaction between the primary and secondary harmonics is an important aspect of the nonlinear interfacial condition. This interaction does not contribute to the amplitude equations away from the interface, since the vertical wavenumber of the second harmonic is not commensurate with the primary harmonic. However, at the interface, where the second harmonic is created, the two harmonics are aligned and their interaction is important.

The second harmonic solution $\zeta ^{(2)}$ is

(3.22)\begin{gather} \zeta^{(2)} = B_1 E^2 \exp({{\rm i} n_{12} z}) + B_1^* {E^*}^2 \exp({-{\rm i} n_{12} z}),\quad z<0, \end{gather}

(3.23)\begin{gather}\zeta^{(2)} = B_2 E^2 \exp({-{\rm i} n_{22} z}) + B_2^* {E^*}^2 \exp({{\rm i} n_{22} z}),\quad z>0, \end{gather}

where

(3.24)\begin{gather} n_{12}^2 = n_1^2 - 3 k^2, \end{gather}

(3.25)\begin{gather}n_{22}^2 = n_2^2 - 3 k^2, \end{gather}

as in McHugh (Reference McHugh2009). The values of $B_1$ and $B_2$ are determined by the quadratic terms in (2.12) and (2.17). Substitute (3.4), (3.5), (3.22) and (3.23) into (2.12), and keep only second harmonic terms, resulting in

(3.26)\begin{equation} [ B_2 - B_1 - {\rm i} n_2 T^2 + {\rm i} n_1 ( I^2 - R^2 )] E^2 + cc = 0 \end{equation}

on $z = 0$.

The dynamic condition (2.17) contains the pressure. A useful relationship between the pressure and wave amplitude can be obtained from the horizontal component of the governing equations (3.11):

(3.27)\begin{equation} p_{xx} = \rho_0 ( \zeta_{ttz}- {Q_1}_x + {Q_5}_{tz} ), \end{equation}

which includes nonlinear contributions. For this second harmonic,

(3.28)\begin{equation} \left.\begin{gathered} Q_1 \approx \alpha^2 [ u^{(1)} u_x^{(1)} + w^{(1)} u_z^{(1)}],\\ Q_5 \approx \alpha^2 [ u^{(1)} \zeta_x^{(1)} + w^{(1)} \zeta_z^{(1)}]. \end{gathered}\right\} \end{equation}

Direct substitution of linear quantities in these expressions shows that $Q_5$ is zero, but $Q_1$ is non-zero for the lower layer. The resulting expressions for the second harmonic of $p$ are

(3.29)\begin{gather} p^{(2)} = \rho_0\,\frac{\sigma^2}{k^2} [{\rm i} n_{12} B_1 \exp({-{\rm i} n_{12} z}) - 2 n_1^2 I ] E^2 + cc,\quad z < 0, \end{gather}

(3.30)\begin{gather}p^{(2)} = {\rm i} \rho_0\,\frac{\sigma^2}{k^2}\,n_{22} B_2 E^2 \exp({-{\rm i} n_{22} z}) + cc,\quad z > 0. \end{gather}

Using (3.29) and (3.30) along with the above expressions for $\zeta$ in (2.17) gives

(3.31)\begin{align} & [n_{22} B_2 - n_{12} B_1 - {\rm i} n_1^2 2 I R \vphantom{\tfrac{1}{2}} -{\rm i}\,\tfrac{1}{2} ( k^2 + n_1^2 ) ( I + R )^2 + {\rm i}\,\tfrac{1}{2} ( k^2 + n_2^2 ) T^2 ] E^2 + cc = 0 \end{align}

on $z = 0$. Finally, (3.26) and (3.31) determine $B_1$ and $B_2$:

(3.32)\begin{align} B_1 &= {\rm i}\,\frac{1}{( n_{22} - n_{12} )}\left[n_1^2 2 I R + \frac{1}{2}\,( k^2 + n_1^2 ) ( I + R )^2 + n_1 n_{22} ( I^2 - R^2 ) \right.\nonumber\\ &\quad \left.{}- \frac{1}{2}\,( k^2 + n_2^2 ) T^2- n_2 n_{22} T^2 \right], \end{align}

(3.33)\begin{align} B_2 &= {\rm i}\,\frac{1}{( n_{22} - n_{12} )} \left[ n_1^2 2 I R + \frac{1}{2}\,( k^2 + n_1^2 ) ( I + R )^2 + n_1 n_{12} ( I^2 - R^2 ) \right.\nonumber\\ &\quad \left.{}-\frac{1}{2}\,( k^2 + n_2^2 ) T^2 - n_2 n_{12} T^2 \right]. \end{align}

3.2.2. The kinematic condition

Now consider the interfacial conditions (2.12) and (2.17) for the primary harmonic. Using (3.1) in (2.12), and keeping only those terms that will contribute to the primary harmonic up to $O(\alpha ^3)$, results in

(3.34) \begin{equation} \alpha \zeta^{(1)} + \alpha^3\,\frac{\partial}{\partial z} [\bar{\zeta} \zeta^{(1)} + \zeta^{(1)} \zeta^{(2)}] + \alpha^3 \frac{1}{2}\,\frac{\partial}{\partial z} [(\zeta^{(1)} )^2 \zeta_z^{(1)} ]|_-^+ = 0 \end{equation}

on $z=0$. Using (3.4), (3.5), (3.22) and (3.23) gives

(3.35)\begin{equation} T - \left( I + R \right) + \alpha^2 Q_k = 0 \end{equation}

on $z=0$, where $Q_k$ is the sum of the nonlinear effects, and

(3.36)\begin{equation} Q_k = M_k + S_k + C_k. \end{equation}

Here, $M_k$ is the nonlinear interaction between the mean and first harmonic, $S_k$ is the interaction between first and second harmonics, and $C_k$ is the cubic nonlinearity of the first harmonic:

(3.37)\begin{gather} M_k = 2 n_1^2 [ I R ( I^* + R^* ) + 3 ( I^2 R^* + I^* R^2 )], \end{gather}

(3.38)\begin{gather}S_k = {\rm i} [ ( n_2 - n_{22} ) T^* ] B_2 - {\rm i} [ ( n_1 +n_{12} ) I^* -( n_1 - n_{12} ) R^* ] B_1, \end{gather}

(3.39)\begin{gather}C_k ={-} \tfrac{1}{2} n_2^2 [ T^2 T^* ] + \tfrac{1}{2} n_1^2 [ I^2 I^* + R^2 R^* + 2 I R ( I^* + R^* ) + 9 ( I^2 R^* + I^* R^2 ) ]. \end{gather}

3.2.3. The dynamic condition

To process the dynamic condition (2.17) in a similar manner, first use (3.27) again to evaluate the pressure, where the nonlinear terms that contribute to the primary harmonic are

(3.40)\begin{gather} Q_1 \approx \alpha^3 ( \bar{u} u_x^{(1)} + \bar{u}_z w^{(1)} + u^{(1)} u_x^{(2)} + u^{(2)} u_x^{(1)} + w^{(1)} u_z^{(2)} + w^{(2)} u_z^{(1)} ), \end{gather}

(3.41)\begin{gather}Q_5 \approx \alpha^3 ( \bar{u} \zeta_x^{(1)} + \bar{\zeta}_z w^{(1)} + u^{(1)} \zeta_x^{(2)} + u^{(2)} \zeta_x^{(1)} + w^{(1)} \zeta_z^{(2)} + w^{(2)} \zeta_z^{(1)} ). \end{gather}

Evaluating pressure in this way, and using (3.1) in the dynamic condition, results in

(3.42) \begin{align} & \left(\frac{1}{k^2} [ \alpha \zeta_{ttz}^{(1)} - \alpha^3 ( {Q_1}_x^{(1)} - {Q_5}_{tz}^{(1)} ) ] \right.\nonumber\\ &\quad \left.+\,\alpha^3 \left( N^2 \bar{\zeta} \zeta^{(1)} + N^2 \zeta^{(1)} \zeta^{(2)} + \frac{1}{6}\, \frac{N^4}{g}\,{\zeta^{(1)}}^3 \right) \right)_-^+ = 0. \end{align}

Slow variables must be added to the linear term $\zeta _{ttz}$ in (2.17):

(3.43) \begin{align} & \left(\frac{1}{k^2} [ \zeta_{ttz}^{(1)} + \epsilon ( \zeta_{ttz_1}^{(1)} + 2 \zeta_{tzt_1}^{(1)} ) \right.\nonumber\\ &\quad \left.{}+\epsilon^2 ( \zeta_{z t_1 t_1}^{(1)} + 2 \zeta_{t t_1 z_1}^{(1)} + \zeta_{ttz_2}^{(1)} + 2 \zeta_{tzt_2}^{(1)} ) - \alpha^2 ( {Q_1}_x^{(1)} - {Q_5}_{tz}^{(1)} ) ]\right.\nonumber\\ &\quad \left.{}+\alpha^2 \left( N^2 \bar{\zeta} \zeta^{(1)} + N^2 \zeta^{(1)} \zeta^{(2)} + \frac{1}{6}\, \frac{N^4}{g}\,{\zeta^{(1)}}^3 \right)\right)_-^+ = 0. \end{align}

The derivatives with respect to $z_1$ may be converted into temporal derivatives using

(3.44) \begin{equation} I_{t_1} + c_g I_{z_1} = 0, \end{equation}

and similar operations on $R,T$. The derivatives with respect to $z_2$ may also be converted into temporal derivatives with the amplitude equations. The nonlinear parts of the amplitude equations do not add further nonlinear terms here, as these terms wind up at higher order. The nonlinear terms are treated by inserting (3.4), (3.5), (3.22), (3.23), (3.8), (3.9) and (3.10), retaining only contributions to the primary harmonic up to third order. The final dynamic condition is

(3.45) \begin{align} & -{\rm i} \left(n_1 I - n_2 T - n_1 R \right)\nonumber\\ &\quad + \frac{1}{\sigma} \left[ \frac{n_1^2 - k^2}{n_1} \right] ( \epsilon I_{t_1} + \epsilon^2 I_{t_2} )+ {\rm i}\,\frac{1}{2 \sigma^2} \left[\frac{k^4 - 5 k^2 n_1^2 - 4 n_1^4}{n_1^3} \right] ( \epsilon^2 I_{t_1 t_1} ) \nonumber\\ &\quad - \frac{1}{\sigma} \left[ \frac{n_2^2 - k^2}{n_2} \right] ( \epsilon T_{t_1} + \epsilon^2 T_{t_2} ) - {\rm i}\,\frac{1}{2 \sigma^2} \left[\frac{k^4 - 5 k^2 n_2^2 - 4 n_2^4}{n_2^3} \right] (\epsilon^2 T_{t_1 t_1})\nonumber\\ &\quad - \frac{1}{\sigma} \left[\frac{n_1^2 - k^2}{n_1} \right] (\epsilon R_{t_1} + \epsilon^2 R_{t_2}) - {\rm i}\,\frac{1}{2 \sigma^2} \left[\frac{k^4 - 5 k^2 n_1^2 - 4 n_1^4}{n_1^3} \right](\epsilon^2 R_{t_1 t_1})\nonumber\\ &\quad + \alpha^2 Q_d \end{align}

on $z=0$, where

(3.46)\begin{equation} Q_d = M_d + S_d + C_d, \end{equation}

(3.47) \begin{align} M_d = &- {\rm i} 2 (2 n_2 ( k^2 + n_2^2 ) T^2 T^* \nonumber\\ &- n_1 ( k^2 + n_1^2 ) [ ( I + R )^2 ( I^* - R^* ) + 2 ( I^2- R^2 )( I^* + R^* ) + 3 I R ( I^* - R^* ) ] \nonumber\\ & - 2 n_1^3 [ I R ( I^* - R^* ) + 3 ( I^* R^2 - I^2 R^* ) ]), \end{align}

(3.48) \begin{align} S_d &= ( ( k^2 + n_2 n_{22} - n_{22}^2 + 3 n_2^2 ) T^* B_2 \nonumber\\ &\quad - [ ( k^2 - n_{12}^2 + 3 n_1^2 ) ( I^* + R^* ) - n_1 n_{12} ( I^* -R^* ) ] B_1), \end{align}

(3.49)\begin{equation} C_d = \frac{1}{2}\,k\,\frac{N_1^2}{g k} \left[ ( k^2 + n_2^2 )\,\frac{N_2^2}{N_1^2}\,T^2 T^* - ( k^2 + n_1^2 ) ( I + R )^2 ( I^* + R^* ) \right]. \end{equation}

Once again, $M_d$ accounts for the interaction between the mean quantities and the primary harmonic, $S_d$ accounts for the interaction between the second harmonic and the primary harmonic, and $C_d$ accounts for the cubic effects of the primary harmonic.

The dimensionless quantity ${N_1^2}/{gk}$, using the definition in (2.6), is

(3.50)\begin{equation} \frac{N_1^2}{g k} ={-} \frac{1}{k}\,\frac{{\rho_0}_z}{\rho_0}, \end{equation}

which is the inverse of the density scale height, made dimensionless with the wavenumber $k$. This quantity is a free parameter in the analysis, but has the generally small value of approximately $1/100$ in the troposphere ($N_1 \approx 0.01\,\textrm {s}^{-1}$, $g \approx 10\,\textrm {m}\,\textrm {s}^{-2}$, $k \approx 1/1000\,\textrm {m}^{-1}$). This suggests that the effect of the cubic term $C_d$ in the dynamic interfacial condition is rather less important.

3.2.4. Final conditions

If slow variables are reverted to original variables in (3.45), then the operator for the linear terms in (3.45) may be defined by

(3.51)\begin{equation} L_m ={-} {\rm i} n_m + a_m\,\frac{\partial}{\partial t} + {\rm i} b_m\,\frac{\partial^2}{\partial t^2}, \end{equation}

where $m = 1,2$, and

(3.52)\begin{equation} \left.\begin{gathered} a_m = \frac{2}{\sigma} \left[ \frac{n_m^2 - k^2}{n_m} \right],\\ b_m = \frac{1}{2 \sigma^2}\left[\frac{k^4 - 5 k^2 n_m^2 - 4 n_m^4}{n_m^3} \right]. \end{gathered}\right\} \end{equation}

The dynamic condition becomes

(3.53)\begin{equation} L_1 I - L_2 T - L_1 R + \alpha^2 Q_d = 0. \end{equation}

This may be combined with the kinematic condition (3.35) to obtain final expressions governing transmitted and reflected waves:

(3.54)\begin{gather} \left( L_1 + L_2 \right) R = \left( L_1 - L_2 \right) I + \alpha^2 Q_R, \end{gather}

(3.55)\begin{gather}\left( L_1 + L_2 \right) T = 2 L_1 I + \alpha^2 Q_T \end{gather}

on $z=0$, where

(3.56)\begin{equation} \left.\begin{gathered} Q_R = Q_d - {\rm i} n_2 Q_k,\\ Q_T = Q_d + {\rm i} n_1 Q_k. \end{gathered}\right\} \end{equation}

4. Linear cases

The interfacial conditions (3.54) and (3.55) are two coupled equations that may be solved for the reflected and transmitted wave amplitudes, $R$ and $T$, at the interface. These equations are coupled since the nonlinear terms contain both $R$ and $T$. Before the calculation of $R$ and $T$ can proceed, the value of the incident wave amplitude $I$ at the interface is required. This value is set by the ascending wave packet.

McHugh (Reference McHugh2015) created a wave packet at the bottom of the computational domain and determined its evolution as the packet approached the interface and interacted with it. The initial packet shape was the ‘raised cosine’. Previous authors (Sutherland Reference Sutherland2001, Reference Sutherland2006; Tabaei & Akylas Reference Tabaei and Akylas2007) chose an initial packet shape with a Gaussian profile. For both profiles, the packet shape evolves at it ascends due to both nonlinearity and dispersion, arriving at the interface with a more complex shape, and continuing to evolve as the packet interacts with the interface. This complexity makes interpretation of the reflection and transmission difficult, and it is instructive to consider $I$ at the interface to be prescribed. Both the ‘raised cosine’ and the Gaussian profile are used here for this purpose.

It is also instructive to consider linear versions of the interfacial conditions. Previous results by McHugh (Reference McHugh2009, Reference McHugh2015), Mathur & Peacock (Reference Mathur and Peacock2009) and Grimshaw & McHugh (Reference Grimshaw and McHugh2013) used linear conditions accurate to $O(\alpha )$ for reflection and transmission of the primary harmonic:

(4.1)\begin{gather} R_0 = \frac{n_1 -n_2}{n_1 + n_2}\,I, \end{gather}

(4.2)\begin{gather}T_0 = \frac{2 n_1}{n_1 + n_2}\,I, \end{gather}

on $z = 0$. These algebraic expressions are identical to Snell's law in optics.

Another set of linear interfacial conditions is obtained from (3.54) and (3.55) by neglecting products of the wave amplitudes (setting $Q_R = Q_T = 0$) but retaining the linear terms containing derivatives of $R$ and $T$. These higher-order linear conditions include second- and third-order effects due to the non-constant wave amplitude.

For this higher-order linear case, the homogeneous part of the solution may be treated separately, which is governed by (3.54) with forcing set to zero:

(4.3)\begin{equation} \left( L_1 + L_2 \right) R = 0 \end{equation}

on $z=0$, which is

(4.4)\begin{equation} {\rm i} \left( b_1 + b_2 \right) R_{tt} + \left( a_1 + a_2 \right) R_t - {\rm i} \left( n_1 + n_2 \right) R = 0. \end{equation}

The homogeneous solutions for $T$ obey this same equation. The exponential function $\textrm {e}^{\lambda t}$ gives

(4.5)\begin{equation} \lambda = {\rm i}\,\frac{1}{2} \left[ \frac{a_1 + a_2}{b_1 + b_2} \pm \frac{\sqrt{(a_1 + a_2)^2 - 4 (b_1 + b_2)(n_1 + n_2)}} {b_1 + b_2} \right]. \end{equation}

If the quantity under the radical is positive, then the two values of $\lambda$ are purely imaginary, implying an oscillatory solution. However, if the argument of the radical is negative, then one of the $\lambda$ values is complex with positive real part, implying exponential growth and instability. Thus the sign of

(4.6)\begin{equation} (a_1 + a_2)^2 - 4 (b_1 + b_2)(n_1 + n_2) \end{equation}

dictates stability. Numerical evaluation shows that this quantity is negative for small $n_1/k$, and positive for larger values, with the critical value depending on $N_2/N_1$. With $N_2/N_1 = 2$,

(4.7)\begin{equation} \left( \frac{n_1}{k} \right)_{cr} \approx 0.335. \end{equation}

Thus an instability exists in this case with $n_1/k < 0.335$. Numerical solutions of the interfacial conditions, in both linear and nonlinear cases, show solutions with exponentially growing amplitude with $n_1/k$ less than the critical value, in agreement with this linear theory. For values of $n_1/k$ larger than the critical value, the homogeneous solutions are purely oscillatory, and determined by initial conditions. If the initial value and first derivative of $R$ and $T$ are zero, then the homogeneous solution is also zero.

The behaviour of $( n_1/k )_{cr}$ with $N_2/N_1$ is shown in figure 1. The value of $( n_1/k )_{cr}$ increases dramatically as $N_2/N_1$ becomes small. Thus as the upper layer becomes only weakly stratified, $N_2/N_1$ is near zero, and the interval of $n_1/k$ where this instability occurs includes all realistic values of $n_1/k$.

Figure 1. Critical value of $n_1/k$ for an interval of $N_2/N_1$. For values of $n_1/k$ less than the critical value, the higher-order linear terms result in instability.

This instability may explain the results of McHugh et al. (Reference McHugh, Dors, Jumper, Roadcap, Murphy and Hahn2008), who launched a sequence of weather balloons over Hawaii, measuring atmospheric velocity, temperature and other quantities. The measurements showed excessively large values of vertical velocity only at the elevation of the tropopause. The above instability showing exponential growth at the interface is likely responsible.

If the packet shape at the interface is chosen to be the ‘raised cosine’

(4.8)\begin{equation} I = \tfrac{1}{2} \left[ 1 - \cos{\omega t} \right] \end{equation}

on $z = 0$, where

(4.9)\begin{equation} \omega = 2 {\rm \pi}\epsilon c_g, \end{equation}

then an exact solution for the forced solution is easily obtained for $R$ and $T$:

(4.10)\begin{align} R &= \frac{1}{2} \left[ \frac{n_1 - n_2}{n_1 + n_2} \right] - \frac{1}{4} \left[\frac{(n_1 - n_2) - \omega (a_1 - a_2) + \omega^2 (b_1 - b_2)}{(n_1 + n_2) - \omega (a_1 + a_2) + \omega^2 (b_1 + b_2)} \right] \exp({\rm i}\omega t)\nonumber\\ &\quad - \frac{1}{4} \left[\frac{(n_1 - n_2) + \omega (a_1 - a_2) + \omega^2 (b_1 - b_2)}{(n_1 + n_2) + \omega (a_1 + a_2) + \omega^2 (b_1 + b_2)} \right] \exp({- {\rm i} \omega t}), \end{align}

(4.11)\begin{align} T &= \frac{1}{2} \left[ \frac{2 n_1}{n_1 + n_2} \right] - \frac{1}{4} \left[\frac{2 (n_1 - \omega a_1+ \omega^2 b_1)} {(n_1 + n_2) - \omega (a_1 + a_2)+ \omega^2 (b_1 + b_2)} \right] \exp({\rm i}\omega t)\nonumber\\ &\quad - \frac{1}{4} \left[\frac{2 (n_1 + \omega a_1+ \omega^2 b_1)} {(n_1 + n_2) + \omega (a_1 + a_2)+ \omega^2 (b_1 + b_2)} \right] \exp({- {\rm i} \omega t}). \end{align}

Since $\omega = O(\epsilon )$, these expressions are relatively minor corrections to $R_0$ and $T_0$ given by (4.1) and (4.2).

Consider numerical results using a prescribed incident wave profile. The ‘raised cosine’ profile discussed above has a discontinuous second derivative at the instant the incident wave packet begins interacting with the interface, and also at the instant the packet ends. For this reason, the numerical results will employ a Gaussian profile to prescribe $I$ at the interface:

(4.12)\begin{equation} I = \exp({- c_g^2 \epsilon^2 t^2}) \end{equation}

at $z = 0$. The numerical techniques used here are discussed in Appendix A.

Linear simulations with higher-order linear terms included and $I$ prescribed by (4.12) show that if $\epsilon \le 0.1$ approximately, then the numerical values of $R$ and $T$ are indistinguishable from $R_0$ and $T_0$. It takes nonlinear effects, discussed below, for $R,T$ to be significantly different from $R_0, T_0$, respectively, for $\epsilon \le 0.1$.

Results of a linear simulation with a short packet length (larger $\epsilon$) are given in figure 2. Figure 2 shows a time history of the incident wave amplitude $I$, and reflected $R$ and transmitted $T$ wave amplitudes determined by including higher-order linear effects but neglecting nonlinear effects. Also shown are linear wave amplitudes $R_0$ and $T_0$, determined using (4.10) and (4.11). With this relatively large value of $\epsilon$, $R$ and $T$ display large oscillations that continue after the incident wave amplitude is negligible. The frequencies match the values in (4.5), thus these oscillations are homogeneous solutions. Yet the exact solution shows that the homogeneous solution should be zero, therefore these oscillations are artefacts of the numerics. Arguably, such short packets are outside the validity of the basic assumption of a long wave packet.

Figure 2. Incident (black), reflected (red) and transmitted (blue) wave amplitudes with $n_1/k = 1$ and $\epsilon = 0.5$. Higher-order linear effects are included, but nonlinear effects are neglected. Dashed lines are the linear amplitudes $R_0$ (red) and $T_0$ (blue). The incident wave packet profile is Gaussian.

One strategy for eliminating these homogeneous oscillations in the numerical results is to add a condition that $R$ and $T$ be zero when the incident wave packet terminates. However, such a condition has proven to be difficult to achieve, partly due to the ill-defined termination of the incident wave packet: a Gaussian wave packet does not end abruptly, and the ‘raised cosine’ packet ends with a discontinuous second derivative. A wave packet that is allowed to evolve as it ascends will be even more problematic than the above prescribed incident wave packets.

A further issue with forcing $R$ and $T$ to zero at the termination of $I$ concerns the nonlinear effects. If $I$ is zero in the interfacial conditions (3.54) and (3.55), then what remains is a pair of coupled equations for $R$ and $T$: e.g. the nonlinear equations do not require $R$ and $T$ to be zero after $I$ has reached zero. It may be possible for the nonlinear effects to cause wave energy to be absorbed by the interface, then released after $I = 0$.

The homogeneous oscillations may be eliminated simply by neglecting the higher-order linear terms in (3.54) and (3.55). As shown by the above exact solution, the effect of the forced part of the higher-order linear terms is small. The resulting nonlinear algebraic conditions must still be treated numerically using the iterative scheme described in Appendix A. Results for short wave packets are obtained below with this approximation.

5. Nonlinear results

Figure 3 shows a time history of the incident, reflected and transmitted wave amplitudes at the interface for the example case $n_1/k = 1$ and $\epsilon = \alpha = 0.05$. All terms in the interfacial conditions (3.54) and (3.55) are included, and the results are obtained numerically. The incident wave amplitude at the interface is prescribed by (4.12). Also shown in figure 3, with dashed lines, are the linear reflected and transmitted amplitudes $R_0$ and $T_0$. Figure 3 shows that the reflected amplitude $R$ is less than $R_0$, and the transmitted amplitude $T$ is greater than $T_0$ for all time. This means that for this case, the nonlinear effects make the interface more transparent to internal waves, transmitting more energy and reflecting less than the linear version.

Figure 3. Incident (black), reflected (red) and transmitted (blue) wave amplitudes with $n_1/k = 1$ and $\epsilon = 0.05$. Dashed lines are the linear amplitudes $R_0$ and $T_0$. The incident wave packet profile is Gaussian.

Figure 4 also shows a time history of the wave amplitudes at $n_1/k = 1$, but with $\epsilon = \alpha = 0.1$, larger than previously. When $I$ is maximum (at $t = 0$), the profile of $R$ is significantly less than the profile of $R_0$, while $T$ is larger than $T_0$, indicating that the interface is more transparent, as before. For times $t > 0$, after the maximum of $I$ has appeared, both $R$ and $T$ in figure 4 show complex behaviour. This complex behaviour is due to the higher-order linear terms, as in the linear results in figure 2. The value of $\epsilon$ for the results in figure 4 is smaller ($\epsilon = 0.1$) than in figure 2 ($\epsilon = 0.5$), and a simulation without the nonlinear terms with $\epsilon = 0.1$ does not show significant homogeneous oscillations. Thus these erroneous homogeneous oscillations are also problematic in the numerical solution to the nonlinear equations, if higher-order linear terms are included.

Figure 4. Incident (black), reflected (red) and transmitted (blue) wave amplitudes with $n_1/k = 1$ and $\epsilon = 0.1$. Dashed lines are the linear amplitudes $R_0$ and $T_0$. The incident wave packet profile is Gaussian.

The values of $|R|$ and $|T|$ in figure 4 are seen to be non-zero beyond the time when $I$ has become negligibly small. As with the linear simulations that include higher-order linear effects, these oscillations in the nonlinear equations appear to continue indefinitely, and would cause waves to radiate away from the interface in both layers. While it may be possible for the interface to absorb energy due to nonlinear effects and then radiate energy after $I$ is zero, such radiation cannot continue indefinitely. Thus physically meaningful solutions with larger values of $\epsilon$ (short wave packets) must have the homogeneous oscillations suppressed, most easily achieved simply by neglecting the higher-order linear terms in (3.54) and (3.55).

Figure 5 shows results for the nonlinear case without higher-order linear effects with $k = 1$ and $\epsilon = 0.1$, the same parameter values as in figure 4. After the incident wave packet in figure 5 is gone, both $R$ and $T$ also are approximately zero. Thus there are no lingering oscillations when the higher-order linear terms are neglected.

Figure 5. Incident (black), reflected (red) and transmitted (blue) wave amplitudes with $n_1/k = 1$ and $\epsilon = 0.1$. Dashed lines are the linear amplitudes $R_0$ and $T_0$. Nonlinear effects are included, but higher-order linear effects are neglected. The incident wave packet profile is Gaussian.

When the higher-order linear terms are neglected, the instability discussed above is not present, and $R$ and $T$ may be determined for any small value of $n_1/k$. However, the nonlinear algebraic interfacial conditions do not have roots if the value of $\alpha$ is too large, depending also on the value of $n_1/k$. Figure 6 shows the value of $\epsilon$ above which the numerical iteration does not find any solutions. Also shown in figure 6 are parameter values for the cases of figures 3, 4 and 5.

Figure 6. Values of $\epsilon$ above which the nonlinear equations do not have any roots. Marker $a$ shows the parameters for the case in figure 3, while marker $b$ shows parameters for figures 4 and 5.

The numerical results show that nonlinear effects cause a reduction in wave reflection and a corresponding increase in wave transmission. This conclusion is true for all parameter values considered. Thus the interface is more transparent than predicted by linear theory. This is indicated in figure 7, which shows the maximum values of $R$ and $T$ for an interval of $n_1/k$ with $\epsilon = 0.05$. The maximum values of $R$ and $T$ occur when $I$ is also maximum, and the data in figure 7 correspond to maximum $I$. For all parameter values, $|R| < |R_0|$ and $|T| > |T_0|$. Furthermore, larger values of $n_1/k$ show a larger increase in $|T|$ and decrease in $|R|$ compared to linear values. This same trend with $n_1/k$ is found with all other values of $\epsilon$.

Figure 7. Maximum reflected (red) and transmitted (blue) wave amplitudes with $\epsilon = 0.05$ over an interval of $n_1/k$. Dashed lines are corresponding maximum values of the linear amplitudes $R_0$ and $T_0$. Nonlinear effects are included, but higher-order linear effects are neglected. The incident wave packet profile is Gaussian.

The interface becomes more transparent to internal waves as $\epsilon$ increases. This is shown in figure 8 for $n_1/k = 0.4$ and in figure 9 for $n_1/k = 1$. Both figures show the maximum values of $|R|$ and $|T|$, as well as the maxima of $|R_0|$ and $|T_0|$. Clearly, $|R| < |R_0|$ and $|T| > |T_0|$ for all values of $\epsilon$ in both figures, indicating a decrease in reflection and increase in transmission with increasing $\epsilon$. There are no results in figure 9 beyond $\epsilon \approx 0.11$. The nonlinear interfacial conditions do not have roots beyond this value, as indicated in figure 6.

Figure 8. Maximum reflected (red) and transmitted (blue) wave amplitudes with $n_1/k = 0.4$ over an interval of $\epsilon$. The maximum value for convergence with $n_1/k = 0.4$ is $\epsilon =0.23$. Dashed lines are corresponding maximum values of the linear amplitudes $R_0$ and $T_0$. Nonlinear effects are included, but higher-order linear effects are neglected. The incident wave packet profile is Gaussian.

Figure 9. Maximum reflected (red) and transmitted (blue) wave amplitudes with $n_1/k = 1$ over an interval of $\epsilon$. The maximum value for convergence with $n_1/k = 1$ is $\epsilon = 0.11$. Dashed lines are corresponding maximum values of the linear amplitudes $R_0$ and $T_0$. Nonlinear effects are included, but higher-order linear effects are neglected. The incident wave packet profile is Gaussian.

Not all nonlinear terms contribute significantly to the reflected and transmitted wave amplitudes. The strongest contributions, $S_d$, $M_d$ and $M_k$, are plotted in figure 10 for the example case $n_1/k = \sqrt {2}$, $\epsilon = 0.05$ and $N_2/N_1 = 2$. The quantity $S_d$ is the interaction of the secondary harmonic with the primary harmonic in the dynamic condition, and can be seen to have the largest value for all times. The value of $S_d$ is largest for $n_1/k > 0.6$, approximately. For smaller values of $n_1/k$, the interaction between the wave-induced mean and the primary harmonic $M_d$ is largest. The other contributions, $S_k$, $C_k$ and $C_d$ (not shown in figure 10) remain smaller than 10 % of the maximum contributor for all $n_1/k$. This suggests a simpler model for the interfacial conditions, where $S_d$ and/or $M_d$ are retained and the other terms are neglected.

Figure 10. Profile of nonlinear contributions $M_k$, $M_d$ and $S_d$ for $n_1/k = \sqrt {2}$, $\epsilon = 0.05$ and $N_2/N_1 = 2$. The incident wave packet profile is Gaussian.