Interaction of mode-one internal solitary waves of opposite polarity in double-pycnocline stratifications

Abstract

Numerical simulations of the interaction of internal solitary waves (ISWs) of opposite polarity are conducted by solving the incompressible Euler equations under the Boussinesq approximation. A double-pycnocline stratification is used. A method to determine when ISWs of both polarities exist is also presented. The simulations confirm previous work that the interaction of waves of the same polarity are soliton-like; however, here it is shown that when a fast ISW with the same polarity as a Korteweg–de Vries (KdV) solitary wave catches up and interacts with a slower ISW of opposite polarity, the interaction can be far from soliton-like. The energy in the fast KdV-polarity wave can increase by more than a factor of 5 while the energy in the slower negative-KdV-polarity wave can decrease by 50 %. Large trailing wave trains may be generated and in some cases multiple ISWs with KdV polarity may be formed by the interaction.

1. Introduction

Internal solitary waves (ISWs) in a density-stratified fluid have been the subject of many theoretical, numerical and experimental studies and their common occurrence in the coastal ocean has motivated numerous field experiments to study their role in mixing, sediment resuspension, transport and other oceanic processes (Shroyer, Moum & Nast 2010; Jackson, Da Silva & Jeans 2012; Boegman & Stastna 2019). They are often modelled with weakly nonlinear models such as the Korteweg–de Vries (KdV) equation, the Gardner equation (often referred to as the extended KdV equation) and the regularized long-wave (RLW) equation (also known as the PBBM, or BBM, equation after Peregrine, Benjamin, Bona and Mahoney). The former two equations belong to the class of completely integrable nonlinear–dispersive wave equations. An important property of these equations is that their solitary wave solutions are solitons: they interact nonlinearly with other waves while preserving their properties. In particular two solitary waves of different amplitudes that interact have the same shapes before and after their interaction, the only evidence of the interaction being a small phase shift (Chow, Grimshaw & Ding 2005). The RLW equation is not completely integrable. It has solitary wave solutions; however, when two of them interact there is a small change in the wave amplitudes (less than 0.1 %) and a small-amplitude (approximately 1 % of the solitary waves) dispersive wave train is generated (Bona, Pritchard & Scott 1980).

The fully nonlinear and dispersive incompressible Euler equations, with and without the Boussinesq approximation, also have ISW solutions which can be numerically obtained by solving the Dubreil–Jacotin–Long (DJL) equation (Turkington, Eydeland & Wang 1991). Numerical investigations using the Boussinesq approximation have shown that these ISWs are not solitons: their interaction results in small changes in wave amplitude (approximately 2 % or less) and in the generation of a small-amplitude dispersive wave train similar to the predictions of the RLW equation (Lamb 1998). Past investigations have considered the interactions of ISWs of the same polarity (e.g. isopycnal displacements in the same direction). Here we extend these investigations to investigate the interaction of ISWs of opposite polarity. Internal solitary waves of opposite polarity only exist for certain types of stratifications, including double-pycnocline stratifications under certain conditions as described below. In this paper I present a method for determining when ISWs of both polarities exist, based on the existence of conjugate flows. Following this, results of numerical simulations of the interaction of ISWs using the incompressible Euler equations are presented. The results show that the interaction of two ISWs of opposite polarity may result in very large changes in wave amplitudes and, in some cases, in the generation of many ISWs.

In § 2 the numerical model used for the simulations is presented and pertinent aspects of relevant theories are given in § 3. In § 4 conditions under which ISW solutions of the DJL equation of both polarities exist are discussed. Simulations of the interaction of ISWs of the same and opposite polarities are presented in § 5. Finally, some conclusions are presented in § 7.

2. Numerical model and model set-up

The governing equations used herein are the incompressible Euler equations with the Boussinesq approximation. A rigid lid is employed and rotational affects are not considered. The stable background densities considered here have the form of two thin pycnoclines separating three layers of constant densities $\rho _3$ in the lower layer, $\rho _2$ in the middle layer and $\rho _1$ in the upper layer with $\rho _3 > \rho _2 > \rho _1$. In dimensional terms the stratifications have the form

(2.1)\begin{align} \bar{\rho}(z) = \rho_3 - \frac{1}{2}\left(\rho_3-\rho_2\right) \left(1 +\tanh\left(\frac{z-z_2}{d_2}\right)\right) - \frac{1}{2}\left(\rho_2-\rho_1\right) \left(1 +\tanh\left( \frac{z-z_1}{d_1}\right)\right), \end{align}

where $z_2< z_1$ are the centres of the two pycnoclines and the constants $d_1$ and $d_2$ determine their thicknesses. The density is non-dimensionalized and scaled via

(2.2)\begin{equation} \rho = \rho_0\left(1 + \Delta\rho\tilde{\rho}\right), \end{equation}

where $\rho _0 = {1}/{2}(\rho _1 + \rho _3)$ is the reference density used in the Boussinesq approximation and $\rho _0\Delta \rho = \rho _3-\rho _1$. The dimensionless density $\tilde {\rho }$ is then approximately equal to $-0.5$ and 0.5 at the top and bottom of the water column respectively if the pycnoclines are thin compared with their distance from the upper and lower boundaries as assumed here. The remaining terms in the governing equations are non-dimensionalized using the constant water depth $H$ as the length scale, $T = \sqrt {H/(g\Delta \rho )}$ as the time scale and $U = H/T$ as the velocity scale. This gives the dimensionless equations (dropping the tilde on the dimensionless density)

(2.3a)$$\begin{gather} u_t+uu_x+wu_z=- p_x, \end{gather}$$

(2.3b)$$\begin{gather}w_t+uw_x+ww_z=- p_z-\rho, \end{gather}$$

(2.3c)$$\begin{gather}\rho_t+u\rho_x+w\rho_z=0, \end{gather}$$

(2.3d)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0 \end{gather}$$

as the governing equations, where $z\in [-1,0]$ is now dimensionless. Here $\boldsymbol {u} = (u,w)$ is the velocity in the vertical $xz$ plane and $p$ is a dimensionless dynamic pressure. The total dimensional pressure is

(2.4)\begin{equation} p^* = \left(\rho_0 \Delta\rho gH\right) p - \rho_0gHz \end{equation}

up to an arbitrary constant. For simulations of interacting ISWs these equations are solved using a two-dimensional non-hydrostatic internal gravity wave model (Lamb 1994, 2007). The model uses Godunov flux limiting which acts as an implicit large-eddy simulation model and a variable time step satisfying a Courant–Friedrichs–Lewy condition. The dimensionless background stratifications are given by

(2.5)\begin{equation} \bar{\rho}(z) =-\frac{1}{2}\Delta\rho_1\tanh\frac{z-z_1}{d_1} - \frac{1}{2}(1-\Delta\rho_1)\tanh\frac{z-z_2}{d_2}, \end{equation}

with $-1 < z_2 < z_1 <0$ so the upper pycnocline is centred at $z_1$ and the lower pycnocline is centred at $z_2$. The upper-, middle- and lower-layer thickness are taken as $h_1 = -z_1$, $h_2 = z_1-z_2$ and $h_3 = 1+z_2$. We only consider thicknesses $d_1=d_2=0.04$. For stratifications with pycnoclines close to the top or bottom boundaries the bottom to surface density difference may be slightly less than 1. Figure 1 shows some sample stratifications. The simulations are initialized with ISW solutions of the DJL equation (see next section) which are interpolated onto the computational grid.

Figure 1. Sample stratifications of (a) $\bar \rho$ and (b) $N^2$ for $(z_1,z_2) = (-0.25, -0.9)$ and $d_1=d_2=0.04$ for $\Delta \rho _1 = 0.55$ (solid) and $0.6$ (dashed).

3. Theoretical background

In this section some relevant theoretical results are briefly presented. We first discuss the Gardner equation which is an approximate model for weakly nonlinear long waves. This equation is often used to model ISWs in a stratified fluid; however, the results presented below show that it is not useful for the interaction of waves of opposite polarity. Following this is a discussion of the DJL equation which gives exact solitary wave solutions of the governing equations (2.3) which are used to initialize the simulations. The DJL solutions show that solitary wave amplitudes are limited and that as the limiting amplitude is approached the waves grow indefinitely in length. The resulting conjugate flow in the horizontally uniform centre of long waves is discussed in the final subsection. The conjugate flow determines limiting wave amplitudes and is used to determine when ISW solutions of the governing equations exist. The procedure for doing so is described in the following section.

3.1. Weakly nonlinear theory

The Gardner equation (Grimshaw, Pelinovsky & Talipova 1997; Talipova et al. 1999; Grimshaw, Slunyaev & Pelinovsky 2010), also referred to as the extended KdV equation, is an extension of the KdV equation which includes a cubic nonlinear term. Like the KdV equation, it is a model for weakly nonlinear long waves. It has the form

(3.1)\begin{equation} \zeta_t + c_0\zeta_x + \alpha \zeta\zeta_x + \alpha_1\zeta^2\zeta_x + \beta\zeta_{xxx} = 0, \end{equation}

where $t$ is time, $x$ is the horizontal coordinate and $\zeta (x,t)$ is the wave shape which could represent an isopycnal displacement, the surface current, the streamfunction at some depth or some other quantity of interest. Parameter $c_0$ is the linear long-wave propagation speed. It and the nonlinear and dispersive coefficients are determined by the stratification and background currents and are given in terms of vertical structure functions (Lamb & Yan 1996). If $\alpha _1=0$ the Gardner equation reduces to the KdV equation and if $\alpha = 0$ it reduces to the modified KdV equation:

(3.2)\begin{equation} \zeta_t + c_0\zeta_x + \alpha_1\zeta^2\zeta_x + \beta\zeta_{xxx} = 0.\end{equation}

The Gardner equation (3.1) has soliton solutions which can be written in several forms. One convenient form (Helfrich & Melville 2006) is

(3.3)\begin{equation} \zeta_{GE} (x,t) = \frac{a}{b+(1-b)\cosh^{2}\theta} = \frac{a\,\mbox{sech}^{2}\theta}{1-b+b\,\mbox{sech}^2\theta}. \end{equation}

Here $a$ is the wave amplitude and the phase is

(3.4)\begin{equation} \theta = \frac{x-Vt}{\lambda}. \end{equation}

The propagation speed $V$, wave amplitude $a$ and wave width $\lambda$ are related via

(3.5)\begin{equation} V-c_0 = \frac{4\beta}{\lambda^2} = \frac{a}{3} \left(\alpha+\frac{1}{2}\alpha_1 a\right)\end{equation}

and the parameter $b$ is related to $a$ via

(3.6)\begin{equation} b =-\frac{\alpha_1 a}{2\alpha+\alpha_1 a}. \end{equation}

When the cubic coefficient $\alpha _1=0$ we have $b=0$ and we recover the well-known KdV soliton solution. Note that for $\zeta _{GE}$ to be bounded we cannot have $1-b+b\,\mbox {sech}^2\theta = 0$ for any value of $\theta$, i.e. $({b-1})/{b} \notin [0,1]$. This means we must have $b<1$.

Figure 2 summarizes the types of soliton solutions for different signs of $\alpha$ and $\alpha _1$ following Grimshaw et al. (1999). For our purposes the important properties of the solution are as follows. For the KdV equation soliton solutions are waves of elevation/depression if $\alpha$ is positive/negative. Mathematically the waves have no amplitude bound although of course physical limitations exist (e.g. isopycnal displacements cannot exceed the water depth) and the perturbation expansion used to derive the equation breaks down when the soliton amplitude becomes too large. As the KdV solitons get larger they become narrower. For the Gardner equation, if the cubic coefficient $\alpha _1<0$ then the soliton solutions have the same polarity (i.e. the same sign of $a$) as those for the KdV equation but there is now a limiting amplitude $a_{lim} = -\alpha /\alpha _1$ as $b\to 1^+$. As this limiting amplitude is approached the waves broaden and become horizontally uniform in their centre. If $\alpha _1>0$, soliton solutions with the same polarity as the KdV solitons exist again but there is now no amplitude limitation. Like KdV solitons, they become narrower as their amplitude increases. A new feature is the existence of a second branch of soliton solutions with the opposite polarity. These have a minimum amplitude of $-2\alpha /\alpha _1$ with no upper bound. Solitons with this limiting amplitude decay algebraically and are called algebraic solitons. In addition, when $\alpha _1>0$ there is a new class of pulsating wave solutions called breathers (Slyunyaev 2001; Lamb et al. 2007; Grimshaw et al. 2010). Both solitons and breathers preserve their identities after interactions with other waves, only undergoing a small phase shift (Chow et al. 2005).

Figure 2. A schematic illustration of the shapes of internal solitary wave (soliton) solutions of the Gardner equation following Grimshaw, Pelinovsky & Talipova (1999). Coefficients $\alpha$ and $\alpha _1$ are the quadratic and cubic nonlinear coefficients. When $\alpha _1>0$ waves of either polarity exist with no maximum amplitude. Minimum wave amplitudes exist for waves of elevation/depression for negative/positive quadratic coefficients. When $\alpha _1<0$ the polarity of the waves is determined by the sign of the quadratic coefficient. There is no minimum amplitude but now a maximum amplitude exists.

3.2. The DJL equation

The Gardner equation is an approximate model for small-amplitude long waves. For finite-amplitude solitary waves we can use the DJL equation (Turkington et al. 1991; Stastna & Lamb 2002). This gives exact solitary wave solutions of the incompressible Euler equations (no longer solitons) but has the shortcoming that it only provides waves of permanent form and is restricted to mode-one waves. Thus it cannot be used to study the evolution of an evolving wave field which is one of the strengths of the Gardner equation. Switching to a reference frame moving with the unknown wave propagation speed $c$ and seeking a steady solution of the governing equations (2.3) leads to the DJL equation:

(3.7)\begin{equation} \nabla^2\eta + \frac{N^2(z-\eta)}{c^2}\eta = 0, \end{equation}

which must be solved subject to the boundary conditions

(3.8)$$\begin{gather} \eta = 0\quad\text{at}\ z =-1,0, \end{gather}$$

(3.9)$$\begin{gather}\eta \to 0 \quad\text{as}\ x \to \pm\infty. \end{gather}$$

Here

(3.10)\begin{equation} N^2(z) =-g\frac{\textrm{d}\bar{\rho}}{\textrm{d}z}(z) \end{equation}

is the square of the buoyancy frequency of the undisturbed dimensionless stratification, $\eta (x,z)$ is the vertical displacement of the streamline passing through $(x,z)$ from its undisturbed height in the far field and $c$ is the propagation speed of the solitary wave, the value of which is to be determined. Once $\eta (x,z)$ are known the density and velocity fields in the moving reference frame are given by

(3.11)$$\begin{gather} \rho(x,z) = \bar{\rho}(z-\eta(x,z)), \end{gather}$$

(3.12)$$\begin{gather}(u,w) = (c\eta_z, -c\eta_x). \end{gather}$$

The DJL equation is numerically solved using an iterative method for which the available potential energy (APE) of the wave is specified (Turkington et al. 1991) and boundary condition (3.9) is replaced with $\eta = 0$ at $x = \pm L$, where $L$ is much larger than the length of the solitary wave. For a version of the DJL equation without using the Boussinesq approximation, see Turkington et al. (1991) and for a version with background currents, see Stastna & Lamb (2002). Here, the APE is that in an infinitely long domain which is computed by using the background density as the reference density (Lamb 2008).

3.3. Conjugate flows

For many stratifications, including those with strong pycnoclines separated from the boundaries as considered here, as the APE of the ISW is increased the ISWs asymptotically approach a limiting finite amplitude and broaden indefinitely as ${\rm APE}\to \infty$ (Lamb & Wan 1998), similar to the solutions of the Gardner equation for $\alpha _1<0$ illustrated in figure 2. In this limit the ISW becomes horizontally uniform in its centre. This horizontally uniform flow is said to be conjugate to the far-field flow (Benjamin 1966). The streamline displacement in the conjugate flow $\eta _{cf}(z)$ is determined from (3.7) by dropping the $x$ dependence which reduces the problem to the nonlinear second-order ordinary differential equation eigenvalue problem:

(3.13)$$\begin{gather} \eta''_{cf}(z) + \frac{N^2(z-\eta_{cf}(z))}{c_{cf}^2}\eta_{cf}(z) = 0, \end{gather}$$

(3.14)$$\begin{gather}\eta_{cf}(-1) = \eta_{cf}(0) = 0. \end{gather}$$

Solutions of this equation can be found for a range of initial values of $\eta '_{cf}(-1)$ and unlike the eigenvalue problem for linear long waves, obtained by replacing $N^2(z-\eta _{cf})$ by $N^2(z)$, these solutions are not simple scalings of one another. Conservation of momentum requires that solutions must also satisfy an auxiliary condition (Benjamin 1966; Lamb & Wan 1998):

(3.15)\begin{equation} T = \int_{-1}^0(\eta_{cf}'(z))^3\textrm{d}z = 0 \end{equation}

which determines the value of $\eta '_{cf}(-1)$. Numerically, (3.13) is solved for the initial condition $\eta _{cf}'(-1) = s$ and then the integral in (3.15) is evaluated to find $T(s)$. A root-finding method is then used to find values of $s$ for which $T=0$. There can be more than one non-zero solution as illustrated in the next section.

4. Theoretical results: When do ISWs of both polarities exist?

The Gardner equation predicts the existence of solitary waves of both polarities when the cubic nonlinear coefficient is positive. The value of $\alpha _1$, including its sign, is, however, not uniquely determined. It depends on the physical interpretation of $\zeta$ in (3.1) (Lamb & Yan 1996). Because of this non-uniqueness, and the fact that the Gardner equation is an approximate evolution equation for small-amplitude waves, we use conjugate flow solutions to determine regions in parameter space where ISWs of both polarities exist. In the following we refer to solitary waves with the polarity predicted by the KdV equation (e.g. depressions (elevations) if $\alpha$ is less than (greater than) 0) as KdV polarity waves and waves of the opposite polarity as negative KdV polarity waves. The KdV polarity waves are always waves of depression for the stratifications considered here.

To illustrate the procedure consider stratifications with $(z_1,z_2) = (-0.25,-0.9)$ (see figure 1). Figure 3 shows $T(s)$ for five values of $\Delta \rho _1$. In all cases $T(0) = 0$. This gives the trivial zero solution and is not of interest (henceforth a conjugate flow solution is assumed to be non-zero). When $\Delta \rho _1=0$ there is a single pycnocline at $z = -0.9$ and there is a single conjugate flow solution which is an elevation because the pycnocline is in the lower half of the water column (not shown). Only ISWs of elevation exist for this stratification. For $\Delta \rho _1 = 0.05$ this is still the case (figure 3): the solid blue curve has a single non-zero root at $s\approx 0.98$. When the upper pycnocline has been strengthened to $\Delta \rho _1 = 0.2$ (and the lower one correspondingly weakened), there are three non-zero roots of $T(s)=0$: two negative roots and one positive root. The positive root is an elevation similar to that for $\Delta \rho _1 = 0.05$. The two new roots at $s = -0.81$ and $-0.47$ (the latter not visible) are depressions. When $\Delta \rho _1 = 0.6$ the upper pycnocline is stronger than the lower pycnocline. There are now two positive roots of $T(s)=0$ and a single negative root. As $\Delta \rho _1$ increases further the smaller positive root increases, the larger root decreases more slowly until the two positive roots merge and disappear as the local maximum of $T(s)$ drops below 0. This is the case for $\Delta \rho _1 = 0.8$ which has a single negative root at $s = -0.58$. When $\Delta \rho _1 = 1$ there is only a single pycnocline at $z = -0.25$ and because it is above the mid-depth there is only a single non-zero conjugate flow solution which is a depression. Now only ISWs of depression exist.

Figure 3. Plots of $T(s)$ for two stratifications using $(z_1,z_2) = (-0.25, -0.9)$ and $d_1=d_2=0.04$. The orange dashed curve for $\Delta \rho _1 = 0.2$ overlays the solid blue curve for $\Delta \rho _1 = 0.05$ for positive values of $s$.

We will focus on stratifications for which there is one negative root and two positive roots. An analogous discussion would apply to the opposite case by symmetry. The existence of three conjugate flow solutions for some stratifications has been noted before for three-layered flows (i.e. in the limit $d_j\to 0$) for some parameter values (Lamb 2000). For the smaller positive root the propagation speed $c_{cf}$ is less than the linear long-wave propagation speed $c_0$ and hence does not appear to have any physical significance, at least for the existence of ISWs. The larger positive root $c_{cf}$ is greater than $c_0$ if $\Delta \rho _1$ is small enough but is less than $c_0$ when $\Delta \rho _1$ is sufficiently large. In the former case ISWs of elevation exist while in the latter case they do not.

Figure 4 shows properties of some conjugate flow solutions of elevation as a function of $\Delta \rho _1$. In these solutions the height of the lower pycnocline is fixed at $z_2 = -0.9$ and three values of $z_1$ are considered, namely $z_1 = -0.25$, $-0.2$ and $-0.15$, focusing on parameter values for which there is a single negative root of $T(s)=0$ and two positive roots. Solutions are shown for the largest positive root. For these three values of $z_1$ the value of $\Delta \rho _1$ is increased from $0.5$ until the solution no longer exists. Figure 4(a) plots $\max \{\eta _{cf}\}$ as a function of $\Delta \rho _1$. As the strength of the upper pycnocline increases the conjugate flow amplitude decreases, and as the upper pycnocline is raised, and hence is further from the lower pycnocline, the amplitude increases. Figure 4(b) plots $c_{cf}$ as a function of $\Delta \rho _1$. Speed $c_{cf}$ decreases as $\Delta \rho _1$ increases and as the upper pycnocline moves towards the upper boundary. Also shown with dashed lines are linear long-wave propagation speeds $c_0$ over a smaller range of $\Delta \rho _1$ values. A striking difference is that these values increase as $\Delta \rho _1$ increases which is not surprising as linear long-wave propagation speeds increase as the stratification shifts towards the mid-depth (e.g. for a two-layer fluid, for a given density jump across the interface it is maximized when the interface is at the mid-depth) and here as $\Delta \rho _1$ increases a larger portion of the density variation is shifted to the pycnocline that is closest to the mid-depth. In each case there is a value of $\Delta \rho _1$ for which $c_0 = c_{cf}$. For smaller values of $\Delta \rho _1$ we have $c_{cf}>c_0$ and solitary wave solutions of the DJL equation can be found for a range of amplitudes with a lower limiting amplitude that is strictly positive. For larger values of $\Delta \rho _1$ no solitary wave solutions of the DJL equation appear to exist although to my knowledge there is no proof of their non-existence.

Figure 4. Conjugate flow solutions as a function of $\Delta \rho _1$ for the largest positive root of $T(s)=0$ for stratifications with $z_2 = -0.9$ and $d_1=d_2=0.04$. (a) Conjugate flow amplitude $\max \{\eta _{cf}\}$ and (b) conjugate flow propagation speed $c_{cf}$. In (b) the dashed curves indicate the linear long-wave propagation speed.

Figure 5 illustrates regions where ISWs of both polarities coexist on the $z_1$–$\Delta \rho _1$ plane for three locations of the lower pycnocline: $z_2 = -0.9$ (figure 5a), $-0.8$ (figure 5b) and $-0.7$ (figure 5c). The dotted lines indicate conjugate flow boundaries: above/below the upper/lower dotted lines there are no conjugate flow elevations/depressions. The upper/lower solid curve denotes conjugate flows of elevation/depression with $c_{cf} = c_0$. Between these two curves (grey shaded region) ISWs of both polarities exist. Along the dashed curve the quadratic nonlinear coefficient $\alpha = 0$. Coefficient $\alpha$ is positive/negative below/above the dashed curve. In the grey region above the dashed line the KdV equation predicts waves of depression. In this region ISWs (DJL solutions) of depression exist with no minimum amplitude and ISWs of elevation exist with a minimum amplitude. The opposite is the case in the grey region below the dashed line. The black circles in figure 5(a) indicate some of the stratifications used for simulations of ISW interactions.

Figure 5. Region where ISWs of both polarities exist (grey shaded regions): (a) $h_3 = 0.1$; (b) $h_3 = 0.2$; and (c) $h_3 = 0.3$. The dotted curves are the conjugate flow boundaries: above/below the upper/lower dotted curves there are no conjugate flow elevations/depressions. The upper/lower solid curve denotes conjugate flows of elevation/depression with $c_{cf} = c_0$. Along the dashed curve the quadratic nonlinear coefficient $\alpha = 0$. Coefficient $\alpha$ is positive/negative below/above the dashed curve. In the grey region above the dashed line the KdV equation predicts waves of depression: ISWs of depression exist with no minimum amplitude and ISWs of elevation exist with a minimum amplitude. The opposite is the case in the grey region below the dashed line. Three of the stratifications used in the interaction simulations are indicated with black circles in (a). Here $h_1 = -z_1$ is the upper-layer thickness and $h_3 = 1+z_2$ is the lower-layer thickness.

As an example, consider the behaviour along the vertical line $z_1 = -0.25$ ($h_1 = 0.25$) in figure 5(a) ($z_2 = -0.9$). Some example $T(s)$ curves for these pycnocline heights are shown in figure 3. Starting at the top ($\Delta \rho _1 = 1$) there is a single pycnocline at $z_1 = -0.25$. There is only one conjugate flow solution. It is a depression (the pycnocline being above the mid-depth) and all solitary waves are waves of depression with no minimum amplitude. As $\Delta \rho _1$ is decreased a lower pycnocline at $z_2 = -0.9$ appears and it increases in strength as $\Delta \rho _1$ continues to decrease (i.e. as the upper pycnocline at $z_1 = -0.25$ weakens). At $\Delta \rho _1 = 0.78$ conjugate flows of elevation appear (the point where $T(s)$ has a positive double root; see figure 3). At this point $c_{cf}$ is less than the linear long-wave propagation speed (figure 4) and there are no ISWs of elevation. When $\Delta \rho _1 = 0.6$ we now have $c_{cf} = c_0$ and below this curve we have a conjugate flow of elevation with $c_{cf}>c_0$ and ISWs of elevation now exist (they can now be computed by solving the DJL equation). There exist ISWs of depression throughout this range of $\Delta \rho _1$ values. The ISWs of elevation initially have a positive minimum amplitude which decreases to zero as the dashed curve is approached along which $\alpha = 0$. Below the dashed curve there are always ISWs of elevation with no minimum amplitude. It is now ISWs of depression which have a minimum amplitude. The ISWs of depression exist until $\Delta \rho _1$ drops below the lower solid curve. The lower pycnocline is now sufficiently strong relative to the upper pycnocline such that only ISWs of elevation exist. The conjugate flow of depression ceases to exist when $\Delta \rho _1$ drops below the lower dashed curve.

Other types of behaviour are illustrated by this figure. For example for $(z_1,z_2) = (-0.45, -0.9)$ ($h_1 = 0.45$) there are never ISWs of both polarities with the waves of depression having a minimum amplitude (no grey shaded region below the dashed curve) and for $(z_1,z_2) = (-0.49,-0.9)$ there is no region where ISWs of both polarities exist. In terms of the $T(s)$ roots, $T(s)$ never has two negative roots for the first example and while $T(s)$ does have two positive roots for a small range of $\Delta \rho _1$ values in the second example, they always correspond to conjugate flow solutions with $c_{cf}< c_0$.

Moving the lower pycnocline up towards the mid-depth appears to reduce the thickness and range of $z_1$ values for which ISWs of both polarities exists as suggested in figure 5(b,c) for $z_2 = -0.8$ and $-0.7$. When $z_2 = -0.7$, for $z_1$ below $-0.375$ ($h_1>0.375$) there are never ISWs of both polarities: for ISWs of both polarities to exist one of the pycnoclines must be sufficiently far from the mid-depth, a distance that depends on the relative strengths of the two pycnoclines.

Figure 6 shows the wave amplitudes and propagation speeds of ISWs obtained by solving the DJL equation for the stratification $(z_1,z_2) = (-0.25, -0.9)$ and $\Delta \rho _1 = 0.55$ (indicated by the rightmost black circle in figure 5a) as a function of the APE. For this stratification ISWs of KdV polarity are waves of depression ($\alpha <0$) and ISWs of elevation exist. The waves of elevation have a minimum amplitude of 0.210 for an APE of 0.013. The black circles in figure 6 indicate waves used in the simulations of ISW interactions for this stratification in the next section. For waves of negative KdV polarity the DJL solutions break down by forming a solitary wave on a pedestal that fills the domain used by the DJL solver.

Figure 6. Internal solitary wave (a) amplitudes and (b) propagation speeds as a function of the wave APE for stratification $S_1$ which has $(z_1,z_2) = (-0.25, -0.9)$ and $\Delta \rho _1 = 0.55$. Solid: values for ISWs of KdV polarity (depressions). Dashed: values for ISWs of negative KdV polarity (elevations). The black circles indicate waves used in simulations of ISW interactions (see table 2). The dark horizontal grey line in (b) at $c= 0.308969$ indicates the linear long-wave propagation speed.

5. Internal solitary wave interactions

We now consider the interaction of two ISWs, both propagating to the right relative to the fluid. The simulations are done in a reference frame moving with the average propagation speed of the two initial ISWs so the wave on the left, which has the largest propagation speed, moves to the right while the wave on the right moves to the left. Table 1 gives the parameters of the various stratifications considered. Results are only shown for a few of these. Cases discussed in some detail are given in table 2. In all cases $\alpha <0$, i.e. KdV polarity waves are waves of depression. The table includes values of the ratio $r$ of kinetic energy (KE) and APE. Turkington et al. (1991) showed that $r>1$. In case $C_6$ the KE is slightly more than 20 % larger than the APE.

Table 1. Stratifications used in ISW interaction simulations.

Table 2. Cases for which results are presented. Here $r = {\rm KE}/{\rm APE}$ is the ratio of the wave KE to its APE. The amplitude $a$ is the extreme isopycnal displacement. In all cases negative $a$ is a KdV polarity wave while a positive amplitude is a wave with the opposite polarity. Speed $c$ is the wave propagation speed. The ‘s’ and ‘f’ in the wave names indicate the slow and fast waves of the given case. A dash indicates the value is the same as in the previous line.

In the isopycnal waterfall plots presented below, the isopycnal at the centre of the upper pycnocline is used. The upper pycnocline is chosen because for all stratifications considered the linear mode-one eigenfunction has its maximum in the upper pycnocline so vertical displacements for linear waves are larger in the upper pycnocline than in the lower pycnocline. Hence, use of this isopycnal, rather than one in the lower pycnocline, makes it easier to see small-amplitude waves generated by the interaction of the two ISWs. In ISWs of depression the upper pycnocline undergoes a larger displacement than the lower pycnocline does but the opposite is true of ISWs of elevation. Hence the amplitude of ISWs of elevation is under-represented by the use of this isopycnal for the waterfall plots. This under-representation is illustrated in the density contour plots presented below.

5.1. Stratification $S_1$

We first consider the stratification $S_1$ for which $(z_1,z_2) = (-0.25, -0.9)$ and $\Delta \rho _1 = 0.55$ (table 1, figure 1). The stratification is indicated by a solid circle in figure 5(a) and the amplitude and propagation speed of ISWs as a function of APE are presented in figure 6. Here ISWs of elevation exist with a minimum amplitude of 0.210.

Figure 7 shows waterfall plots for four cases using this stratification. In figure 7(a) the interaction of two ISWs of KdV polarity (Case $C_1$, table 2) is shown. Their initial APEs are $0.04$ and $0.01$. The initial wave amplitudes (extreme isopycnal displacements) are $-0.140$ and $-0.231$ (figure 6) and these are also the initial displacements of the isopycnal at the centre of the upper pycnocline. At the scale of this waterfall plot there is no evidence of any additional waves being generated by the interaction; however, some small waves are generated with amplitude 0.0008 which is 0.6 % of the amplitude of the smaller initial ISW (see figure 8). At $t=7000$ the extreme isopycnal displacements in the two ISWs are $-0.139$ and $-0.230$, i.e. a decrease in amplitude by 0.7 % and 0.4 %. For comparison, in a simulation with only the single large wave its amplitude decreased by 0.2 % during the same time period due to numerical dissipation and a slight adjustment of the initial DJL ISWs. The interaction is similar to that of KdV solitary waves with amplitudes that are not too different: during the interaction the larger wave decreases in amplitude and slows down while the wave ahead of it grows and speeds up. There are always two local minimums in the isopycnal displacements. The interaction of these two solitary waves is very soliton-like, similar to the interactions of ISWs in an exponential stratification reported by Lamb (1998).

Figure 7. Waterfall plots showing the time evolution of the isopycnal at the centre of the upper pycnocline. Stratification $S_1$ with $(z_1,z_2) = (-0.25, -0.9)$ and $\Delta \rho _1 = 0.55$. (a) Case $C_1$: interaction of two solitary waves of KdV polarity with APEs of 0.01 (amplitude $-0.140$) and 0.04 (amplitude $-0.231$). (b) Case $C_2$: interaction of two solitary waves of negative KdV polarity with APEs of 0.02 (amplitude $0.258$) and 0.05 (amplitude $0.351$). (c) Case $C_3$: interaction of a fast wave of negative KdV polarity with APE of 0.05 (amplitude 0.351) with a slow wave of KdV polarity with APE of 0.002 (amplitude $-0.059$). (d) Case $C_4$: interaction of a fast wave of KdV polarity with APE of 0.003 (amplitude $-0.075$) with a slow wave of negative KdV polarity with APE of 0.03 (amplitude 0.301). Isopycnal displacements are scaled identically in all cases. The amplitude of the negative polarity waves (waves of elevation) is under-represented because of the use of the upper pycnocline.

Figure 8. Zoom in on the isopycnal at the centre of the upper pycnocline at $t=4500$ for case $C_1$ (same case shown in figure 7a). Small-amplitude linear dispersive waves generated by the interaction lie between $x=-60$ and $-30$.

Figure 7(b) shows the interaction of two ISWs of negative KdV polarity for the same stratification (Case $C_2$). Their initial APEs are $0.02$ and $0.05$ with amplitudes of 0.258 and 0.351 (figure 6, table 2). The interaction generates a much larger trailing wave train and a small-amplitude solitary wave both of which can be seen in the waterfall plot. The displacement of the centre of the upper pycnocline decreases by 0.5 % and 1.9 % for the large and small ISWs respectively.

The interaction of two solitary waves of opposite polarity (Case $C_3$; see figure 6, table 2) with the fast wave being a wave of negative KdV polarity is shown in figure 7(c). The negative KdV polarity wave has an APE of 0.05 (amplitude 0.351) while the KdV polarity wave has an APE of 0.002 (amplitude $-0059$). A few small waves are generated by the interaction which is again soliton-like. The displacement of the upper pycnocline decreases by about 8 % and 0.9 % for the KdV and negative KdV polarity waves respectively.

When the slow and fast waves are negative KdV and KdV polarity waves, respectively, the results are quite different. An example, Case $C_4$ (see figure 6, table 2), is shown in figure 7(d). The initial faster-propagating wave of depression centred at $x = 30$ is a small-amplitude wave of KdV polarity (APE 0.003, amplitude $-0.0746$). The wave of elevation ahead of it at $x = 60$ has negative KdV polarity and is larger (APE 0.03, amplitude $0.3008$). The interaction results in the generation of a much larger-amplitude wave train than in the previous three cases. A striking result of the interaction is that the KdV polarity ISW undergoes a large increase in amplitude. It increases its amplitude by 113.1 % while the negative polarity ISW decreases its amplitude by 14.0 %. The energy of the KdV polarity wave (KE plus APE) increases by a factor of 4.43 from 0.0062 to 0.0273 while the energy of the negative polarity wave decreases by a factor of 0.53 from 0.0646 to 0.0340. The combined energy in the two ISWs has decreased by about 13.4 %. Approximately 9 % of the decrease is due to transfer of energy to the smaller trailing waves and the remaining 4.4 % is loss due to numerical dissipation (next section). An accompanying change in the propagation speeds of the two ISWs is evident. This change is exaggerated in the waterfall plot because of the use of a moving reference frame. Relative to the fluid the propagation speed of the KdV polarity wave increases by 5 % while the propagation speed of the negative KdV polarity wave decreases by 4 %.

Figure 9 shows contour plots of the density field before and after the interaction for Case $C_4$. This plot illustrates how displacements of the upper pycnocline can greatly under-represent the amplitude of the waves of elevation. It also illustrates that isopycnal displacements in the trailing wave train are larger in the upper pycnocline.

Figure 9. Density contour plot of the interaction of two solitary waves of opposite polarities. Case $C_4$ (see figure 7d). (a) Shortly before the interaction at $t=4000$. (b) After the interaction at $t=6000$.

5.2. Stratification $S_2$

The second stratification considered has $(z_1,z_2) = (-0.15,-0.9)$ and $\Delta \rho _1 = 0.6$ (table 1, figure 5a). Compared with stratification $S_1$, the upper pycnocline has been raised and slightly strengthened with a corresponding decrease in the strength of the lower pycnocline. This stratification is indicated by a solid circle in figure 5(a). As for stratification $S_1$, KdV polarity waves are waves of depression and negative polarity waves of elevation exist, this time with a minimum amplitude of 0.144.

The interaction of two ISWs of KdV polarity or two waves of negative KdV polarity is similar to those for the previous stratification. Here we consider the interaction of two ISWs of opposite polarity (Case $C_5$, table 2) with a KdV polarity ISW (${\rm APE} = 0.01$, $a=-0.174$, $c=0.341$) initially trailing a slower-propagating negative KdV polarity ISW (${\rm APE} = 0.03$, $a=0.328$, $c=0.303$) as in Case $C_4$. A waterfall plot illustrating the interaction is presented in figure 10(a). The interaction is much more complicated than it was in Case $C_4$. As in Case $C_4$, the KdV polarity ISW increases its amplitude and propagation speed during the interaction while the negative KdV polarity wave undergoes a reduction in amplitude and propagation speed. The interaction generates several other waves including a KdV polarity solitary wave of depression which initially trails the negative KdV polarity ISW. Its propagation speed is negative in this reference frame. Relative to the fluid it is travelling faster than the negative KdV polarity wave so it catches up and interacts with it. During the interaction its amplitude increases and it then propagates to the right. A new packet of small-amplitude waves and a new KdV polarity ISW are produced. This cycle repeats several times. After each interaction the negative KdV polarity wave decreases in amplitude and each subsequent new KdV polarity ISW is smaller than the preceding one. After the third interaction of ISWs of opposite polarity (at $t\approx 1800$) the KdV polarity ISW is drifting to the left as its propagation speed relative to the fluid is smaller than the speed of the moving reference frame. Contour plots of the density field at three different times are presented in figure 11. A movie for this case is provided in the supplementary movie (movie 1) available at https://doi.org/10.1017/jfm.2023.284.

Figure 10. Waterfall plots showing the time evolution of the isopycnal at the centre of the upper pycnocline. (a) Case $C_5$. Stratification $S_2$ with $(z_1,z_2) = (-0.15, -0.9)$ and $\Delta \rho _1 = 0.6$. (b) Case $C_6$. Interaction of two solitary waves of opposite polarities for the symmetric stratification $S_3$ with $(z_1,z_2) = (-0.1, -0.9)$ and $\Delta \rho _1 = 0.5$. Waves have APEs of $0.005$ (amplitude 0.1602) and $0.02$ (amplitude $-0.2754$).

Figure 11. Density contour plot of the interaction of two solitary waves of opposite polarities. Case $C_5$ (see figure 10a). (a) Before the first interaction at $t=600$. (b) Second interaction at $t=1200$. (c) Third interaction at $t=1800$.

5.3. Stratification 3

The third stratification considered has $(z_1,z_2) = (-0.1,-0.9)$ and $\Delta \rho _1 = 0.5$ (table 1, figure 5a). This stratification is symmetric; hence the quadratic coefficient $\alpha =0$ and the cubic coefficient $\alpha _1$ is uniquely determined and is positive. The interaction of two ISWs of opposite polarity is shown in figure 10(b). The initial ISW of elevation has APE of 0.005 and amplitude of $0.1602$ while the larger initial ISW of depression has APE of 0.02 and amplitude of $-0.2754$. The interaction results in the generation of a packet of trailing small-amplitude waves. This is typical of the results obtained using other symmetric stratifications (e.g. stratifications $S_4$ and $S_6$, results not shown). These trailing waves are much smaller than those produced in the previous two examples of the interaction of ISWs of opposite polarity. After the interaction the energy (KE plus APE) in the small ISW has decreased by 24.0 % while that in the large wave has increased by 3.3 %. The total energy in the two waves has decreased by 2.2 %.

6. Energy considerations

The use of a second-order finite-volume method to solve the governing equations implies the presence of numerical damping and a gradual loss of energy from the wave field. This is illustrated in figure 12 where the time evolution of the total APE, KE and half the total energy in the computational domain is shown for Case $C_3$. The total energy decreases linearly in time until $t\approx 6000$ when waves in the wave train generated by the interaction (which occurs at about $t = 5200$) start to leave the domain (see figure 7c). During the interaction the KE undergoes a sharp, transient decrease while the APE has a corresponding sharp, transient increase. They recover somewhat but the total KE and APE undergo a permanent decrease and increase respectively after the interaction. The total energy shows no such behaviour, continuing to decrease at the same rate during the interaction.

Figure 12. Time evolution of energy in the computational domain for Case $C_3$. Shown are the APE (red), KE (blue) and half the total energy (green).

7. Discussion and summary

There are two main results of this paper. The first is the presentation of a method for determining where in parameter space ISWs of both polarities exist. This method is somewhat tedious as it requires the determination of where the conjugate flow speed is equal to the linear long-wave propagation speed. While a mathematical proof that this method is correct has not been provided, the search for DJL solutions of both polarities strongly supports this contention. The second main result is that numerical simulations using the incompressible Euler equations show that the interaction of ISWs of opposite polarity in sufficiently asymmetric stratifications are generally far from soliton-like if the faster trailing wave has KdV polarity: the interactions can result in large changes in the amplitudes of the two ISWs and a significant transfer of energy to additional waves, including other solitary waves.

Figure 13 shows the energy changes resulting from the interaction of two ISWs of opposite polarity in four series of simulations: two series using the symmetric stratification $S_3$ for two different slow waves (note by symmetry it does not matter which waves are waves of elevation or depression) and one series for each of the asymmetric stratifications $S_1$ and $S_5$ (see table 1). For the asymmetric stratifications the initial leading (slow) wave has negative KdV polarity and the initial faster trailing waves have KdV polarity. The energy changes are plotted as a function of $\Delta c$, the difference in the initial propagation speed of the two waves, which is a proxy for the interaction time of the two waves. Figure 13(a,d) shows the ratio $Ea_s/Eb_s$, where $Eb_s$ and $Ea_s$ are the total energies (APE plus KE) in the slow ISW before (at $t=0$) and after the interaction. Figure 13(b,e) shows the corresponding ratio $Ea_f/Eb_f$ for the fast ISW. The ISW energies after the interaction are computed at the earliest time that the final ISWs have separated from the trailing waves generated by the interaction. For all series the initial slow/fast ISW loses/gains energy during the interaction with the change increasing monotonically as $\Delta c$ decreases, i.e. as the time the two waves interact increases. Large energy changes are possible, especially in series $S_1$ in which the fast wave increases its energy by up to a factor of five after the interaction while the energy in the slow wave decreases by slightly more than half. The total fractional change $(E_a-E_b)/E_b$ in the combined energy of the two ISWs is shown in Figure 13(c,f). Here $E_b = Eb_s+Eb_f$ and $E_a = Ea_s+Ea_f$ are the combined energies in the two ISWs before and after the interaction (only cases with two ISWs after the interaction are considered here).

Figure 13. Energy changes after ISW interactions as a function of the difference in propagation speed of the two initial ISWs. (a–c) Red (circles) is for stratification $S_3$ with small wave with ${\rm APE} =0.005$ and larger waves with ${\rm APE} = 0.008$, 0.01, 0.015, 0.02 and 0.04; orange (squares) is for stratification $S_3$ with small wave with ${\rm APE} =0.01$ and larger waves with ${\rm APE} = 0.015$, 0.02, 0.04 and 0.06. (d–f) Blue (asterisks) is for stratification $S_1$ with a slow negative KdV polarity wave with ${\rm APE} =0.03$ and fast KdV waves with ${\rm APE} = 0.0035$, 0.003, 0.004, 0.005, 0.01, 0.015 and 0.02; green (triangles) is for stratification $S_5$ with a slow anti-KdV polarity wave with ${\rm APE} =0.02$ and fast KdV waves with ${\rm APE}= 0.006$, 0.008, 0.01, 0.12 and 0.02. (a,d) Fraction of initial energy in the slow ISW after the interaction. (b,e) Fraction of initial energy in the fast wave after the interaction. (c,f) Change in total energy in the two ISWs after the interaction as a fraction of the initial total energy (solid curves) and the adjusted values taking into account energy loss due to numerical dissipation (dotted curves).

The effects of numerical dissipation are accounted for by assuming numerical dissipation removes energy at a constant rate throughout the simulation as suggested by figure 12, that is, we compute the rate $r$ at which energy is lost before the interaction, which has a typical value of $-3\times 10^{-7}$, and use a corrected final energy value $\tilde {E}_a = E_a - r T_a$, where $T_a$ is the time at which $E_a$ was calculated. We also assume that energy is lost from each ISW via numerical dissipation at the same rate. The resulting values, corrected for numerical dissipation, are indicated by the dotted curves in figure 13. Note that accounting for energy dissipation increases the energy increase of the fast wave while decreasing the energy decrease of the slow waves. The change in the total energy provides an estimate of the energy lost to the trailing waves. It is as large as 10 % of the initial wave energy in series $S_1$.

Only a few simulations with two ISWs of opposite polarity were done with the initial fast trailing wave having negative KdV polarity (it being far easier to find DJL solutions with faster KdV polarity waves). The interactions all resulted in small changes in the amplitudes of the two ISWs, as in figure 7(c). When the fast trailing wave has KdV polarity the interactions result in large changes when the difference in propagation speed is small (figure 13). Note that $C_3$ and $C_4$ (see figure 7c,d), which have the KdV and negative KdV polarity waves in reversed order, have similar $\Delta c$ but the results of the interactions are very different.

This work has implications for the use of the Gardner equation with positive coefficient to model ‘soliton-turbulence’ in the context of internal waves (e.g. Didenkulova 2019), because it demonstrates that for internal waves the Gardner equation does a poor job of modelling interactions of waves of opposite polarity. Hence, better simplified models are needed; however, it may be difficult to replace full two-dimensional numerical simulations. The interactions of two mode-two internal solitary-like waves, one on each of two well-separated pycnoclines ($H/\lambda \gg 1$, where $H$ is the pycnocline separation distance and $\lambda$ is a typical wavelength), have been studied in the laboratory (Weidman & Johnson 1982) and theoretically (Liu, Pereira & Ko 1982; Nitsche et al. 2010; Weidman, Nitsche & Howard 2011). (Solitary-like because these waves are not exact ISWs. They are in general accompanied by a small mode-one dispersive wave train whose phase speed is equal to the propagation speed of the mode-two solitary-like wave. Because their group velocity is smaller that their phase speed this means there is a continual drain of energy from the mode-two solitary-like wave which therefore slowly decays in amplitude (Akylas & Grimshaw 1992).) The waves were modelled using coupled intermediate-depth equations (Joseph 1977; Kubota, Ko & Dobbs 1978). Here a reverse energy transfer occurs with the leading faster-propagating solitary-like wave transferring energy to the trailing solitary-like wave. In this process the trailing wave grows, speeds up and passes the leading wave which has decreased in amplitude. This ‘leap-frogging’ process then repeats. Each interaction results in the loss of energy to trailing small-amplitude dispersive waves and eventually the leap-frogging stops. Gear & Grimshaw (1984) considered the interaction of internal solitary waves of different modes in cases where the total water depth is small compared with the wavelength of the waves. In this case KdV-type equations are obtained. For cases where the linear long-wave speeds of the two modes differ by an $O(1)$ amount the coupling is weak; however, when the linear long-wave speeds differ by a small amount the evolution is governed by a set of nonlinearly and dispersively coupled KdV-type equations. Gear & Grimshaw (1984) demonstrated the existence of stratifications for which the linear long-wave speeds of the two modes are similar. These stratifications consisted of upper and lower layers of constant buoyancy frequency with the middle layer unstratified. In the simulations discussed here $H/\lambda$ varies between about 0.4 and 1 (based on the half-width of the wave-induced surface currents). None of the existing theories apply for the situation considered here, namely the interaction of two ISWs of the same mode with opposite polarities; however, as the waves of depression and elevation have very different vertical structures there is some similarity to the cases considered by Gear & Grimshaw (1984) suggesting a future avenue of exploration. It is important to emphasize, however, that the mode-one solitary waves of elevation and depression have vastly different vertical structures. In particular, the negative KdV polarity wave has a vertical structure that is not a small correction to the linear mode-one modal function and hence it is not a weakly nonlinear wave.