2.1. Governing equation

Let $\boldsymbol {u}$ be the velocity field, $p'$ the pressure perturbation from hydrostatic, $\rho '$ the density perturbation from the background stratification, $\rho _{0}$, and $e_i$ the unit vector along the i axis. Then, the three nonlinear governing equations are the conservation of momentum (Euler equation),

(2.2)\begin{equation} \rho_{00} \left( \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} \right) = - \boldsymbol{\nabla} p' - \rho' g \boldsymbol{e}_z , \end{equation}

the conservation of mass,

(2.3)\begin{equation} \frac{\partial \rho'}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} ( \rho_0 + \rho' ) = 0 , \end{equation}

and the conservation of volume,

(2.4)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 . \end{equation}

We re-express these equations in terms of the buoyancy, $b = - {g \rho '}/{\rho _{00}}$, and the streamfunction, $\psi$, which is defined by $\boldsymbol {u} = \boldsymbol {\nabla } \times ( \psi \boldsymbol {e}_y ) = ( -{\partial \psi }/{\partial z} , {\partial \psi }/{\partial x} )$ and automatically satisfies volume conservation (2.4), thereby reducing the number of simultaneous scalar equations to solve from four to three.

Taking the curl of the momentum equation (2.2) and defining the Jacobian determinant of two scalars $\alpha$ and $\beta$,

(2.5)\begin{equation} \left| \frac{\partial ( \alpha, \beta )}{\partial ( x, z )} \right| = \frac{\partial \alpha}{\partial x} \frac{\partial \beta}{\partial z} - \frac{\partial \alpha}{\partial z} \frac{\partial \beta}{\partial x} , \end{equation}

which has the form and algebraic properties of the Poisson bracket in classical Hamiltonian dynamics, yields the vorticity equation,

(2.6)\begin{equation} \frac{\partial}{\partial t} \nabla^2 \psi + \left| \frac{\partial \left( \psi , \nabla^2 \psi \right)}{\partial ( x, z )} \right| = \frac{\partial b}{\partial x} . \end{equation}

The quantity $-\nabla ^2 \psi$ is the vorticity, which points in the $y$ direction for two-dimensional flows.

The conservation of mass (2.3) is transformed by simple substitution of variables,

(2.7)\begin{equation} \frac{\partial b}{\partial t} + \left| \frac{\partial ( \psi , b )}{\partial ( x, z )} \right| = - N^2 \frac{\partial \psi}{\partial x} . \end{equation}

This formulation explicitly shows the buoyancy frequency, $N$, is intrinsic to the flows in a stratified fluid. All of the nonlinear terms are now contained in the two Jacobian determinants, which are the transformation of the advection operator, $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla }$.

We expand the streamfunction, $\psi$, in powers of the small dimensionless amplitude, $a$,

(2.8)\begin{equation} \psi = a \psi_1 + a^2 \psi_2 + a^3 \psi_3 + \dots = \sum_{n=1}^{\infty} {a^n \psi_n} , \end{equation}

and similarly for the buoyancy, $b$. For the following analysis, we assume that the coefficient functions $\psi _n ( x, z, t )$ are no greater than $\textrm {ord} ( 1 )$, which is required for the sum to converge. Substitution of these expansions into the vorticity equation (2.6) gives

(2.9)\begin{equation} \frac{\partial}{\partial t} \nabla^2 \left( \sum_{n=1}^{\infty} a^n \psi_n \right) + \left| \frac{\partial \left( \sum_{p=1}^{\infty} {a^p \psi_p} , \nabla^2 \left( \sum_{q=1}^{\infty} a^q \psi_q \right) \right)}{\partial ( x, z )} \right| = \frac{\partial}{\partial x} \left( \sum_{n=1}^{\infty} {a^n b_n} \right) . \end{equation}

Setting $n = p + q$ in the Jacobian term and noting that $1 \leqslant p = n - q \leqslant n - 1$, so that the summation is now over $n$ and $p$, enables factorisation to yield an outer sum in terms of powers of $a$,

(2.10)\begin{equation} \sum_{n=1}^{\infty} { a^n \left\{ \frac{\partial}{\partial t} \nabla^2 \psi_n + \sum_{p=1}^{n-1}{\left| \frac{\partial \left( \psi_p , \nabla^2 \psi_{n-p} \right)}{\partial ( x, z )} \right|} \right\} } = \sum_{n=1}^{\infty} {a^n \frac{\partial b_n}{\partial x}} . \end{equation}

Similarly, inserting the perturbation expansion (2.8) into the equation of conservation of mass (2.7) and summing over powers of $a$ gives

(2.11)\begin{equation} \sum_{n=1}^{\infty} {a^n \left\{ N^2 \frac{\partial \psi_n}{\partial x} + \sum_{p=1}^{n-1}{\left| \frac{\partial ( \psi_p , b_{n-p} )}{\partial ( x, z )} \right|} \right\} } = - \sum_{n=1}^{\infty} {a^n \frac{\partial b_n}{\partial t}} . \end{equation}

Comparing coefficients of powers of $a$ in the expanded governing equations (2.10) and (2.11) gives the two equations at $\textrm {ord} ( a^n )$,

(2.12a)\begin{gather} \frac{\partial}{\partial t} \nabla^2 \psi_n + \sum_{p=1}^{n-1}{\left| \frac{\partial \left( \psi_p , \nabla^2 \psi_{n-p} \right)}{\partial ( x, z )} \right|} = \frac{\partial b_n}{\partial x} , \end{gather}

(2.12b)\begin{gather}N^2 \frac{\partial \psi_n}{\partial x} + \sum_{p=1}^{n-1}{\left| \frac{\partial \left( \psi_p , b_{n-p} \right)}{\partial ( x, z )} \right|} = - \frac{\partial b_n}{\partial t} . \end{gather}

The buoyancy at each order, $b_n$, can be eliminated by differentiating (2.12a) with respect to $t$ and (2.12b) with respect to $x$ and then adding the resulting equations to give the inhomogeneous internal wave equation for $\psi _n$,

(2.13)\begin{equation} \frac{\partial^2}{\partial t^2} \nabla^2 \psi_n + N^2 \frac{\partial^2 \psi_n}{\partial x^2} = - \sum_{p=1}^{n-1}{\left\{ \frac{\partial}{\partial t} \left| \frac{\partial \left( \psi_p , \nabla^2 \psi_{n-p} \right)}{\partial ( x, z )} \right| + \frac{\partial}{\partial x} \left| \frac{\partial \left( \psi_p , b_{n-p} \right)}{\partial ( x, z )} \right| \right\}}. \end{equation}

The homogeneous part of this equation consists of the sum of two temporal derivatives and two spatial derivatives, so forms a wave equation. Its spatially anisotropic structure yields the unusual properties of internal waves. At first order $( n = 1 )$, the summation vanishes, leaving just the equation for linear internal gravity waves. For all higher-order contributions to the streamfunction, the internal wave equation is inhomogeneous, but all terms in the summation arise from lower orders. Consequently, we can inductively evaluate all orders. This set of equations governs all weakly nonlinear wave–wave interactions in free space. However, especially in the case of a flow driven by a moving material surface, such as of our wave maker, it is necessary to consider in detail the role of boundary conditions.

2.2. Boundary conditions

The kinematic boundary condition on the wave maker is of no penetration. Since it is assumed inviscid, the fluid may slip along the boundary. Because the actuating rods of the wave maker move vertically, the velocity of its surface is in the vertical direction,

(2.14)\begin{equation} \boldsymbol{U} ( x, t ) = \frac{\partial h}{\partial t} \boldsymbol{e}_z . \end{equation}

No penetration of the boundary requires that the normal velocity of the fluid, in the direction of unit vector $\boldsymbol {n}$, matches that of the surface of the wave maker at $z = h ( x, t )$,

(2.15)\begin{equation} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{n} . \end{equation}

Let $\alpha$ be the angle the local tangent to the flexible wave maker surface makes with the horizontal, so that $\tan {\alpha } = {\partial h}/{\partial x}$, then using trigonometry, the normal vector pointing into the fluid can be expressed as $\boldsymbol {n} = ( -\sin {\alpha } , \cos {\alpha } )$. Extracting a common factor of $\cos {\alpha }$, substituting into the boundary condition (2.15), expressing $\boldsymbol {u}$ in terms of $\psi$ and $\boldsymbol {U}$ following (2.14), then written in the order $n_z u_z + n_x u_x = n_z U_z$, we obtain

(2.16)\begin{equation} \left. \frac{\partial \psi}{\partial x} \right|_{z = h} + \frac{\partial h}{\partial x} \left. \frac{\partial \psi}{\partial z} \right|_{z = h} = \frac{\partial h}{\partial t} . \end{equation}

A physical interpretation of this equation is that there is no penetration of the fluid material surfaces, $u_z = {\textrm {D}h}/{\textrm {D}t}$, where we have used the material (total) time derivative.

Solving partial differential equations on domains with time-varying, curved boundaries $( z = h )$ is usually analytically intractable and here is no exception. Instead, under the low steepness approximation, $\left | a \right | \ll 1$, we Taylor expand the streamfunction about $z = 0$ with the summation variable $q$,

(2.17)\begin{equation} \sum_{q = 0}^{\infty} \frac{h^{q}}{q!} \left( \frac{\partial}{\partial x} + \frac{\partial h}{\partial x} \frac{\partial}{\partial z} \right) \left. \frac{\partial^{q} \psi}{\partial z^{q}} \right|_{z=0} = \frac{\partial h}{\partial t} . \end{equation}

This expansion will be specialised and fully expanded in powers of $a$ in §§ 3.1 and 5.2 according to the configuration under consideration.

3.1. Boundary conditions

In the weakly nonlinear regime, we assume that the prescribed boundary displacement, $h$, is no greater than $\textrm {ord} ( a )$, so we define $\hat {h} = \textrm {ord} ( 1 )$ such that $h = a \hat {h}$. In addition, we assume that the characteristics of the internal waves each only intersect the wave maker once, which requires $\max{\left | {\partial h}/{\partial x} \right |} < \min{\cot {\varTheta }}$, where, as we will see in § 3.2, $\varTheta$ is the angle of the direction of energy propagation of one such internal wave to the vertical.

Then, we expand the streamfunction (2.8) in the kinematic boundary condition (2.16), collect terms of equal order (powers of $a$) and express as a double sum,

(3.1)\begin{equation} \sum_{q = 0}^{\infty} \sum_{r = 1}^{\infty} \frac{a^{q + r}}{q !} \hat{h}^{q} \left( \frac{\partial}{\partial x} + a \frac{\partial \hat{h}}{\partial x} \frac{\partial}{\partial z} \right) \left. \frac{\partial^{q} \psi_{r}}{\partial z^{q}} \right|_{z=0} = a \frac{\partial \hat{h}}{\partial t} . \end{equation}

Now all quantities are either $\textrm {ord} ( 1 )$ or are powers of $a$, so this has been fully expanded. It is convenient to factor out powers of $a$ and sum over them with summation variable $n$ and use separate inner summations over $q$ for each term. After adjusting the summation limits accordingly, we have

(3.2)\begin{equation} \sum_{n = 1}^{\infty} a^{n} \left\{ \sum_{q = 0}^{n - 1} \frac{\hat{h}^q}{q !} \left. \frac{\partial^{q + 1} \psi_{n - q}}{\partial x \, \partial z^{q}} \right|_{z=0} + \sum_{q = 0}^{n - 2} \frac{\hat{h}^q}{q !} \frac{\partial \hat{h}}{\partial x} \left. \frac{\partial^{q + 1} \psi_{n - q - 1}}{\partial z^{q + 1}} \right|_{z=0} \right\} = a \frac{\partial \hat{h}}{\partial t} . \end{equation}

The system is forced only at $\textrm {ord} ( a )$, so the right-hand side contains no contributions for $n \geqslant 2$, and terms of those orders on the left-hand side must themselves balance.

At $\textrm {ord} ( a^{n} )$, the $q = 0$ term in the first summation reduces to $\left . {\partial \psi _{n}}/{\partial x} \right |_{z=0}$. The remaining terms in the first $q$ summation all arise from the Taylor's expansion that extrapolates evaluation of the vertical fluid velocity from $z = 0$ to the material surface at $z = h$. The second $q$ summation contains corrections due to the variations of the surface normal, $\boldsymbol {n}$, about the vertical and these unavoidably contain the horizontal fluid velocity, which we also Taylor expand to extrapolate from $z = 0$ to $z = h$. Except for the $q = 0$ term in the first summation, which yields $\psi _{n}$, all terms that appear at $\textrm {ord} ( a^{n} )$ are contributions from lower orders. The forcing of the governing equation for $\psi _{n}$ (2.13) also only depends on lower orders. Hence, as shown in figure 1, there exists a unidirectional cascade of dependence from lower- to higher-order streamfunction contributions.

Figure 1. Graph of dependencies of contributions to the streamfunction at each order. The black triangular arrows indicate vertical $( z )$ derivatives and grey triangular arrows indicate horizontal $( x )$ derivatives. The other arrows only show the direction of dependence; no other operations occur.

In addition to the kinematic boundary condition (3.2), the solution must satisfy causality: the time-averaged energy flux must be directed away from $z \leqslant h$ for all components of the generated flow. For internal waves, this is equivalent to saying the group velocity has a positive vertical component. Let the time average over one period of oscillation be denoted by angle brackets $\left \langle \cdot \right \rangle$. Then, causality requires $\left \langle\,p' w \right \rangle \geqslant 0$ for all linearly independent components of the flow (derived, for example, in Dobra Reference Dobra2018, pp. 143–144).

3.2. D'Alembert's solution for arbitrary boundary displacements

Setting $n=1$ in the expansion of the governing equation (2.13) yields the first-order contribution to the streamfunction, $\psi _1$,

(3.3)\begin{equation} \frac{\partial^2}{\partial t^2} \nabla^2 \psi_1 + N^2 \frac{\partial^2 \psi_1}{\partial x^2} = 0 . \end{equation}

This is the linearised form of the wave equation for internal waves. As noted earlier, it has an anisotropic structure, and here we use a method of characteristics that generalises d'Alembert's solution to the classical wave equation (d'Alembert Reference d'Alembert1747). While a Fourier transform could be performed to obtain a dispersion relation directly, in general a Fourier approach cannot be used at higher orders containing quadratic Jacobian determinants, except for cases exhibiting a special symmetry, which we discovered and will report in § 3.3. Furthermore, our approach identifies the hyperbolic structure and the geometry of the characteristics, which we have shown in Dobra (Reference Dobra2018) is an important consideration for wave–wave interactions. The algebra given here is a preparatory step for extension to higher-order harmonics from a monochromatic boundary displacement, which we will discuss in § 3.4.

The linear kinematic boundary condition is given by taking all terms of $\textrm {ord} ( a )$ in the expansion (3.2) $( n = 1, q = 0 )$,

(3.4)\begin{equation} \left. \frac{\partial \psi_1}{\partial x} \right|_{z = 0} = \frac{\partial \hat{h}}{\partial t} . \end{equation}

We integrate this with respect to $x$,

(3.5)\begin{equation} \left. \psi_1 \right|_{z = 0} = \int{\frac{\partial \hat{h}}{\partial t} \, \textrm{d} x} , \end{equation}

where the arbitrary constant of integration will be chosen such that $\psi$ represents the perturbation to the (constant) background streamfunction with no net volume flux through $z = 0$; in other words, $\left \langle \left . \psi _1 \right |_{z = 0} \right \rangle = 0$.

Any smooth boundary displacement profile can be expressed as a real Fourier transform,

(3.6)\begin{equation} \hat{h} = \iint{A ( k , \omega ) \sin{( kx - \omega t )} + B ( k , \omega ) \cos{( kx - \omega t )} \textrm{d} \omega \, \textrm{d} k} , \end{equation}

where the functions $A$ and $B$ of $k$ and $\omega$ are the Fourier coefficients. Substituting this form into the kinematic boundary condition (3.5) gives

(3.7)\begin{equation} \left. \psi_1 \right|_{z = 0} = - \frac{\omega}{k} \iint{A ( k , \omega ) \sin{( kx - \omega t )} + B ( k , \omega ) \cos{( kx - \omega t )} \,\textrm{d} \omega \, \textrm{d} k} . \end{equation}

Since the operation of integration, the governing equation (3.3) and the boundary conditions are all linear, we will consider each term independently for a particular $( k , \omega )$ and then integrate over these contributions to recover the full streamfunction field.

Taking only the terms at a particular frequency $\omega$, which we denote as $\psi _{\omega }$, we seek a wave-like sinusoidal solution, so the linear internal wave equation reduces to the two-dimensional partial differential equation

(3.8)\begin{equation} - \omega^2 \nabla^2 \psi_{\omega} + N^2 \frac{\partial^2 \psi_{\omega}}{\partial x^2} = 0 , \end{equation}

which readily rearranges into the form of the classical wave equation,

(3.9)\begin{equation} \left( \frac{N^2}{\omega^2} - 1 \right) \frac{\partial^2 \psi_{\omega}}{\partial x^2} - \frac{\partial^2 \psi_{\omega}}{\partial z^2} = 0 . \end{equation}

In the case $\omega > N$, this is an elliptic equation, so does not admit propagating wave solutions, but instead evanescent waves form, which we will discuss later in this section. Internal waves are the solutions that occur along characteristics when $\omega < N$ and thus the system is hyperbolic. Although elliptic equations can often be solved more readily than hyperbolic equations (for example, Hurley (Reference Hurley1972) used analytic continuation to extend an elliptic solution to propagating internal waves), here we specialise d'Alembert's direct approach for the solution of hyperbolic forms (d'Alembert Reference d'Alembert1747) to linear internal waves. Solutions are projected along the characteristics, so satisfying the boundary condition at $z = 0$ provides a streamfunction everywhere in the fluid interior.

Factorising the hyperbolic differential operator yields

(3.10)\begin{equation} \left( \sqrt{\frac{N^2}{\omega^2} - 1} \frac{\partial}{\partial x} + \frac{\partial}{\partial z} \right) \left( \sqrt{ \frac{N^2}{\omega^2} - 1} \frac{\partial}{\partial x} - \frac{\partial}{\partial z} \right) \psi_{\omega} = 0 . \end{equation}

This form clearly shows the fundamental property of internal waves (when $\omega < N$) that the characteristics of the streamfunction, which are also the streamlines, are parallel at a constant angle to the vertical. Let $\varTheta _1$ be the angle these make with the vertical, where $0 < \varTheta _1 < {{\rm \pi} }/{2}$, and $\eta _{1\pm }$ be the normalised characteristic variables,

(3.11)\begin{equation} \eta_{1\pm} = x \cos{\varTheta_1} \pm z \sin{\varTheta_1} . \end{equation}

The difference in $\eta$ between two parallel characteristics is the perpendicular distance between them in $( x , z )$ space. The derivatives with respect to $\eta _{1\pm }$ are found using the chain rule,

(3.12)\begin{equation} \frac{\partial}{\partial \eta_{1\pm}} = \sec{\varTheta_1} \frac{\partial}{\partial x} \pm \textrm{cosec}{\varTheta_1} \frac{\partial}{\partial z} = \textrm{cosec}{\varTheta_1} \left( \tan{\varTheta_1} \frac{\partial}{\partial x} \pm \frac{\partial}{\partial z} \right) . \end{equation}

Comparing this with the factorised form of the wave equation (3.10) shows that

(3.13)\begin{equation} \frac{\partial^2 \psi_{\omega}}{\partial \eta_{1+} \partial \eta_{1-}} = 0 \end{equation}

and $\tan {\varTheta _1} = \sqrt{(N/\omega)^2 - 1}$, so

(3.14)\begin{equation} \omega = N \cos{\varTheta_1} , \end{equation}

which we identify as the dispersion relation for linear internal waves. Therefore, the characteristics are parallel to the group velocity. Although the tangent function could take either sign, we take $\tan {\varTheta _1}$ to be positive throughout this paper, because it represents a positive square root. The general solution of the transformed equation (3.13) is the sum of two arbitrary functions each of one variable,

(3.15)\begin{equation} \psi_{\omega} = f ( \eta_{1+} ) + g ( \eta_{1-} ) . \end{equation}

Applying the boundary condition (3.7) at this chosen frequency, $\omega$, to the general solution implies that this contribution to the streamfunction is of the form

(3.16)\begin{align} \psi_{\omega} &= - \frac{\omega}{k} \int CA \sin{[ k ( x + z \tan{\varTheta_1} ) - \omega t ]} + ( 1 - C ) A \sin{[ k ( x - z \tan{\varTheta_1} ) - \omega t ]}\nonumber\\ & \quad + DB \cos{[ k ( x + z \tan{\varTheta_1} ) - \omega t ]} + ( 1 - D ) B \cos{[ k ( x - z \tan{\varTheta_1} ) - \omega t ]} \, \textrm{d} k , \end{align}

where $C$ and $D$ are constants to be determined from the causality condition, $\left \langle \kern0.05em p_{\omega } w_{\omega } \right \rangle \geqslant 0$ on $z = 0$. In fact, due to the characteristic nature of the system, this condition holds everywhere in the fluid domain, $z \geqslant 0$. The vertical velocity component $w_{\omega }$ is given by ${\partial \psi _{\omega }}/{\partial x}$, and we find the corresponding pressure perturbation by integrating the linearised horizontal momentum equation (2.2) with respect to $x$,

(3.17)\begin{equation} p_{\omega} = \rho_{00} \int{\frac{\partial^2 \psi_{\omega}}{\partial t \partial z} \, \textrm{d} x} , \end{equation}

where the integration constant will be set to zero to ensure zero time-averaged perturbation. Alternatively, one could derive this by considering the force balance on a fluid parcel; see Dobra (Reference Dobra2018, pp. 144–145) for details. We now consider each sinusoid in turn, noting that time averages of cross terms equal zero. For the first sinusoid,

(3.18)\begin{align} \left\langle \kern0.05em p_{\omega} w_{\omega} \right\rangle = C^2 A^2 \left\langle - \rho_{00} \frac{\omega k \tan{\varTheta_1}}{k} k \cos^2{[ k ( x + z \tan{\varTheta_1} ) - \omega t ]} \right\rangle = - \frac{1}{2} C^2 A^2 \rho_{00} k \omega \tan{\varTheta_1} , \end{align}

where we have used that the mean square of a sinusoid is half its amplitude. In preparation for considering phase-locked waves in § 3.3, we define $c_x = {\omega }/{k}$ to be the horizontal phase velocity, and so

(3.19)\begin{equation} \left\langle \kern0.05em p_{\omega} w_{\omega} \right\rangle = - \tfrac{1}{2} C^2 A^2 \rho_{00} k^2 c_x \tan{\varTheta_1} . \end{equation}

Causality is only satisfied for waves generated by the lower boundary when this quantity is positive, so $\sin {[ k ( x + z \tan {\varTheta _1} ) - \omega t ]}$ is physical only when $c_x < 0$ (i.e. $k$ and $\omega$ have opposite signs). Conversely, for the second sinusoid, the same method shows that $\left \langle \kern0.05em p_{\omega } w_{\omega } \right \rangle = \frac {1}{2} ( 1 - C )^2 A^2 \rho _{00} k^2 c_x \tan {\varTheta _1}$, so is only causal when $c_x > 0$. These properties hold for the third and fourth sinusoids, when the sine is replaced by a cosine, because it is simply a phase shift. Therefore, the coefficients $C$ and $D$ are either zero or one, according to the sign of the horizontal phase velocity. We succinctly express this using the sign function,

(3.20)\begin{align} \psi_{\omega} &= - \frac{\omega}{k} \int A ( k , \omega ) \sin{[ k ( x - \textrm{sgn}{( k \omega )} z \tan{\varTheta_1} ) - \omega t ]} \notag\\ & \quad + B ( k , \omega ) \cos{[ k ( x - \textrm{sgn}{( k \omega )} z \tan{\varTheta_1} ) - \omega t ]} \,\textrm{d} k . \end{align}

Instead, if $\omega > N$, linear internal waves cannot propagate and are evanescent. Furthermore, the spatial equation (3.9) becomes elliptic, meaning that there are no real characteristics and information at one point propagates throughout the whole domain. We seek a separable solution, which we will denote $\psi _e$, that is harmonic in $x$, so must be exponential in $z$ with growth/decay rate $k \sqrt {1 - (N/\omega)^2}$. In order to satisfy causality, the disturbance decays into the fluid domain, so the contribution to the streamfunction in the evanescent case is

(3.21)\begin{align} \psi_e &= - \frac{\omega}{k} \int A ( k , \omega ) \exp\left(-kz \sqrt{1-\left(\frac{N}{\omega}\right)^2}\right) \sin{( kx - \omega t )} \nonumber\\ &\quad + B ( k , \omega ) \exp\left(-kz \sqrt{1-\left(\frac{N}{\omega}\right)^2}\right) \cos{( kx - \omega t )} \,\textrm{d} k . \end{align}

This is simply the (unstratified) potential flow response, but with a rescaled vertical coordinate, $z \mapsto z \sqrt {1 - (N/\omega)^2}$; potential flow is smoothly recovered in the unstratified limit, $N \to 0$. These forced oscillations are in phase with the boundary forcing. The disturbance extends further up into the fluid as the strength of the stratification increases and does not decay at all at the point where internal waves start to propagate, $N = \omega$. Unlike propagating internal waves, evanescent waves are reversible in time, meaning that it would not be possible to determine if a video of one is being played backwards. Thus, steady evanescent waves do not transport any energy.

Assembling the propagating (3.20) and evanescent (3.21) wave solutions gives the linear contribution to the streamfunction generated by an arbitrary boundary displacement expressed in the form (3.6),

(3.22)\begin{equation} \psi_1 = \int_{-\infty}^{-N}{\psi_e \, \textrm{d} \omega} + \int_{-N}^{N}{\psi_{\omega} \, \textrm{d} \omega} + \int_{N}^{\infty}{\psi_e \, \textrm{d} \omega} . \end{equation}

3.3. Symmetries of phase-locked internal waves

Here, we derive a symmetry of phase-locked internal waves, both propagating and evanescent, which have the same horizontal phase velocity $c_x$. Such propagating waves may have an arbitrary amplitude spectrum according to

(3.23a)\begin{align} \psi &= \int A ( k ) \sin{[ k ( x - \textrm{sgn}{( c_x )} z \tan{\varTheta} - c_x t ) ]} \nonumber\\ &\quad + B ( k ) \cos{[ k ( x - \textrm{sgn}{( c_x )} z \tan{\varTheta} - c_x t ) ]} \,\textrm{d} k , \end{align}

where, from the dispersion relation (3.14), the angle $\varTheta = \cos ^{-1}({k c_x}/{N})$ and thus depends on $k$. The corresponding form of evanescent waves is

(3.23b)\begin{align} \psi = \int \exp\left(-kz \sqrt{1 - \left(\frac{N}{\omega}\right)^2}\right) \left\{A ( k ) \sin{[ k ( x - c_x t ) ]} + B ( k ) \cos{[ k ( x - c_x t ) ]} \right\}\textrm{d} k . \end{align}

This is a general description of travelling wavepackets of both finite and infinite extent along a material surface, such as the surface of the wave maker. This includes classes of problem such as atmospheric lee waves (e.g. Scorer Reference Scorer1949; Dalziel et al. Reference Dalziel, Patterson, Caulfield and Le Brun2011; Dobra et al. Reference Dobra, Lawrie and Dalziel2019), though excludes cases such as standing waves because they are superpositions of waves of opposing phase velocities. For such a propagating wave spectrum, the Jacobian terms (which correspond to the advection terms, $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla }$, of the vorticity equation (2.6) and conservation of mass (2.7)) vanish. This important symmetry shows not only that resonant interactions in the domain interior are not admissible but in fact that all second-order interactions between waves arising from a horizontally phase-locked spectrum are inadmissible. We also note that, although linear, such a spectrum fully satisfies the nonlinear governing equations (2.2)–(2.4) at all amplitudes, which is a remarkable generalisation of this property observed for monochromatic plane waves by McEwan (Reference McEwan1973) and Tabaei & Akylas (Reference Tabaei and Akylas2003).

We now derive this symmetry by first differentiating the phase-locked form of the propagating streamfunction (3.23a) to obtain the negative of the vorticity,

(3.24)\begin{align} \nabla^2 \psi_1 &= - \int \left( k^2 + k^2 \tan^2{\varTheta} \right) A ( k ) \sin{[ k ( x - \textrm{sgn}{( c_x )} z \tan{\varTheta} - c_x t ) ]} \nonumber\\ & \quad + \left( k^2 + k^2 \tan^2{\varTheta} \right) B ( k ) \cos{[ k ( x - \textrm{sgn}{( c_x )} z \tan{\varTheta} - c_x t ) ]} \,\textrm{d} k . \end{align}

Using trigonometry and the dispersion relation (3.14) gives

(3.25)\begin{equation} k^2 + k^2 \tan^2{\varTheta} = k^2 \sec^2{\varTheta} = \frac{N^2}{c_x^2} , \end{equation}

which is a constant and so can be factored out of the integral. Therefore, the vorticity is proportional to the streamfunction with the constant of proportionality depending only on the buoyancy frequency, $N$, and the horizontal phase velocity, $c_x$,

(3.26)\begin{equation} \nabla^2 \psi = - \frac{N^2}{c_x^2} \psi . \end{equation}

Applying the Laplacian to the evanescent form of the streamfunction (3.23b) gives the same result. From the linear and antisymmetric properties of Jacobians, we now show that the Jacobian corresponding to the vorticity is zero,

(3.27)\begin{equation} \left| \frac{\partial \left( \psi , \nabla^2 \psi \right)}{\partial ( x, z )} \right| = \left| \frac{\partial \left( \psi , - \dfrac{N^2}{c_x^2} \psi \right)}{\partial ( x, z )} \right| = - \frac{N^2}{c_x^2} \left| \frac{\partial ( \psi , \psi )}{\partial ( x, z )} \right| = 0 . \end{equation}

Similarly considering the buoyancy, $b$, each term in the integrand is a plane internal wave, which satisfies the linear internal wave equation (3.3), so we calculate the buoyancy for each component separately, denoted by a prime, using the linearised conservation of mass (2.6),

(3.28)\begin{equation} \frac{\partial b'}{\partial t} = - N^2 \frac{\partial \psi'}{\partial x} , \end{equation}

and then integrate over the resulting contributions. For both propagating and evanescent waves, integrating mass conservation with respect to time gives $b' = ({N^2}/{c_x} )\psi '$, where the constant of integration has been set to zero to enforce zero average perturbation. Like the vorticity, the constant of proportionality is independent of the horizontal wavenumber, $k$, so may be factored out of the integral, yielding $b = ({N^2}/{c_x} )\psi$. Therefore, the Jacobian containing the buoyancy is also zero,

(3.29)\begin{equation} \left| \frac{\partial ( \psi , b )}{\partial ( x, z )} \right| = \frac{N^2}{c_x} \left| \frac{\partial ( \psi , \psi )}{\partial ( x, z )} \right| = 0 , \end{equation}

and the symmetry of phase-locked internal waves is proven.

Consequently, for a phase-locked boundary displacement $\hat {h}$, the streamfunction contribution, $\psi _n$, is also phase locked and thus is generated solely at the wave maker surface at all orders. We will now prove this using the strong principle of induction by first assuming that $\psi _q$ is phase locked with horizontal phase velocity $c_x$ for all $q < n$. Then, the Jacobian terms in the expanded internal wave equation (2.13) at $\textrm {ord} ( a^n )$, which depend only on the lower, phase-locked orders, are all zero, so no wave–wave interactions can occur and the fluid response, $\psi _n$, is generated solely at the surface of the wave maker. The kinematic boundary condition (3.2) at $\textrm {ord} ( a^n )$ consists of ${\partial \psi _n}/{\partial x} |_{z=0}$ and terms proportional to

(3.30a,b)\begin{equation} \hat{h}^{q} \left. \frac{\partial^{q + 1} \psi_{n - q}}{\partial x \partial z^{q}} \right|_{z=0} \quad \textrm{and} \quad \hat{h}^{q} \frac{\partial \hat{h}}{\partial x} \left. \frac{\partial^{q + 1} \psi_{n - q - 1}}{\partial z^{q + 1}} \right|_{z=0} , \end{equation}

which all sum to zero for $n \geqslant 2$, or ${\partial \hat {h}}/{\partial t}$ when $n = 1$. By the induction assumption, each of these terms is an integral over products of sines and cosines with uniform horizontal phase velocity $c_x$. The product of a pair of such sinusoids also has phase velocity $c_x$, because, for example,

(3.31)\begin{align} &\cos{[ A ( x - c_x t ) ]} \cos{[ B ( x - c_x t ) ]}\nonumber\\ &\quad = \tfrac{1}{2} \left( \cos{[ ( A + B ) ( x - c_x t ) ]} + \cos{[ ( A - B ) ( x - c_x t ) ]} \right) , \end{align}

where $A$ and $B$ are arbitrary constants, and thus all of the product terms in the kinematic boundary condition have horizontal phase velocity $c_x$. Therefore, integrating the boundary condition with respect to $x$ and setting the integration constant to zero to ensure no net flux through $z = 0$ gives that $\psi _n |_{z = 0}$ is also phase locked. Since the Jacobian determinants are zero, the internal wave equation (2.13) at $\textrm {ord} ( a^n )$ reduces to the linear internal wave equation (3.3), and so, from the linear solution (3.20), the streamfunction contribution $\psi _n$ is phase locked everywhere in the domain. Finally, we already know from the linear solution that $\psi _1$ is phase locked. Therefore, by induction, the streamfunction is phase locked at all orders.

In general, the Jacobian determinant gives the area of the image of a unit element having undergone a coordinate transformation. Here, these zero Jacobian determinants indicate that the transformations into the two-dimensional streamfunction–buoyancy and streamfunction–vorticity spaces are singular for arbitrary superpositions of phase-locked internal waves, namely that all points map onto straight lines through the origin of gradients ${N^2}/{c_x}$ and ${N^2}/{c_x^2}$ respectively. Conversely, the image space remains two-dimensional for an unconstrained superposition of internal waves and other flows.

3.4. Algorithmic evaluation of higher-order contributions for monochromatic boundary displacement

We now present the process by which we obtain contributions to monochromatic boundary displacements for arbitrary order. The steps in this process we divide into a pair of interconnected algorithms 1 and 2, then for convenience we illustrate their use by explicitly calculating key expressions at first, second and third orders in tables 1–3.

Table 1. Kinematic boundary condition at the first three orders for a monochromatic boundary displacement.

In § 3.3, we showed that the expansion of the internal wave equation (2.13) is linear at all orders for a phase-locked boundary displacement. A special case is of a monochromatic sinusoid travelling to the right, which is infinite in extent, $h = A \sin {( kx - \omega t )}$, where we use the convention $k , \omega > 0$. Defining the dimensionless amplitude as $a = Ak$, we have $\hat {h} = ({1}/{k} )\sin {( kx - \omega t )}$. For this case, we may derive analytic expressions for the produced spectrum of harmonics at each of the first three orders by noting that the flow at each order is only generated at the boundary. The expansion of the kinematic boundary condition (3.2) becomes

(3.32)\begin{align} \sum_{n = 1}^{\infty} a^{n} \left\{ \sum_{q = 0}^{n - 1} \frac{\sin^{q}{\phi}}{q ! k^{q}} \left. \frac{\partial^{q + 1} \psi_{n - q}}{\partial x \, \partial z^{q}} \right|_{z=0} + \sum_{q = 0}^{n - 2} \frac{\sin^{q}{\phi} \, \cos{\phi}}{q ! k^{q}} \left. \frac{\partial^{q + 1} \psi_{n - q - 1}}{\partial z^{q + 1}} \right|_{z=0} \right\} = - a \frac{\omega}{k} \cos{\phi} . \end{align}

Since this condition at $\textrm {ord} ( a^n )$ depends on all of the lower orders, the contribution to the streamfunction at each order is evaluated in turn, according to algorithm 1.

Algorithm 1 Calculation of streamfunction $\psi$

To calculate the contribution to the streamfunction at $\textrm {ord} ( a^n )$, denoted by $\psi _n$, firstly we take the first $n$ terms of the outer summation in (3.32). These are shown for the first three orders in table 1. The boundary condition at all orders contains ${\partial \psi _n}/{\partial x} |_{z=0}$, which is the vertical velocity at $\textrm {ord} ( a^n )$. Higher orders also include derivatives of lower-order contributions to the streamfunction, and these are multiplied by sines and cosines of integer multiples of the horizontal phase, $\phi = kx - \omega t$. All derivatives are evaluated at the equilibrium height of the wave maker, $z = 0$. Secondly, we evaluate and substitute for the derivatives of the lower-order contributions to the streamfunction. For example, the required derivatives of $\psi _1$ follow the pattern

(3.33a)\begin{equation} \left. \frac{\partial^{q + 1} \psi_1}{\partial x \, \partial z^{q}} \right|_{z=0} = \begin{cases} ( -1 )^{({q + 2})/{2}} \omega k^{q - 1} \tan^{q}{\varTheta_1} \cos{\phi} & \textrm{for even } q,\\ ( -1 )^{({q + 1})/{2}} \omega k^{q - 1} \tan^{q}{\varTheta_1} \sin{\phi} & \textrm{for odd } q, \end{cases}\end{equation}

and

(3.33b)\begin{equation} \left. \frac{\partial^{q + 1} \psi_1}{\partial z^{q + 1}} \right|_{z=0} = \tan{\varTheta_1} \left. \frac{\partial^{q + 1} \psi_1}{\partial x \partial z^{q}} \right|_{z=0} . \end{equation}

We are left with a product of sines and cosines of several multiples of $\phi$ for each term in the boundary condition. The next stage is to simplify these as sums of harmonics, or equivalently express a Fourier series, using the formulae derived using standard methods in appendix A. We first expand all of the higher harmonic terms into powers of trigonometric functions of the fundamental using, for $p \in \mathbb {Z}$,

(3.34a)\begin{equation} \cos{(\kern0.05em p \phi )} = \sum_{\beta = 0}^{{p}/{2}}{( -1 )^{\beta} \binom{p}{2 \beta} \cos^{p - 2 \beta}{\phi} \sin^{2 \beta}{\phi}} \end{equation}

and

(3.34b)\begin{equation} \sin{(\kern0.05em p \phi )} = \sum_{\beta = 0}^{({p - 1})/{2}}{( -1 )^{\beta} \binom{p}{2 \beta + 1} \cos^{p - 2 \beta - 1}{\phi} \sin^{2 \beta + 1}{\phi}} , \end{equation}

where $\binom {n}{r}$ is the binomial coefficient. Then, we collect the terms into a single product and expand as a series of harmonics using formula (A 15), which we re-express for convenience as algorithm 2. Collecting like terms shows that $\left . {\partial \psi _n}/{\partial x} \right |_{z=0}$ is equal to a sum of harmonics with constant amplitudes. Moreover, we find that these harmonics need to be represented by cosine functions in order to match the symmetry of the boundary displacement about $x = 0$, and up to the $n{\textrm {th}}$ harmonic, denoted by $n \phi$, is included. For odd $n$, all and only the odd-numbered harmonics are present (up to the $n{\textrm {th}}$ harmonic); conversely, for even $n$, we have all and only even-numbered harmonics.

Algorithm 2 Expressing $\cos^{\alpha}{\phi} \sin^{\beta}{\phi}$ as a sum of harmonics.

Such a form is readily integrated with respect to $x$ to give the contribution to the streamfunction, $\psi _n$, evaluated at $z = 0$. We find that it is equal to a sum of sinusoids of phase $n \phi$. The integration constant is set to zero to enforce that the equilibrium height of the wave maker is at $z = 0$.

Since the streamfunction is a discrete sum of linearly independent temporal (and spatial) harmonics along the boundary and it satisfies the linear internal wave equation (3.3), we project each harmonic with a frequency less than the buoyancy frequency along the corresponding characteristics, which are at angle $\varTheta _n$ to the vertical, given by the dispersion relation (3.14). Then, each harmonic takes the form of the linear solution (3.20). The harmonics above the buoyancy frequency generate evanescent waves, whose contribution takes the form of the linear evanescent waves (3.21), with $k$ and $\omega$ multiplied by the appropriate value of $n$. Finally, the contribution to the streamfunction, $\psi _n$, is given by the linear superposition of these propagating and evanescent internal waves, even though the solution is nonlinear. Provided the third harmonic is not evanescent, these contributions are listed in table 2 and are plotted in figure 2 together with a phase plot in physical space. We note that in this case all harmonics are in phase with the boundary displacement.

Figure 2. Predictions for the first three harmonics generated by an input sinusoid of frequency $\omega = 0.4 = 0.25N$ and wavenumber $k = 0.4$, where units for frequency and wavenumber are freely chosen provided they are self-consistent: (a) vertical displacement amplitude; (b) time-averaged energy flux, which has a similar profile to the energy density; and (c) phase profile showing the characteristics where $\psi$ decreases through zero, equivalently where the vertical displacement increases through zero. The expansion is valid for $Ak < 1$, since thereafter some of the characteristics will intersect the flexible boundary more than once.

Table 2. Contributions to the streamfunction at the first three orders, provided $3 \omega < N$.

The leading-order contribution to the $n{\textrm {th}}$ harmonic comes from $n{\textrm {th}}$ order and so grows as $a^n$. Higher-order corrections to this harmonic arise at orders $( n + 2 )$, $( n + 4 )$, $( n + 6 )$, $\dots$, but these corrections become decreasingly significant as $a$ reduces. All of the corrections to the lower harmonics are also given by sine functions, thereby ensuring odd symmetry about $x = 0$. Considering the expression $4 \tan {\varTheta _2} - 3 \tan {\varTheta _1}$ as a function of $\omega$ shows that the third-order correction to the first harmonic reinforces its amplitude for $0 < \omega < \frac {N}{2}$, and this reinforcement is more pronounced for smaller $\omega$. Superlinear growth has been observed previously in experiments (Ermanyuk, Shmakova & Flór Reference Ermanyuk, Shmakova and Flór2017); here, in figure 2(a), we find it appearing on all three propagating harmonics within and beyond the domain of applicability $( kA < \cot {\varTheta _1} )$ in which each internal wave characteristic intersects the sinusoidal boundary exactly once.

Since each mode satisfies the linear equation (3.3), the energy density of each mode is proportional to the square of the amplitude with uniform constant of proportionality $\frac {1}{2} \rho _{00} N^2$ for all harmonics, and their time-averaged energy fluxes $( \langle \kern0.05em p' | \boldsymbol {u} | \rangle )$ are equal to their energy densities multiplied by their group velocities, which are given by $c_x \sin {\varTheta _n}$. Although product terms are developed for all pairs of harmonics in the series, by orthogonality, the time averages of the cross terms are zero, leaving only the linear contributions for each mode. Thus, the energy density and the energy flux have very similar profiles, and as an illustration, we show the energy flux in figure 2(b). We see that for a given monochromatic input, the total energy flux is greater than that contained in the single wave beam predicted by the linear theory. The increased energy flux is not a violation of causality, because the power that the flexible boundary transmits to the fluid is not specified, only the position of its surface.

On the other hand, if at least one of the first three harmonics were evanescent $( 3 \omega > N )$, some of the tangent functions would be replaced by explicit square roots, $\tan {\varTheta _n} = \sqrt {( {N}/{n \omega } )^2 - 1} \mapsto \sqrt {1 - ( {N}/{n \omega } )^2}$, as can be seen in the linear evanescent solution (3.21). Moreover, the $z$ derivatives of an evanescent wave have a different phase to those of the corresponding propagating wave, so if the $m{\textrm {th}}$ harmonic is the lowest evanescent one, all contributions at $( m + 1 ){\textrm {th}}$ and higher orders are shifted in phase relative to the boundary displacement. For any given order of perturbation expansion, the explicit form of the solution depends on the number of propagating harmonics, and we provide up to the third-order contributions in table 3 for when only the first harmonic is propagating. In this case, the first three harmonics again grow superlinearly.

Table 3. Contributions to the streamfunction at the first three orders when the first harmonic is propagating but the second is evanescent, $\omega < N < 2 \omega$. The tangent functions are replaced by explicit square roots when the corresponding frequency is above the buoyancy frequency.

4.1. Magic carpet

The Arbitrary Spectrum Wave Maker (ASWaM, Dobra et al. Reference Dobra, Lawrie and Dalziel2019) is a 1 m flexible section in the base of a tank that is 11 m long, 0.255 m wide and 0.48 m deep. The wave maker's shape is controlled by an array of $100$ Portescap 26DBM10D1B-L linear stepper motors positioned at a pitch of 10 mm along the flexible section, each of which has a vertical resolution of 0.0127 mm and a stroke of 48 mm.

For generating the digital input signals to these stepper motors, we constructed a coupled set of Texas Instruments Beaglebone Blacks (revision C). Each Beaglebone contains a processor where every instruction takes exactly 5 ns, on which we deploy an efficient assembly-language algorithm to issue signals to motor drivers. The signal timings are compiled from analytic functions specified in a text file. The waveforms produced for this paper have a temporal resolution of 30 ns.

The surface of the wave maker is a nylon-faced neoprene foam sheet that is 3 mm thick (similar to that used for wetsuits). At zero displacement, the neoprene surface is flush with the base of the tank, but in operation is deformed by $100$ horizontal rods, each spanning the width of the tank and driven by one of the stepper motors. The lengthwise edges of the sheet are not sealed to the tank wall, and there is a cavity of fluid 80 mm deep beneath the neoprene with both sides of the sheet wetted. However, there is almost no pressure gradient to drive a leakage flow from the underlying cavity into the working section of the tank, provided the chosen waveform conserves volume. To leading order, three-dimensional effects are limited to wall boundary layers.

The neoprene attaches to sleeves around the horizontal rods using hook-and-loop fasteners. The material has some resistance to bending, and conveniently the sleeves can rotate about the rods, minimising the tensile stress in the sheet and the bending moments on the actuators. Our modelling (Dobra et al. Reference Dobra, Lawrie and Dalziel2019) indicates that this produces $C^2$-continuous profiles, despite being specified by a discrete set of actuation rods. We find that the wave maker can reliably produce sinusoids of steepness $\left | {\partial h}/{\partial x} \right | \leqslant 0.6$ without the motors stalling or neoprene detaching.

5.1. Approach

We saw in § 3.4 that a monochromatic boundary forcing produces a full spectrum of internal wave harmonics. However, to study the free-space dynamics of internal waves experimentally, such as for the interaction of two incident wave beams (Smith & Crockett Reference Smith and Crockett2014), wave fields without the extra harmonics are desirable.

One approach is to modify the wave maker so that it is mounted perpendicular to the characteristics of the intended internal wave, thus at angle $\varTheta$ to the horizontal. Then, the velocity of the surface of the modified wave maker, $\boldsymbol {U}$, is always perpendicular to its equilibrium plane, thus has the same direction as the characteristics of the internal wave and hence the fluid velocity, $\boldsymbol {u}$. Therefore, the kinematic boundary condition (2.15) implies that $\boldsymbol {U} = \boldsymbol {u}$ on the wave maker surface, and the monochromatic response exactly satisfies the nonlinear boundary condition, so no additional harmonics are generated. Moreover, a monochromatic sinusoidal internal wave of any amplitude satisfies the linear internal wave equation (3.3) (McEwan Reference McEwan1973), so the response is monochromatic even in the strongly nonlinear regime, provided there is no overturning or shear instability. In particular, we note that for our unmodified horizontal wave maker, critically evanescent internal waves $( \omega = N )$ have vertical characteristics, which are perpendicular to our wave maker, so these are the only monochromatic fluid oscillations for which our horizontal wave maker can eliminate harmonics entirely.

Alternatively, we can use the unique ability of our wave maker to choose a polychromatic input waveform that generates a monochromatic wave at some other angle to the vertical, $\varTheta$. As an example, suppose we wish to solve the inverse problem of constructing the input waveform, $h$, that produces exactly the internal wave field (3.20) of the linear solution in § 3.2, then we would have

(5.1)\begin{equation} \psi = a \hat{\psi} = - \frac{a \omega}{k^2} \sin{[ k ( x - z \tan{\varTheta} ) - \omega t ] } . \end{equation}

We know from § 3 that the wave maker profile

(5.2)\begin{equation} h = a h_1 = A \sin{( kx - \omega t )} = \frac{a}{k} \sin{\phi} \end{equation}

is the leading-order (linear) input required to generate $\psi$, but it also produces higher harmonics that are in this case unwanted. Thus, seeking a solution valid in the weakly nonlinear regime $( | a | \ll 1 )$, this time we expand $h$ and seek a series solution of the form,

(5.3)\begin{equation} h = \sum_{n = 1}^{\infty}{a^n h_n} ,\end{equation}

that generates the monochromatic internal wave field (5.1).

At $\textrm {ord} ( a^2 )$, the second harmonic that is generated by the linear forcing (5.2) is given by $\psi _2$, which is stated in table 2. We can cancel this second harmonic by superposing a corresponding correction, $a^2 h_2$, on the wave maker. Since the linearised kinematic boundary condition (3.4) is ${\partial h}/{\partial t} = \left . {\partial \psi }/{\partial x} \right |_{z=0}$ and it needs to negate $\psi _2$, we deduce that

(5.4)\begin{equation} h_2 = - \int{\left. \frac{\partial \psi_2}{\partial x} \right|_{z=0} \, \textrm{d} t} = - \frac{1}{2 k} \tan{\varTheta} \sin{[ 2 ( k x - \omega t ) ]} ,\end{equation}

with the constant of integration set to zero so that $\left \langle h \right \rangle = 0$. It then follows that

(5.5)\begin{equation} h = \frac{a}{k} \sin{\phi} - \frac{a^2}{2 k} \tan{\varTheta} \sin{( 2 \phi )} + O ( a^3 ) .\end{equation}

However, since the input waveform has now been modified, $h \neq A \sin {\phi }$, the expansion of the kinematic boundary condition (2.16) needs to be recalculated to obtain the internal wave field at $\textrm {ord} ( a^3 )$, $\psi _3$, which would then lead to further such corrections.

5.2. Kinematic boundary condition

Such approaches rapidly become unwieldy. Instead, we take the approach that our wave field is entirely specified by $\psi = a \hat {\psi }$ and by this definition one cannot make higher-order corrections to $\psi$. Instead, we choose to expand the dependent function, $h$, using (5.3) in the Taylor-expanded kinematic boundary condition (2.17),

(5.6)\begin{equation} \sum_{q = 0}^{\infty} \frac{h^{q}}{q!} \left( \frac{\partial}{\partial x} + \frac{\partial h}{\partial x} \frac{\partial}{\partial z} \right) \left. \frac{\partial^{q} \psi}{\partial z^{q}} \right|_{z=0} = \frac{\partial h}{\partial t} , \end{equation}

which gives

(5.7)\begin{equation} \sum_{q = 0}^{\infty} \frac{1}{q!} \left( \sum_{s = 1}^{\infty}{a^{s} h_{s}} \right)^{q} \left( \frac{\partial}{\partial x} + \left( \sum_{r = 1}^{\infty}{a^{r} \frac{\partial h_{r}}{\partial x}} \right) \frac{\partial}{\partial z} \right) \left. \frac{\partial^{q} \left( a \hat{\psi} \right)}{\partial z^{q}} \right|_{z=0} = \sum_{n = 1}^{\infty}{a^{n} \frac{\partial h_{n}}{\partial t}} . \end{equation}

Although a truncation of this expansion of $h$ may generate evanescent harmonics, this possibility does not need to be considered here, because the only fluid flows are those specified in $\hat {\psi }$, which can consist of arbitrary non-internal wave motions.

Next, we manipulate this expansion to factor out all the powers of $a$. Firstly, we re-express the infinite sum raised to an arbitrary finite integer $q$ as a new power series,

(5.8)\begin{equation} \left( \sum_{s = 1}^{\infty}{a^{s} h_{s}} \right)^{q} = a^{q} \left( \sum_{s = 0}^{\infty}{a^{s} h_{s + 1}} \right)^{q} = a^{q} \sum_{s = 0}^{\infty}{a^{s} c_{s}} , \end{equation}

where the coefficients $c_{s}$ are given by the recurrence relation,

(5.9)\begin{equation} c_{s + 1} ( q ) = \frac{1}{( s + 1 ) h_1} \sum_{p = 0}^{s}{[ q ( s + 1 ) - p ( q + 1 ) ] c_{p} h_{s - p + 2}} , \end{equation}

and $c_0 = h_1^{q}$. While aspects of this formula are standard material (see, for example, Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014), appendix B contains our derivation, from which we obtain the next three coefficients,

(5.10a)\begin{gather} c_1 = q h_1^{q - 1} h_2 , \end{gather}

(5.10b)\begin{gather}c_2 = \frac{1}{2 h_1} [ 2 q h_3 c_0 + ( q - 1 ) h_2 c_1 ] = q h_1^{q - 1} h_3 + \frac{1}{2} q ( q - 1 ) h_1^{q - 2} h_2^2 , \\[-2.2pc]\nonumber\end{gather}

(5.10c)\begin{align} c_3 &= \frac{1}{3 h_1} [ 3 q h_4 c_0 + ( 2 q - 1 ) h_3 c_1 + ( q - 2 ) h_2 c_2 ] \nonumber\\ &= q h_1^{q - 1} h_4 + q ( q - 1 ) h_1^{q - 2} h_2 h_3 + \frac{1}{6} q ( q - 1 ) ( q - 2 ) h_1^{q - 3} h_2^3 . \end{align}

Then, the kinematic boundary condition becomes

(5.11)\begin{equation} \sum_{q = 0}^{\infty} \frac{a^{q + 1}}{q!} \left( \sum_{s = 0}^{\infty}{a^{s} c_{s} \left( q \right)} \right) \left( \left. \frac{\partial^{q + 1} \hat{\psi}}{\partial x \partial z^{q}} \right|_{z=0} + \left. \frac{\partial^{q + 1} \hat{\psi}}{\partial z^{q + 1}} \right|_{z=0} \sum_{r = 1}^{\infty}{a^{r} \frac{\partial h_{r}}{\partial x}} \right) = \sum_{n = 1}^{\infty}{a^{n} \frac{\partial h_{n}}{\partial t}} . \end{equation}

The first term can be straightforwardly rearranged to isolate powers of $a$. The second term requires the Cauchy product (B 5), which evaluates the product of two summations, before the reorganisation in powers of $a$,

(5.12)\begin{equation} \sum_{q = 0}^{\infty} \sum_{s = 0}^{\infty} \frac{a^{q + s + 1}}{q!} \left( c_{s} ( q ) \left. \frac{\partial^{q + 1} \hat{\psi}}{\partial x \partial z^{q}} \right|_{z=0} + \left. \frac{\partial^{q + 1} \hat{\psi}}{\partial z^{q + 1}} \right|_{z=0} \sum_{r = 1}^{s}{c_{s - r} ( q ) \frac{\partial h_{r}}{\partial x}} \right) = \sum_{n = 1}^{\infty}{a^{n} \frac{\partial h_{n}}{\partial t}} . \end{equation}

Finally, letting $n = q + s + 1$ and adjusting the summation limits accordingly gives the expansion of the kinematic boundary condition factorised into powers of $a$,

(5.13)\begin{align} \sum_{n = 1}^{\infty} \sum_{q = 0}^{n - 1} \frac{a^{n}}{q!} \left( c_{n - q - 1} ( q ) \left. \frac{\partial^{q + 1} \hat{\psi}}{\partial x \partial z^{q}} \right|_{z=0} + \left. \frac{\partial^{q + 1} \hat{\psi}}{\partial z^{q + 1}} \right|_{z=0} \sum_{r = 1}^{n - q - 1}{c_{n - q - r - 1} ( q ) \frac{\partial h_{r}}{\partial x}} \right) = \sum_{n = 1}^{\infty}{a^{n} \frac{\partial h_{n}}{\partial t}} . \end{align}

In this structure, at $\textrm {ord} ( a^{n} )$, $h_{n}$ only appears on the right-hand side. On the left-hand side, the $c_{s}$ terms produce orders of $h$ up to $s + 1$, but $c_0 ( 0 ) = 1$ and $c_{s} ( 0 ) = 0$ for $s \geqslant 1$, so the highest order appearing is $h_{n - 1}$. Thus, $h_{n}$ depends only on lower-order contributions to the solution, and we obtain a similar hierarchy of dependencies to that found for $\psi _{n}$ in § 2.2 (depicted in figure 1).

This boundary condition (5.13) holds for any fluid flow, $\hat {\psi }$, in the weakly nonlinear regime, which need not consist of internal waves. Since it is derived only from the kinematic boundary condition (2.15) of no penetration and thus is only evaluated at the boundary, this equation is independent of the fluid dynamics in the interior of the domain, provided the flow is inviscid and incompressible, and holds for arbitrary density stratifications, or indeed no stratification at all. As a result, not only can we prescribe the wave maker displacement, $h ( x, t )$, required for any arbitrary flow field, but we can also solve the inverse problem of deducing a suitable displacement on the wave maker that will fully absorb any incoming waves: a non-reflecting boundary condition for internal waves. Furthermore, given sufficiently many spatially separate measurements of velocity distant from $z = 0$, the Taylor's expansion at $z = 0$ can be computed and thus the spectrum of the source may be inferred.

5.3. Algorithmic calculation of boundary displacement for a single spectrum of internal wave harmonics

We consider a single spectrum of harmonics to be one arising from a common fundamental, so have frequencies $n \omega$ that are integer multiples of the fundamental and have a common horizontal phase velocity, $c_x$, which restricts the wavevectors to be $\boldsymbol {k}_n = ( nk , -nk \tan {\varTheta _n} )$. This is sufficiently general to admit a polychromatic spectrum constructed with arbitrary amplitudes of such harmonics to form a Fourier series and thus may represent arbitrary translating periodic shapes. In this section, we present a procedure to explicitly calculate order-by-order the boundary displacement, $h$, required to generate a single spectrum of internal wave harmonics with streamfunction $\psi = a \hat {\psi }$; this is summarised in algorithm 3. As an example, we illustrate how a polychromatic spectrum of three harmonics would be expanded to obtain $h$ correct to second order, with related expressions listed in tables 4 and 5. We then specialise to a monochromatic wave and give the corresponding boundary displacement in table 6.

Table 4. Kinematic boundary condition at the first three orders, after the $c_s$ coefficients have been expanded in terms of $h_q$.

Table 5. Contributions to the boundary displacement at the first two orders that generate three in-phase internal wave harmonics (5.14).

Table 6. Contributions to the boundary displacement at the first three orders that generate a monochromatic internal wave (5.1).

Algorithm 3 Calculation of boundary displacement, h, to obtain a single set of internal wave harmonics with a common phase angle.

To calculate $h_n$, first we take all of the terms at $\textrm {ord} ( a^n )$ in the kinematic boundary condition (5.13). We note that the linear condition is the same as for the forwards problem (3.4). Second, we evaluate the coefficients $c_s$ in terms of $h_q$ and substitute these into the boundary condition; these are listed for the first three orders in table 4. Third, we evaluate and substitute all of the required derivatives and boundary displacement contributions. The expansion is now a sum of products of sines and cosines with phases of the form $\alpha \phi$, where $\alpha \in \mathbb {Z}$. Exactly as in § 3.4, we re-express these as a sum of terms of the form $\sin ^{\alpha }{\phi } \cos ^{\beta }{\phi }$, where $\alpha , \beta \in \mathbb {Z}$, using the general compound angle formulae (3.34), and then convert them to a sum of harmonics using algorithm 2 (see appendix A for derivations of these formulae). After simplification, we are left with ${\partial h_n}/{\partial t}$ equal to a sum of harmonics with fundamental phase $\phi$. We integrate this with respect to time, $t$, setting the integration constant to zero to enforce no net displacement. This yields the contribution to $h$ at $n{\textrm {th}}$ order.

For example, the contributions to the boundary displacement correct to second order, $h = a h_1 + a^2 h_2$, for a polychromatic internal wave field consisting of three harmonics that are in phase at $z = 0$,

(5.14)\begin{align} \psi = a \hat{\psi} &= A_1 \sin{[ k ( x - z \tan{\varTheta_1} ) - \omega t ]} + A_2 \sin{ [ 2 k ( x - z \tan{\varTheta_2} ) - 2 \omega t ]} \nonumber\\ & \quad + A_3 \sin{[ 3 k ( x - z \tan{\varTheta_3} ) - 3 \omega t ]} , \end{align}

are listed in table 5. In line with § 3, we define $a$ to be the characteristic steepness of the first harmonic of $h$ predicted by linear theory, $a = {A_1 k^2}/{\omega }$. Expanding to third order would introduce up to five harmonics along the boundary, but if expanded to all orders, these components would cancel to produce only three internal wave harmonics.

Specialising further to generate a monochromatic internal wave field (5.1) with $A_1 = - {a \omega }/{k^2}$ and $A_2 = A_3 = 0$, we note that ${\partial }/{\partial z} = - \tan {\varTheta } ({\partial }/{\partial x})$ due to the characteristic structure, so the kinematic boundary condition (5.13) specialises to

(5.15a)\begin{equation} \sum_{n = 1}^{\infty} \sum_{q = 0}^{n - 1} \frac{a^{n}}{q!} \left. \frac{\partial^{q + 1} \hat{\psi}}{\partial x \, \partial z^{q}} \right|_{z=0} \left( c_{n - q - 1} - \tan{\varTheta} \sum_{r = 1}^{n - q - 1}{c_{n - q - r - 1} \frac{\partial h_{r}}{\partial x}} \right) = \sum_{n = 1}^{\infty}{a^{n} \frac{\partial h_{n}}{\partial t}} . \end{equation}

Recalling that $\phi = kx - \omega t$, the derivatives of the streamfunction, following formula (3.33a), are given by

(5.15b)\begin{equation} \left. \frac{\partial^{q + 1} \hat{\psi}}{\partial x \partial z^{q}} \right|_{z=0} = \begin{cases} ( -1 )^{({q + 2})/{2}} \omega k^{q - 1} \tan^{q}{\varTheta} \cos{\phi} & \textrm{for even } q \\ ( -1 )^{({q + 1})/{2}} \omega k^{q - 1} \tan^{q}{\varTheta} \sin{\phi} & \textrm{for odd } q.\end{cases} \end{equation}

The contributions to the boundary displacement at the first three orders that generate a monochromatic internal wave are listed in table 6. The first two orders agree with the solution for $h$ (5.5) inferred from the internal wave field generated by a monochromatic boundary forcing in § 5.1. As we would expect from the forwards problem at third order, a third harmonic is required on the boundary to eliminate the third harmonic internal wave that would be generated by a monochromatic boundary displacement. However, this is not simply the negative of the third-order wave field, $\psi _3$ (as listed in table 2), generated by a monochromatic forcing along the boundary. Nonetheless, it does exhibit a third-order reduction that is cubic in $a$ to the amplitude fundamental frequency along the wave maker. This qualitatively agrees with the observation in § 3.4 that there is a cubically increasing response in the fundamental frequency internal wave due to a monochromatic forcing, so we expect a cubically decreasing input to counteract this and generate an internal wave field of a given amplitude.

We remark that we could have alternatively derived the expanded kinematic boundary condition for a monochromatic internal wave (5.15) by directly considering the fluid velocities projected onto the direction of motion of the wave maker. Doing so for arbitrarily large amplitudes produces physical inconsistencies, because our wave maker cannot take multiple values of $h$ at any value of $x$. However, within the single-valued constraint, it is possible to compute an $h ( x, t )$ that matches a wave of arbitrary amplitude. One obtains a strongly nonlinear equation where the dependent variable appears both inside and outside a trigonometric function. This can be resolved by Taylor expanding on those trigonometric functions and this leads to an expansion in $h$ that is identical to equation (5.15). The details of this calculation can be found in appendix C.

5.4. Experimental verification

We experimentally tested the predictions for a single spectrum of harmonics in § 5.3 using the apparatus and method described in §§ 4.1 and 4.2. For these experiments, the tank contained a nearly linear stratification of buoyancy frequency $N = {1.4}\ \textrm {rad}\ \textrm {s}^{-1}$.

Initially, we displaced the magic carpet with a right-travelling monochromatic sinusoid of frequency $\omega = {0.3}\ \textrm {rad}\ \textrm {s}^{-1} = 0.21 N$, wavenumber $k = {40}\ \textrm {rad}\ \textrm {m}^{-1}$ and steady amplitude $A = {4}\ \textrm {mm}$, giving $Ak = 0.16$; the resulting wave field is shown in figure 6(a). As expected from § 3, there is a dominant first harmonic plus a visible second harmonic, but negligible third harmonic. In contrast, we applied the corresponding second-order correction of table 6 in figure 6(b) to almost eliminate the second harmonic but consequently generated a significant third harmonic. We were unable to completely remove the second harmonic using our theoretical waveform because of the nonlinear stratification and flow in the boundary layer highlighted in § 4.3, which cannot be accommodated in this solution. Nevertheless, we have demonstrated a useful technique in the experimental study of internal waves: the substantial attenuation of an unwanted harmonic, which allows a clearer view of the desired fundamental wave beam.

Figure 6. Vertical gradient of the normalised density perturbation $({1}/{\rho _{00}}) ({\partial \rho '}/{\partial z})$ for: (a) a monochromatic sinusoid of amplitude 4 mm, frequency ${0.3}\ \textrm {rad}\ \textrm {s}^{-1}$ and wavenumber ${40}\ \textrm {rad}\ \textrm {m}^{-1}$ in a stratification where the harmonic sequence decays in amplitude; and (b) the corresponding polychromatic input to remove the second-order contributions to the second harmonic, which in this configuration generates a significant third harmonic, as expected. Harmonic analysis confirms that wavy perturbations in phase lines are not intrinsic to the first harmonic.

To test the polychromatic expansion given in table 5, we estimated the amplitudes of the three internal wave harmonics in figure 6(b) using the method in § 4.2 and then reconstructed the corresponding theoretical boundary displacement correct to second order. We found that the second and third harmonics were in antiphase relative to the first harmonic, so we multiplied the corresponding amplitudes in the model by $-1$. Figure 7 compares the inferred displacement, shown with a solid line, with the actual input along the ‘magic carpet’ linearly scaled by a factor of $0.19$, shown with a dashed line. A pure sinusoid is also drawn in dots to demonstrate the modulation of a sinusoid introduced by our expansion (5.13). The very similar shapes of the inferred and input waveforms, except at phases corresponding to distance 0 m along the wave maker, confirm that the second-order correction accurately determines the amplitude of the second harmonic relative to the first harmonic. The small disagreement between the two curves arises principally from an overestimation of the third harmonic. This is partially due to already identified difficulty in measuring the amplitudes of weak harmonics but also due to the boundary layer around the wave maker. The calculated profile is in fact for a material surface just outside the boundary layer. Despite this small error, we have successfully calculated the boundary displacement required to produce an observed spectrum of waves, with a superior accuracy than would be given by a linear model.

Figure 7. Vertical displacement profile calculated from the experiment in figure 6(b) showing a good match to the input waveform, scaled down linearly to match the amplitudes. Also shown, for reference, is a monochromatic sinusoid.