Linear Non-Modal Growth of Planar Perturbations in a Layered Couette Flow

: Layered ﬂows that are commonly observed in stratiﬁed turbulence are susceptible to the Taylor–Caulﬁeld Instability. While the modal stability properties of layered shear ﬂows have been examined, the non-modal growth of perturbations has not been investigated. In this work, the tools of Generalized Stability Theory are utilized to study linear transient growth within a ﬁnite time interval of two-dimensional perturbations in an inviscid, three-layer constant shear ﬂow under the Boussinesq approximation. It is found that, for low optimization times, small-scale perturbations utilize the Orr mechanism and achieve growth equal to that in the case of an unstratiﬁed ﬂow. For larger optimization times, transient growth is much larger compared to growth for an unstratiﬁed ﬂow as the Kelvin–Orr waves comprising the continuous spectrum of the dynamical operator and the gravity edge-waves comprising the discrete spectrum interact synergistically. Maximum growth is obtained for perturbations with scales within the region of instability, but signiﬁcant growth is maintained for modally stable perturbations as well. For perturbations with scales within the unstable region, the unstable normal modes are excited at high amplitude by their bi-orthogonals. For perturbations with modally stable scales, the Orr mechanism is utilized to excite at high amplitude neutral propagating waves resembling the neutral Taylor–Caulﬁeld modes.


Introduction
Stratified shear instability is ubiquitous in the atmosphere and the ocean and contributes to turbulence production and, more importantly, to the mixing of momentum and trace constituents (see the recent review by Caulfield [1]). A vast amount of research for over a century now-since the seminal works of Helmholtz [2], Lord Kelvin [3] and Lord Rayleigh [4]-has been devoted to the understanding of the inner workings and the evolution of the instability as well as to the quantification of its mixing efficiency.
Depending on the shear and the density stratification of the flow, there are three main types of instabilities. The first and most well studied is Kelvin-Helmholtz instability that is very common in cases in which density varies on a depth that is large compared to the depth of the shear layer. This is an unstratified instability arising from opposite gradients of vorticity at the edges of the shear layer in the flow and is partially stabilized by the stratification of the fluid [5]. The non-linear evolution of the instability produces elliptical billows [6] that are susceptible under certain conditions to secondary instabilities [7][8][9] and collapse, producing efficient mixing.
In the other two types of instabilities, stratification plays a crucial role. If the density varies substantially on a scale that is much shorter compared to the depth of the shear layer, Holmboe instability arises [10]. This instability is markedly different compared to the Kelvin-Helmholtz instability regarding its onset, its characteristics and its non-linear evolution. Instead of a single, stationary unstable mode that grows in the case of Kelvin-Helmholtz instability, there are two unstable propagating waves with slower growth rates compared to the Kelvin-Helmholtz modes. The phase speed of the waves are opposite in the case of symmetric flow, that is, in the case in which the density interface is located at the midpoint of the shear layer. In the case of asymmetry, the phase speeds as well as the growth rates differ [11]. Experimental studies [12][13][14], as well as numerical studies [15][16][17], have examined the non-linear evolution of the instability and showed that the unstable modes evolve into propagating vorticity patches above/below the density interfaces and form cusps with sporadic turbulence and mixing produced as fluid is being pulled from the tops of the cusps.
In the case in which there are multiple layers of the fluid with sharp density interfaces, then even for inflectionless shear flows, a third type of instability can arise. This instability was first noted by Taylor [18] for a three-layer fluid and was later thoroughly investigated by Caulfield [19]. It is thus referred to as Taylor-Caulfield Instability (TCI) and is the focal point of this work. Its non-linear evolution was investigated in experiments [20] as well as in numerical simulations [21][22][23][24]. It was found that the instability manifests as elliptical vortices, which resemble the Kelvin-Helmholtz billows and are situated between the density interfaces. An interesting feature observed in all numerical studies is the existence of 'parasitic' waves on the sheared density interfaces in the flow that can locally produce mixing or disrupt the structure of the billows.
The TCI has recently gained interest for two reasons. The first is that it is the least studied among the three. The second is that the fluxes in stratified turbulence are expected to amplify any fluctuations in a uniform density gradient [25]. This anti-diffusive mechanism leads to an instability of the statistical state dynamics and forms well-mixed layers separated by thin layers of sharp density gradients [26], a result that is verified by numerical simulations of stratified turbulence [27,28]. As these layers are under shear, the TCI might play a role in the rich phenomenology observed with layer merging/splitting, as well as the re-homogenization of the flow [29].
The linear instability characteristics of TCI have been investigated for an inviscid flow [19], for a viscous flow [23] and for an arbitrary number of interfaces [30]. The mechanism underlying the instability can be interpreted as arising from the interaction of gravity edge-waves. According to this idea by Taylor [18] and Goldstein [31], each density interface supports two localized gravity waves propagating upstream and downstream relative to the mean flow velocity at the position of the interface. The gravity waves supported at the different interfaces interact through action at a distance. That is, they impose velocities at the other interface and can enhance/reduce the amplitude as well as affect the propagation speed of the waves at the other interface. Instability arises when the waves supported at the different interfaces phase-lock in a mutually growing configuration. Attempting to provide a mechanistic understanding of baroclinic instability, Bretherton [32] expanded this idea to the interaction of edge-waves supported at vorticity interfaces (the waves are termed in this case as vorticity edge-waves or counter-propagating Rossby waves). This led to the interpretation of the other two main shear instabilities through this mechanistic approach [33]. Kelvin-Helmholtz instability occurs due to the resonance between two vorticity edge-waves [34] while Holmboe instability occurs due to the interaction between a vorticity edge-wave and a gravity edge-wave [35]. It was also shown that, although the analytic treatment of the edge-wave interactions is possible in the case of piecewise-constant density and vorticity profiles, the analysis and the results qualitatively hold even for smooth profiles [15,36].
Apart from modal exponential growth, shear flows also support non-modal growth even in cases in which modal instability is absent due to the non-normality of the linear operator governing the evolution of small perturbations [37][38][39]. Such transient growth can occur either due to the interaction of the edge-waves that comprise the discrete, modal part of the spectrum of the linear operator [40,41], or due to the interactions of the members of the continuous spectrum of the operator [42,43]. A paradigm for the dynamics of the latter is transient growth in parallel Couette flow, which lacks any discrete modes and exhibits only a continuum of neutral singular modes. Lord Kelvin [44] and Orr [45] showed, through plane-wave solutions with time-dependent wavenumber, that the energy of perturbations grows transiently and that the energy growth has no bound for inviscid flows. This so-called Orr mechanism was found to produce robust linear growth in many stratified two-dimensional shear flows [46][47][48][49][50][51][52] as the dynamics of the Kelvin-Orr waves do not depend on the details of the shear flow and stratification was not found to decrease the obtained growth dramatically. In addition, these waves can interact with the vorticity or the gravity edge-waves and deposit their energy to the modal part of the spectrum, thereby exciting at high amplitude the unstable or the stable analytic modes [38,48,[53][54][55]. As a result, recent studies addressing transient growth in the non-linear regime of perturbation evolution, showed that the linear transient growth is robust and capable of leading to mixing of a stratified fluid even in the absence of modal instabilities [56,57].
However, there is no study investigating non-modal growth for multi-layered stratified flows as in the TCI. In this work, we undertake this task and utilize the tools of Generalized Stability Theory [38] to investigate linear transient growth of perturbations in a three-layer constant shear flow under the Boussinesq approximation. This flow lacks any inflectional points. Therefore the TCI can be investigated in isolation as the Kelvin-Helmholtz instability is absent. We also consider the perturbations to be two-dimensional and the fluid to be inviscid. While viscosity affects the modal stability properties of the flow [23] and three-dimensional structures might play a role on the late stages of evolution for these layered flows [58], these two assumptions were made as a first step of probing the transient growth dynamics in the simplest setting. However, some of the potential effects of viscosity and three-dimensional dynamics on the obtained results will be briefly discussed. This paper is organized as follows. In Section 2, we present the linearized Boussinesq equations governing the evolution of small two-dimensional perturbations in the inviscid three-layer stratified shear flow. In Section 3 we discuss the modal characteristics of TCI, while in Section 4 we investigate the non-modal growth of perturbations. Specifically, we calculate the initial perturbations leading to the maximum energy growth over a specified time interval. We also analyze the growth achieved, the perturbation evolution as well as the mechanisms underlying this growth. Finally, we end with a discussion of our results in Section 5, their relevance in the dynamics and evolution of layered flows as well as a brief discussion on the effect of viscosity and three-dimensionality of perturbations on our results.

Linear Equations for the Evolution of Small Perturbations
We consider the parallel shear flow of a stably stratified fluid in a channel of width 2h (z * ∈ [−h, h]) with mean velocity: and mean density, that is hydrostatically balanced by the mean pressure field p * 0 (z * ). Asterisks denote dimensional variables, while s is the shear rate. The fluid consists of three layers with densities ρ m − ∆ρ, ρ m and ρ m + ∆ρ and approximate width 2h/3. The sharpness of the transition between the three layers is measured by the parameter R.
We examine the stability of this flow under two-dimensional, small-amplitude perturbations u * (x * , z * , t * ), ρ * (x * , z * , t * ), p * (x * , z * , t * ). The linearized, non-dimensional equations governing the evolution of the perturbations under the Boussinesq approximation are: where length, time, velocities, density and pressure have been non-dimensionalized by the channel half-width h, the shear rate 1/s, the velocity scale sh, the density difference ∆ρ and the pressure scale ρ m s 2 h 2 respectively and J is the bulk Richardson number J = g∆ρ/ρ m s 2 h. The non-dimensional mean velocity and density profiles are given by: and are shown in Figure 1 for R = 30. The non-dimensional boundary conditions are: Taking the curl of (3) to eliminate pressure, introducing a streamfunction ψ through the relation u = ∇ × (ψ k) and taking harmonic solutions of the form [ψ , ρ ] = [ψ(z, t),ρ(z, t)]e ikx , we rewrite the system of Equations (3)-(5) in the compact form: where is the dynamical operator,x is the state vectorx = ψ ,ρ T , ∆ = ∂ 2 /∂z 2 − k 2 is the Laplacian, and ∆ −1 is its inverse that is rendered unique by the boundary conditions:

Modal Characteristics of the TCI
The modal growth of perturbations is revealed by eigenanalysis of the dynamical operator A. We discretize the streamfunction and the density in a staggered vertical grid with N points. The state vector then becomes a (2N + 1) × 1 vector and the dynamical operator becomes a (2N + 1) × (2N + 1) matrix, whose eigenvalues λ i and eigenvectors w i are calculated numerically. All calculations were performed with N = 1600 points and numerical convergence was tested by doubling the grid points and checking on the convergence of the results.
The growth rate σ = (λ) and phase speed c r = − (λ/k) of the normal modes are shown in Figure 2 as a function of the wavenumber k and the bulk Richardson number J. There are two branches of exponential instability. The first branch at lower values of J is the TCI branch of stationary modes. As was shown by Eaves and Caulfield [23], in the limiting case of R → ∞ the discontinuity of density supports four counter-propagating gravity edge-waves (two at each interface of rapid density changes) that can interact resonantly. The phase speed of the normal modes at resonance was calculated to be where The growth rates in the TCI branch for the smooth density profile correspond closely to −ik (c − t ), therefore the TCI branch results from these resonant interactions. It is also worth noting that the unstable region is centered around the line J = 2k/9 for which resonance of the gravity edge-waves in the long-wave limit is achieved. The same holds for the maximum growth rate that is achieved for k max = 1.9 and J max = 0.5. The structure of the most unstable mode is shown in Figure 3a,b and exhibits the familiar tilt of the streamfunction against the shear, while there is no discontinuity of the streamfunction at the critical layer z = 0.  The second branch for higher values of J has modes with growth rates that are approximately one order of magnitude smaller than the growth rates of the modes in the TCI branch and also have non-zero phase speeds. The two complex conjugate eigenvalues for each mode in this branch, correspond to two counter-propagating waves that are similar to the higher-order Holmboe instabilities found by previous studies [36,59]. Their phase speeds shown in Figure 2b are very close to the phase speeds c − t for J 1 and c + t for J 2 of the neutral modes in the R → ∞ limit. However, these modes are rendered unstable by the smooth density gradients and were referred to as "overtones" as they are caused by resonances of higher order gravity edge-wave harmonics [30]. The structure of one of these modes is shown in Figure 3c,d. We observe that it is markedly different from the mode in the TCI branch, as it is centered around one of the regions of large density gradients (the oppositely propagating wave is centered around the other). We also note that similar to the TCI modes, there is no discontinuity of streamfunction at the critical layer but this layer becomes evident by the concentration of the density perturbation around it. In addition, according to the diagnostic developed by Eaves and Balmforth [24], which is based on the contribution of vorticity and density to the pseudomomentum of the modes, these modes are characterized as Taylor-Caulfield modes since vorticity does not contribute to their pseudomomentum. We therefore interpret them as overtone TCI modes, even though they are propagating and the streamfunction as well as the vorticity of the modes appears to be concentrated around one of the two density interfaces.

Non-Modal Growth of Perturbations
Apart from modal exponential growth, the non-normality of the dynamical operator can give rise to transient non-modal growth of perturbations that can in turn lead to turbulence and mixing. It is therefore of interest to calculate the initial conditions yielding the largest energy growth over a specified time interval T opt , as well as to quantify the amount of growth achieved. The time interval T opt is physically limited by the time scale of the mean flow changes and is typically in the order of 1-100 advective time units. The amount of growth is quantified by specifying a norm for the perturbation amplitude. In this work, we follow previous studies [46,48] and consider the perturbation energy density, where the overbar denotes an average over a horizontal wavelength, as the positive definite norm quantifying growth. We also note that, while a different choice of norm can be made, E has the advantage that it is conserved in the absence of a mean flow.
To address transient growth, we rewrite the energy in terms of the state vector as is the energy norm, δz is the grid spacing and the dagger denotes the Hermitian transpose. We transform the dynamics into energy coordinatesŷ = M 1/2x , where M 1/2 is the Cholesky factorization of the energy norm M = M 1/2 M 1/2 † . In these coordinates, the energy is given by the Euclidean norm of the state vector and the dynamics is governed by the operator A M = M −1/2 AM 1/2 . The optimal perturbations giving rise to the largest energy growth, over an interval T opt are then calculated through a singular value decomposition of the propagator e A M T opt = UΣV † , where U, V are unitary matrices and Σ is a diagonal matrix with the singular values along its diagonal. The optimal growth is the square of the largest singular value of the propagator, while the optimal perturbation is the corresponding column of V [38]. The main physical mechanism underlying transient growth in planar flows is the Orr mechanism [45], which is exemplified by considering an unstratified, unbounded flow. In this case the dynamics reduce to the conservation of vorticity: The solution obtained by Kelvin [60] and Orr [45] is given in terms of plane waves with time dependent vertical wavenumber: The constant phase lines of these Kelvin-Orr waves, rotate clockwise as the vertical wavenumber decreases linearly with time due to shearing of the perturbation. In addition, shearing deforms the perturbation and constantly increases the distance between the constant phase lines for perturbations leaning against the shear (m > 0). To compensate for the decreasing distance and conserve vorticity, the amplitude of the streamfunction increases (e.g., Equation (16)) leading to energy growth. The energy amplification continues until the phase lines of the waves become vertical at t = m/k. Subsequently, their energy decays as the time dependent wavenumber becomes negative and shearing decreases the distance between the phase lines. Therefore, the optimal perturbations have an initial tilt m/k = T opt , such that their phase lines become vertical at T opt and achieve a growth equal to G Orr (T opt ) = 1 + T 2 opt [61]. From the geometric point of view of non-normality, this growth arises due to the properties of the linear operator A M , which has a continuous spectrum of singular modes with eigenvalues kc and eigenfunctionŝ centered around z = c. As these modes are non-orthogonal in the energy norm, an initial perturbation can have a very large projection on some of these modes. This large projection is initially masked by the projections onto other modes, thereby yielding an order one initial norm. However, as the perturbation evolves, the projection is unmasked yielding a large norm and producing growth [38]. In the presence of boundaries, the growth achieved as a function of wavenumber is shown in Figure 4a. We observe that, for large wavenumbers, G 0 asymptotically reaches the value of growth for the unbounded flow. The reason is that for large k, most of the singular modes have the form of (17) rapidly decaying away from z = c and are unaffected by the boundaries at |z| = 1. In contrast, singular modes with lower wavenumbers have a slower decay with height and are significantly affected by the boundaries producing lower growth.  For J = 0, the energy growth achieved for low optimization times (T opt = 5) is shown in Figure 4b as a function of wavenumber k and bulk Richardson number J. We observe that, for such low optimization times, the growth achieved is very weakly dependent on the bulk Richardson number for J 2 and increases with wavenumber, asymptoting to a constant value. This value is approximately equal to G 0 as revealed in Figure 4c illustrating the ratio of the optimal growth to the corresponding growth for an unstratified flow. Therefore for early times, the kinetic energy growth produced by the Orr mechanism underlies the transient growth even for the stratified three-layer flow. Figure 5 shows the finite time Lyapunov exponent ln(G(T opt ))/2T opt for larger optimization times. For low wavenumbers, we again observe lower values of optimal growth along with a weak dependence on the Richardson number for J 2. Maximum growth is achieved either at large wavenumbers for T opt = 10 or within the TCI branch for larger T opt . Note however, that the optimal growth rate in the TCI branch is systematically larger that the modal growth rate shown in Figure 2. The reason is that the optimal perturbation achieving this growth is the initial perturbation with the maximal projection on the unstable mode. For a non-normal operator, this perturbation is not the mode w i itself, but its bi-orthogonal w ad i [38]. To find this projection, we write an initial perturbation x(0) as the sum of the eigenmodes of A M : The projection coefficients α i can be calculated by utilizing the orthogonality between w i and w ad i , as for a non-normal operator the orthogonal to the mode w i with eigenvalue λ i is the bi-orthogonal mode w ad i of the adjoint operator A † M corresponding to the eigenvalue λ * i . As a result, the coefficients α i are given by: The maximum projection on the unstable mode is thus found by substituting the bi-orthogonal as the initial perturbation to yield: for normalized bi-orthogonal vectors. Figure 6a shows the maximum projection p i for the unstable modes in the TCI branch, which is the amplified amplitude with which these modes are excited. We observe that the amplified amplitude is of the order of 50 for the most unstable wavenumbers but can reach very high values for wavenumbers close to the cut-off scale of modal instability. To illustrate the process of the excitation of the unstable mode with increased amplitude, we plot in Figure 6b the energy growth for the bi-orthogonal and compare it with the energy growth of the mode itself. We can see that the energy of both initial perturbations grows at large times exponentially with the same rate 2λ i , but during the early stages of evolution there is an initial boost for the energy of the bi-orthogonal. Figure 7 shows snapshots of the streamfunction during this early stage of evolution of the bi-orthogonal. We observe that the bi-orthogonal shown in Figure 7a has a higher tilt compared to the unstable mode shown in Figure 3a. As it is sheared over (Figure 7b), it is able to achieve rapid transient growth for early times through the Orr mechanism until it assumes the form of the normal mode (Figure 7c). For wavenumbers close to the cut-off scale, the initial growth phase lasts longer yielding larger growth and larger projections. Asymptotically close to the modal marginal stability, the duration of the growth phase tends to infinity along with the growth due to the Orr mechanism and leads to infinite values for the projection (cf. Figure 6a). Therefore, we conclude that, for wavenumbers in the TCI branch, the bi-orthogonals utilize the Orr mechanism to excite the unstable normal modes at large amplitude. More importantly though, we observe in Figure 5 that there is transient growth even for exponentially stable modes. Such growth is lower compared to the corresponding growth within the TCI branch, but is nevertheless significant. Comparison with the optimal growth for an unstratified flow revealed that it can be up to five times larger than G 0 for T opt = 50 and up to 30 times larger for T opt = 100, showing that the perturbations take advantage of the interactions between the Kelvin-Orr waves and the gravity edge-waves to produce more growth compared to the unstratified Orr dynamics.
This synergistic action is evident in the evolution of the optimal perturbation in this case which is shown in Figure 8a. We observe a large energy amplification until T opt as would be expected by the Orr dynamics with an equipartition of the energy between kinetic and potential forms (not shown). However, instead of an energy decay after t = T opt we observe a large energy oscillation which is also present at lower amplitude in the amplifying stage. This oscillation is produced by the excitation of propagating waves that are shown in Figure 8b where the streamfunction at the upper density interface (z = 1/3) is plotted as a function of x and time t. The tilted alternating bands in streamfunction indicate a wave propagating upstream with a specific phase speed. At the lower density interface, a counterpart to this wave propagating downstream with the same phase speed is also observed (not shown). To examine these waves into detail, we project the optimal perturbation v opt onto the modes of A M as in (19). Figure 8c shows the projection coefficients α i as a function of the phase speed of the modes w i . We have maximum projection near the phase speeds c ± t of the neutral waves in the R → ∞ limit, showing that the optimal perturbation is mainly composed of these neutral modes and a collection of continuous spectrum modes with similar phase speeds. The structure of these waves is similar to the unstable modes in the 'overtone' branch (not shown) with maximum values of the streamfunction at the two density interfaces and the density perturbation centered around the level z = c − t . We therefore conclude that the initial perturbation excites at very large amplitude neutral waves that are supported by the interaction of the gravity edge-waves at the two density interfaces.
We investigate the excitation of these neutral waves by the Kelvin-Orr waves and plot in Figure 9 the evolution of the optimal perturbation and its discrete Fourier decomposition in z. The optimal perturbation that is shown in Figure 9a has a very high tilt against the shear and it comprises of five bands of large wavenumbers as shown in Figure 9b. The snapshots of the evolution of the streamfunction (Figure 9c-e), as well as the Fourier decomposition of the streamfunction at later times ( Figure 9f) reveals that as the perturbation is sheared over, its wavenumbers decrease linearly with time according to the Orr mechanism. For example, as it is shown in Figure 9f after t = 10 advective time units, the first band of wavenumbers for the initial perturbation that was centered around wavenumber m = 10 has decreased to zero. As this first band of wavenumbers goes to zero and the perturbation starts to tilt vertically (Figure 9d), the next band of wavenumbers is still in the growing phase and the streamfunction tilts again against the shear after a few advection times (Figure 9e). Therefore the perturbation undergoes another period of growth. During these successive growth periods the amplitude of the neutral waves increases as they are gradually excited. That is, contours of [ψ(z = 1/3, t)e ikx ] as a function of x and t for the optimal perturbation. Also shown (dashed line) is the phase speed of the mode onto which the optimal perturbation projects (see panel (c)). (c) Projection coefficients of the optimal perturbation for k = 3 on the normal modes and the continuous spectrum modes. Also shown (dashed lines) are the c ± t phase speeds of the neutral normal modes in the R → ∞ limit.
To understand the spacing of the five bands of wavenumbers that constitute the optimal perturbation, we investigate the late stage of the evolution of the optimal perturbation that occurs after t = T opt and is shown in Figure 10. We observe that the neutral waves at the two interfaces constantly repeat a cycle in which they come in phase ( Figure 10a) and out of phase ( Figure 10c) as they counter-propagate. In between, as the waves come partially in and out of phase, the streamfunction tilts with the shear (Figure 10b) and against the shear (Figure 10d) leading to energy decay and to energy growth of the perturbation respectively. The result is the large energy oscillation observed in Figure 8. During the excitation phase for these waves, the energy oscillation induced by the neutral waves and the energy amplification brought by the Kelvin-Orr waves should synchronize in order to achieve maximum growth. To achieve this, the spacing of the five bands of wavenumbers coincides with the period of the excited neutral waves as shown in Figure 9b. The same phenomenology of neutral wave excitation holds for all other stable regions in the k-J parameter space. The only difference is that when we approach the TCI branch, the phase speed of the neutral modes goes to zero. Therefore, the period of the oscillation increases as shown in Figure 8a illustrating the energy evolution for the optimal perturbations for two other wavenumbers approaching the stability boundary. More importantly, the amplitude of the oscillation increases rapidly.

Discussion
In this work, linear non-modal growth of planar perturbations in a constant shear flow with three density layers was examined. Robust growth at a finite time interval was found with maximum growth obtained for modally unstable perturbations and suboptimal but large values of growth obtained for modally stable perturbations. It was also shown that the optimal perturbations utilize the Orr mechanism to either excite the unstable normal modes or a pair of neutral counter-propagating Taylor-Caulfield modes at high amplitude.
This dynamics might play a role in the evolution of the flow even in the non-linear stage. Numerical simulations of TCI have shown that in many cases the elliptical vortices formed by the primary instability are eventually sheared out and decay [23,24]. After this period, there are persistent oscillations in energy associated with the appearance of parasitic waves that are characterized by increased shear and mixing in one of the interfaces. While these waves were interpreted as Holmboe waves, we note that the neutral Taylor-Caulfield modes as well as the unstable overtone modes share similar characteristics. Taking into account the fact that the flow is modified by the elliptical vortices in many of the cases described above and that the growth of the parasitic waves occurs at a finite time interval, we believe that the Orr mechanism that does not depend sensitively on the details of the flow and produces rapid growth at finite time intervals, can play a significant role in the excitation of these waves, whether they are unstable modes of the flow at the given time or not.
A possible avenue for future work would be to extend the analysis to three-dimensional perturbations in a viscous flow. In full three-dimensional dynamics there is an additional transient growth mechanism, the lift-up effect, which involves the tilt-up of the background spanwise vorticity by the cross-stream perturbation velocity to produce a large amplification of the streamwise velocity [62,63]. In addition, three-dimensional unstable modes will appear at least for Richardson numbers approximately lower than J max . The reason is that Yih's theorem [64] establishes that the most unstable modes are two-dimensional only if the maximally growing perturbation at fixed bulk-Richardson number grows slower than the square root of the Richardson number. That is if f σ = max k (σ)/ √ J is a decreasing function of J. Calculation of f σ reveals that the criterion is satisfied for J 0.4, which means that the maximally growing modal instability is three-dimensional for low Richardson numbers.
Regarding the effect of viscosity on non-modal growth, it is known that for unbounded, unstratified flows, transient growth is bounded. For streamwise independent perturbations, growth is maximized at times of O(Re) and is of the order of O(Re 2 ), while for planar perturbations growth is maximized at times of O(Re 1/3 ) and is of the order of O(Re 2/3 ) [39]. We therefore expect that stable or unstable modes including those in the new branch that appears in the presence of viscosity [23], would be excited at amplitudes of that magnitude. However, the exact interplay of the Orr and the lift-up mechanisms in producing growth and in the excitation of normal modes remains to be seen in a future work.

Conflicts of Interest:
The authors declare no conflict of interest.