3.1 The pseudomomentum of a linearised disturbance

The inviscid linearised equations for perturbations $\check{\unicode[STIX]{x1D713}}$ , $\check{\unicode[STIX]{x1D701}}$ and $\check{\unicode[STIX]{x1D70C}}$ about the base flow profiles $U(y)$ and $\bar{\unicode[STIX]{x1D70C}}(y)$ are

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \check{\unicode[STIX]{x1D701}}_{t}+U\check{\unicode[STIX]{x1D701}}_{x}-U^{\prime \prime }\check{\unicode[STIX]{x1D713}}_{x}+J\check{\unicode[STIX]{x1D70C}}_{x}=0, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \check{\unicode[STIX]{x1D70C}}_{t}+U\check{\unicode[STIX]{x1D70C}}_{x}+\bar{\unicode[STIX]{x1D70C}}^{\prime }\check{\unicode[STIX]{x1D713}}_{x}=0. & \displaystyle\end{eqnarray}$$

Formulating the combination

(3.3)

where $\langle \cdots \rangle$ denotes an integral over the area of the spatial domain and $c_{\ast }$ is an arbitrary constant velocity, we arrive at a conservation law:

(3.4) $$\begin{eqnarray}\frac{\text{d}{\mathcal{E}}}{\text{d}t}=0,\quad {\mathcal{E}}=\left\langle \frac{1}{2}(\check{\unicode[STIX]{x1D713}}_{x}^{2}+\check{\unicode[STIX]{x1D713}}_{y}^{2})+\frac{J\check{\unicode[STIX]{x1D70C}}^{2}}{2\bar{\unicode[STIX]{x1D70C}}^{\prime }}-(U-c_{\ast })U^{\prime \prime }\frac{\check{\unicode[STIX]{x1D70C}}^{2}}{2\bar{\unicode[STIX]{x1D70C}}^{\prime 2}}-(U-c_{\ast })\frac{\check{\unicode[STIX]{x1D70C}}\check{\unicode[STIX]{x1D701}}}{\bar{\unicode[STIX]{x1D70C}}^{\prime }}\right\rangle\end{eqnarray}$$

(cf. (2.28) of Abarbanel et al. Reference Abarbanel, Holm, Marsden and Ratiu1986).

For normal modes, we have $\check{\unicode[STIX]{x1D713}}(x,y,t)=\hat{\unicode[STIX]{x1D713}}(y)\text{e}^{\text{i}k(x-ct)}+\text{c.c}.$ , where $k$ is the wavenumber and $c=c_{r}+ic_{i}$ is the complex phase speed. We may further exploit the normal-mode equations, $\hat{\unicode[STIX]{x1D701}}=(U^{\prime \prime }\hat{\unicode[STIX]{x1D713}}-J\hat{\unicode[STIX]{x1D70C}})/(U-c)$ and $\hat{\unicode[STIX]{x1D70C}}=-\bar{\unicode[STIX]{x1D70C}}^{\prime }\hat{\unicode[STIX]{x1D713}}/(U-c)$ , to arrive at the conserved quantity,

(3.5) $$\begin{eqnarray}M=\frac{{\mathcal{E}}}{(c_{r}-c_{\ast })}=\left\langle \frac{1}{2}\check{\unicode[STIX]{x1D713}}^{2}\left[\frac{U^{\prime \prime }}{|U-c|^{2}}+\frac{2J(U-c_{r})\bar{\unicode[STIX]{x1D70C}}^{\prime }}{|U-c|^{4}}\right]\right\rangle =\langle {\mathcal{M}}_{V}+{\mathcal{M}}_{B}\rangle ,\end{eqnarray}$$

with

(3.6a,b ) $$\begin{eqnarray}{\mathcal{M}}_{V}=\frac{1}{2}\frac{U^{\prime \prime }\check{\unicode[STIX]{x1D713}}^{2}}{|U-c|^{2}}\quad \text{and}\quad {\mathcal{M}}_{B}=\frac{J(U-c_{r})\bar{\unicode[STIX]{x1D70C}}^{\prime }\check{\unicode[STIX]{x1D713}}^{2}}{|U-c|^{4}}.\end{eqnarray}$$

The quantity $M$ represents the net pseudomomentum of a linear disturbance; the time rate of change of its density ${\mathcal{M}}_{V}+{\mathcal{M}}_{B}$ represents the local acceleration of the mean flow induced by the disturbance (Bühler Reference Bühler2014). Moreover, because $M$ is conserved, the exponential growth of any linear perturbation with normal-mode form is forbidden unless $M=0$ for that mode, as noted originally by Synge (Reference Synge1933) and discussed by Miles (Reference Miles1961) and Howard (Reference Howard1961).

For unstratified flow with $\bar{\unicode[STIX]{x1D70C}}^{\prime }=0$ , the condition $M\equiv \langle {\mathcal{M}}_{V}\rangle =0$ leads to the inflexion-point criterion. With stratification, the fact that the contributions to $M$ must cancel for an unstable normal mode constrains the types of wave interactions that can lead to instability: gravity waves localised at a density interface have a sign for the pseudomomentum density ${\mathcal{M}}_{B}$ given by $-\text{sgn}(U-c_{r})$ (given $\bar{\unicode[STIX]{x1D70C}}^{\prime }<0$ for stable stratification). Thus, waves riding on two density interfaces with phase speeds between the flow speeds of the interfaces have opposite signs of ${\mathcal{M}}_{B}$ , and set the scene for TCI. In particular, if $U^{\prime \prime }=0$ , TCI modes must have cancelling interfacial contributions to the pseudomomentum to give $M\equiv \langle {\mathcal{M}}_{B}\rangle =0$ .

Likewise, HWI can only arise when a wave on a density interface, with a pseudomomentum contribution through $\langle {\mathcal{M}}_{B}\rangle$ , interacts with a wave on a vorticity interface, with an oppositely signed contribution from $\langle {\mathcal{M}}_{V}\rangle$ . This demands that $U^{\prime \prime }(U-c_{r})>0$ , a Fjortoft-like condition that rules out instability for certain flow profiles which nevertheless allow wave resonances (cf. Carpenter et al. Reference Carpenter, Balmforth and Lawrence2010). For such stable interfaces, we also see that the criterion $M=0$ is stronger than the Miles–Howard criterion, as the local Richardson number vanishes throughout the bulk of the flow.

Some examples of the pseudomomentum contributions of unstable modes are shown in figure 1. These modes are computed from the linear eigenvalue problem for the basic profiles:

(3.7) $$\begin{eqnarray}\displaystyle \text{KHI}\quad U(y) & = & \displaystyle \tanh (10y),\nonumber\\ \displaystyle \bar{\unicode[STIX]{x1D70C}}(y) & = & \displaystyle -\text{tanh}(30y),\quad J=0.05,\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle \text{HWI}\quad U(y) & = & \displaystyle 1+{\textstyle \frac{4}{135}}(30y-30-\log \{\cosh [30(y-1/8)]/\cosh (105/4)\}),\nonumber\\ \displaystyle \bar{\unicode[STIX]{x1D70C}}(y) & = & \displaystyle -\text{tanh}[30(y+1/8)],\quad J=0.2,\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle \text{TCI}\quad U(y) & = & \displaystyle y,\qquad \nonumber\\ \displaystyle \bar{\unicode[STIX]{x1D70C}}(y) & = & \displaystyle -{\textstyle \frac{1}{2}}\{\tanh [30(y-1/4)]+\tanh [30(y+1/4)]\},\quad J=0.14.\end{eqnarray}$$

These profiles have relatively sharp but smooth interfaces in the background density and vorticity and represent ‘typical’ examples that lead to KHI, HWI and TCI (respectively). Figure 1 also displays the base profiles and the vorticity and density disturbances associated with their modal eigenfunctions. For the KHI case, ${\mathcal{M}}_{B}$ is negligible and there are cancelling contributions through ${\mathcal{M}}_{V}$ from the two vorticity interfaces. For HWI, there are equal and opposite contributions to $M=0$ from ${\mathcal{M}}_{V}$ and ${\mathcal{M}}_{B}$ that are centred at the vorticity and density interfaces, respectively. For TCI, ${\mathcal{M}}_{V}\equiv 0$ and the contributions of each density interface cancel.

Figure 1. Columns from left to right: sample background profiles $U(y)$ (blue solid) and $\bar{\unicode[STIX]{x1D70C}}(y)$ (red dot-dashed) for (a) KHI, (b) HWI and (c) TCI. Absolute value of vorticity (blue) and density (red dot-dashed) for the associated growing mode. ${\mathcal{M}}_{V}$ (blue solid, when non-zero) and ${\mathcal{M}}_{B}$ (red dot-dashed) for the background profiles and the associated modes.

3.2 Pseudomomentum-based mode characterisation

For arbitrary spatial profiles of vorticity and density, one can exploit the decomposition into ${\mathcal{M}}_{V}$ and ${\mathcal{M}}_{B}$ to classify the instabilities: first, the relative sizes of $|{\mathcal{M}}_{B}|$ and $|{\mathcal{M}}_{V}|$ indicate whether an interfacial region acts like a density or vorticity interface. Given this ‘flavour’ of the interface, one can then pick out the two largest cancelling contributions to $M$ , and thereby classify the instability.

To illustrate this scheme, we solve the linear stability problem analytically for two piecewise linear profiles. In the first of these examples, a transition occurs from KHI to TCI (Taylor Reference Taylor1931) as the flow parameters are varied, and the distinction between the two becomes ambiguous. Similarly, the second example exhibits a transition between KHI and HWI (Holmboe Reference Holmboe1962), and is the flow adopted by Carpenter et al. (Reference Carpenter, Balmforth and Lawrence2010) in applying their diagnostic for distinguishing between these two instabilities.

The first example is Taylor’s three-layer flow profile,

(3.10a,b ) $$\begin{eqnarray}U(y)=\left\{\begin{array}{@{}ll@{}}-1,\quad \quad & \text{if }y\leqslant -1\\ y,\quad \quad & \text{if }-1<y<1\\ 1,\quad \quad & \text{if }y\geqslant 1\end{array}\right.\quad \text{and}\quad \bar{\unicode[STIX]{x1D70C}}(y)=\left\{\begin{array}{@{}ll@{}}1,\quad \quad & \text{if }y\leqslant -1,\\ 0,\quad \quad & \text{if }-1<y\leqslant 1,\\ -1,\quad \quad & \text{if }y>1,\end{array}\right.\end{eqnarray}$$

in which there are two, coincident vorticity and density interfaces. The dispersion relation for normal modes is

(3.11) $$\begin{eqnarray}\displaystyle & & \displaystyle c^{4}+c^{2}\left[\frac{e^{-4k}-(2k-1)^{2}}{4k^{2}}-1-\frac{J}{k}\right]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{J^{2}}{4k^{2}}\left(1-e^{-4k}\right)-\frac{J}{k}\left(\frac{2k-1+e^{-4k}}{2k}\right)-\left(\frac{e^{-4k}-(2k-1)^{2}}{4k^{2}}\right)=0.\end{eqnarray}$$

The pseudomomentum reduces to

(3.12) $$\begin{eqnarray}\displaystyle M & = & \displaystyle L_{x}\left[|\hat{\unicode[STIX]{x1D713}}(1)|^{2}\left(\frac{1}{|1-c|^{2}}-\frac{2J}{|1-c|^{4}}\right)+|\hat{\unicode[STIX]{x1D713}}(-1)|^{2}\left(-\frac{1}{|1+c|^{2}}+\frac{2J}{|1+c|^{4}}\right)\right]\qquad\end{eqnarray}$$

(3.13) $$\begin{eqnarray}\displaystyle & \equiv & \displaystyle L_{x}\left[|\hat{\unicode[STIX]{x1D713}}(1)|^{2}\left(M_{V}^{+}-M_{B}^{+}\right)+|\hat{\unicode[STIX]{x1D713}}(-1)|^{2}\left(-M_{V}^{-}+M_{B}^{-}\right)\right],\end{eqnarray}$$

where $M_{V}^{\pm },M_{B}^{\pm }>0$ are the magnitudes of the contributions from the upper and lower vorticity and density interfaces, respectively. The unstable normal modes have $|\hat{\unicode[STIX]{x1D713}}(1)|^{2}=|\hat{\unicode[STIX]{x1D713}}(-1)|^{2}$ , and $c_{r}=0$ so that $M_{V}^{+}=M_{V}^{-}\equiv M_{V}$ and $M_{B}^{+}=M_{B}^{-}\equiv M_{B}$ , allowing us to define the relative contribution of the vorticity interfaces as

(3.14) $$\begin{eqnarray}R_{1}=\frac{M_{V}}{M_{V}+M_{B}}=\frac{1+c_{i}^{2}}{1+c_{i}^{2}+2J}.\end{eqnarray}$$

Then, $R_{1}=1$ corresponds to pure KHI in the sense that only the vorticity interfaces contribute to the dynamics. Likewise, $R_{1}=0$ corresponds to pure TCI dynamics, given that only the density interfaces then contribute. Intermediate values of $R_{1}$ indicate that all four interfaces play a role in the dynamics. Figure 2(a) shows contours of constant growth rate $kc_{i}$ on the $k-J$ plane, superposed on a density plot of $R_{1}$ . For $J=0$ , the density variations play no role, implying pure KHI dynamics, and we see that $R_{1}\rightarrow 1$ . For $J\gg 1$ , $R_{1}\rightarrow 0$ , indicating an approach to pure TCI. Between these two limits there is a smooth transition in which all four interfaces play varying roles in the dynamics; the contour $R_{1}=1/2$ , which occurs when $J\approx 0.6$ , suggests a convenient distinction between KHI and TCI. However, we emphasise that this distinction is a relative one, and does not imply that only the density interfaces are active when $R_{1}<1/2$ , or that the instability is a pure vorticity interface interaction when $R_{1}>1/2$ .

Figure 2. (a) Contours of constant growth rate of the unstable mode (black lines) for the linear profile in (3.10) showing the transition between KHI and TCI. The colour shading shows a density plot of $R_{1}=M_{V}/(M_{V}+M_{B})$ . (b) Contours of constant growth rate of the rightward propagating mode (black lines) for the piecewise linear profile in (3.15) showing the transition between KHI and HWI. The colour shading shows a density plot of $R_{2}=1-|M_{B}/M_{V}^{+}|$ , where $M_{V}^{+}$ is the contribution to $M_{V}$ due to the upper vorticity interface alone.

The second example is Holmboe’s flow profile,

(3.15) $$\begin{eqnarray}\displaystyle U(y)=\left\{\begin{array}{@{}ll@{}}-1,\quad \quad & \text{if }y\leqslant -1\\ y,\quad \quad & \text{if }-1<y<1\\ 1,\quad \quad & \text{if }y\geqslant 1\end{array}\right.\quad \text{and}\quad \bar{\unicode[STIX]{x1D70C}}(y)=\left\{\begin{array}{@{}ll@{}}1,\quad \quad & \text{if }y\leqslant 0,\\ 0,\quad \quad & \text{if }y>0,\end{array}\right. & & \displaystyle\end{eqnarray}$$

with dispersion relation,

(3.16) $$\begin{eqnarray}c^{4}+c^{2}\left[\frac{e^{-4k}-(2k-1)^{2}}{4k^{2}}-\frac{J}{k}\right]+\frac{J}{k}\left[\frac{e^{-2k}+(2k-1)}{2k}\right]^{2}=0.\end{eqnarray}$$

The pseudomomentum is given by

(3.17) $$\begin{eqnarray}M=L_{x}\left[\frac{|\hat{\unicode[STIX]{x1D713}}(-1)|^{2}}{2|1+c|^{2}}-\frac{|\hat{\unicode[STIX]{x1D713}}(1)|^{2}}{2|1-c|^{2}}+\frac{2Jc_{r}|\unicode[STIX]{x1D713}(0)|^{2}}{|c|^{4}}\right]\equiv M_{V}^{-}+M_{V}^{+}+M_{B}.\end{eqnarray}$$

Figure 2(b) shows contours of constant growth rate $kc_{i}$ on the $k$ – $J$ plane for the rightward propagating unstable mode. These contours are superposed on a density plot of $R_{2}=1-|M_{B}/M_{V}^{+}|$ . For large wavenumber $k$ , the modal eigenfunctions decay exponentially away from the interfaces, leading to mode interactions that are characterised by cancelling pseudomomentum contributions from just two of $M_{V}^{\pm }$ and $M_{B}$ . For the rightward propagating modes, $R_{2}\rightarrow 0$ indicates that the instability is generated by the interaction of the upper vorticity interface with the density interface, and has the character of HWI. In view of the symmetry of the profiles, simultaneously there is a leftward travelling HWI with $1-|M_{B}/M_{V}^{-}|\rightarrow 0$ . On the other hand, when $R_{2}\rightarrow 1$ , the density interface cannot contribute to the interaction and so the mode must take a KHI character. Indeed, for these flow profiles, any mode with zero phase speed $c_{r}=0$ has this feature, since $M_{B}$ then vanishes identically. Intermediate values of $R_{2}$ indicate balancing contributions to $M=0$ from all three of $M_{V}^{\pm }$ and $M_{B}$ . The contour $R_{2}=1/2$ again suggests a convenient distinction between KHI-like and HWI-like interactions. As shown in figure 2(b), the $R_{2}$ diagnostic cleanly identifies the one-sided rightward travelling HWI at higher wavenumber and bulk Richardson number. Moreover, there is an area with $1/2<R_{2}<1$ just above the transition from stationary to propagating waves (highlighted by the sharp changes in the growth-rate contours). Here, the lower vorticity interface plays a role in the mode dynamics, and these propagating instabilities inherit a KHI-like character. This result complements the findings of Carpenter et al. (Reference Carpenter, Balmforth and Lawrence2010), but the extent of the region that they identify as KHI is slightly different.

4.1 Numerical methods and initial conditions

To simulate TCI, we consider initial-value problems in which small perturbations are added to the base profiles,

(4.1) $$\begin{eqnarray}U(y)=y,\quad \bar{\unicode[STIX]{x1D70C}}(y)=-{\textstyle \frac{1}{2}}\left\{\tanh [30(y-1/4)]+\tanh [30(y+1/4)]\right\}.\end{eqnarray}$$

The perturbations take the form of the unstable mode of the diffusive linear problem associated with these profiles, with an amplitude corresponding to a perturbation energy of $E(0)\approx 10^{-7}$ . We solve (2.1)–(2.3) in the domain $y\in [-1,1]$ and $x\in [0,L_{x}=2\unicode[STIX]{x03C0}/k]$ , with the horizontal wavenumber $k$ as a parameter.

Table 1 lists the parameter values of the simulations conducted. We consider four different pairings of the Reynolds and Prandtl numbers, $(Re,Pr)=(300\,000,0.6)$ , (180 000,1), (60 000,3) and (20 000,10), chosen so that the Péclet number is nearly the same ( $Pe=180,000$ for the first three pairings and 200 000 for the last). The simulations form six groups in the table, so that each group has the same bulk Richardson number $J$ and horizontal wavenumber $k$ , and consists of a simulation from each of the four Reynolds/Prandtl pairings. In addition, group ‘X’ represents an additional set of two TCI simulations at small Reynolds number (although still high Péclet number) at the two extremes of wavelength, plus a pair of representative KHI and HWI computations at comparable Péclet numbers. For the KHI and HWI cases, the initial base states are given by (3.7) and (3.8), respectively.

Table 1. Parameter values of simulations. Group X includes additional TCI simulations at low $Re$ and representative KHI and HWI simulations.

Figure 3. Stability boundaries in the $(k,J)$ plane for TCI with background profiles and boundary conditions given by an infinitely sharp interface version of (4.1), along with values of $(k,J)$ for the six groups of simulations listed in table 1. Colours indicate the observed primary dynamics of the saturated instability.

The locations of the simulations on the $(k,J)$ -plane are shown in figure 3. The comparison of the small wavenumber, small bulk Richardson number simulations of group 1 with the high wavenumber, high bulk Richardson number of groups 5 and 6 offers insight into how the mechanism of saturation varies across the strip of instability, which is also drawn in figure 3 (taking the interfaces to be sharp and modes to be inviscid). Groups 2–4, which share the wavenumber $k=4/3$ but have bulk Richardson numbers that span those of group 1 and group 5, allow one to disentangle any combined effect of varying $k$ and $J$ . Note that the values of $k$ were chosen so that simulations could be compared at different pairings of $J$ and $k$ , whilst trying to ensure that there was at most one unstable mode in the domain to avoid any complication due to mode interactions. The choices do not correspond to the fastest growing mode over all wavenumbers at each fixed $J$ , which is a common strategy for constraining parameters in numerical studies. However, the nonlinear dynamics we find for each $J$ appears to be insensitive to the precise choice for $k$ , at least in two dimensions (it is possible that in three dimensions, the secondary instabilities are more sensitive to the choice of $k$ ).

The simulations were conducted using Diablo, a parallel Fortran-based numerical integration program developed by Bewley and Taylor (Taylor Reference Taylor2008), which uses a second-order finite difference scheme in the $y$ -direction and a de-aliased Fourier decomposition in the $x$ -direction. Time-stepping is implemented by a combined implicit–explicit third-order Runge–Kutta–Wray Crank–Nicolson scheme. The grid resolution for the geometry of groups 1–4 is $2048\times 2048$ and for groups 5 and 6 it is $1024\times 2048$ . Movies illustrating the simulations are provided as supplementary information available at https://doi.org/10.1017/jfm.2018.827.

4.2 Diagnostics

To compare the simulations with one another, we make use of a number of global measures: first, the average perturbation energy $E(t)$ (potential plus kinetic) is defined as

(4.2) $$\begin{eqnarray}E(t)=\frac{k}{8\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}/k}\int _{-1}^{1}\left\{[u-U(y)]^{2}+v^{2}+J[\unicode[STIX]{x1D70C}-\bar{\unicode[STIX]{x1D70C}}(y)]^{2}\right\}\,\text{d}y\,\text{d}x\end{eqnarray}$$

(an expansion to second order in mode amplitude connects this quantity to the pseudo-energy ${\mathcal{E}}$ defined in § 3.1; Bühler (Reference Bühler2014)). Second, following Peltier & Caulfield (Reference Peltier and Caulfield2003), a monotonic adiabatic rearrangement of the current density field $\unicode[STIX]{x1D70C}(x,y,t)$ furnishes the background stratification $\unicode[STIX]{x1D70C}_{\ast }(y,t)$ , from which we obtain the background potential energy $P_{B}(t)$ ,

(4.3) $$\begin{eqnarray}P_{B}(t)=-\frac{1}{2}\int _{-1}^{1}Jy\unicode[STIX]{x1D70C}_{\ast }\,\text{d}y.\end{eqnarray}$$

The instantaneous mixing rate $m(t)$ can then found from the time rate of change of $P_{B}(t)$ , correcting for the natural diffusive increase of $P_{B}$ due to the Péclet number:

(4.4) $$\begin{eqnarray}m(t)=\frac{\text{d}P_{B}}{\text{d}t}-\frac{J}{Pe}.\end{eqnarray}$$

Third, the mixing efficiency $\unicode[STIX]{x1D702}(t)$ is given by

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D702}(t)=\frac{m(t)}{m(t)+D(t)},\end{eqnarray}$$

where the dissipation rate is

(4.6) $$\begin{eqnarray}D(t)=\frac{k}{4\unicode[STIX]{x03C0}Re}\int _{0}^{2\unicode[STIX]{x03C0}/k}\int _{-1}^{1}\left[\left(\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}x}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}y}\right)^{2}\right]\,\text{d}y\,\text{d}x.\end{eqnarray}$$

Figure 4. (a) Instantaneous mixing efficiency $\unicode[STIX]{x1D702}(t)$ and (b) instantaneous mixing rate $m(t)$ for representative KHI (blue solid) and HWI (red dashed, for which we plot $10m(t)$ ) simulations.

It should be emphasised that for the two-dimensional simulations presented here, the mixing efficiencies are expected to be significantly higher than typical values for three-dimensional flow; two-dimensional flows dissipate far less than their three-dimensional counterparts (e.g. Peltier & Caulfield Reference Peltier and Caulfield2003). Nevertheless, we find the mixing efficiency to be a useful indicator to distinguish between our simulations in addition to the perturbation energy density $E(t)$ since it is possible to have energetic flows which do little mixing and vice versa. The final scene is set for our discussion of the mixing character of TCI by figure 4, which shows $\unicode[STIX]{x1D702}(t)$ and $m(t)$ for the representative KHI and HWI computations of table 1, and illustrates the vigorous, but short-lived mixing of KHI and the lower, but more long-lived action of HWI.

4.3 Small $k$ , small $J$ dynamics (group 1)

Figure 5. Time evolution of (a) the perturbation energy $E(t)$ , (b) mixing efficiency $\unicode[STIX]{x1D702}(t)$ and (c) mixing rate $m(t)$ for the simulations in group 1, with $(Re,Pr)=(300\,000,0.6)$ (blue solid), (180 000, 1) (red dashed), (60 000, 3) (yellow dot-dashed) and (20 000, 10) (purple dotted).

4.3.1 Energetic characteristics

Figure 5 shows the time evolution of $E(t)$ , $\unicode[STIX]{x1D702}(t)$ and $m(t)$ for group 1, for each of the four $(Re,Pr)$ pairings. The initial evolution of the perturbation energy $E(t)$ collapses onto the same curve during the initial linear growth and peaks around $t=100$ for all four simulations. Subsequently, $E(t)$ settles into a decaying oscillation for the $(Re,Pr)=(300\,000,0.6)$ and (180 000, 1) simulations; the (60 000, 3) simulation levels off somewhat, whereas the (20 000, 10) simulation transitions to a much higher amplitude, chaotic state.

During the linear growth, the mixing efficiency $\unicode[STIX]{x1D702}(t)$ and rate $m(t)$ remain small, but increase substantially after the first maximum in $E(t)$ . The four simulations reach different peak efficiencies, mainly because of the enhanced dissipation at lower $Re$ ; the time evolution of the mixing rate $m(t)$ through the first peak is more comparable between the simulations. Evidently, the simulations with higher Prandtl numbers ( $Pr=3$ and 10) experience higher levels of sustained mixing than in the lower Prandtl number cases ( $Pr=0.6$ and 1).

Figure 6. Snapshots of the density $\unicode[STIX]{x1D70C}(x,y,t)$ (left) and vorticity $\unicode[STIX]{x1D701}(x,y,t)+1$ (right) for four different times and parameter values in group 1: (a)  $(Re,Pr)=(300\,000,0.6)$ at $t=110$ , (b)  $(Re,Pr)=(180\,000,1)$ at $t=141$ , (c)  $(Re,Pr)=(20\,000,10)$ at $t=188$ and (d)  $(Re,Pr)=(20\,000,10)$ at $t=375$ . Colour bars show density (left) and vorticity (right).

4.3.3 Billow cascade as a secondary TCI

Figure 7 shows the gradient Richardson number $-J\unicode[STIX]{x1D70C}_{y}/u_{y}^{2}$ for vertical profiles of $\unicode[STIX]{x1D70C}$ and $u$ , horizontally averaged over a region $2.4<x<4.2$ outside the primary billow in the top row of images of figure 6. Also shown is the corresponding initial condition. The region outside the primary billow retains two relatively strong density interfaces, but their vertical separation has changed, which suggests that the region may still be susceptible to TCI. Indeed, one can use the horizontal averages of $\unicode[STIX]{x1D70C}$ and $u$ over this region, together with an adjusted wavenumber of $k=3$ consistent with its rough horizontal extent, to conduct a second linear stability analysis. This calculation reveals a TCI with a growth rate of approximately $0.1$ . By contrast, the initial basic state is not unstable at $k=3$ . Thus, the nonlinear distortion of the mean profile by the primary instability generates a secondary TCI, and occurs because the original instability moves the density interfaces closer together to reduce the local length scale. Moreover, the density interfaces remain sharp through the combined action of the relatively high Péclet number limiting the diffusive effects and a favourable straining flow induced by the primary billow. Evidently, this mechanism continues to operate at higher Prandtl numbers, yielding further secondary instabilities and generating the cascade described above.

Figure 7. Gradient Richardson number $-J\unicode[STIX]{x1D70C}_{y}/u_{y}^{2}$ for the horizontally averaged buoyancy and velocity profiles in the region between the billow $2.4<x<4.2$ for the $(Re,Pr)=(300\,000,0.6)$ group 1 simulation at $t=110$ (blue), along with the gradient Richardson number for the initial condition (grey, dotted). Also shown is the value $Ri=0.25$ (black dashed).

4.4 Large $k$ , large $J$ dynamics (groups 5–6)

4.4.1 Energetic characteristics

Figure 8. Time evolution of (a,d) the perturbation energy $E(t)$ , (b,e) mixing efficiency $\unicode[STIX]{x1D702}(t)$ and (c,f) mixing rate $m(t)$ for the simulations in (a–c) group 5 and (d–f) group 6, with $(Re,Pr)=(300\,000,0.6)$ (blue thick solid), (180 000, 1) (red dashed), (60 000, 3) (yellow dot-dashed) and (20 000, 10) (purple thin solid).

Figure 8 shows the time evolution of $E(t)$ , $\unicode[STIX]{x1D702}(t)$ and $m(t)$ for the $(k,J)=(2,0.3)$ simulations of group 5 and the $(k,J)=(4,0.5)$ simulations of group 6. Once more, the simulations show the same initial linear growth of the perturbation energy $E(t)$ . The peak in $E(t)$ is again followed by decaying oscillations about a slowly diffusing nonlinear state at lower $Pr$ , or sustained higher-frequency oscillations symptomatic of emergent parasitic Holmboe waves at larger $Pr$ . As for the group 1 simulations, the mixing efficiency $\unicode[STIX]{x1D702}(t)$ and rate $m(t)$ increase to significant levels only after the end of the phase of linear growth. The maximum values reached, however, are substantially higher and are comparable to those observed for two-dimensional KHI simulations (see Peltier & Caulfield Reference Peltier and Caulfield2003, and figure 4). The implied enhanced mixing events have only previously been associated with KHI; we uncover the origin of this presently.

4.4.2 Billow break-up

Figure 9. Snapshots of the density $\unicode[STIX]{x1D70C}(x,y,t)$ (left) and vorticity $\unicode[STIX]{x1D701}(x,y,t)+1$ (right) for six different times and parameter values in group 5:  $(Re,Pr)=(60\,000,3)$ at (a)  $t=56$ , (b) $t=66$ , (c)  $t=68$ and (d)  $t=86$ ; (e)  $(Re,Pr)=(20\,000,10)$ at $t=240$ , (f)  $(Re,Pr)=(300\,000,0.6)$ at $t=250$ .

Figure 9 shows snapshots of the density $\unicode[STIX]{x1D70C}(x,y,t)$ and vorticity $\unicode[STIX]{x1D701}(x,y,t)$ at a number of key times in the simulations of group 5. Figure 9(a,b) shows the $(Re,Pr)=(60\,000,3)$ simulation; once again, during the linear growth of the TCI, cusp-like deformations appear on each density interface and are swept sideways by the mean flow. This time, however, the resulting nonlinear travelling waves do not collide and merge to form a billow. Instead, they retain their identity and generate more baroclinic vorticity due to the higher bulk Richardson number; these alternative primary nonlinear structures correspond to the cuspy wave observed by Eaves & Caulfield (Reference Eaves and Caulfield2017). Moreover, the waves interact strongly to draw out large-amplitude filaments of vorticity as time progresses. These filaments subsequently roll up through localised KHI-type secondary instabilities, creating vortices throughout the middle density layer and precipitating a vigorous mixing event, as illustrated by figure 9(c,d) from the $(Re,Pr)=(60\,000,3)$ simulation. Eventually, the finer-scale vortices viscously decay, coarse graining a residual nonlinear structure (see figure 9 f). At long times, a parasitic Holmboe wave emerges for $Pr=10$ (see figure 9 e).

Figure 10. Snapshots of the density $\unicode[STIX]{x1D70C}(x,y,t)$ (left) and vorticity $\unicode[STIX]{x1D701}(x,y,t)+1$ (right) for nine different times and parameter values in group 6: $(Re,Pr)=(300\,000,0.6)$ at (a)  $t=55$ , (b) $t=62$ , (c)  $t=74$ , (d)  $t=82$ , (e)  $t=88$ , (f)  $t=95$ , (g)  $t=124$ and (h)  $t=250$ ; (i)  $(Re,Pr)=(20\,000,10)$ at $t=200$ .

Figure 10 collects snapshots of $\unicode[STIX]{x1D70C}(x,y,t)$ and $\unicode[STIX]{x1D701}(x,y,t)$ from simulations of group 6. The horizontal extent of the domain is now sufficiently short that the cuspy nonlinear travelling waves no longer lock into place when they approach, but pass by one another to create complex vorticity dynamics in the middle density layer; sharp filaments of vorticity once more appear and subsequently roll up due to secondary KHI (see figure 10 d–f). This fills the central density layer with finer-scale vortical structures that again viscously decay to a smoother state but for the emergence of parasitic Holmboe waves at $Pr=10$ , as seen in figure 10(g–i).

Figure 11. Time evolution of (a,d,g) the perturbation energy $E(t)$ , (b,e,h) mixing efficiency $\unicode[STIX]{x1D702}(t)$ and (c,f,i) mixing rate $m(t)$ for the simulations in (a–c) group 2, (d–f) group 3 and (g–i) group 4, with $(Re,Pr)=(300\,000,0.6)$ (blue thick solid), (180 000, 1) (red dashed), (60 000, 3) (yellow dot-dashed) and (20 000, 10) (purple thin solid).

The roll-up of baroclinically generated vorticity filaments via secondary KHI rationalises the heightened mixing efficiencies and rates in the simulations of groups 5 and 6; the primary TCI itself generates very little mixing. Because the secondary roll-ups span the entire middle density layer, the mixing efficiency reaches values more typically associated with a primary KHI. The baroclinic generation of vorticity arises through the term $J\unicode[STIX]{x1D70C}_{x}$ in (2.3), and so can be associated with either large values of $J$ or small horizontal length scales. The simulations of groups 5 and 6 possess both attributes in comparison to group 1, and so it remains unclear whether the secondary behaviour is due to a larger $J$ or smaller horizontal wavenumber $k$ . The simulations of groups 2, 3 and 4, described below, indicate that the bulk Richardson number $J$ is chiefly responsible.

Unlike the cascade of Taylor–Caulfield billows for weaker stratification in larger domains, the destruction of saturated TCI states by secondary KHI roll-up is not robust to changes in Reynolds number. This is indicated by the $(J,k,Re,Pr)=(0.5,4,600,300)$ TCI simulation which shows the emergence of large-amplitude cuspy nonlinear travelling waves on each density interface. These waves interact relatively weakly, however, and propagate past each other without forming any small-scale vorticity filaments due to the higher viscosity. The waves thereby slowly decay without any small-scale KHI roll-up. Such $Re$ -dependent secondary KHI roll-up also distinguishes the two simulations presented by Eaves & Caulfield (Reference Eaves and Caulfield2017).

4.5 Intermediate $k$ , varying $J$ dynamics (groups 2–4); KHI roll-up versus TCI cascade

Figure 11 shows the time evolution of $E(t)$ , $\unicode[STIX]{x1D702}(t)$ and $m(t)$ for groups 2, 3 and 4. The plots of mixing efficiency $\unicode[STIX]{x1D702}(t)$ and rate $m(t)$ illustrate how the small $J=0.14$ simulations of group 2 resemble those of group 1, whereas the large $J=0.3$ simulations of group 4 are more similar to group 5. The evolution of $\unicode[STIX]{x1D702}(t)$ and $m(t)$ for the intermediate $J=0.23$ simulations of group 3 appears to blend the behaviours of groups 2 and 4.

Figure 12. Snapshots of the density $\unicode[STIX]{x1D70C}(x,y,t)$ (left) and vorticity $\unicode[STIX]{x1D701}(x,y,t)+1$ (right) for three different times and parameter values in groups 2, 3 and 4:  $(Re,Pr)=(20\,000,10)$ for (a) $J=0.14$ at $t=169$ , (b)  $J=0.23$ at $t=104$ and (c)  $J=0.3$ at $t=93$ .

Figure 12 collects snapshots from simulations of groups 2, 3 and 4. Figure 12(a), from the $(Re,Pr)=(20\,000,10)$ simulation with $J=0.14$ , illustrates how secondary and tertiary billows again appear in addition to the primary billow, much as in the simulations of group 1 but now over a shorter horizontal domain. In figure 12(b), taken from the $(Re,Pr)=(20\,000,10)$ simulation with $J=0.23$ , a cascade of billows almost forms, but at the expense of large baroclinic vorticity generation in the region outside the primary billow. Secondary KHIs then occur over that region, but the primary billow is shielded from their destructive effect and persists to the end of the simulation. Finally, figure 12(c) shows the cascade of billows forming in the $(Re,Pr)=(20\,000,10)$ simulation with $J=0.3$ . This time, the billows all suffer secondary KHI due to the elevated baroclinic vorticity generation; small-scale vortices then fill the middle layer, as in the simulations of groups 5 and 6.

4.6 Nonlinear TCI phenomenology

To summarise, the nonlinear saturation of TCI creates a cascade of billow-like structures at small $J$ , cuspy travelling waves that draw out tight vorticity filaments and induce vigorous mixing via small-scale KHI at large $J$ , and parasitic nonlinear Holmboe waves at large $Pr$ . The occurrence in the $(k,J)$ plane of the former two primary saturation dynamics are shown in figure 3. The two latter features are captured and summarised in figure 13. Figure 13(a) plots the maximum mixing rate, $\max [m(t)]$ , for all TCI simulations in groups 1–6 against the corresponding value of $J$ . The maximum mixing rate increases almost exponentially with $J$ , quantifying the mixing associated with the observed breakdown to small-scale KHI. Figure 13(b) plots a measure of the vertical excursions of the two density interfaces, $\overline{h_{\unicode[STIX]{x1D70C}}}$ , against Prandtl number. This quantity is computed from the contours along which $\unicode[STIX]{x1D70C}(x,y_{\pm },t)=\pm 1/2$ : we first determine the instantaneous excursions, $h_{\unicode[STIX]{x1D70C}}^{\pm }(t)\equiv \max _{x}(y_{\pm })-\min _{x}(y_{\pm })$ , and then average $h_{\unicode[STIX]{x1D70C}}=[h_{\unicode[STIX]{x1D70C}}^{+}(t)+h_{\unicode[STIX]{x1D70C}}^{-}(t)]/2$ over late times to obtain $\overline{h_{\unicode[STIX]{x1D70C}}}$ . Parasitic Holmboe waves generate localised large-amplitude vertical motions, and are therefore identified by elevated levels of $h_{\unicode[STIX]{x1D70C}}$ and its standard deviation, as in the $Pr=10$ simulations.

Figure 13. (a) Peak mixing rate as a function of $J$ for all simulations. Groups 1 and 5–6 (circles) and groups 2–4 (triangles) for $(Re,Pr)=(300\,000,0.6)$ (blue), (180 000, 1) (red), (60 000, 3) (yellow) and (20 000, 10) (purple). (b) Late-time, time-averaged vertical excursion $\overline{h_{\unicode[STIX]{x1D70C}}}$ of the $\unicode[STIX]{x1D70C}=\pm 1/2$ density contours, normalised by $L_{x}=2\unicode[STIX]{x03C0}/k$ and plotted against Prandtl number. Groups 1 (dark blue), 2 (red), 3 (yellow), 4 (purple), 5 (green) and 6 (light blue). The size of the points indicates the standard deviation of $h_{\unicode[STIX]{x1D70C}}$ , again normalised by $L_{x}$ .

To address the question of how robust these finding are when we progress from two to three spatial dimensions, we conducted a single simulation in three dimensions corresponding to the $(Re,Pr)=(20\,000,10)$ simulation of group 6. Even for the relatively small domain of this simulation, the resolution requirements at $Pe=200\,000$ are substantial. In order to make the three-dimensional (3-D) computation manageable, we use the 2-D simulation at $t=50$ as an initial condition with spanwise wavenumber $2\unicode[STIX]{x03C0}$ , adding extra low-level noise to force 3-D dynamics, and then continue the simulation up to $t=100$ . The added noise changed the perturbation kinetic energy from $4.28\times 10^{-4}$ to $5.15\times 10^{-4}$ . Figure 14 compares $E(t)$ for the 2-D and 3-D simulations. The 3-D simulation shadows the 2-D simulation for much of its evolution, but begins to deviate at $t=93$ . Nevertheless, the 3-D simulation follows the 2-D simulation for sufficiently long to exhibit the same breakdown to small-scale KHI roll-ups. In figure 15(a) we show surfaces of the two density interfaces $\unicode[STIX]{x1D70C}=\pm 1/2$ for the 3-D simulation at $t=100$ . Three-dimensionality is clearly visible. Figure 15(b–d) also compares a 2-D slice through the full 3-D density field with the density field of the 2-D computation. As highlighted by the detachment of a feature from a density interface, the 2-D solution at $t=102$ shows most similarity with the 3-D slice at $t=100$ , indicating a slight lengthening of the time scales in three dimensions. Thus, at least for the addition of weak 3-D noise, it appears that the breakdown to small-scale secondary KHI is a robust phenomenon for TCI at large $k$ and $J$ .

Figure 14. Evolution of the perturbation kinetic energy for the 2-D simulation with $(Re,Pr)=(20\,000,10)$ in group 6 (thin blue) and a 3-D simulation initialised by randomly perturbing the 2-D simulation at $t=50$ (thick red).

Figure 15. (a) Surfaces of $\unicode[STIX]{x1D70C}=-1/2$ (blue) and $\unicode[STIX]{x1D70C}=1/2$ (red) at $t=100$ for a $(Re,Pr)=(20\,000,10)$ 3-D simulation simulation of group 6. (b) Vertical slice through the 3-D density field at $t=100$ , and (c,d) the 2-D density fields at $t=100$ and 102. White boxes highlight a feature detaching from the upper density interface. Density colour bar as in figures 6, 9, 10 and 12.

5.1 Streamfunction relations of the saturated states

Without dissipation (i.e. $Pe,Re\rightarrow \infty$ ), the equations of motion admit steady nonlinear wave solutions (Long Reference Long1953) that satisfy, in the frame of the wave, the relations

(5.1a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713},\quad \unicode[STIX]{x1D70C}={\mathcal{R}}(\unicode[STIX]{x1D713}),\quad \unicode[STIX]{x1D701}=Jy\frac{\text{d}{\mathcal{R}}}{\text{d}\unicode[STIX]{x1D713}}+{\mathcal{L}}(\unicode[STIX]{x1D713}),\end{eqnarray}$$

for two arbitrary functions, ${\mathcal{R}}(\unicode[STIX]{x1D713})$ and ${\mathcal{L}}(\unicode[STIX]{x1D713})$ . These streamfunction relations can be used to nonlinearly characterise the states which result after the saturation of KHI, HWI and TCI.

Figure 16. (a) Density field (left) and vorticity field (right) of the final state of the representative KHI simulation. Functional relationships (b)  ${\mathcal{R}}$ and (c)  ${\mathcal{L}}$ found from averaging over streamlines shown with black dots. Light grey lines show functional relations of the initial condition. Vertical lines indicate streamfunction values corresponding to the streamlines overlaid on the density and vorticity fields. Coloured lines show approximations to the functional relationships ${\mathcal{R}}(\unicode[STIX]{x1D713})$ and ${\mathcal{L}}(\unicode[STIX]{x1D713})$ from individual sparse 1-D profiles of density $\unicode[STIX]{x1D70C}$ , horizontal velocity $u$ and dissipation $\unicode[STIX]{x1D716}$ , where the streamfunction $\unicode[STIX]{x1D713}\approx \int u\,\text{d}y$ and the vorticity $\unicode[STIX]{x1D701}\approx -\sqrt{\unicode[STIX]{x1D716}}$ at $x=2.36$ (blue) and $x=2.66$ (red). Red and blue lines in (a) show these $x$ -locations in the spatial density and vorticity fields. Grey dots show relations taken from scattered values of $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D701}$ at $1/100\text{th}$ of the total grid points, where $\text{d}{\mathcal{R}}/\text{d}\unicode[STIX]{x1D713}$ is taken from the smooth black data.

To illustrate the construction, we use simulations in which secondary instabilities do not immediately destroy the nonlinear states reached at saturation, but viscous diffusion smooths the structures to furnish slowly evolving nonlinear waves. Unfortunately, scatter plots of the solution at each point of the computational domain at a representative time provide poor representations of the effective streamfunction relations because of the small residual time-dependence of the solutions and the need to differentiate the density relation. We therefore perform averages of $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D701}$ and $y$ over streamlines within the domain for a flow snapshot towards the end of each simulation, and then use these averages and their derivatives to construct ${\mathcal{R}}(\unicode[STIX]{x1D713})$ and ${\mathcal{L}}(\unicode[STIX]{x1D713})$ . More specifically, from the representative snapshot of the solution, we identify the streamlines from the contours of $\unicode[STIX]{x1D713}(x,y,t)=c$ , for a set of constant values $c$ , and then average the density and vorticity fields and $y$ over these curves. The former provides the relation $\unicode[STIX]{x1D70C}={\mathcal{R}}(c)={\mathcal{R}}(\unicode[STIX]{x1D713})$ ; this function can be differentiated numerically in $c$ and combined with the averages of $y$ , denoted by $Y(c)$ , and of $\unicode[STIX]{x1D701}$ , denoted by $Z(c)$ , to arrive at ${\mathcal{L}}(\unicode[STIX]{x1D713})={\mathcal{L}}(c)=Z(c)-J\;Y(c){\mathcal{R}}^{\prime }(c)$ .

The results of this averaging procedure are shown in figures 16–18 for the representative KHI and HWI computations and the TCI simulation with $(Re,Pr)=(180\,000,1)$ from group 1. The figures display the final density and vorticity fields used for the analysis, and the resulting streamfunction relations. Also included are scatter plots of $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D701}-Jy{\mathcal{R}}^{\prime }(\unicode[STIX]{x1D713})$ against $\unicode[STIX]{x1D713}$ taken from a subset of the grid points, using the streamline-averaged ${\mathcal{R}}(\unicode[STIX]{x1D713})$ function. The variability of these data about the streamline average curves is mostly indicative of residual time-dependence. Note that in the case of HWI, the non-zero phase speed of the final nonlinear state was subtracted from the velocity field prior to calculating the streamfunction in order to shift into the wave frame.

The KHI solution in figure 16 displays a clear three-branch structure for ${\mathcal{R}}$ , reflective of the existence of three density layers. The middle density layer with $\unicode[STIX]{x1D70C}\approx 0$ is not present in the initial condition but corresponds to the well-mixed region forming at the core of the overturning KH billow, and forms the long left-hand tail of the $\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D713}$ relation. The plot of ${\mathcal{L}}$ is reminiscent of Stuart’s cat’s eye vortex solution for unstratified shear flow (Stuart Reference Stuart1967), for which $\unicode[STIX]{x1D701}\propto \exp (-2\unicode[STIX]{x1D713})$ . The only significant deviation from such a smooth curve is the slight kink that arises as one moves across the border of the KH billow and the baroclinic term $Jy\,\text{d}{\mathcal{R}}/\text{d}\unicode[STIX]{x1D713}$ enters the fray.

Figure 17. A similar plot to figure 16, but for the representative HWI simulation. The approximate relations built from 1-D profiles are taken at $x=0.78$ (blue) and $x=1.09$ (red).

Figure 18. A similar plot to figure 16, but for the TCI simulation with $(Re,Pr)=(180\,000,1)$ from group 1. The approximate relations built from 1-D profiles are taken at $x=0.1$ (blue) and $x=1$ (red).

The HWI solution in figure 17 again has a three-branch structure for ${\mathcal{R}}$ , even though there are only two distinct density layers. The complication arises because two pieces of the ${\mathcal{R}}$ function have $\unicode[STIX]{x1D70C}\approx -1$ , one of which corresponds to the upper density layer whereas the other (left-hand) part characterises the vortex propagating above the density interface. That vortex draws a filament of dense fluid away from the density interface to generate the familiar cusp-like form of the Holmboe wave; the entrainment of this fluid into the vortex slightly elevates its density above the original level of $-1$ . The plot of ${\mathcal{L}}$ is also fairly complex, with three distinct branches corresponding to the two density layers and the vortex.

The TCI relations of figure 18 preserve the three-branch structure for ${\mathcal{R}}$ of the initial condition. Unlike the corresponding KHI relation, the tail of the $\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D713}$ relation in the central density layer $\unicode[STIX]{x1D70C}=0$ is very short, indicating a much weaker circulation than for KHI. This feature distinguishes the Taylor–Caulfield billow from the strong vorticity-driven overturning of KHI. Further distinction with KHI is made when examining ${\mathcal{L}}$ , which has a rather different, non-monotonic form for TCI inherited from the initial condition.

The streamfunction relations can also be built approximately using vertical profiles cutting through the nonlinear structures, such as might be furnished by observational measurements. Such observations do not readily access the streamfunction or vorticity. Instead, to mirror a reconstruction from observational data, we first extract the streamfunction by vertically integrating the horizontal velocity. Then, inspired by the fact that the total dissipation rate is equal to the integral of the square of the vorticity, we approximate the vorticity using the negative square root of the dissipation-rate profile. We also substantially coarsen the computational data, using a vertical grid with spacing $0.01$ that is approximately ten times larger than the grid spacing in the simulation. Figures (16–18) include the approximate streamfunction relations built in this fashion from two particular vertical cuts. The relations are reproduced qualitatively, if not quantitatively in all cases.

5.2 A comparison with reality

Guided by the considerations above, we examine some existing oceanographic and laboratory measurements of KHI, HWI and TCI in order to investigate the feasibility of distinguishing their respective nonlinear states using the functions ${\mathcal{R}}$ and ${\mathcal{L}}$ . More specifically, KHI is distinguished by a monotonic ${\mathcal{L}}$ function in which ${\mathcal{R}}^{\prime }$ plays little role with the density field acting almost like a passive tracer. HWI is distinguished by a three-branch structure for ${\mathcal{R}}$ in which two of the branches share nearly the same value of $\unicode[STIX]{x1D70C}$ , and ${\mathcal{L}}$ has two branches away from the region associated with the vortex core. Finally, TCI inherently contains three layers in ${\mathcal{R}}$ , and has a sharply peaked hump in ${\mathcal{L}}$ , which is otherwise reasonably flat. Nevertheless, real evolving flows involve three-dimensional time-dependent dynamics on all scales down to the dissipation scales, which our nonlinear diagnostic is unlikely to capture in any detail. The most we may hope for is that the framework provides a ‘flavour’ of the associated nonlinear dynamics.

Figure 19. Approximations to the functional relationships ${\mathcal{R}}(\unicode[STIX]{x1D713})$ and ${\mathcal{L}}(\unicode[STIX]{x1D713})$ for approximated profiles taken from (a) the observational data of van Haren et al. (Reference van Haren, Gostiaux, Morozov and Tarakanov2014), and experimental data of (b) Tedford et al. (Reference Tedford, Pieters and Lawrence2009) and (c) Caulfield et al. (Reference Caulfield, Peltier, Yoshida and Ohtani1995). Relations plotted for fits to the data; see appendix A for details of the fits.

For KHI, we construct approximate streamfunction relations using oceanographic observations of a train of KH-like billows by van Haren et al. (Reference van Haren, Gostiaux, Morozov and Tarakanov2014); figure 19(a) shows the results, exploiting fits of their 1-D data profiles. The details of the fits are provided in appendix A. A large-scale layering of the density field is evident in ${\mathcal{R}}$ , and the functional form of ${\mathcal{L}}(\unicode[STIX]{x1D713})$ is clearly of KH type. Note that a suitable choice for the phase speed of the observed structures is crucial in the construction; in producing figure 19(a), we adjusted the supposed phase speed to increase the degree of overlap between the two branches evident in ${\mathcal{L}}(\unicode[STIX]{x1D713})$ .

Figure 19 also plots approximate relationships for observations of HWI (Tedford et al. Reference Tedford, Pieters and Lawrence2009) and TCI (Caulfield et al. Reference Caulfield, Peltier, Yoshida and Ohtani1995) from wave-tank experiments. The fits to the observed 1-D profiles used for the constructions are again given in appendix A. For the data of Tedford et al. (Reference Tedford, Pieters and Lawrence2009), the three-branch structure of ${\mathcal{R}}$ is evident, as are the two distinct branches of ${\mathcal{L}}$ connected by a turning point at the minimum value of $\unicode[STIX]{x1D713}$ . The form of ${\mathcal{L}}$ is much smoother here because the mean profiles provided do not contain much vertical structure. Nevertheless, both ${\mathcal{R}}$ and ${\mathcal{L}}$ are similar to the theoretical relations for HWI. For the data of Caulfield et al. (Reference Caulfield, Peltier, Yoshida and Ohtani1995), we see a clear three-layer structure in ${\mathcal{R}}$ , and a characteristic humped structure in ${\mathcal{L}}$ , as expected for a nonlinearly saturated TCI.