2.1 Equations of motion

The fluid motion is governed by the incompressible Navier–Stokes system in the Boussinesq approximation, written formally, in non-dimensional form as

for the velocity $\boldsymbol{u}(x,z,t)$ , the density $s(x,z,t)$ and the pressure field $p(x,z,t)$ . The solution vector is $\boldsymbol{q}=(\boldsymbol{u},p,s)^{\text{T}}$ and $\boldsymbol{E}$ is the projection operator onto the velocity and density fields such that $\boldsymbol{E}\boldsymbol{q}=(u,v,w,0,s)^{\text{T}}$ . The variables have been scaled using the reduced gravity $g^{\prime }=g(\unicode[STIX]{x1D70C}_{b}-\unicode[STIX]{x1D70C}_{0})/\unicode[STIX]{x1D70C}_{0}$ , where $g$ is the acceleration of gravity, $\unicode[STIX]{x1D70C}_{b}$ and $\unicode[STIX]{x1D70C}_{0}$ are, respectively, the densities at the bottom and top of the domain of depth $H$ .

The non-dimensional quantities are given by

where the $\ast$ denotes a dimensional quantity. Note that $\unicode[STIX]{x1D70C}^{\ast }$ is the density and $s$ is a scaled departure from $\unicode[STIX]{x1D70C}_{0}$ . The Reynolds number $Re=\sqrt{g^{\prime }H^{3}}/\unicode[STIX]{x1D708}$ and $Sc=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ is the Schmidt number. Here $\unicode[STIX]{x1D708}$ is the kinematic viscosity of the fluid and $\unicode[STIX]{x1D705}$ is the diffusivity of the stratifying agent.

The two-dimensional flow domain is $-L\leqslant x\leqslant L$ , $0\leqslant z\leqslant 1$ . Adiabatic and free-slip conditions, $\left[\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z,w,\unicode[STIX]{x2202}s/\unicode[STIX]{x2202}z\right]=0$ , are imposed on $z=0$ and $1$ and the lateral boundary condition at $x=\pm L$ are taken as a specified upstream inflow and, for the numerical solutions, a downstream outflow with an advective-type condition

is employed. Here $\hat{u}$ is the fluid velocity just upstream of the outflow boundary.

2.2 DJL theory for internal solitary waves

The finite-amplitude internal solitary waves that form the base state for the stability analysis are obtained from the Dubreil-Jacotin–Long (DJL) equation (Dubriel-Jacotin Reference Dubriel-Jacotin1934; Long Reference Long1953). It is an exact reduction of the steady, two-dimensional Euler equations (equation (2.1) with $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D705}=0$ ). For a Boussinesq fluid of depth $H$ , the DJL equation is, in dimensional variables, (cf. Stastna & Lamb Reference Stastna and Lamb2002)

with boundary conditions

Here $c$ is the wave phase speed, $\unicode[STIX]{x1D702}(x,z)$ is the displacement of an isopycnal from its upstream resting position. The buoyancy frequency $\bar{N}$ of the resting density profile $\bar{\unicode[STIX]{x1D70C}}(z)$ is given by

The scaled background density profile $\bar{S}(z)=[0,1]$ . In the non-dimensionalization of (2.2) with $\unicode[STIX]{x1D702}$ scaled by $H$ , equations (2.4) and (2.5) are unchanged except that $\bar{N}^{2}=-\text{d}\bar{S}/\text{d}z$ . In the frame moving with the wave the streamfunction $\unicode[STIX]{x1D6F9}=c(\unicode[STIX]{x1D702}-z)$ , the velocities $(U,W)=(\unicode[STIX]{x1D6F9}_{z},-\unicode[STIX]{x1D6F9}_{x})$ , and the density field is $S(x,z)=\bar{S}(z-\unicode[STIX]{x1D702}(x,z))$ . Upper case symbols are used for the ISW fields to distinguish them from the perturbation variables introduced in the next section.

Given $\bar{S}(z)$ , a family DJL solutions that branch from the linear long wave with phase speed $c_{0}$ are obtained for increasing values of $c$ $({>}c_{0})$ using Newton–Raphson iterations, where $c$ is added as a parameter in a pseudo-arclength continuation method (Luzzatto-Fegiz & Helfrich Reference Luzzatto-Fegiz and Helfrich2014). Standard second-order finite differences are used for the Laplacian operator. Iterations are continued until the $L_{2}$ norm of the corrections is less than $10^{-10}$ . Since solitary waves (see below) are symmetric about the wave crest, we take $\unicode[STIX]{x2202}\unicode[STIX]{x1D702}/\unicode[STIX]{x2202}x=0$ at $x=0$ and reduce the domain to $0<x<L$ and take $\unicode[STIX]{x1D702}(L,z)=0$ . $L\geqslant 6$ is made large enough that the ISW solutions are not affected by the finite size. The calculation is started from a weakly nonlinear solitary wave solution to the Korteweg–de Vries (KdV) equation (Helfrich & Melville Reference Helfrich and Melville2006). Typically $250$ cells in $z$ and $500$ in $x$ are used. However, when a particular wave is needed for a calculation, the solution is interpolated onto the desired fine grid and adjusted to convergence by additional Newton–Raphson iterations.

2.3 DJL solitary waves

A two-layer stratification is given by

where $z_{0}$ is the location of the interface with thickness scale $\unicode[STIX]{x1D706}^{-1}$ . In what follows $z_{0}=0.85$ and $\unicode[STIX]{x1D706}=80$ are used to produce a thin pycnocline close to the upper boundary. Note that because the interface is in the upper half of the domain, the solitary waves are waves of depression with negative wave amplitudes. However, the amplitude $\unicode[STIX]{x1D702}_{MAX}=|\text{min}[\unicode[STIX]{x1D702}(0,z)]|$ is defined as the magnitude for convenience. This definition is close to, but not precisely, the maximum displacement of the $S=0.5$ isopycnal.

Figure 1.

Three DJL internal solitary waves for the density profile (2.7) with $z_{0}=0.85$ and $\unicode[STIX]{x1D706}=80$ . Values for $c$ , $Ri_{min}$ and $L_{Ri}$ are indicated. The thin lines are the $S=[0.1,0.5,0.9]$ isolines and the heavy line shows the $Ri=0.25$ contour. Only a portion of the full domain is shown.

Figure 2.

(a) $\unicode[STIX]{x1D702}_{MAX}$ and (b) $L_{Ri}$ (solid) and $L_{Ri}/\unicode[STIX]{x1D709}$ (dashed) versus $c$ for the DJL solitary wave solutions with $z_{0}=0.85$ and $\unicode[STIX]{x1D706}=80$ . The square (circle) indicates $Ri_{min}=0.25$ ( $0.1$ ), and the diamond indicates the maximum conjugate state wave.

Three example DJL solutions are shown in figure 1 for increasing $c$ . The figures show the isopycnals $S(x,z)=[0.1,0.5,0.9]$ and the boundary of the $Ri<0.25$ region. Each of these waves has $Ri_{min}<0.25$ , which always occurs at the wave crest, and decreases as $c$ increases. The half-length of this region, $L_{Ri}$ , amplitude $\unicode[STIX]{x1D702}_{MAX}$ , and the wave width, $\unicode[STIX]{x1D709}$ , all increase with $c$ (the latter only for this range of $c$ ).

The relationship between $c$ and $\unicode[STIX]{x1D702}_{MAX}$ for the full family of ISWs for this stratification is shown in figure 2(a). The solutions branch from infinitesimal linear long waves at $c_{0}=0.3485$ and end at the conjugate state wave (cf. Lamb & Wan Reference Lamb and Wan1998) with $c_{cs}=0.4930$ and $\unicode[STIX]{x1D702}_{MAX}=0.3565$ . The flat-crested, infinitely broad conjugate state is found from a one-dimensional version of the DJL model following Lamb & Wan (Reference Lamb and Wan1998). Figure 2(b) shows the behaviour of $L_{Ri}$ and $L_{Ri}/\unicode[STIX]{x1D709}$ with $c$ . Waves with $Ri_{min}<0.25$ are found for $c>0.4389$ and $Ri_{min}=0.1$ at $c=0.4792$ where $L_{Ri}=0.807$ and $L_{Ri}/\unicode[STIX]{x1D709}=0.928$ .

In the following, we consider the dynamics of infinitesimal perturbations, governed by (2.1), to ISWs from the family of DJL solutions shown in figures 1 and 2. These waves are representative of internal solitary waves on similar thin-interface background stratifications. The aim here is to compare the most amplified transient dynamics of two-dimensional wave packets and compare their evolution with the amplification provided by a local asymptotic expansion in terms of unstable normal modes where a separation of scales between the DJL wave and the wave packet is assumed following Camassa & Viotti (Reference Camassa and Viotti2012).

2.4 Viscous adjusted steady states

The question naturally arises whether viscosity modifies the DJL waves sufficiently, at least for laboratory scales, such that the transient growth results are significantly affected. However, viscosity does not allow for an equation of the form of (2.4) and viscous DJL waves will be unsteady, although the temporal changes will be slow for large $Re$ . See, for example, Grimshaw, Pelinovsky & Talipova (Reference Grimshaw, Pelinovsky and Talipova2003) for an adiabatic approximation for decaying KdV solitary waves. Quasi-steady, viscously adjusted DJL waves can be found using the selective frequency damping (SFD) method (Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006) where the Navier–Stokes equations (2.1) are coupled to a low-pass temporal filtered solution $\bar{\boldsymbol{q}}$ (a sliding average of the state solution $\boldsymbol{q}$ ) of the form

where $\unicode[STIX]{x1D712}^{-1}$ is the width of the temporal filter, and $\unicode[STIX]{x1D701}$ is the amplitude of the damping applied to the Navier–Stokes equations. When initialized with a DJL wave, this approach rapidly relaxes the inviscid wave to a viscously adjusted wave as (2.8) is integrated in time. The optimal cutoff frequency $\unicode[STIX]{x1D712}=12.7423$ and damping factor $\unicode[STIX]{x1D701}=0.0534$ at $Re=10^{5}$ were computed using the approach of Cunha, Passaggia & Lazareff (Reference Cunha, Passaggia and Lazareff2015).

3.1 Direct and adjoint problem formulation

The solution is now decomposed between a base flow, the solitary wave and a perturbation such that

where $\boldsymbol{Q}=(\boldsymbol{U},P,S)^{\text{T}}$ is the DJL solution or viscous adjusted DJL wave and $\boldsymbol{q}^{\prime }=(\boldsymbol{u}^{\prime },p^{\prime },s^{\prime })^{\text{T}}$ contains the perturbation velocity field $\boldsymbol{u}^{\prime }$ , the perturbation pressure $p^{\prime }$ and the density $s^{\prime }$ .

The evolution of infinitesimal-amplitude perturbations are governed by the linearized Navier–Stokes system, where the second-order nonlinear terms ( $\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime },\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\prime }$ ) are negligible,

The sponge layer,

is added to the advection and production terms and is closest to zero near the inlet and the outlet. It allows the direct solutions to vanish as they approach the lateral boundaries of the domain far from the shear-induced regions of the ISW and without reflection. The length of the sponge layer, $x_{0}=0.73L$ , was chosen such that the maximum distance a perturbation can travel is at least twice as long as the unstable region ( $Ri<1/4$ ) of the larger-amplitude DJL wave.

In the following optimization problem, we seek to maximize an objective function $G(\boldsymbol{q}^{\prime }(\boldsymbol{x},T))$ which is a measure of the energy of the perturbation $\boldsymbol{q}^{\prime }$ at a finite time $T$ , normalized by the initial energy $E_{0}$ at time $t=0$ and is given by

where $\unicode[STIX]{x1D6FA}$ denotes the computational domain. In the present linear analysis $E_{0}=1$ . Here $N_{\ast }^{2}$ is the maximum Brunt-Väisälä frequency of the base flow $S$ , chosen to avoid division by zero in regions where $N^{2}(\boldsymbol{x})$ approaches zero.

This optimization problem is constrained by the linearized Navier–Stokes equations (3.2), which can be solved by introducing Lagrange multipliers $\boldsymbol{q}^{+}$ for the solution vector, $\boldsymbol{q}_{0}^{+}$ for the initial condition and $E_{0}^{+}$ for the initial energy. (See Joslin et al. (Reference Joslin, Gunzburger, Nicolaides, Erlebacher and Hussaini1997) for a step by step derivation in the non-stratified case and Passaggia & Ehrenstein (Reference Passaggia and Ehrenstein2013) for guidelines on the present notations.) Hence the Lagrangian is given by

and is to be rendered stationary. The scalar product $\langle \boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}\rangle$ is defined by the spatial integral $\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}\,\text{d}\boldsymbol{x}$ whereas $\langle \langle \boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}\rangle \rangle$ is defined by the double integral $\int _{0}^{T}\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}\,\text{d}\boldsymbol{x}\,\text{d}t$ where the optimization window is taken in the time interval $[0,T]$ . The constraints for the initial state $\boldsymbol{g}(\boldsymbol{u}^{\prime },s^{\prime },\boldsymbol{u}_{0}^{\prime },s_{0}^{\prime })$ and the initial amplitude $h(\boldsymbol{u}_{0}^{\prime },s_{0}^{\prime },E_{0})$ associated with $\boldsymbol{u}_{0}^{\prime }$ and $E_{0}$ respectively are

Taking variations of the Lagrangian (3.5) with respect to the state variable $\boldsymbol{q}$ and setting the result equal to zero, the adjoint system is equivalent to the one derived by Kaminski, Caulfield & Taylor (Reference Kaminski, Caulfield and Taylor2014) and is

The boundary conditions of the direct–adjoint equations are determined by the boundary terms $B$ remaining from the integrations by parts leading to (3.7) and are

where $A$ denotes the boundary part of the computational domain $\unicode[STIX]{x1D6FA}$ . These boundary terms have to be cancelled in order to insure compatibility between the direct system and its adjoint (Luchini & Bottaro Reference Luchini and Bottaro2014). Choosing the boundary conditions for the perturbation such that

leads to the condition $B=0$ with $B$ defined in (3.8) provided that the adjoint boundary conditions

are imposed. The direct flow field $(\boldsymbol{u}^{\prime },s^{\prime })$ is obtained through time marching from $0$ to $T$ and it enters the adjoint system at time $T$ with

which is to be solved backward in time from $T$ to $0$ in a direct–adjoint-loop (DAL) procedure (Schmid Reference Schmid2007).

Taking variations of the Lagrangian (3.5) with respect to the initial solution $\boldsymbol{q}_{0}$ and using the last integral in (3.8), the expression for the gradient of the objective function is

with

where the expression for $E_{0}^{+}$ is found assuming $\unicode[STIX]{x1D735}_{\boldsymbol{q}_{0}}G=0$ .

3.2 Numerical methods

Two different numerical methods were used to compute the optimal perturbations. In the case of the viscous adjusted ISWs, both the base flow, through (2.8), and the optimal perturbation were computed using the same grid which is a mixed spectral/finite differences-type discretization. The streamwise direction ( $x$ ) with $L=3$ is discretized using finite differences with $2001$ discretization points, the vertical direction ( $z$ ) uses 200 Chebyshev polynomials and the pressure is solved using a pressure projection method. The temporal integration is semi-implicit and is performed using a second-order backward Euler-type scheme for the diffusion part whereas an explicit second-order Adams–Bashforth scheme is used for the advection terms. Further details about the numerical discretization procedure can be found in Marquillie & Ehrenstein (Reference Marquillie and Ehrenstein2002) and Marquillie & Ehrenstein (Reference Marquillie and Ehrenstein2003).

Solutions of the SFD system (2.8) were considered to be converged to a viscous quasi-steady state when the $L_{\infty }$ -norm of the temporal derivative of each component for the solution vector became constant (i.e. where the relaxation imposed by the SFD is only due to viscous effects). Initial transients adjustments are rapidly damped by the SFD and within $6$ – $10$ time units depending on the value of $c$ where the viscous stabilized DJL solutions are found to decay at a rate lower than $10^{-3}$ for the velocity and the density. The convergence of the algorithm is illustrated in figure 3 for a simulation initialized with a DJL wave computed at $c=0.4681$ , $Re=10^{5}$ and $Sc=1$ . The adjusted wave is then used as a steady base state for an optimal perturbation calculation under the assumption that the time scale for the transient dynamics, ${\sim}2L_{Ri}/c$ , is much faster than the subsequent changes to the viscously adjusted DJL wave. An adjusted wave develops only slight asymmetry about its crest and was found to propagate at a speed that is only ${\approx}4\,\%$ slower than the inviscid wave. Similar to the procedure used to solve the Navier–Stokes system (2.1) and the stabilized system (2.8), the linearized perturbation dynamics and the adjoint system are solved using the same projection method to recover a divergence-free velocity field.

Figure 3.

Example of the convergence of the infinity norm $\Vert \unicode[STIX]{x2202}_{t}u\Vert _{\infty }$ (- - -), $\Vert \unicode[STIX]{x2202}_{t}w\Vert _{\infty }$ ( $\cdots$ ),  $\Vert \unicode[STIX]{x2202}_{t}s\Vert _{\infty }$ (——) using the SFD algorithm (2.8), initialized with the DJL ISW computed at $c=0.4681$ and for the control parameters $Re=10^{5}$ and $Pr=1$ . The vertical line shows when the flow was considered to be at a quasi-steady state.

Because of the shorter vertical scale of the interfacial region in an inviscid DJL wave, it proved advantageous to switch to a scheme with a uniform grid of $513$ points in the vertical. Additionally, both the forward and adjoint equations were cast in streamfunction-vorticity form. Spatial derivatives in both $x$ and $z$ were computed with sixth-order compact finite differences (Lele Reference Lele1992) and the temporal integration uses a third-order Runge–Kutta method.

We note that for both the inviscid and viscously adjusted DJL waves the domain with total length $2L=6$ was adequate to eliminate effects of the boundaries because of both the use of the windowing function (3.3) and, as shown below, because the optimal disturbances were localized near the $Ri<0.25$ zone. As a consequence, many of the inviscid DJL wave cases were conducted with no windowing, $\unicode[STIX]{x1D6F7}(x)=1$ , and periodic conditions in $x$ with no discernible effect. Tests with $L=4$ also showed no influence of the domain length.

Finally, the optimization DAL procedure was considered converged when the $L_{2}$ norm of the energy difference between two iterates was smaller than $10^{-3}$ and when the $L_{2}$ norm between two iterations of the gradient (3.12) was smaller than $10^{-2}$ which were typically achieved after $12$ – $15$ iterations.

3.3 Optimal transient growth results

Optimal perturbations were computed from the direct–adjoint system (3.2)–(3.7) for six DJL internal solitary waves from the solution family in figure 2. The wave properties $c$ , $\unicode[STIX]{x1D702}_{MAX}$ , $Ri_{min}$ , $L_{Ri}$ and $L_{Ri}/\unicode[STIX]{x1D709}$ are given in table 1. Figure 4 shows the dependence of the natural logarithm of optimal gain

versus the integration period $T$ for several of these waves. The perturbations were computed for $Re=10^{5}$ and $Sc=1$ . This Reynolds number is representative of laboratory experiments in water ( $\unicode[STIX]{x1D708}\approx 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}=10^{-2}$ and $H\approx 0.5$  m) close to those by Carr et al. (Reference Carr, Fructus, Grue, Jensen and Davies2008, Reference Carr, Franklin, King, Davies, Grue and Dritschel2017). The choice $Sc=1$ is not correct, but used for numerical convenience. The figure shows $\ln [G(T)]$ for both the inviscid and viscously adjusted ISW base states at the same $c$ . In all cases the gain was found to reach a maximal value, $\ln (G_{MAX})$ , at a time $T_{MAX}$ . While viscous effects on the primary solitary wave lead to reduced gains and slightly different optimal times, the overall behaviour is unchanged. The maximal linear gain grows dramatically with $c$ , with $\ln (G_{MAX})\approx 40$ at $c=0.4925$ for the inviscid DJL wave and ${>}30$ for its viscously adjusted counterpart. These correspond to energy gains of up to $10^{17}$ , indicating just how unstable these waves can be. It also suggests that the nonlinear evolution of the optimal perturbations needs to be assessed. The optimal time scaled with $T_{MAX}\approx 2L_{Ri}/c$ , the time for the perturbation to travel through the potentially unstable zone. This is consistent with Camassa & Viotti (Reference Camassa and Viotti2012) who argued that thin-interface ISWs are convectively, but not globally, unstable.

Figure 4.

Linear transient optimal growth gain $G$ versus the integration time $T$ for $c$ as indicated. The gain is found for $Re=10^{5}$ and $Sc=1$ . The solid (dashed) line is for the inviscid ( $Re=10^{5}$ adjusted) DJL wave base state.

Table 1.

Properties of DJL solitary waves used in the optimal perturbation calculations.

Figure 5 shows the effect of $Re$ on the perturbation energy gain, $G$ , at $T=4.43$ for the $c=0.4810$ inviscid DJL. The gain jumps approximately one order of magnitude from $Re=10^{5}$ to $10^{7}$ , although the change is small for $Re>10^{6}$ , suggesting an inviscid asymptote of $\ln (G)\sim 17$ . The effect of the varying the Schmidt number in the range $Sc=[1,12]$ , appropriate for temperature effects in water, with $Re=10^{5}$ , is shown in figure 5(b). Increasing $Sc$ also increases the gain, although the growth is modest and slows with increasing $Sc$ . In what follows we fix $Sc=1$ for computational expediency.

Figure 5.

$G$ at $T=4.43$ versus (a) $Re$ at $Sc=1$ and (b) $Sc$ at $Re=10^{5}$ for $c=0.4810$ and the inviscid DJL wave.

Figure 6.

$G$ at $T=6$ versus $\unicode[STIX]{x1D6FD}$ , the transverse wavenumber, for the $Re=10^{5}$ adjusted DJL wave with $c=0.4925$ .

We also investigated the energy growth obtained by two-dimensional against three-dimensional DAL perturbations since there is no equivalent to Squire’s theorem (Squire Reference Squire1933; Schmid & Henningson Reference Schmid and Henningson2001) for stratified shear flows. We employed a Fourier decomposition in the transverse ( $y$ ) direction and three-dimensional perturbation of the form $\boldsymbol{q}^{\prime }=\tilde{\boldsymbol{q}}(x,z)\text{e}^{\text{i}\unicode[STIX]{x1D6FD}y}$ , where $\unicode[STIX]{x1D6FD}$ is the transverse wavenumber considered. Linear gains were computed for $T=6$ with the $c=0.4925$ , $Re=10^{5}$ adjusted DJL wave (see figure 6). As expected, the optimal gain $G$ was maximum at $\unicode[STIX]{x1D6FD}=0$ , although the change is minimal for $\unicode[STIX]{x1D6FD}\leqslant 2$ . It is worth noting that a large gain, $\ln (G)\approx 20$ , can still be achieved for $\unicode[STIX]{x1D6FD}\approx 20$ . As shown below, this transverse scale is comparable to the longitudinal scale of the optimal disturbance.

The structure and evolution of the optimal perturbation ( $T_{MAX}=4.43$ ) at $Re=10^{5}$ for the $c=0.4810$ inviscid DJL wave is illustrated in figure 7. The figure shows the perturbation density $s^{\prime }$ and streamfunction $\unicode[STIX]{x1D713}$ (defined by $\unicode[STIX]{x1D713}_{z}=u^{\prime }$ and $\unicode[STIX]{x1D713}_{x}=-w^{\prime }$ ) in the frame of a solitary wave that is propagating to the right ( $c>0$ ). The fields have been normalized by their respective maximum values for clarity. The initial ( $t=0$ ) optimal disturbance is a localized wave packet just upstream of the $Ri<0.25$ region (see figure 7 a,c). Both $s^{\prime }$ and $\unicode[STIX]{x1D713}$ are tilted into the ISW-induced shear. The disturbance is dominated by its total kinetic energy which is ${\approx}3\times 10^{3}$ times larger than the potential energy. The tilt and dominance of the velocity field over the density field suggests disturbance amplification through the non-normal Orr mechanism (Orr Reference Orr1907). This will be explored further in § 4.3.

As the optimal perturbation travels through the ISW, the wave packet is amplified and its structure changes. Figure 7(b,d) shows the packet at $t=2.2$ when it is located at the wave crest. The $\unicode[STIX]{x1D713}$ field is still tilted into the shear, but extends further away from the interfacial region. The $s^{\prime }$ field is more confined and tilted with the shear. This structure of both fields and the ratio of total potential to kinetic energies, $0.238$ , are consistent with a standard unstable K–H normal mode.

Figure 7.

(a,c) The structure of the optimal linear perturbation for the $c=0.4810$ DJL wave base state with $Re=10^{5}$ (at $T_{MAX}=4.43$ ). (a,b) The structure of the perturbation density field $s^{\prime }$ and (c,d) the streamfunction $\unicode[STIX]{x1D713}$ . The solid lines are the $S=0.5$ and the $Ri=0.25$ contours. (b,d) The structure of the optimal linear at $t=2.2$ from a forward linear calculation. In both (a,c) and (b,d) the perturbation fields have been normalized by their maximum values.

Figure 8.

An $x$ – $t$ plot of $\unicode[STIX]{x1D713}$ , normalized to a maximum of one at each time, on the $S=0.5$ contour from the forward linear calculation in figure 7(b,d). The dashed lines are at $|x|=L_{Ri}$ .

The temporal evolution of the optimal wave packet further illustrated in figure 8 by means of an $x-t$ diagram of $\unicode[STIX]{x1D713}$ along the $S=0.5$ isopycnal. Again, $\unicode[STIX]{x1D713}$ is scaled by the maximum value at each time. The wave packet remains localized and travels at a quasi-constant group velocity through the solitary wave. The phase and packet group velocities are nearly identical. The carrier, or central, frequency of the waves in the packet is $\unicode[STIX]{x1D714}=-14.85$ . The negative sign is used since $\unicode[STIX]{x1D714}$ is a Doppler shifted value in the frame moving with the solitary wave and we will take all real wavenumbers $k>0$ . The carrier frequency of the wave packet is found using the discrete Hilbert transform in $t$ of $\unicode[STIX]{x1D713}(x,t)$ . At fixed $x$ , this gives $Z(t)=H(\unicode[STIX]{x1D713}(t))=A\exp (\text{i}\unicode[STIX]{x1D703})$ , where $H(\cdot )$ is the Hilbert transform. Then $A(t)=|Z(t)|$ is the packet envelop and $\text{d}\unicode[STIX]{x1D703}/\text{d}t$ is the instantaneous frequency (Oppenheim, Schafer & Buck Reference Oppenheim, Schafer and Buck1999). The carrier frequency is then defined as the average of $\text{d}\unicode[STIX]{x1D703}/dt$ for $A/A_{0}>0.7$ , where $A_{0}=\max (A)$ . These values are then averaged over a range of $x$ around $x=0$ to obtain the reported frequency, $\unicode[STIX]{x1D714}$ . Standard Fourier analysis gives estimates generally within $\pm 0.1$ ; however, the Hilbert analysis is used since it also gives the envelop characteristics and allows $\unicode[STIX]{x1D714}$ to be defined in a way that avoids potential dispersive effects at the edges of a packet. A similar analysis in $x$ at fixed $t$ gives the central wavenumber of the packet (see the discussion in § 4).

The structure and behaviour of optimal disturbances for all the other ISWs (i.e. $c$ ) examined are consistent with figures 7 and 8. The optimal packets are initially situated just upstream of the $Ri<0.25$ zone, remain coherent and compact as they move through the ISW and have well-defined wavenumbers and frequencies. As a consequence, a local approach to transient growth might be expected to be relevant.

4.1 Local stability analysis

Normal modes are sought of the form

where the (real) frequency is $\unicode[STIX]{x1D714}$ and the (complex) wavenumber $k=k_{r}+\text{i}k_{i}$ is the spatial eigenvalue. Here $k_{r}$ is the horizontal wavenumber and $k_{i}$ is the spatial growth rate. Linearizing the Navier–Stokes system for the parallel flow profiles $(U(z),0,S(z))$ shown in figure 9(a,b), taking the curl of (3.2) in the $(x{-}z)$ plane and combining with the density equation, the Taylor–Goldstein (T–G) equation reads

and $N^{2}(z)=-S_{z}$ (in non-dimensional variables). The eigenvalue problem (4.2) is solved for specified background state using MATLAB’s bvp5c adaptive boundary value problem solver (Matlab 2014).

Figure 9.

Vertical profiles of (a) $U(z)$ in the wave frame, (b) $S(z)$ and (c) $Ri(z)$ for the $c=0.4810$ DJL wave at $|x|=0,0.2,0.4,0.6$ and $0.8$ .

Vertical profiles of velocity, density and Richardson number at several positions within $-L_{Ri}<x<L_{Ri}$ are shown in figure 9(a) for the $c=0.4810$ ISW. Note that the velocities are in the frame moving with the solitary wave and that the structure is symmetric about $x=0$ . As $|x|$ decreases, the minimum Richardson number decreases to $Ri_{min}=0.096$ at $x=0$ . The spatial growth rate $k_{i}(x;\unicode[STIX]{x1D714})$ ( ${>}0$ for left-going disturbances) found from (4.2) with profiles at $|x|=[0:0.1:0.8]$ for a range of $\unicode[STIX]{x1D714}$ is plotted in figure 10(a). For fixed $\unicode[STIX]{x1D714}$ , the growth rates increase smoothly as $|x|\rightarrow 0$ and the frequency of maximum growth rate, $\unicode[STIX]{x1D714}\approx -15$ , depends only weakly on $x$ .

Figure 10.

(a) The spatial growth rate $k_{i}$ versus frequency $\unicode[STIX]{x1D714}$ from the solution of the Taylor–Goldstein equation (4.2) for the $c=0.4810$ DJL wave. The curves are shown at $|x|=[0:0.1:0.8]$ from top to bottom. (b) The total gain $G$ versus $\unicode[STIX]{x1D714}$ of an unstable perturbation moving through the $Ri<0.25$ zone ( $L_{Ri}=0.846$ ) from (4.6) using $k_{i}(x;\unicode[STIX]{x1D714})$ in (a).

4.2 WKB approximation to transient growth

The local stability properties are linked to the spatio-temporal growth of the perturbation by integrating $k_{i}(x;\unicode[STIX]{x1D714})$ from $x=L_{Ri}$ to $x=-L_{Ri}$ . The evolution of the amplitude $A(x,t)$ of a small-amplitude perturbation is therefore given by

where $k^{-}$ indicates the waves that grow while propagating to the left ( $k_{i}>0$ ) and $A(L_{Ri})$ is the amplitude at the onset of the region $Ri<1/4$ . The spatial growth of a perturbation with fixed $\unicode[STIX]{x1D714}$ is given by the real part of (4.3)

For these linear disturbances the corresponding energy gain is

and the total gain after passage of the disturbance through the solitary wave is

The dependence of the total gain $G$ from (4.6) as a function of $\unicode[STIX]{x1D714}$ is plotted in figure 10(b) for the results in figure 10(a). The maximum gain $\ln (G_{MAX})=12.48$ occurs at $\unicode[STIX]{x1D714}_{MAX}=-14.74$ . Recall that Troy & Koseff (Reference Troy and Koseff2005) and Barad & Fringer (Reference Barad and Fringer2010) found observable instabilities required $\bar{\unicode[STIX]{x1D714}}_{I}T_{W}>5$ which gives $\ln (G)>10$ . Similarly, Lamb & Farmer (Reference Lamb and Farmer2011) found $2\bar{k}_{i}L_{Ri}>4$ , or $\ln (G)>8$ . Figure 10(b) shows that $\ln (G)>8$ occurs over a significant frequency range $-24<\unicode[STIX]{x1D714}<-7.5$ .

Figure 11 shows a comparison of $\ln (G_{MAX})$ as a function of the ISW phase speed $c$ for the optimal perturbations from the inviscid DJL waves, the viscously adjusted ( $Re=10^{5}$ ) DJL waves and the WKB analysis described above. The energy gain is always largest for the optimal perturbations on the inviscid DJL wave. Recall that both optimal disturbance calculations were based on $Re=10^{5}$ , while the WKB results are for inviscid disturbances. For purely inviscid flows the optimal disturbance growth would be even greater (cf. figure 5 a). The difference between the optimal DAL gain and the WKB estimate becomes more evident as $c$ , and therefore $L_{Ri}$ , decreases. However, the difference in $\ln (G_{MAX})$ is only weakly dependent on $c$ , decreasing from $2.9$ to $2.0$ between $c=0.4442$ and $0.4925$ .

The frequencies for maximal WKB growth and the carrier frequencies of the optimal perturbations from the DAL calculation for the inviscid DJL waves are shown in figure 12. The agreement between the two approaches is similarly quite good, but does degrade as $c$ decreases, where the difference in $\ln (G_{MAX})$ is also greatest.

Figure 11.

The maximum gain, $G_{MAX}$ , versus $c$ . The solid diamonds (solid squares) are for the DAL optimal perturbations to the inviscid ( $Re=10^{5}$ adjusted) DJL wave base state. The circles are the maximal gains from the spatial WKB analysis. The triangles are for the $k_{+}$ free wave packets with carrier frequencies equal to those of the optimal disturbances (see § 5).

Figure 12.

The frequency for maximum gain, $\unicode[STIX]{x1D714}_{MAX}$ , from the DAL optimal disturbances (solid diamonds) and the spatial WKB analysis (circles) versus $c$ .

Figure 13.

$k_{i}$ , $k_{r}$ , $t=\int c_{g}^{-1}\,\text{d}x$ and $G$ versus $x$ through the $Ri<0.25$ zones for inviscid DJL waves with (a) $c=0.4579$ , (b) $c=0.4810$ and (c) $c=0.4925$ . The WKB results (solid) are computed for $\unicode[STIX]{x1D714}=-14.74$ , $-14.85$ and $-16.79$ , respectively, the carrier frequencies of the optimal transient growth wave packets (see figure 12). The dashed curves show the corresponding characteristics from forward linear integrations of the optimal transient growth disturbances with $Re=10^{5}$ .

A detailed comparison between the disturbance properties through the unstable zone from the linear optimal perturbations, from the DAL method and those from the WKB approach is given in figure 13. The dashed lines in the lower three panels of each column give the wavenumber in the centre of the packet, $k_{r}(x)$ , the packet envelop peak position, $x(t)$ (or $t(x)$ ) and the energy gain, $G(x)$ . These were obtained from the forward linear calculations initiated with the optimal DAL disturbances for $c=0.4579$ , $c=0.4810$ and $c=0.4925$ . As discussed earlier, both $k_{r}(x)$ and $x(t)$ were found using Hilbert transforms in $x$ of disturbance $\unicode[STIX]{x1D713}(x,t)$ on the $S=0.5$ isopycnal (see figure 8). The solid lines show $k_{i}(x)$ , $k_{r}(x)$ , the disturbance $x-t$ relation,

and $G(x)$ from the WKB approach. The frequency of the optimal perturbations from the DAL analysis (figure 12) were used in the WKB calculations. The agreement between the central (real) wavenumbers $k_{r}(x)$ is quite good as is the agreement between the WKB group speed $c_{g}(x)$ (the slopes of the optimal disturbance trajectories $x(t)$ ). The differences in the energy growth $G(x)$ curves are due almost entirely to an initial growth phase in $x>L_{Ri}$ for the optimal disturbance and accounts for the difference in $G_{MAX}$ in figure 11. The difference in $\ln (G_{MAX})$ is largely independent of the DJL wave speed $c$ , only varying from $2.9$ to $2$ between $c=0.4442$ and $0.4925$ . Since vertical shear is present for $x>L_{Ri}$ where $Ri>1/4$ , the inviscid Orr mechanism (Orr Reference Orr1907) may be responsible for this initial transient growth and is addressed next.

4.3 Non-modal transient growth

As mentioned above, the WKB analysis does not account for the non-normality associated with the Taylor–Goldstein operator (4.2). In regions where $Ri>1/4$ , shear can still play a destabilizing role. In his original work, Orr (Reference Orr1907) showed that in the case of a simple inviscid parallel shear flow, perturbations with a non-zero streamwise wavenumber $k$ could produce transient growth through the kinematic deformation of the perturbation vorticity by the base flow advection and shear. Later, Farrell & Ioannou (Reference Farrell and Ioannou1993a ) derived an analytic solution for the Orr temporal growth rate,

for two-dimensional perturbations in a constant, unstratified linear shear $U_{c}^{\prime }$ over the optimization time $\tilde{T}$ . Perturbations with small values of $k$ will need more time to accomplish the Orr-induced rotation than perturbations with large values of $k$ . Hence the optimization time $\tilde{T}$ in (4.8) is rescaled to take into account this such effect and it is linked to the optimization time $T$ through the transformation $T=2\unicode[STIX]{x03C0}c\tilde{T}/k$ . Farrell & Ioannou (Reference Farrell and Ioannou1993b ) also derived an approximation for the Orr gain in the case of an unbounded stratified shear flow where both the shear and the stratification are linear, but in the large $Ri$ limit. Note that in these idealized cases, $\unicode[STIX]{x1D70E}_{Orr}$ is independent of the streamwise wavenumber $k$ .

In the following we examine the possibilities for Orr-type transient growth in regions of the flow just upstream of the $Ri<0.25$ zone using the Taylor–Goldstein equation (4.2) and maximizing the disturbance energy for short times $T$ . Therefore we are seeking a non-modal approach to the optimization problem for the gain

with $\tilde{\boldsymbol{q}}=[\unicode[STIX]{x1D713},s^{\prime }]$ . The energy norm $\Vert \cdot \Vert ^{2}$ is defined equivalently to (3.4) which in matrix form is given by $\Vert \tilde{\boldsymbol{q}}(T)\Vert ^{2}=\tilde{\boldsymbol{q}}^{\text{T}}(T)\unicode[STIX]{x1D648}\tilde{\boldsymbol{q}}(T)$ where $\unicode[STIX]{x1D648}$ is the matrix containing the weights of integration. Such a problem can be solved considering the initial value problem associated with the parallel flow Taylor–Goldstein equation (4.2) using the local buoyancy and velocity profiles. In streamfunction/density perturbation formulation (4.2) becomes

Recasting (4.10) in matrix form

the initial value problem becomes

where the initial condition $\tilde{\boldsymbol{q}}_{0}=[\unicode[STIX]{x1D713}_{0},s_{0}]$ is to be optimized to maximize (4.9) at a given time $T$ . The system (4.12) can be solved efficiently by performing a singular value decomposition (Schmid & Henningson Reference Schmid and Henningson2001). The system (4.10) is first recast into a generalized eigenvalue problem of the form

where the matrices $\unicode[STIX]{x1D726}$ and $\unicode[STIX]{x1D63E}$ contain the eigenvalues and eigenvectors, respectively. The energy norm $\Vert \cdot \Vert ^{2}$ is then converted to the $L_{2}$ -norm (Euclidian norm) $(\cdot )_{2}^{2}$ through the transformation $\Vert \unicode[STIX]{x1D63E}\Vert ^{2}=\unicode[STIX]{x1D63E}^{H}\unicode[STIX]{x1D648}\unicode[STIX]{x1D63E}=\boldsymbol{F}^{H}\boldsymbol{F}=(\boldsymbol{F})_{2}^{2}$ , where $\boldsymbol{F}^{H}$ denotes the Hermitian conjugate of $\boldsymbol{F}$ . The optimal Orr gain is finally given by

where $\unicode[STIX]{x1D70E}_{1}$ is the first singular value.

Figure 14.

The maximum gain of the Orr mechanism, $\ln (G_{Orr})$ , versus the optimization time $T$ . The solid lines are the optimal transient growth using the Taylor–Goldstein equation (4.2) for the stable profiles at $x=[0.88,0.96,1.12,1.36]$ ( $Ri_{min}=[0.26,0.33,0.59,1.82]$ ) for the inviscid DJL wave at $c=0.4810$ . The dashed line is the prediction for homogeneous shear from (4.8) with $U_{c}^{\prime }=7.5$ (i.e. the maximum shear for $Ri_{min}=0.26$ ).

As shown in figures 11 and 13, the WKB approach differs by a shift in the energy gain $\ln (G)\approx 2{-}3$ over an initial advection time scale $t\approx 0.5$ . Using the $c=0.4810$ case as an example, the optimal gain (4.14) is computed for the ISW $U(z)$ and $S(z)$ profiles at $x=[0.88,1.12,1,24,1.36]$ which correspond to $Ri_{min}=[0.26,0.33,0.59,1.82]$ respectively. The results for $k=36.6$ , the average initial packet wavenumber, are shown in figure 14. The Orr mechanism contributes to the optimal gain $\ln (G)$ by amounts of $6{-}8.5$ and appears to be only weakly affected by changes in $Ri$ . The optimal time $T\approx 1$ is twice the time scale in figure 13. However the gain difference between $T=0.5$ and $T=1$ is minor, and for $T>0.5$ the K–H instability provides a faster growth than the Orr mechanism.

Also included in figure 14 is the gain, $\exp (2\unicode[STIX]{x1D70E}_{Orr}\tilde{T})$ , found using (4.8) with $U_{c}^{\prime }=7.5$ . This value is equal to the maximum shear at $x=0.88$ . Over times $T<1$ , this simple unstratified model provides a relatively good estimate for $\ln (G_{Orr})$ from (4.10) and (4.14) and shows that for small optimization times, stratification does not play a major role in the non-normal growth process (Farrell & Ioannou Reference Farrell and Ioannou1993b ). As expected for longer optimization times $T$ , the stratification damps transient growth.

6.1 Long-time evolution of ISWs forced by optimal perturbations

The evolution of base states subjected to continuous, periodic forcing with spatial structure given by the optimal perturbation has been considered by Brandt et al. (Reference Brandt, Sipp, Pralits and Marquet2011) and Sipp & Marquet (Reference Sipp and Marquet2013). However, a simpler, approximate approach that follows from models of adiabatic decay of KdV solitary waves (Grimshaw et al. Reference Grimshaw, Pelinovsky and Talipova2003) is possible. From the previous section, each encounter of a disturbance packet with a wave extracts a small fraction of energy from the ISW given by $\unicode[STIX]{x1D6FE}(c)$ . Assuming a slow, adiabatic adjustment to each encounter, the ISW energy evolves as

where $c(E_{ISW})$ and $\unicode[STIX]{x1D702}_{MAX}(E_{ISW})$ are found from DJL solution family, and $N_{p}$ is the number of disturbance packets per unit length. The disturbance packets have a length ${\approx}0.5$ (see figure 7), so that $N_{p}=2$ approximates continuous encounters and thus the maximum rate of decay. Cases with $N_{p}<2$ can be obtained by a simple rescaling of  $x$ .

Figure 21 shows solutions for the maximal decay, $N_{p}=2$ , of an initial wave with $\unicode[STIX]{x1D702}_{MAX}=0.3025$ ( $c=0.49$ ) subject to both optimal disturbance and free wave packets. For the optimal disturbance the decay 90 % of the way to the largest stable wave occurs on a spatial scale $x_{diss}\approx 100$ and time scale $t_{diss}\approx 200$ . The free wave decays occurs slightly sooner, although the final waves are quite different. As a consequence, optimal disturbances can result in substantially more total loss of energy from the initial solitary wave. In this example optimal disturbances extract $82\,\%$ of the initial ISW energy compared to the $52\,\%$ for the free wave packets. Taking $H=100$  m and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=2~\text{kg}~\text{m}^{-3}$ as representative of coastal settings, the initial wave has a total energy of $1.71\times 10^{6}~\text{J}~\text{m}^{-1}$ and the dimensional decay scales become $x_{diss}\approx 10$ km and $t_{diss}\approx 4$  h. The rate of dissipation per unit width, $cN_{p}\unicode[STIX]{x1D6FE}E_{ISW}$ , falls from ${\approx}1100~\text{W}~\text{m}^{-1}$ at $x=0$ , to ${\approx}40~\text{W}~\text{m}^{-1}$ at $x=x_{diss}/2$ . This initial rate is extremely large, although the later rate is comparable to observational estimates of $10{-}50~\text{W}~\text{m}^{-1}$ (Moum et al. Reference Moum, Farmer, Smyth, Armi and Vagle2007; Shroyer et al. Reference Shroyer, Moum and Nash2010a ).

Figure 21.

Slowly varying estimates of $\unicode[STIX]{x1D702}_{MAX}(x)$ from (6.3) with $\unicode[STIX]{x1D702}_{MAX}(0)=0.49$ ( $c=0.49$ ),  $N_{p}=2$ and $\unicode[STIX]{x1D6FE}(c)$ from the optimal (solid line) or free wave (dashed line) disturbances. The dash-dot line shows the largest wave stable to optimal disturbances at $\unicode[STIX]{x1D702}_{MAX}=0.158$ and $c=0.451$ .