1. Introduction
Buoyancy budgets of the global oceans suggest that turbulence, on scales too small to simulate directly in computational models, is an important element both to dissipate energy and to close the budget (Wunsch & Ferrari Reference Wunsch and Ferrari2004; Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018). Observations (Baker & Gibson Reference Baker and Gibson1987; Alford & Pinkel Reference Alford and Pinkel2000) have shown that such turbulence is intermittent and localised, and the nonlinear breaking of internal gravity waves has been proposed as a likely candidate for the main source of the turbulence (MacKinnon et al. Reference MacKinnon2017).
Such wave breaking could be caused by a number of mechanisms. Lombard & Riley (Reference Lombard and Riley1996) used linear stability analysis to show that instabilities on an internal wave are strongly dependent upon both the amplitude and the propagation angle of the wave, with strongly three-dimensional and two-dimensional modes being dominant in different cases. One important effect, not taken into account in these analyses, is the amplification of waves as they approach critical layers within the flow (Booker & Bretherton Reference Booker and Bretherton1967), where the flow velocity matches the phase speed. Motivated by the observation in Alford & Pinkel (Reference Alford and Pinkel2000) of coinciding shear and large amplitude waves, Howland, Taylor & Caulfield (Reference Howland, Taylor and Caulfield2021), henceforth denoted HTC21, numerically simulated the idealised flow arising from a superposition of a plane internal wave in a uniform density gradient, and a simple sinusoidal shear profile. They observed a clear three-dimensional, convective-like structure, with a definite spanwise wavelength, very different from the primarily two-dimensional Kelvin–Helmholtz billows and other shear instabilities often thought to dominate breakdown to turbulence. We exactly recreate the basic flow used in that work, but study it from a more theoretical viewpoint to analyse the mechanisms responsible.
Since the flow profile adopted by HTC21 is not steady, a traditional linear stability analysis is unsuitable for determining the structures which are important to its breakdown. One approach would be to perform linear stability analyses on ‘frozen’ background flows at different times, which has been done extensively, for example, on Kelvin–Helmholtz billows (Klaassen & Peltier Reference Klaassen and Peltier1985; Caulfield & Peltier Reference Caulfield and Peltier2000; Mashayek & Peltier Reference Mashayek and Peltier2012, to name but a few). This is a valid strategy for slowly varying background flows, or for quickly growing instabilities, but otherwise just gives a hint on the possible nonlinear behaviour.
The approach we take here is to ask, over a fixed finite time, which initial, infinitesimal perturbation is amplified by the greatest amount. This is still an entirely linear approach, but typically requires a lot more computation than traditional linear stability analyses. Indeed, even for a steady background flow, the finite time ‘optimal growth’ is still an interesting problem, since for non-normal linear operators such as in the Orr–Sommerfeld equations, the most unstable normal mode is not necessarily the structure that grows the most over a finite time interval (Schmid Reference Schmid2007).
The method we employ, direct–adjoint looping (DAL) (Corbett & Bottaro Reference Corbett and Bottaro2000; Luchini Reference Luchini2000), is essentially equivalent to that used by Arratia, Caulfield & Chomaz (Reference Arratia, Caulfield and Chomaz2013) to study optimal growth on a two-dimensional (2-D) time-evolving Kelvin–Helmholtz billow. It is an iterative method, with each solution of the Navier–Stokes equations followed by a solution of the corresponding adjoint equations, which gives the flow sensitivity with respect to a given quantity of interest. Optimising, for example, the energy at time 
$T$ of a infinitesimal perturbation at time 
$t=0$ gives a particularly simple formulation.
For the non-parallel, non-steady flow studied herein, we have no reason a priori to assume that the fastest growing disturbance is a 2-D one, since recourse cannot be made to Yih's theorem (Yih Reference Yih1955). However, since the background flow is two-dimensional, the lack of nonlinearity makes it sufficient to study individual Fourier modes in the (third) spanwise direction, for which we introduce a method using a fully 3-D numerical simulation code. Our aims in this paper are thus twofold. First, we wish to determine whether the optimal perturbations predicted by our DAL calculations are qualitatively or even quantitatively similar (for example in terms of the predicted spanwise wavelength) to the structures observed to trigger breakdown in the fully nonlinear simulations presented in HTC21. Second, we wish to determine the relative importance of shear-driven and convective growth mechanisms in the amplification of these optimal perturbations as they grow on the evolving background flow.
To address these twin aims, the remaining three sections of the paper are as follows: § 2 gives the precise flow we are considering, and gives the derivation and implementation details of the DAL algorithm. § 3 presents our results for different target times and discusses in detail two different cases, with subsections considering the influence of Prandtl number and Richardson number, and § 4 gives concluding remarks.
2. Methods
The Boussinesq equations consist of the Navier–Stokes equation, the advection–diffusion equation for buoyancy and the incompressiblity condition, governing the evolution of velocity 
$\boldsymbol{u}$, pressure p and buoyancy b
$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} ={-}\boldsymbol{\nabla} p + Ri_b b \,\boldsymbol{e}_z + \frac{1}{Re}\nabla^{2} \boldsymbol{u}, \end{gather}$$
$$\begin{gather}\frac{\partial b}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} b = \frac{1}{Re Pr}\nabla^{2} b, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0. \end{gather}$$
These equations have been non-dimensionalised using a typical length scale 
$L$, velocity 
$U$, gravitational acceleration 
$g$, density 
$\rho$, density gradient 
$\rho _z$, kinematic viscosity 
$\nu$ and density diffusion coefficient 
$\kappa$ to give the non-dimensional Reynolds number 
$Re=U L/\nu$, Prandtl number 
$Pr=\nu /\kappa$ and bulk Richardson number 
$Ri_b=g \rho _z L^{2}/\rho U^{2}$.
Following HTC21, we consider an internal gravity wave with wavevector 
$\boldsymbol {k}=(k_1,0,k_3)$ (and define 
$k=\|\boldsymbol {k}\|$) and ‘wave steepness’ 
$s$ incident on a background flow that is uniformly stratified and has a sinusoidal velocity profile, which gives a particularly simple and periodic model of a stratified shear flow. At time 
$t=0$, then,
$$\begin{gather} u = \sin{z} + \frac{s \omega}{k_1}\sin{\left( \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}\right)}, \end{gather}$$
$$\begin{gather}w ={-}\frac{s \omega}{k_3}\sin{\left( \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}\right)}, \end{gather}$$
$$\begin{gather}b = z + \frac{s}{k_3}\cos{\left( \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}\right)}, \end{gather}$$
where 
$\omega ^{2}=Ri_b ({k_1^{2}}/{k^{2}})$ is the (squared) frequency of the internal wave, so that the phase speed of this wave in isolation is given by 
$\boldsymbol {k}\omega /k^{2}$. The wave steepness 
$s$ is defined such that 
$s>1$ produces a region of negative buoyancy gradient where the wave overturns. The evolution of this 2-D background flow is complex, as will be seen in § 3. The initial evolution of the wave can be characterised as refraction of the wave by the mean shear flow. We restrict this study to that of a 2-D base flow, since motion out of the plane of the wave cannot affect this linear refraction process. However, we cannot preclude the possibility that a 3-D base flow modifies the subsequent nonlinear breakdown of the wave that we analyse here. Further research is needed to isolate how such 3-D base flows may impact internal wave breaking.
As discussed in more detail by HTC21 the idealised flow arising from (2.1) and (2.2) necessarily neglects a range of important processes that can lead to internal wave breaking in the ocean. For example, we neglect rotation in (2.1) by assuming that the Coriolis frequency 
$f$ is significantly smaller than the buoyancy frequency 
$N=-g\rho _z/\rho$. Although rotation is absent in our system, one could associate the background shear in (2.2) with rotationally dominated processes on a larger horizontal scale such as eddies or near-inertial waves. The use of an initial value problem to model wave breaking can also be questioned. However, the alternative setup of continuously forcing a stratified flow at large scales can exhibit results that depend non-trivially on the details of the forcing (Howland, Taylor & Caulfield Reference Howland, Taylor and Caulfield2020).
An infinitesimal (now three-dimensional) perturbation to (2.1) satisfies the linear equations (the primes denote the perturbation)
$$\begin{gather} \frac{\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}' + \boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} ={-}\boldsymbol{\nabla} p' + Ri_b b' \boldsymbol{e}_z + \frac{1}{Re}\nabla^{2} \boldsymbol{u}', \end{gather}$$
$$\begin{gather}\frac{\partial b'}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} b' + \boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla} b = \frac{1}{Re Pr}\nabla^{2} b', \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}' = 0. \end{gather}$$
In these equations, the 2-D background fields 
$\boldsymbol {u}(x,z,t)$ and 
$b(x,z,t)$ evolve with time according to (2.1) from initial conditions (2.2) parameterised by 
$s$, 
$Ri_b$ and 
$\boldsymbol {k}$. Since the background fields are purely two-dimensional and the perturbation is infinitesimal, there is no nonlinearity in the spanwise 
$y$ direction, and therefore it is sufficient to consider individual spanwise Fourier modes, and then vary the domain size to investigate different wavelengths. We consider the evolution of perturbations following (2.3), starting from time 
$t=0$. The choice of when to perturb the background flow in its evolution is somewhat arbitrary, but 
$t=0$ is the most obvious and agrees with HTC21.
2.1. Direct–adjoint looping
The power iteration DAL algorithm, as described by Schmid (Reference Schmid2007) and Arratia et al. (Reference Arratia, Caulfield and Chomaz2013), on the evolving background flow, allows us to find optimal perturbations at a target time 
$T$ with respect to the perturbation energy
\begin{equation} E:=\frac{1}{2L_xL_yL_z}\int (\boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{u}'+Ri_b {b^{\prime}}^{2})\,\mathrm{d}V, \end{equation}where the integration is performed over the full periodic domain.
Consider the space of state vectors 
$X=(\boldsymbol {u}',b')$ satisfying 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}'=0$ (
$p'$ can be determined from these by solving a Poisson equation). Let us define a linear operator 
$\varPhi _T$ acting on this space, defined as the solution of (2.3a)–(2.3c) up to time 
$t=T$. Further, we define an inner product
\begin{equation} \left\langle X,Y\right\rangle := \frac{1}{L_x L_y L_z}\int \left(\boldsymbol{u}'_X \boldsymbol{\cdot} \boldsymbol{u}'_Y + Ri_b b'_X b'_Y \right)\,\mathrm{d}V, \end{equation}
so that an energy for the perturbation 
$X$ is given by 
$\left \langle X,X\right \rangle /2$. We wish to find the maximum possible energy growth of a perturbation of fixed energy over a time 
$T$, i.e. to maximise the Lagrangian
\begin{equation} \mathcal{L} =
\tfrac{1}{2}\left\langle X_T,X_T\right\rangle +
\lambda\left(\tfrac{1}{2}\left\langle X_0,X_0\right\rangle -
\tfrac{1}{2}\right) +\langle \tilde{X},X_T-\varPhi_T
X_0\rangle. \end{equation}
Here, 
$\lambda$ is a Lagrange multiplier that enforces the normalisation of the initial state 
$X_0$. The precise choice of this initial energy is irrelevant, since the system is linear and we are only interested in the energy gain, i.e. the ratio of final energy to initial energy, but for a well-posed problem we nevertheless must constrain the initial energy. 
$\tilde {X}$ is a Lagrange multiplier state we call the adjoint state – for reasons which will become clear below – that enforces that 
$X_T=\varPhi _T X_0$. For an optimal perturbation, all variations of the Lagrangian vanish, so that
$$\begin{gather} 0 = \frac{\delta\mathcal{L}}{\delta\lambda} = \frac{1}{2}\left\langle X_0,X_0\right\rangle - \frac{1}{2}, \end{gather}$$
$$\begin{gather}0 = \frac{\delta\mathcal{L}}{\delta\tilde{X}} = X_T - \varPhi_T X_0, \end{gather}$$
$$\begin{gather}0 = \frac{\delta\mathcal{L}}{\delta X_0} = \lambda X_0 - \varPhi_T^{{{\dagger}}} \tilde{X}, \end{gather}$$
$$\begin{gather}0 = \frac{\delta\mathcal{L}}{\delta X_T} = X_T + \tilde{X} , \end{gather}$$
where 
$\varPhi _T^{{\dagger} }$ is the adjoint operator to 
$\varPhi _T$ with respect to our inner product. The precise definition of 
$\varPhi _T^{{\dagger} }$ is the solution of the so-called adjoint equations
$$\begin{gather} -\frac{\partial \tilde{\boldsymbol{u}}}{\partial t} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{u}} + \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \left(\boldsymbol{\nabla} \boldsymbol{u}\right)^{T} +Ri_b \tilde{b} \boldsymbol{\nabla} b = \boldsymbol{\nabla} \tilde{p} + \frac{1}{Re}\nabla^{2} \tilde{\boldsymbol{u}}, \end{gather}$$
$$\begin{gather}-\frac{\partial \tilde{b}}{\partial t} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{b} = \tilde{w} + \frac{1}{Re}{Pr} \nabla^{2} \tilde{b}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{\boldsymbol{u}} = 0, \end{gather}$$
integrated backwards in time from 
$t=T$ to 
$t=0$. The derivation of these is given in Appendix A. Equations (2.7) can be solved to give
\begin{equation} X_0 = \frac{\varPhi_T^{{{\dagger}}} \varPhi_T X_0}{\sqrt{\langle \varPhi_T^{{{\dagger}}} \varPhi_T X_0,\varPhi_T^{{{\dagger}}} \varPhi_T X_0\rangle}} \end{equation}at an optimal, which suggests the iterative method
\begin{equation} X_{n+1} = \frac{\varPhi_T^{{{\dagger}}} \varPhi_T X_n}{\sqrt{\langle \varPhi_T^{{{\dagger}}} \varPhi_T X_n,\varPhi_T^{{{\dagger}}} \varPhi_T X_n\rangle}}. \end{equation}
This is in fact precisely the power iteration eigenvalue algorithm to find the eigenvalue of greatest modulus of the linear operator 
$\varPhi _T^{{\dagger} } \varPhi _T$, and so will converge given a unique such eigenvalue. Since 
$\varPhi _T^{{\dagger} } \varPhi _T$ is self-adjoint, the eigenvalue will be real. The value of the eigenvalue is given by
\begin{equation} \lim_{n\to\infty} \langle X_n,\varPhi_T^{{{\dagger}}} \varPhi_T X_n\rangle = \lim_{n\to\infty} \left\langle \varPhi_T X_n, \varPhi_T X_n\right\rangle, \end{equation}which is exactly (twice) the energy growth we wish to maximise. Therefore, so long as the ‘initial guess’ state has a component in the direction of the optimal, (2.10) will find the maximum energy growth and the state needed to excite it.
2.2. Algorithm implementation
Because of the large storage requirements of the 2-D background state 
$(\boldsymbol {u}, b)$, which must be known at every point in time, the following algorithm is used, which employs ‘checkpointing’:
(i) The 2-D background state is evolved according to (2.1) from
$t=0$ to $t=T$. Every $100$ timesteps, the state is stored to disk.(ii) An initial perturbation state
$X_1$ is generated randomly. Set $n=1$.(iii) The value of
$X_n$ is scaled to have unit energy (required by (2.7a)).(iv) The perturbation state X n is evolved from
$t=0$ to $t=T$ in blocks of $100$ timesteps (required by (2.7b)):
(a) The background state is loaded at the start of the block, and evolved by
$100$ timesteps according to (2.1), with results at each timestep stored in memory.(b) The perturbation state is evolved from the start to the end of the block according to (2.3), using the background states stored in memory.
(v) The adjoint state
$\tilde {X}_n$ is initialised as the negative of the result of step 4 at $t=T$ (required by (2.7d)).(vi) The adjoint state is evolved from
$t=T$ to $t=0$ in blocks of $100$ timesteps (required by (2.7c)):
(a) The background state is loaded at the start of the block, and evolved by
$100$ timesteps according to (2.1), with results at each timestep stored in memory.(b) The adjoint state is evolved from the end to the start of the block according to (2.8), using the background states stored in memory (in reverse order).
(vii) The next state
$X_{n+1}$ is initialised to be the negative of the result of step 5 at $t=0$.(viii) Repeat from step 3 with
$n\to n+1$ until the (normalised) residual
(2.12)\begin{equation} \left\langle X_{n+1}-X_n,X_{n+1}-X_n\right\rangle/\left\langle X_n,X_n\right\rangle < 10^{{-}5}. \end{equation}
The algorithm was started with noisy states created by randomly exciting the first six of the Fourier modes in the vertical and streamwise directions. The precise choice of initial condition does not affect the results. Convergence was found to require 10–100 iterations before the residual was sufficiently small, with those converging to entirely 2-D results requiring more iterations.
3. Results
The (2.1), (2.3) and (2.8) are solved on a triply periodic grid with a pseudo-spectral method, using a modified version of the code developed for Parker, Caulfield & Kerswell (Reference Parker, Caulfield and Kerswell2019). We use 
$2048$ gridpoints in the streamwise (
$x$) direction, with a domain length 
$L_x=8{\rm \pi}$, and 
$512$ gridpoints in the vertical (
$z$) direction, with a domain height of 
$L_z=2{\rm \pi}$. The resolutions in these directions match those employed in the non-turbulent phase of the flow evolution by HTC21. As no nonlinearity is present in the spanwise (
$y$) direction, only the first two Fourier modes are evaluated, allowing mode-0 (i.e. spanwise independent, exactly two-dimensional) disturbances, and mode-1 disturbances, with a wavelength that matches the domain depth 
$L_y$. We vary 
$L_y$ to determine the spanwise wavelength of the fastest growing disturbance. This is a straightforward way of simulating a single Fourier mode with a full pseudo-spectral DNS code by using only three gridpoints in this direction (with dealiasing deactivated), and has the added benefit of determining for which wavelengths the spanwise independent optimal perturbation grows faster than the mode-1 optimal perturbation.
All calculations employ 
$Re=5000$, the lowest used by HTC21 for computational efficiency and direct comparison. This is too low to be realistic for oceanic flows, but is sufficiently high to capture the initial behaviour studied in HTC21, before the breakdown to turbulence. Initially we used 
$Pr=1$ and 
$Ri_b=1$; variation of these parameters is discussed below. In the initial conditions (2.2) we took the wave steepness to be 
$s=0.75$ and the wave vector to be 
$\boldsymbol {k}=({1}/{4},0,3)$ as in HTC21. With this choice, one wavelength in the streamwise direction and three wavelengths in the vertical direction fit within the periodic box. We consider only a single choice of 
$Re$ and 
$s$, rather than vary them as in HTC21, and choose instead to examine the effects of other 
$Pr$ and 
$Ri_b$. Figure 1 shows the complex, nonlinear evolution of this two-dimensional flow. By tracing ray paths for internal waves through the sinusoidal shear, HTC21 predict a cluster of critical layers close to the central height 
$z={\rm \pi}$ (shown in figure 2 of that paper). A clear amplification of the wave near these central critical layers is apparent, as predicted by classical wave theory. Regions with negative vertical buoyancy gradient are visible near the critical layers after approximately 
$t=5$.
Figure 1. The complex 2-D evolution of the background flow, a superposition of an internal gravity wave and a sinusoidal shear. (a,c,e,g) Vorticity 
$\partial u / \partial z - \partial w/\partial x$. (b,d,f,h) Buoyancy gradient 
$\partial b/\partial z$. The black contour surrounds regions with negative buoyancy gradient.
Figure 2. Perturbation vorticity (a,c,e,g) and buoyancy (b,d,f,h) of the nonlinear simulation from HTC21 whose parameters match those considered here, i.e. 
$Re=5000$, 
$s=0.75$. The black contour on the right surrounds those regions of the background flow for which the buoyancy gradient is statically unstable, which are strongly correlated with increased perturbation growth.
Figure 2 shows the nonlinear evolution of a finite amplitude, 3-D perturbation to this background flow. See HTC21 for details of this simulation. Up to 
$t=20$, the behaviour is simple, with a clear concentration of activity in those regions with negative background buoyancy gradient (indicated by a black contour in panels b,d,f,h). Subsequently, turbulence develops, and we should not expect to capture this behaviour in the present study.
We perform DAL with target times 
$T\in \{5,10,20,30\}$, chosen to capture the time horizon of the initial behaviour of the simulations of HTC21. The spanwise domain size was varied over a range 
$L_y\in [0.1,1.6]$ (initially with increments of 
$0.2$, with additional calculations where necessary to smooth the curves), which is sufficiently large to capture mode-1 structures in the 
${\rm \pi} /2$ domain of HTC21. The results are shown in figure 3(c). In each case, when the spanwise domain size 
$L_y$ is sufficiently small, the optimal structures become entirely two-dimensional, and are independent of 
$L_y$. These results are shown as dashed horizontal lines. HTC21 used a periodic domain of size 
$L_y = {\rm \pi}/2$, so assuming a normal mode structure in this direction, wavelengths 
${\rm \pi} /2$, 
${\rm \pi} /4$, 
${\rm \pi} /6$, etc. are permissible, as well as purely 2-D structures. The first six of these possible wavelengths are shown as vertical lines on the figure. Note that the results of figure 3(c) show no evidence of the ‘ultraviolet catastrophe’ found when performing stability analyses of frozen background profiles, for example in Salehipour, Peltier & Mashayek (Reference Salehipour, Peltier and Mashayek2015). This is perhaps unsurprising due to the inherent time-dependent evolution of the background flow.
Figure 3. Maximum possible final energy of a linear perturbation with initial energy 
${1}/{2}$, as the spanwise wavelength varies. The vertical grid lines mark the wavelengths of discrete modes which could be supported in the finite-sized box employed in the nonlinear DNS to which we are comparing. The wavelength of the disturbance which was actually observed, 
${\rm \pi} /8$, is marked in red. The horizontal dashed lines show the 2-D results, which are independent of 
$L_y$. The left panels show the streamwise vorticity in a 
$yz$-plane at 
$t=20$ from 
$(a)$ the DNS and 
$(b)$ the optimal perturbation calculated for a spanwise wavelength of 
$0.4$ and a target time of 
$T=30$. A limited range of 
$z$ values is shown – outside this central region there is very little contribution.
Figure 3(a) shows a 
$x$–
$y$ slice through the simulation of figure 2, which shows a clear, if noisy, mode-4 structure. The wavelength of this spanwise mode is marked by the vertical red line in figure 3(c). There is a strong agreement here with our results, as the wavelength measured from this DNS corresponds very closely to the wavelength of maximum growth from our analysis, at both 
$T=20$ and 
$T=30$.
Figures 4 and 5 show the development of the optimal perturbation for target time 
$T=30$, in respectively the 2-D case (which was found by the 3-D computations for 
$L_y\leq 0.1$), and the maximal growth case with 
$L_y=0.4$, for which there is a simple normal mode structure in the spanwise direction, shown in figure 3(b). The two figures are typical of the 2-D and 3-D mechanisms, respectively, which are qualitatively completely different from one another.
Figure 4. The 
$x$–
$z$ plane slices of the perturbation spanwise vorticity field 
$\partial u' / \partial z - \partial w'/\partial x$ (a,c,e,g) and buoyancy 
$b'$ (b,d,f,h) for the 2-D optimal perturbation (calculated using 
$L_y=0.1$) with 
$T=30$. Alternating spanwise vortices are tilted and distorted by the background shear, as is typical of the Orr mechanism. The black contour on the right surrounds regions of negative (and hence statically unstable) background buoyancy gradient.
Figure 5. Slices of the perturbation vorticity (a,c,e,g) and buoyancy field (b,d,f,h) for 
$L_y=0.4$ with 
$T=30$. The streamwise-aligned vortices are greatly amplified as they are advected. The black contour on the right surrounds statically unstable regions for which the background buoyancy gradient is negative.
The 2-D optimal perturbations, exemplified by figure 4, exploit the Orr mechanism, the transient amplification of elongated spanwise vortices as they are rotated ‘tilted over’ and distorted by a shear (Orr Reference Orr1907), which is commonly found in transient growth analyses of shear flows (Arratia et al. Reference Arratia, Caulfield and Chomaz2013; Kaminski, Caulfield & Taylor Reference Kaminski, Caulfield and Taylor2014). In this case, a patch of alternating-signed vortices is visible, at locations of high shear within the background flow. This grows in both spatial extent and amplitude as it is sheared and advected by the background. Compared with the 2-D optimal perturbation for lower target times, the vortex pattern visible is of particularly high wavenumber, which allows the Orr mechanism more time to amplify the disturbance. There does not appear to be any component, in these optimal perturbations, of a Kelvin–Helmholtz-type shear instability – which would manifest as spanwise vortices of only one sign which are not significantly distorted and sheared as the flow evolves – as opposed to the transient Orr process. The patch in figure 4 consists of vertically counterrotating vortices, either side of the black contour which surrounds the region of negative buoyancy gradient in the background flow. The optimal perturbation is therefore exploiting the unstable stratification for energy growth via spanwise counterrotating convective rolls.
Figure 6(a) shows the buoyancy flux, defined as
\begin{equation} \mathcal{J} = \frac{1}{L_x L_y L_z} \int -Ri_b b' u'_i\partial_i(b-z) \,\mathrm{d}V, \end{equation}and the shear production
\begin{equation} \mathcal{P}_S = \frac{1}{L_x L_y L_z} \int -u'_iu'_j\partial_i u_j \,\mathrm{d}V, \end{equation}
for this 2-D optimal perturbation. These terms from the perturbation energy equation are derived in Appendix B. While the buoyancy flux monotonically increases exponentially, as the convection is increasingly exploited towards the end of the time window, the shear production, orders of magnitude less important, shows two small peaks at 
$t\approx 5$ and 
$t \approx 25$ associated with Orr-like transient processes, before becoming strongly negative at the end of the time window. This suggests that the transient Orr mechanisms, while present, are not the dominant process.
Figure 6. Components of the energy budget for the 
$T=30$ optimal perturbations. Blue: buoyancy flux. Green: shear production. (a) Two-dimensional optimal. (b) Three-dimensional optimal.
The 3-D optimal perturbations, exemplified by figure 5, show no evidence of any Orr mechanism, and instead take the form of a single patch of quasi-streamwise-independent, counterrotating, streamwise-aligned vortices. As the flow evolves, this patch is advected and significantly amplified. The patch exactly aligns with one of the regions of negative buoyancy gradient in figure 1, strongly suggesting that these are indeed convective rolls, being amplified by the statically unstable stratification, though it is likely that the lift-up mechanism, a viscous algebraic instability of shear flows (Landahl Reference Landahl1980) is also being exploited. Figure 6(b) shows the buoyancy flux and shear production for this optimal perturbation. In this case, both components are roughly equally important, and both grow monotonically and roughly exponentially throughout the time window. This is in stark contrast to the 2-D results in figure 6(a), and suggests that in this case, the result is not merely transient growth but a genuine instability.
Figure 7 shows the energy growth for both of these 
$T=30$ optimal perturbations. Both of these, after some initial waviness, show apparently exponential energy growth, suggesting the dominant mechanism in each case is a convective instability, rather than the transient Orr mechanism or the algebraic lift-up mechanism. The 3-D optimal perturbation is many orders of magnitude more energetic. The figure additionally shows the energies of the 
$T=20$ optimal perturbations, when continued up to 
$t=30$, which is the target time which corresponds most closely to the results of HTC21, as a breakdown to turbulence was observed around 
$t=20$. The 3-D result for 
$T=20$ appears almost identical to that with 
$T=30$, but the 2-D result is markedly different. At the larger target time, the 2-D optimal perturbation has much smaller initial growth, but ends with exponential growth, whereas at smaller target times, there is no exponential final region, but a much larger initial growth, suggesting stronger exploitation of the Orr mechanism.
Figure 7. The evolution of the energy for the 
$T=30$ (solid) and 
$T=20$ (dashed) optimal perturbations (see figure 3). Blue: the 3-D results computed with 
$L_y=0.4$. Green: the 2-D results, using 
$L_y=0.1$.
3.1. The effects of Prandtl number
The results have thus far been restricted to the case of 
$Pr=1$, which is consistent with HTC21. However, this value is inappropriate for the oceans, as the ratio of the thermal diffusivity and kinematic viscosity of water has a typical value of 
$Pr\approx 7$. Numerous computational studies have noted that nonlinear behaviour is strongly dependent on 
$Pr$ (Smyth & Moum Reference Smyth and Moum2000; Brucker & Sarkar Reference Brucker and Sarkar2007; Mashayek & Peltier Reference Mashayek and Peltier2011; Salehipour et al. Reference Salehipour, Peltier and Mashayek2015; Parker, Caulfield & Kerswell Reference Parker, Caulfield and Kerswell2021). We therefore repeat our 
$T=30$ study at 
$Pr=7$, with an increased resolution of 
$3072$ streamwise gridpoints and 
$768$ vertically.
In this case it was found that the energy growth was approximately maximised at 
$L_y=0.3$ rather than 
$L_y = 0.4$, but the form of the optimal perturbation, shown in figure 8, is qualitatively very similar to that found for 
$Pr=1$, with the same mechanisms at work. The final energy of the perturbation in this case was 
$6.21\times 10^{13}$, significantly greater than the 
$6.77\times 10^{10}$ which was found for 
$Pr=1$. This difference is despite the fact that the background flow fields have very similar evolution and energy for flows with 
$Pr=1$ and 
$Pr=7$, as the value of 
$Pr$ has little effect on the evolution of the internal wave as it is distorted by the shear layer.
Figure 8. Initial conditions of vorticity (a) and buoyancy (b) for the optimal perturbation with 
$T=30$ when 
$Pr=7$, 
$Ri_b=1$, which is a simple normal mode in the spanwise direction with period 
$L_y=0.3$. It is very similar to the 
$Pr=1$ optimal perturbation in figure 5, despite the fact that the gain is orders of magnitude higher in this case.
3.2. The effects of Richardson number
We also repeated the 
$T=30$ calculations at 
$Ri_b=0.1$ and 
$Pr=1$, as opposed to 
$Ri_b=1$. For domain sizes 
$L_y\in \{0.1,0.4,2.0\}$, the zeroth Fourier mode was always preferred, so that 2-D spanwise-independent optimal perturbations very similar to those found at low 
$L_y$ for 
$Ri_b=1$ give more growth than any spanwise-varying perturbations. This is not particularly surprising, since with a decreased influence of stratification, 2-D shear mechanisms are expected to be strengthened, whereas 3-D convective instabilities will be damped.
4. Conclusion
The linear optimal energy growth analysis we have employed has provided an explanation of the early phase of the simulations described by Howland et al. (Reference Howland, Taylor and Caulfield2021, herein referred to as HTC21). For target time 
$T=20$, around which time the DNS of HTC21 begins to break down into turbulence for the matching parameter values, we see a clear optimal wavelength which matches that found in the DNS (figure 3). The magnitude of the energy amplification found (
$> 10^{6}$ at 
$t=20$) indicates why streamwise rolls are apparent in the simulations: if the initial disturbance has any component of this wavenumber, it is so massively amplified it will necessarily be visible as the simulation progresses. As turbulence develops as a result of this energy growth, the rolls break down and are no longer visible.
More generally, we have shown that both spanwise-independent and fully 3-D perturbations are able to exploit negative (statically unstable) buoyancy gradients which arise in the background flow as the internal wave is amplified near the critical layer, despite the high value of 
$Ri_b$ which was used. For longer times, the fully 3-D perturbations can experience orders of magnitude more energy growth than the spanwise-independent perturbations. For example, the 3-D optimal perturbation at the ideal spanwise wavelength of around 
$0.4$ (with 
$T=30$) was found to grow by a factor 
$O(10^{5})$ larger than the equivalent 2-D optimal perturbation. Understanding where these very large energy growth factors come from and how they vary with all the parameters of the problem involves unravelling what is happening in the complicated time-dependent underlying 2-D flow. Figure 1 shows a series of thin statically unstable layers (where 
$\textrm {d}b/\textrm {d}z<0$) concentrated over an 
$O(1)$ width for the parameters considered here. The optimal 3-D disturbance is focused on one of these layers (see figure 5) and possesses a spanwise length scale apparently comparable to the vertical (height) scale of the layer. This height scaling must also set the growth factor and so teasing out how this depends on all the parameters of the problem is an important challenge for the future.
This study was motivated by the hope that linear energy growth computations, which include only a single Fourier mode in the spanwise direction, could capture some essence of the expensive DNS results reported in HTC21. The fact that the optimal wavenumber which emerges in the former resembles that found significant in the latter (admittedly in only one point of comparison) augurs well for a systematic exploration of parameter space which would be impractical using DNS. HTC21 investigated different 
$s$, and found viscous effects to hinder the development of turbulence at smaller values. The method presented here could add considerably more detail on how low-steepness waves break at higher, more geophysically relevant values of 
$Re$, further investigating the subtle interplay between shear-driven and convective processes. We have already performed a brief analysis of the effects of 
$Pr$ and 
$Ri_b$, but this could be expanded greatly, to determine at what 
$Ri_b$ the 2-D mechanisms become dominant, for example. The efficiency of this method means it may be possible, in reasonable computation times, to examine 
$Pr$ up to 
$700$, characteristic of salt-stratified flows, which requires a very fine numerical grid. This study has focused on the case of a shear and wave aligned in the same 2-D plane, but for a wave coming in obliquely we may well have a qualitatively different background flow evolution, and thus the optimal perturbations could be quite different. In this oblique case, shear instabilities in particular would be altered. However, this case would require fully 3-D computations. There is clearly much to explore.
Acknowledgements
The source code used in this work is available at https://github.com/Jezz0r/Stratiflow/tree/triplyperiodic. The data behind figures 3 and 7 can be found at https://doi.org/10.17863/CAM.71802.
Funding
J.P.P. was supported by an EPSRC DTA studentship (number 1940773). The work of C.J.H. was supported by the Natural Environment Research Council through the Cambridge Earth System Science DTP (grant number NE/L002507/1).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the adjoint equations
Recalling the definitions of the state vectors and operators, (2.6) may be rewritten as
$$\begin{gather} \mathcal{L} = \frac{1}{2} \frac{1}{L_x L_y L_z}\int \left(\boldsymbol{u}'_T\boldsymbol{\cdot}\boldsymbol{u}'_T + Ri_b b'_T b'_T \right) \,\mathrm{d}V + \lambda \left(\frac{1}{2} \frac{1}{L_x L_y L_z}\int \left( \boldsymbol{u}'_0\boldsymbol{\cdot}\boldsymbol{u}'_0 + Ri_b b'_0 b'_0 \right) \,\mathrm{d}V - \frac{1}{2}\right)\nonumber\\ + \int_0^{T} \,\mathrm{d}t \frac{1}{L_x L_y L_z}\int \,\mathrm{d}V \left[ \tilde{\boldsymbol{u}}\boldsymbol{\cdot}\left(\frac{\partial \boldsymbol{u}'}{\partial t}+ \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}' + \boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} +\boldsymbol{\nabla} p' - Ri_b b' \boldsymbol{e}_z - \frac{1}{Re}\nabla^{2} \boldsymbol{u}'\right) \right.\nonumber\\ \left.+\, Ri_b \tilde{b}\left(\frac{\partial b'}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} b' + \boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla} b - \frac{1}{Re Pr}\nabla^{2} b'\right) + \tilde{p}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'\right)\right], \end{gather}$$
where the adjoint state 
$(\tilde {\boldsymbol {u}},\tilde {b})$ is now treated as a time-varying Lagrange multiplier state which enforces the evolution of 
$(\boldsymbol {u}',b')$ from 
$(\boldsymbol {u}'_0,b'_0)$ to 
$(\boldsymbol {u}'_T,b'_T)$ according to (2.3), including the addition of an adjoint pressure 
$\tilde {p}$ to enforce incompressibility.
Computing the variations of (A1) with respect to 
$\boldsymbol {u}'$, 
$b'$ and 
$p'$ we are now able to derive the adjoint equations
$$\begin{gather} \frac{\delta \mathcal{L}}{\delta \boldsymbol{u}'} ={-}\frac{\partial \tilde{\boldsymbol{u}}}{\partial t} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{u}} + \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \left(\boldsymbol{\nabla} \boldsymbol{u}\right)^{T} - \frac{1}{Re}\nabla^{2} \tilde{\boldsymbol{u}} +Ri_b \tilde{b} \boldsymbol{\nabla} b- \boldsymbol{\nabla} \tilde{p}, \end{gather}$$
$$\begin{gather}\frac{\delta \mathcal{L}}{\delta b'} ={-}Ri_b \tilde{w} -Ri_b\frac{\partial \tilde{b}}{\partial t} - Ri_b\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{b} - Ri_b\frac{1}{Re}{Pr} \nabla^{2} \tilde{b}, \end{gather}$$
$$\begin{gather}\frac{\delta \mathcal{L}}{\delta p'} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}. \end{gather}$$Equating these expressions with zero yields (2.8).
Appendix B. Derivation of buoyancy flux and shear production
Since the processes we observe are assumed to be inviscid, we start with (2.3) with the viscous terms ignored
$$\begin{gather} \frac{\partial u_i'}{\partial t} + u_j\partial_j u_i' + u_j'\partial_j u_i ={-}\partial_i p' + Ri_b b' \delta_{i3}, \end{gather}$$
$$\begin{gather}\frac{\partial b'}{\partial t} + u_j\partial_j b' + u_j'\partial_j b = 0, \end{gather}$$
$$\begin{gather}\partial_j u_j' = 0. \end{gather}$$
Taking the product of the first of these with 
$u_i'$ (applying the Einstein summation convention) and the product of the second with 
$b'$ we see
$$\begin{gather} \tfrac{\partial}{\partial t}\left(\tfrac{1}{2}u_i'u_i'\right) + \partial_j\left(\tfrac{1}{2}u_ju_i'u_i'\right) + u_i'u_j'\partial_ju_i ={-}\partial_i\left(u_i' p'\right) + Ri_b b' w', \end{gather}$$
$$\begin{gather} \tfrac{\partial}{\partial t}\left(\tfrac{1}{2}b'^{2}\right) + \partial_j\left(\tfrac{1}{2}u_jb'^{2}\right) + b'u_j'\partial_j b = 0, \end{gather}$$
where we have used the incompressibility condition to simplify some of the terms. Then integrating the sum of these two expressions (with an appropriate scaling factor of 
$Ri_b$ for the second expression) gives the equation for energy change
\begin{equation} \frac{\textrm{d}}{\textrm{d} t}\int \left(\frac{1}{2}\boldsymbol{u}'^{2} + \frac{Ri_b}{2}b'^{2}\right)\,\mathrm{d}V = \int\left(Ri_bb'w'-Ri_bb'u_j'\partial_jb -u_i'u_j'\partial_ju_i\right)\,\mathrm{d}V, \end{equation}using the divergence theorem over our periodic domain to eliminate most of the terms. The terms remaining on the right-hand side then give the expressions for the buoyancy flux (3.1) and the shear production (3.2).
Open access
