*Article* **Non-Modal Three-Dimensional Optimal Perturbation Growth in Thermally Stratified Mixing Layers**

**Helena Vitoshkin 1,\* and Alexander Gelfgat <sup>2</sup>**


**Abstract:** A non-modal transient disturbances growth in a stably stratified mixing layer flow is studied numerically. The model accounts for a density gradient within a shear region, implying a heavier layer at the bottom. Numerical analysis of non-modal stability is followed by a full threedimensional direct numerical simulation (DNS) with the optimally perturbed base flow. It is found that the transient growth of two-dimensional disturbances diminishes with the strengthening of stratification, while three-dimensional disturbances cause significant non-modal growth, even for a strong, stable stratification. This non-modal growth is governed mainly by the Holmboe modes and does not necessarily weaken with the increase of the Richardson number. The optimal perturbation consists of two waves traveling in opposite directions. Compared to the two-dimensional transient growth, the three-dimensional growth is found to be larger, taking place at shorter times. The nonmodal growth is observed in linearly stable regimes and, in slightly linearly supercritical regimes, is steeper than that defined by the most unstable eigenmode. The DNS analysis confirms the presence of the structures determined by the transient growth analysis.


**Citation:** Vitoshkin, H.; Gelfgat, A. Non-Modal Three-Dimensional Optimal Perturbation Growth in Thermally Stratified Mixing Layers. *Fluids* **2021**, *6*, 37. https://doi.org/ 10.3390/fluids6010037

Received: 11 December 2020 Accepted: 4 January 2021 Published: 11 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Keywords:** stratified mixing layer; non-modal instability; Kelvin-Helmholtz instability; Holmboe instability

## **1. Introduction**

The mixing layer flow is the simplest configuration allowing for the well-known Kelvin-Helmholtz (or KH) instability that takes place when two parallel flows having different velocities meet. The instability develops as a wave (or KH mode) in the shear layer separating two uniform flows. Studies of this phenomenon started in the early works of Lord Kelvin [1] and Strutt and Lord Rayleigh [2]. The phenomenon appears to be so complicated and to pose so many questions that in spite of hundreds of studies published, it remains the subject of many current researches [3–6]. For the details, the reader is referred to review papers [7–9].

The isothermal mixing layer is known to be linearly unstable either in the inviscid limit or starting from relatively low Reynolds numbers within the Newtonian viscous model [10]. It is also known that the flow remains linearly stable for the streamwise dimensionless wavenumber larger than unity [10,11]. At the same time, studies [12,13] showed that temporal non-modal disturbances growth can take place in the isothermal inviscid and viscous mixing layer flows, respectively. A similar mechanism was discovered recently in round jets [5].

The classical result on the linear stability of an inviscid mixing layer flow stably stratified by density (or temperature) stays that the flow becomes linearly stable for the gradient Richardson number exceeding 0.25 [10–15]. Later studies showed that along with the monotonic Kelvin-Helmholtz instability, an oscillatory instability is also possible [16]. The latter is the Holmboe instability that develops as two waves traveling in opposite

ance.

mary and conclusions are derived in Section 5.

mal expansion acceleration,

notes the Laplacian operator.

**2. Formulation of the Problem and Numerical Techniques**

directions [10,17]. Furthermore, numerical calculations [18–20] supported by experiments of [21] showed that a three-dimensional mode traveling at an angle to the base velocity may attain the largest growth rate. For example, a shear layer characterized by a relatively thin region of stably stratified fluid can exhibit Holmboe instability [16,21,22], whose growth rate increases with increasing stratification. Despite the valid Squire transformation proved in [23], the flow is predicted to be linearly unstable to a three-dimensional disturbance. [21] showed that a three-dimensional mode traveling at an angle to the base velocity may attain the largest growth rate. For example, a shear layer characterized by a relatively thin region of stably stratified fluid can exhibit Holmboe instability [16,21,22], whose growth rate increases with increasing stratification. Despite the valid Squire transformation proved in [23], the flow is predicted to be linearly unstable to a three-dimensional disturbance. It was shown in [10] that the instability at small Richardson numbers is governed by *Fluids* **2021**, *6*, x FOR PEER REVIEW 2 of 17 latter is the Holmboe instability that develops as two waves traveling in opposite direc-

latter is the Holmboe instability that develops as two waves traveling in opposite directions [10,17]. Furthermore, numerical calculations [18–20] supported by experiments of

*Fluids* **2021**, *6*, x FOR PEER REVIEW 2 of 17

It was shown in [10] that the instability at small Richardson numbers is governed by two monotonic KH modes, one of which can become unstable. With further increase of the Richardson number, the pair of monotonic modes are replaced by a pair corresponding to complex conjugated eigenvalues related to the Holmboe instability. It was shown recently that KH instability plays an important role in the formation of interfaces between coating and substrate materials and metallic materials [24,25]. The non-linear evolution of KH and Holmboe instabilities was studied in quite a large number of articles (e.g., [17,26–33]). A difference in the nonlinear disturbances growth compared with the predictions of linear stability theory was reported in [31]. Surprisingly, the issue of non-modal growth of disturbances at short times was addressed only by [14] for non-stratified mixing layers and by [32] for a stratified inviscid model. The study [32] described the non-modal amplification of the stratified mixing layer in terms of the interaction of gravity wave and vortex energy of the optimal perturbation for large wavenumbers. It was shown that the amplified transient growth can be three-dimensional and is associated with wave generation due to a vertical motion at large wavenumbers. two monotonic KH modes, one of which can become unstable. With further increase of the Richardson number, the pair of monotonic modes are replaced by a pair corresponding to complex conjugated eigenvalues related to the Holmboe instability. It was shown recently that KH instability plays an important role in the formation of interfaces between coating and substrate materials and metallic materials [24,25]. The non-linear evolution of KH and Holmboe instabilities was studied in quite a large number of articles (e.g., [17,26–33]). A difference in the nonlinear disturbances growth compared with the predictions of linear stability theory was reported in [31]. Surprisingly, the issue of non-modal growth of disturbances at short times was addressed only by [14] for non-stratified mixing layers and by [32] for a stratified inviscid model. The study [32] described the non-modal amplification of the stratified mixing layer in terms of the interaction of gravity wave and vortex energy of the optimal perturbation for large wavenumbers. It was shown that the amplified transient growth can be three-dimensional and is associated with wave generation due to a vertical motion at large wavenumbers. tions [10,17]. Furthermore, numerical calculations [18–20] supported by experiments of [21] showed that a three-dimensional mode traveling at an angle to the base velocity may attain the largest growth rate. For example, a shear layer characterized by a relatively thin region of stably stratified fluid can exhibit Holmboe instability [16,21,22], whose growth rate increases with increasing stratification. Despite the valid Squire transformation proved in [23], the flow is predicted to be linearly unstable to a three-dimensional disturb-It was shown in [10] that the instability at small Richardson numbers is governed by two monotonic KH modes, one of which can become unstable. With further increase of the Richardson number, the pair of monotonic modes are replaced by a pair corresponding to complex conjugated eigenvalues related to the Holmboe instability. It was shown recently that KH instability plays an important role in the formation of interfaces between coating and substrate materials and metallic materials [24,25]. The non-linear evolution of KH and Holmboe instabilities was studied in quite a large number of articles (e.g., [17,26–33]). A difference in the nonlinear disturbances growth compared with the predictions of linear

This study extends the previous results of [12] and [29] to viscous thermally stratified mixing layers. We show that the non-modal growth in this flow is governed mainly by the Holmboe-type modes. In the following, we examine their dynamics and interaction via consideration of a fully nonlinear three-dimensional model. This study extends the previous results of [12] and [29] to viscous thermally stratified mixing layers. We show that the non-modal growth in this flow is governed mainly by the Holmboe-type modes. In the following, we examine their dynamics and interaction via consideration of a fully nonlinear three-dimensional model. stability theory was reported in [31]. Surprisingly, the issue of non-modal growth of disturbances at short times was addressed only by [14] for non-stratified mixing layers and by [32] for a stratified inviscid model. The study [32] described the non-modal amplification of the stratified mixing layer in terms of the interaction of gravity wave and vortex

The paper is organized in the following way. In Section 2, we outline the fundamental details of the instability analysis and give a formulation of the optimal growth methodology, including the derivation of an adjoint operator. Section 3 discusses the implementation of different techniques for the calculation of the largest possible non-modal growth. Section 4 presents non-modal instability properties as a function of the Richardson number and evolution of the optimal disturbances within three-dimensional models. The summary and conclusions are derived in Section 5. The paper is organized in the following way. In Section 2, we outline the fundamental details of the instability analysis and give a formulation of the optimal growth methodology, including the derivation of an adjoint operator. Section 3 discusses the implementation of different techniques for the calculation of the largest possible non-modal growth. Section 4 presents non-modal instability properties as a function of the Richardson number and evolution of the optimal disturbances within three-dimensional models. The summary and conclusions are derived in Section 5. energy of the optimal perturbation for large wavenumbers. It was shown that the amplified transient growth can be three-dimensional and is associated with wave generation due to a vertical motion at large wavenumbers. This study extends the previous results of [12] and [29] to viscous thermally stratified mixing layers. We show that the non-modal growth in this flow is governed mainly by the Holmboe-type modes. In the following, we examine their dynamics and interaction via consideration of a fully nonlinear three-dimensional model.

#### **2. Formulation of the Problem and Numerical Techniques 2. Formulation of the Problem and Numerical Techniques** The paper is organized in the following way. In Section 2, we outline the fundamental

We consider a flow of an incompressible Newtonian fluid in a thermally stratified mixing layer. After the Boussinesq approximation is applied, the flow is governed by the momentum, continuity, and energy equations We consider a flow of an incompressible Newtonian fluid in a thermally stratified mixing layer. After the Boussinesq approximation is applied, the flow is governed by the momentum, continuity, and energy equations details of the instability analysis and give a formulation of the optimal growth methodology, including the derivation of an adjoint operator. Section 3 discusses the implementation of different techniques for the calculation of the largest possible non-modal growth. Section 4 presents non-modal instability properties as a function of the Richardson num-

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho} \nabla p + \mathbf{v} \Delta \mathbf{u} + \theta \gamma (T - \overline{T}) \mathbf{e}\_z \tag{1a}$$

$$\nabla \cdot \mathbf{u} = 0 \tag{1b}$$

$$\frac{\partial T}{\partial t} + (\mathbf{u} \cdot \nabla)T = \kappa \Delta T \tag{1c}$$
 
$$\text{Then } \mathbf{u} = \begin{pmatrix} u & u & u \end{pmatrix} \text{ is the polaristic weight component is the transmission (v) corresponding to the } \mathbf{u} \text{ and } \mathbf{u} \text{ coordinates.}$$

Here, = (, , ) is the velocity with components in the streamwise ( ), spanwise () and vertical () directions; and are the pressure and the temperature, is the density, ν is the kinematic viscosity, ℊ is gravitational acceleration, is the thermal expansion acceleration, is the unit vector in the -direction (vertical direction), ̅ Here, **u** = *ux*, *uy*, *u<sup>z</sup>* is the velocity with components in the streamwise (*x*), spanwise (*y*) and vertical (*z*) directions; *p* and *T* are the pressure and the temperature, *ρ* is the density, ν is the kinematic viscosity, momentum, continuity, and energy equations <sup>+</sup> ( <sup>∙</sup> ) <sup>=</sup> <sup>−</sup> 1 + ν∆ + ℊ( − ̅) (1a) ∇ ∙ = 0 (1b) is gravitational acceleration, *γ* is the thermal expansion acceleration, *e<sup>z</sup>* is the unit vector in the *z*-direction (vertical direction), *T* is the mean temperature, which is defined below, *κ* is the thermal diffusivity, and ∆ denotes the Laplacian operator.

is the mean temperature, which is defined below, is the thermal diffusivity, and ∆ denotes the Laplacian operator. The stratified mixing layer flow is modeled as two horizontal fluid layers moving in opposite directions with velocities ±*Umax* when the colder fluid layer is located under the

<sup>+</sup> ( <sup>∙</sup> ) <sup>=</sup> ∆ (1c)

Here, = (, , ) is the velocity with components in the streamwise ( ), spanwise () and vertical () directions; and are the pressure and the temperature, is the density, ν is the kinematic viscosity, ℊ is gravitational acceleration, is the ther-

is the mean temperature, which is defined below, is the thermal diffusivity, and ∆ de-

warmer one. It is assumed that the base flow and temperature profiles can be described by hyperbolic tangent profiles

$$\mathcal{U}(z) = \mathcal{U}\_{\text{max}} \tanh\left(\frac{z}{\delta\_{\upsilon}}\right),\\T(z) = T\_1 + \overline{T} \left[1 + \tanh\left(\frac{z}{\delta\_{T}}\right)\right] \tag{2}$$

where *T*<sup>1</sup> and *T*<sup>2</sup> are temperatures of lower and upper fluid layers, respectively, (*T*<sup>2</sup> > *T*1), *T* = 0.5(*T*<sup>2</sup> − *T*1), and *δ<sup>v</sup>* and *δ<sup>T</sup>* are the thicknesses of the velocity and temperature layers (note *δ<sup>v</sup>* > *δ<sup>T</sup>* for a Prandtl number larger than unity). The thickness of the velocity layer is defined following assumptions derived in [30] where the authors showed that presentation of base flow based on the "*tanh*" profile provides a better fit to the experimental data and is proportional to momentum thickness. *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8)

> Seeking a solution as a linear combination of the steady base flow (*U*, *T*) and an infinitesimal disturbance where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by ∗ ∗ ∗ Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 = (¯ **v**, *θ*), we arrive at the dimensionless linearized equations for the velocity, ¯ **v**, and temperature, *θ*, perturbations:

M [

]

,

,

]

〈 ,

〈<sup>1</sup> ,

So that (, , ) = ‖‖

tion, is defined as ([38–40])

So that (,, ) = ‖‖

(ℒ)(0) = ()(0), yields

definite Gram matrix =

definite Gram matrix =

optimal initial vector [43].

optimal initial vector [43].

〈()(0), ()(0)〉 = max

(ℒ)(0) = ()(0), yields

〈()(0), ()(0)〉 = max

‖‖ <sup>2</sup> =

‖‖ <sup>2</sup> =

as the value of the 2-norm of the matrix ()

as the value of the 2-norm of the matrix ()

a given from the eigenvalues of the linearized problem.

a given from the eigenvalues of the linearized problem.

tion, is defined as ([38–40])

2

2

(0)≠0

〈(0), <sup>∗</sup>

〈(0), <sup>∗</sup>

(0)≠0

〈 ,

〈 ,

(, , ) = max

(0)≠0

(,,) = max

max 

where <sup>∗</sup>

(0)≠0

max 

where <sup>∗</sup>

〉 = 〈, M〉 = ∫ [

$$\frac{\partial \mathbf{\bar{v}}}{\partial t} + \left(\mathbf{\bar{v}} \cdot \nabla\right) \mathcal{U} + (\mathcal{U} \cdot \nabla) \mathbf{\bar{v}} = -\nabla p + \mathcal{R}e^{-1} \Delta \mathbf{\bar{v}} + \mathcal{R}i\overline{\theta}e\_2 \tag{3a}$$
 
$$\Box$$

$$
\nabla \cdot \mathbf{v} = 0 \tag{3b}
$$

] (9)

(10)

(12)

∗

 (0)

)] (8)

)] (8)

] (9)

(10)

()()

(12)

()()‖ =

1<sup>2</sup> ∗

] (9)

(10)

()()

()()‖ (11)

(12)

()()‖ =

$$\frac{\partial \overline{\theta}}{\partial t} + \left(\overline{\mathbf{v} \cdot \nabla}\right) T + (\mathcal{U} \cdot \nabla) \overline{\theta} = P e^{-1} \Delta \overline{\theta} \tag{3c}$$

 (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = Equations (3a–3c) are rendered dimensionless using the scales *δv*, *Umax*, *δv*/*Umax*, *ρU*<sup>2</sup> *max*, and (*T*<sup>2</sup> − *T*1) for length, velocity, time, pressure, and temperature, respectively. The Reynolds number is defined by *Re* = *Umaxδv*/*ν*, *Pr* = *ν*/*κ* is the Prandtl number, *Pe* <sup>=</sup> *RePr* <sup>=</sup> *<sup>U</sup>maxδv*/*<sup>κ</sup>* is the Péclet number, *Ri* <sup>=</sup> *<sup>γ</sup>g*(*T*<sup>2</sup> <sup>−</sup> *<sup>T</sup>*1)*δv*/*U*<sup>2</sup> *max* is the bulk Richardson number, and *δ* = *δv*/*δ<sup>T</sup>* is the velocity and temperature thicknesses ratio. *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〉 = ∫ [ 1 (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + 1<sup>2</sup> ∗ )] (8) *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〈<sup>1</sup> , 〉 = ∫ [ 1 (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + 1<sup>2</sup> ∗ )] (8) *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value de-For the modal analysis, we assume an exponential dependence of the perturbations on time and apply the Fourier integral expansion in the infinite *x-* and *y*-direction, i.e., e 2(<sup>2</sup> + 2) (0) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 ] (9) = ˇ 2(<sup>2</sup> + 2) (0) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 ] (9) (*z*) exp[*i*(*αx* + *βy*) + *λt*]. Subsequently, the linearized continuity and momentum equations can be rearranged into a system of Orr-Sommerfeld and linearized energy equations (see, e.g., [28]) that read 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> +

$$
\lambda \tilde{\mathbb{A}}(z) = \mathcal{L} \tilde{\mathbb{A}}(z), \tag{4a}
$$

where \* denotes the complex conjugate. In the following, the inner product of two com-

$$\mathcal{L} = \begin{pmatrix} \mathcal{L}\_{\text{OS}} & \mathcal{L}\_{\text{w}\theta} \\ \mathcal{L}\_{\text{Tw}} & \mathcal{L}\_{\text{S}q\_{\theta}} \end{pmatrix} \tag{4b}$$

where \* denotes the complex conjugate. In the following, the inner product of two com-

$$\mathcal{L}\_{\rm OS} = \Delta^{-1} \left( \dot{\mathbf{a}} \left( -\mathbf{U} \Delta + \mathbf{U}\_{\rm zz} \right) + \mathrm{Re}^{-1} \Delta^2 \right), \tag{4c}$$
 
$$\Delta^{-1} \quad \Box$$

$$\mathcal{L}\_{w\theta} = -\text{Ri}\left(a^2 + \beta^2\right)\Delta^{-1} \tag{4d}$$

$$\mathcal{L}\_{\theta} = \text{er} \tag{4d}$$

$$\mathcal{L}\_{\rm Tw} = -T\_{z\nu}^{\prime} \tag{4e}$$

(12)

<sup>2</sup> = ∑

,

composition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector

> ∗ 〈 ∗ , ,

> > <sup>∗</sup> ∗

<sup>∗</sup> = ∑

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem.The second approach is the iterative forward/backward time integration applying <sup>a</sup> random initial perturbation. The evolution forward is governed by the operator ℒ and is

(12)

∗ 〈 ∗ , ,

composition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector

〉

ℳ = ∑

<sup>∗</sup> ∗

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

, as follows

〉

,

, the maximum energy growth at time can be calculated

<sup>∗</sup>

ℳ = ∑

, as follows

,

, the maximum energy growth at time can be calculated

<sup>∗</sup>

−1 and is equal to the square of its

∗

. The latter is

∗

. The latter is

−1 and is equal to the square of its

$$\mathcal{L}\_{S\eta\_{\theta}} = Pe^{-1}\Delta - \text{ind}\,\mathrm{I}.\tag{46}$$

on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent ap- (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent ap-Here, *λ* is a complex time increment, *α* and *β* are real wavenumbers in the *x*− and *y*− directions respectively, L is the linearized generalized operator, L*OS* is the Orr-Sommerfeld operator that determines the evolution of perturbation of the vertical velocity component *w*; L*Tw* and L*w<sup>θ</sup>* are operators coupling the velocity perturbation *w* with temperature perturbations *θ*; L*Sq<sup>θ</sup>* is the Squire operator that determines the evolution of perturbation of temperature *θ*; *T<sup>z</sup>* and *Uzz* are the first and second derivatives of temperature and velocity profiles (Equation (2)), respectively. Note that the Laplacian operator reduces to ∆ = *<sup>∂</sup>* 2 *∂z* <sup>2</sup> − *α* <sup>2</sup> + *β* 2 , and ∆ <sup>−</sup><sup>1</sup> denotes the inverse Laplacian operator.

<sup>∗</sup>

,

decomposition of the matrix : ‖‖<sup>2</sup>

that involves forward/backward time integrations yields the correct results.

‖‖ <sup>2</sup> =

, the maximum energy growth at time can be calculated

as the value of the 2-norm of the matrix ()

,

ℳ = ∑

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

definite Gram matrix =

optimal initial vector [43].

, the maximum energy growth at time can be calculated

decomposition of the matrix : ‖‖<sup>2</sup>

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

> ‖‖ <sup>2</sup> =

as the value of the 2-norm of the matrix ()

<sup>∗</sup>

−1 and is equal to the square of its

−1 and is equal to the square of its

<sup>∗</sup> = ∑

<sup>2</sup> = ∑

,

∗

a given from the eigenvalues of the linearized problem.

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

. The latter is

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

followed by the integration backward governed by the adjoint operator, ℒ

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

. The latter is

∗

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

definite Gram matrix =

<sup>∗</sup> ∗

,

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

optimal initial vector [43].

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

optimal initial vector [43].

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

,

<sup>∗</sup> = ∑

‖‖

〈

〈 ,

Equations (4a–4f) define an eigenproblem for eigenvalues *λ* and eigenvectors e plex vectors and associates with the energy norm (7), (8) and is given by 〈 , 〉 = 〈, M〉 = ∫ [ Ω <sup>Ω</sup> with M = [ So that (, , ) = ‖‖ tion, is defined as ([38–40]) . The existence of an eigenvalue *λ* having a positive real part is a sufficient condition for linear instability. Therefore, the linear stability problem reduces to finding a Reynolds number at which the real part of the leading eigenvalue (i.e., the eigenvalue with the largest real part) changes its sign from negative to positive and then minimizing *Re* over all possible values of *Ri*, *α*, and *β*. An extended analysis of Squire's theorem [27,30] demonstrated that a three-dimensional unstable mode can become most unstable in strongly stratified viscous flow. *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

#### *2.1. Non-Modal Analysis*

(0)≠0

〈()(0), ()(0)〉 = max

(ℒ)(0) = ()(0), yields

(, , ) = max

*Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

〈(0), <sup>∗</sup>

〈 ,

(0)≠0

〈(0), <sup>∗</sup>

(0)≠0

<sup>2</sup> = ∑

max 

∗ 〈 ∗ , ,

<sup>2</sup> = ∑

∗

max 

(0)≠0

(, , ) = max

(0)

∗

(0)≠0

〈<sup>1</sup> ,

(ℒ)(0) = ()(0), yields

(0)≠0

max 

〈<sup>1</sup> ,

〈

〈()(0), ()(0)〉 = max

(ℒ)(0) = ()(0), yields

〈()(0), ()(0)〉 = max

tion, is defined as ([38–40])

decomposition of the matrix : ‖‖<sup>2</sup>

So that (, , ) = ‖‖

definite Gram matrix =

decomposition of the matrix : ‖‖<sup>2</sup>

2

definite Gram matrix =

(ℒ)(0) = ()(0), yields

optimal initial vector [43].

definite Gram matrix =

(0)≠0

(, , ) = max

max 

where <sup>∗</sup>

max 

(0)≠0

where <sup>∗</sup>

(, , ) = max

optimal initial vector [43].

optimal initial vector [43].

decomposition of the matrix : ‖‖<sup>2</sup>

〈 ,

〈 ,

optimal initial vector [43].

as the value of the 2-norm of the matrix ()

(ℒ)(0) = ()(0), yields

a given from the eigenvalues of the linearized problem.

decomposition of the matrix : ‖‖<sup>2</sup>

decomposition of the matrix : ‖‖<sup>2</sup>

definite Gram matrix =

definite Gram matrix =

optimal initial vector [43].

optimal initial vector [43].

‖‖ <sup>2</sup> =

tion, is defined as ([38–40])

〈 ,

〈 ,

(, , ) = max

(, , ) = max

(0)≠0

max 

where <sup>∗</sup>

(0)≠0

(, , ) = max

max 

where <sup>∗</sup>

〈 ,

(0)≠0

max 

where <sup>∗</sup>

(, , ) = max

Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent ap-Considering the non-modal disturbances growth, we are interested in a short-time behavior of small finite-amplitude perturbation governed by the linearized equations. We exploit the homogeneity of the *x*- and *y*-direction by assuming solutions in the form ˆ 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 0 0 ⁄(0) ] (9) So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) () = 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 0 0 ⁄(0) ] (9) So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) (*z*, *t*) exp(*iαx* + *iβy*). It can be shown that by applying the continuity equation and excluding the pressure, an arbitrary three-dimensional infinitesimal perturbation is described by 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by , 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 0 0 ⁄(0) ] (9) So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth func-= *ux*, *uy*, *uz*, *θ* . Furthermore, the two velocity components *u<sup>x</sup>* and *u<sup>y</sup>* can be replaced by the vertical vorticity component *η* = *∂uy*/*∂x* − *∂ux*/*∂y*. Thus, consideration of the perturbation vector (*w*, *η*, *θ*) *<sup>T</sup>* where *u<sup>z</sup>* = *w* is the vertical velocity component yields an equivalent formulation. The resulting linearized dynamical matrix operator L is expressed as *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〈1, 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) where \* denotes the complex conjugate. In the following, the inner product of two com-*Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + where \* denotes the complex conjugate. In the following, the inner product of two com-

∗

∗

$$\frac{\partial \mathfrak{A}}{\partial t} = \mathcal{L}\mathfrak{A}\_{\prime} \tag{5a}$$

$$\mathfrak{A}\_{\prime} \triangleq \begin{array}{c} \mathfrak{A} \\ \mathfrak{A} \end{array} \tag{5b}$$

*Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

1 2(<sup>2</sup> + 2)

> ∗ , ∗ , ∗ ] M [ ]

2

(0)≠0

1<sup>2</sup> ∗

 (0) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> 

where \* denotes the complex conjugate. In the following, the inner product of two com-

(, , ) = max

〈(0), <sup>∗</sup>

)] (8)

1<sup>2</sup> ∗

<sup>2</sup> = ∑

] (9)

,

<sup>∗</sup> = ∑

∗ 〈 ∗ , ,

] (9)

)] (8)

<sup>∗</sup> ∗

(10)

()()

(12)

()()‖ (11)

()()

()()‖ (11)

()()‖ =

(10)

(12)

()()

()()‖ =

()()

()()‖ =

(10)

(10)

∗

<sup>∗</sup>

<sup>∗</sup>

∗

. The latter is

∗

∗

. The latter is

. The latter is

−1 and is equal to the square of its

, the maximum energy growth at time can be calculated

<sup>∗</sup>

1<sup>2</sup> ∗

<sup>∗</sup> +

 (0)

 (0)

−Δ 0 0 0 1 0 0 0 ⁄(0)

−Δ 0 0 0 1 0 0 0 ⁄(0)

<sup>∗</sup> +

. The latter is

. The latter is

(12)

(12)

(12)

∗

−1 and is equal to the square of its

−1 and is equal to the square of its

()()‖ =

] (9)

()()‖ (11)

1<sup>2</sup> ∗

)] (8)

1<sup>2</sup> ∗ (0)≠0

() (0)

()()(0)〉 = ‖<sup>∗</sup>

〉

(10)

()()

)] (8)

] (9)

] (9)

)] (8)

()()‖ =

ℳ = ∑

, as follows

,

<sup>∗</sup>

−1 and is equal to the square of its

∗

. The latter is

<sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup>

<sup>∗</sup> +

is the energy amplification at time = , or growth func-

 (0)

−Δ 0 0 0 1 0 0 0 ⁄(0)

1<sup>2</sup> ∗

()()‖ (11)

)] (8)

] (9)

(10)

()()

(12)

()()‖ =

〈<sup>1</sup> ,

〉 = ∫ [

$$
\mathcal{L} = \begin{pmatrix}
\mathcal{L}\_{\mathcal{S}\mathcal{S}} & 0 & \mathcal{L}\_{\pi\mathcal{S}} \\
\mathcal{L}\_{\mathcal{S}\mathcal{S}} & \mathcal{L}\_{\mathcal{S}\mathcal{S}} & 0 \\
\mathcal{L}\_{\mathcal{S}\mathcal{S}} & 0 & \mathcal{L}\_{\mathcal{S}\mathcal{S}}
\end{pmatrix} \tag{5b}
$$

$$\stackrel{\circ}{\mathcal{L}}\_{\eta w} = -i\beta l \mathcal{L}\_z \tag{5c}$$

$$\mathcal{L}\_{S\eta} = R e^{-\frac{1}{2}} \Delta - \text{in} \, \mathcal{U} \tag{5d}$$

which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent ap-Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> (0)≠0 (0) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> In this case, the coupling operator L*η<sup>w</sup>* and the Squire operator governing the vorticity perturbation L*Sq* are added to the equations system (4a–4f). The non-modal analysis examines the maximum energy growth at a particular time period of the flow disturbance, 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 ] (9) (0), over all possible initial conditions in terms of energy norm, produced by a sum of kinetic and potential energies, and assuming that the temperature perturbation alters the latter. The squared energy norm is represented by the inner product (to be defined below): *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〈<sup>1</sup> , 〉 =∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup>

inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than can be connected with the Euclidean norm of the vector () defined by the eigenvector 〉 , as follows can be connected with the Euclidean norm of the vector () defined by the eigenvector 〈 ∗ , ,〉 , as follows proaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> = is the largest eigenvalue that corresponds to the maximum possible relative is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> = is the largest eigenvalue that corresponds to the maximum possible relative 0 0 ⁄(0) is the energy amplification at time = , or growth func-*E*(*t*) = plex vectors and associates with the energy norm (7), (8) and is given by 2 *<sup>E</sup>* = plex vectors and associates with the energy norm (7), (8) and is given by , plex vectors and associates with the energy norm (7), (8) and is given by *E* (6)

where \* denotes the complex conjugate. In the following, the inner product of two com-

where \* denotes the complex conjugate. In the following, the inner product of two com-

where \* denotes the complex conjugate. In the following, the inner product of two com-

a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent apgrowth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. tion, is defined as ([38–40]) (, , ) = max () (10) 〈 ,〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 0 0 ⁄(0) 〈 ,〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ 〈 , 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ For fixed wavenumbers, the integrals in *x*- and *y*- coordinates scale out of the problem, and the energy norm in terms of *w*, *η,* and *θ* becomes

∗

The second approach is the iterative forward/backward time integration applying a

, as follows

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

Substitution for () in (10), with assuming (0) = 1, and recalling that () =

decomposition of the matrix : ‖‖<sup>2</sup>

(0)≠0

, the maximum energy growth at time can be calculated

Substitution for () in (10), with assuming (0) = 1, and recalling that () =

,

optimal initial vector [43].

 = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is

To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector

() is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup>

definite Gram matrix =

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

ℳ = ∑

(, , ) = max

followed by the integration backward governed by the adjoint operator, ℒ derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive , the maximum energy growth at time can be calculated where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive , the maximum energy growth at time can be calculated ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> , ℳ = ∑ <sup>∗</sup> , (12) where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, proaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector ∗ To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector (0)≠0 (0) Substitution for () in (10), with assuming (0) = 1, and recalling that () = So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) () So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) () So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) () *E*(*t*, *α*, *β*) = Z " 1 2(*α* <sup>2</sup> + *β* 2) |*w*| <sup>2</sup> + *dw dz* 2 + |*η*| <sup>2</sup> <sup>+</sup> *<sup>C</sup>*|*θ*<sup>|</sup> 2 !#*dz* (7) *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

optimal initial vector [43]. To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results. *2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations* as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the specdecomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ 〈∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ , where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = = is the largest eigenvalue that corresponds to the maximum possible relative (, , ) = max (0)≠0 (0) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) (, , ) = max (0)≠0 (0) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> (, , ) = max (0)≠0 (0) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> where the coefficient *C* is a positive definite scalar. The choice of *C* will not change the norm qualitatively [34,35]. It was shown in [36] that the perturbation energy is associated with the Brunt-Väisäla frequency and can be expressed in terms of the local Richardson number *Ri<sup>l</sup>* = *Ri Tz* (see also [37]). Following this result, we define *C* = *Ril*(*z* = 0) = *Ri Tz* (*z*=0) . It is easy to see that the above energy norm is produced by the inner product 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) 〈1, 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8)

the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a trum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent apwhere <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative h where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by −Δ 0 0 1, where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by −Δ 0 0 <sup>2</sup>i = Z 1 2(*α* <sup>2</sup> + *β* 2) *w*1*w* ∗ <sup>2</sup> + *dw*<sup>1</sup> *dz dw*∗ 2 *dz* <sup>+</sup> *<sup>η</sup>*1*<sup>η</sup>* ∗ <sup>2</sup> + *Ri Tz*(0) *θ*1*θ* ∗ 2 *dz* (8)

a given from the eigenvalues of the linearized problem.

() (0)

−1 and is equal to the square of its

decomposition of the matrix : ‖‖22 = ∑

decomposition of the matrix : ‖‖<sup>2</sup>

<sup>∗</sup>

() (0)

(0)≠0

‖‖ <sup>2</sup> =

()()(0)〉 = ‖<sup>∗</sup>

()()(0)〉 = ‖<sup>∗</sup>

as the value of the 2-norm of the matrix ()

definite Gram matrix =

a given from the eigenvalues of the linearized problem.

, the maximum energy growth at time can be calculated

ℳ = ∑

〉

, as follows

,

ℳ = ∑∗ ,

, the maximum energy growth at time can be calculated

∗

a given from the eigenvalues of the linearized problem.

‖‖ <sup>2</sup> =

. The latter is

0 1 0 0 0 ⁄(0)

<sup>2</sup> = ∑

,

. The latter is

on a separated part of the spectrum. Following the arguments of [13], only a discrete part of

growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is

] (9)

inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

,

()()‖ (11)

<sup>∗</sup> = ∑

on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spec-

growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is

inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

proaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

∗

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

∗ 〈 ∗ , ,

> <sup>∗</sup> ∗

()()‖ (11)

<sup>∗</sup> = ∑

(12)

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

that involves forward/backward time integrations yields the correct results.

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

−1 and is equal to the square of its

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

∗

∗

. The latter is

. The latter is

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

−1 and is equal to the square of its

followed by the integration backward governed by the adjoint operator, ℒ

followed by the integration backward governed by the adjoint operator, ℒ

. The latter is

To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector

To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector

growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is

] (9)

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

∗ 〈 ∗ , ,

ℳ = ∑

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

()()‖ =

()()

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

<sup>∗</sup> ∗

∗ 〈 ∗ , ,

(10)

<sup>∗</sup> ∗

, as follows

(10)

〉

,

, the maximum energy growth at time can be calculated

()()

<sup>∗</sup>

ℳ = ∑

ℳ = ∑

〉

, as follows

, as follows

,

,

, the maximum energy growth at time can be calculated

〉

<sup>2</sup> = ∑

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

()()‖ =

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

(12)

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

(12)

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

,

∗

the disturbance which yields this maximal growth value.

0 1 0 0 0 ⁄(0)

followed by the integration backward governed by the adjoint operator, ℒ

‖‖ <sup>2</sup> =

<sup>∗</sup> = ∑

that involves forward/backward time integrations yields the correct results.

as the value of the 2-norm of the matrix ()

as the value of the 2-norm of the matrix ()

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

. The latter is

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

∗

, as follows

followed by the integration backward governed by the adjoint operator, ℒ

a given from the eigenvalues of the linearized problem.

<sup>∗</sup>

a given from the eigenvalues of the linearized problem.

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

that involves forward/backward time integrations yields the correct results.

〈(0), <sup>∗</sup>

is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup>

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

(, , ) = max

〉

 ] Ω <sup>Ω</sup> with M = [

 ] Ω <sup>Ω</sup> with M = [

optimal initial vector [43].

is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup>

〈(0), <sup>∗</sup>

definite Gram matrix =

() is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup>

 = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is

To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector

> ∗ 〈 ∗ , ,

∗ 〈 ∗ , ,

> <sup>∗</sup> ∗

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

<sup>2</sup> = ∑

<sup>∗</sup> ∗

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

,

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

random initial perturbation. The evolution forward is governed by the operator ℒ and is

random initial perturbation. The evolution forward is governed by the operator ℒ and is

proaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation ()

∗ , ∗ , ∗ ] M [

followed by the integration backward governed by the adjoint operator, ℒ

2

a given from the eigenvalues of the linearized problem.

∗ , ∗ , ∗ ] M [

<sup>2</sup> = ∑

2

that involves forward/backward time integrations yields the correct results.

,

the spectrum and results in the growth function value at a given time and the corresponding

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

∗ 〈 ∗

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

<sup>∗</sup> ∗

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

(0)≠0

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

followed by the integration backward governed by the adjoint operator, ℒ

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

<sup>∗</sup> = ∑

optimal initial vector [43].

<sup>2</sup> = ∑

,

<sup>∗</sup> = ∑

optimal initial vector [43].

optimal initial vector [43].

〉

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

(0)≠0

followed by the integration backward governed by the adjoint operator, ℒ

〉 = 〈, M〉 = ∫ [

So that (, , ) = ‖‖

〉 = 〈, M〉 = ∫ [

So that (, , ) = ‖‖

tion, is defined as ([38–40])

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

followed by the integration backward governed by the adjoint operator, ℒ

the disturbance which yields this maximal growth value.

the disturbance which yields this maximal growth value.

that involves forward/backward time integrations yields the correct results.

as the value of the 2-norm of the matrix ()

as the value of the 2-norm of the matrix ()

a given from the eigenvalues of the linearized problem.

a given from the eigenvalues of the linearized problem.

‖‖ <sup>2</sup> =

‖‖ <sup>2</sup> =

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

<sup>∗</sup> = ∑

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

〈()(0), ()(0)〉 = max

〈()(0), ()(0)〉 = max

(ℒ)(0) = ()(0), yields

where \* denotes the complex conjugate. In the following, the inner product of two complex vectors where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by *<sup>i</sup>* and where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by *<sup>j</sup>* associates with the energy norm (7), (8) and is given by , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> ∗ )] (8) *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> <sup>∗</sup>)] (8)

<sup>∗</sup> +

<sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup>

 (0)

1<sup>2</sup> ∗

 (0)

<sup>∗</sup> +

)] (8)

proaches. The first one is the Gram matrix factorization, followed by the singular value de-

()()‖ =

,

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

. The latter is

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of

perturbed base flow as the initial condition. This validation shows that the second approach

optimal initial vector [43].

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

∗

. The latter is

that involves forward/backward time integrations yields the correct results.

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

∗

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

a given from the eigenvalues of the linearized problem.

)] (8)

(10)

 (0)

∗ , ∗ , ∗ ] M [ ] 

2

1

(10)

1<sup>2</sup> ∗

> (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup>

)] (8)

] (9)

(10)

() (0)

()()(0)〉 = ‖<sup>∗</sup>

(10)

<sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup>

<sup>∗</sup> +

is the energy amplification at time = , or growth func-

 (0)

−Δ 0 0 0 1 0 0 0 ⁄(0)

1<sup>2</sup> ∗

()()‖ (11)

)] (8)

] (9)

(10)

()()

(12)

()()‖ =

()()

] (9)

(0)≠0

)] (8)

(12)

()()‖ =

()()

〉

(12)

ℳ = ∑

, as follows

,

<sup>∗</sup>

−1 and is equal to the square of its

∗

. The latter is

()()‖ =

(, , ) = max

〈(0), <sup>∗</sup>

1<sup>2</sup> ∗

<sup>2</sup> = ∑

()()‖ (11)

,

<sup>∗</sup> = ∑

∗ 〈 ∗ , ,

> <sup>∗</sup> ∗

(12)

<sup>∗</sup>

()()‖ =

()

(0)≠0

(12)

∗

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

∗

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

. The latter is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

. The latter is

random initial perturbation. The evolution forward is governed by the operator ℒ and is

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

1<sup>2</sup> ∗

<sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup>

*Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

(1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> 

()

()

*Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

1 2(<sup>2</sup> + 2)

〉 = ∫ [

(0)≠0

1 2(<sup>2</sup> + 2)

〉 = ∫ [

(0)≠0

(1<sup>2</sup>

1

(1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> 

1 2(<sup>2</sup> + 2)

*Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

〈<sup>1</sup> ,

〈

〈<sup>1</sup>

〈 ,

〈,

(, , ) = max

(0)≠0

max 

(, , ) = max

where <sup>∗</sup>

max 

(, , ) = max

(, , ) = max

〈 ,

(0)≠0

(, , ) = max

(, , ) = max

max 

where <sup>∗</sup>

(0)≠0

max 

where <sup>∗</sup>

(0)≠0

where <sup>∗</sup>

〈 ,

(0)≠0

tion, is defined as ([38–40])

tion, is defined as ([38–40])

〈<sup>1</sup> ,

〈

max 

(, , ) =max

max 

(0)≠0

where <sup>∗</sup>

definite Gram matrix =

max 

optimal initial vector [43].

definite Gram matrix =

‖‖ <sup>2</sup> =

‖‖

‖‖ <sup>2</sup> =

〉 = ∫ [

〈<sup>1</sup> ,

tion, is defined as ([38–40])

〈 ,

tion, is defined as ([38–40])

〈<sup>1</sup> ,

〉 = ∫ [

$$\begin{aligned} \langle \Psi\_i, \Psi\_j \rangle\_{\Xi} &= \langle \Psi\_i, \Psi \Omega \rangle = \left\langle \begin{bmatrix} m\_i^\*, \psi\_i^\*, \psi\_i^\* \end{bmatrix} \right\rangle \mathbf{M} \begin{bmatrix} m\_i \\ \psi\_i \\ \psi\_i \end{bmatrix}\_{\downarrow} \text{d}\Omega \text{ with } \mathbf{M} = \begin{bmatrix} -\Delta & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \Omega I(T\_2(0)) \end{bmatrix} \end{aligned} \tag{9}$$

tion, is defined as ([38–40]) (, , ) = max () tion, is defined as ([38–40]) (,, ) = max () 0 0 ⁄(0) So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth func- 0 0 ⁄(0) So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth func- 0 0 ⁄(0) So that (, , ) = ‖‖ 2 00 ⁄(0) So that (, ,) = ‖‖ 2 So that *E*(*t*, *α*, *β*) = plex vectors and associates with the energy norm (7), (8) and is given by −Δ 0 0 2 *E* is the energy amplification at time *t* = *τ*, or growth function, is defined as ([38–40]) 〈<sup>1</sup> ,〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 0 0 ⁄(0) So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth func-

$$\text{a. } \text{norm} \text{ as } (\text{"} \land \text{"} \land \text{"} \land \text{"})$$

$$G(\mathbf{r}, a, \theta) = \max\_{\mathbf{Q}\_{(0)} \neq \mathbf{0}} \frac{E(\mathbf{r})}{E(\mathbf{0})} \tag{10}$$

 (ℒ)(0) = ()(0), yields 〈()(0), ()(0)〉 = max 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) (ℒ)(0) = ()(0), yields 〈()(0), ()(0)〉 = max 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) (, , ) = max (0)≠0 (0) (10) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (, , ) = max (0)≠0 (0) (10) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (, , ) = max (0)≠0 (0) (10) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (, , ) = max (0)≠0 (0) (10) Substitution for () in (10), with assuming (0) = 1, and recalling that () = So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) 〈 , 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ −Δ 0 0 0 1 0 0 0 ⁄(0) Substitution for *E*(*τ*) in (10), with assuming *E*(0) = 1, and recalling that plex vectors and associates with the energy norm (7), (8) and is given by (*τ*) = *exp*(L*τ*) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by (0) = Φ(*τ*) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by (0), yields (0)≠0 (0) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields *Fluids* **2021**, *6*, x FOR PEER REVIEW 5 of 17

$$\mathcal{G}(\mathbf{r}, \mathfrak{a}, \boldsymbol{\beta}) = \max\_{\mathfrak{q}(0) \neq 0} \langle \boldsymbol{\Phi}(\mathbf{r}) \mathfrak{q}(0), \boldsymbol{\Phi}(\mathbf{r}) \mathfrak{q}(0) \rangle = \max\_{\mathfrak{q}(0) \neq \mathfrak{q}} \langle \mathfrak{q}(0), \boldsymbol{\Phi}'(\mathbf{r}) \mathfrak{q}(\mathbf{r}) \mathfrak{q}(0) \rangle = \| \boldsymbol{\Phi}'(\mathbf{r}) \boldsymbol{\Phi}(\mathbf{r}) \| \tag{11}$$

()

()

 = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value de- = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value de-(0)≠0 () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is (0)≠0 () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is (0)≠0 where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is (0)≠0 (0)≠0 where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = (, , ) = max (0)≠0 () (0) Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields So that (, , ) = ‖‖ tion, is defined as ([38–40]) Substitution for () in (10), with assuming (0) = 1, and recalling that () = So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) (, , ) = max (0)≠0 () (0) (10) So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) (, , ) = max (0)≠0 () (0) (10) where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent apwhere Φ∗ (*τ*) is the adjoint of the operator Φ(*τ*) (see Section 2.2). The operator Φ∗ (*τ*)Φ(*τ*) is self-adjoint and normal, its eigenvalues are real and non-negative, and ||Φ<sup>∗</sup> (*τ*)Φ(*τ*)|| = max *i σ<sup>i</sup>* = *σ<sup>m</sup>* is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector 〈<sup>1</sup> , 〉 = ∫ [ 1 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> <sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by −Δ0 0 *<sup>m</sup>* corresponding to the eigenvalue *σ<sup>m</sup>* is the disturbance which yields this maximal growth value.

composition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) composition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ 〈 ∗ , , 〉 , as follows the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖22 = ∑ ∗ 〈 ∗ , , 〉 , as follows the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value de-(, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() Substitution for () in (10), with assuming (0) = 1, and recalling that () = (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0),()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> ()()‖ (11) where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() proaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ ℳ = ∑ <sup>∗</sup> 〈 , 〉 = 〈, M〉 = ∫ [ ∗ , ∗ , ∗ ] M [ ] Ω <sup>Ω</sup> with M = [ 0 1 0 0 0 ⁄(0) So that (, , ) = ‖‖ 2 is the energy amplification at time = , or growth function, is defined as ([38–40]) To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation *q*(*t*) can be connected with the Euclidean norm of the vector *k*(*t*) defined by the eigenvector decomposition of the matrix <sup>L</sup>: ||*q*||<sup>2</sup> <sup>2</sup> = ∑*i*,*<sup>j</sup> k* ∗ *i* h*p* ∗ *j* , *pi*i*k<sup>j</sup>* , as follows

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive , the maximum energy growth at time can be calculated where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive , the maximum energy growth at time can be calculated <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup>∗ , ℳ = ∑ <sup>∗</sup> , (12) composition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent ap- = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated (, , ) = max (0)≠0 (0) Substitution for () in (10), with assuming (0) = 1, and recalling that () = ||*q*||<sup>2</sup> *<sup>E</sup>* = *q* <sup>∗</sup>M*<sup>q</sup>* = ∑ *i*,*j k* ∗ *i p* ∗ *<sup>j</sup>* M*ji <sup>p</sup>ik<sup>j</sup>* = ∑ *i*,*j k* ∗ *i Sk<sup>j</sup>* (12)

,

max

as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. where is the mass matrix, *p* is the eigenvector, and =〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude composition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> ‖‖ <sup>2</sup> = where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is derived in the next section. This approach involves both discrete and continuous parts of (ℒ)(0) = ()(0), yields (, , ) = max (0)≠0 〈()(0), ()(0)〉 = max (0)≠0 〈(0), <sup>∗</sup> ()()(0)〉 = ‖<sup>∗</sup> where <sup>∗</sup> () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () where M is the mass matrix, *p* is the eigenvector, and *S* = h*p*,M*p*i is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix *S*= *F <sup>H</sup>F*, the maximum energy growth at time *τ* can be calculated as the value of the 2-norm of the matrix *Fexp*(Λ*τ*)*F* <sup>−</sup><sup>1</sup> and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given *<sup>ε</sup>* from the eigenvalues of the linearized problem.The second approach is the iterative forward/backward time integration applying a

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43]. To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results. derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43]. To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43]. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43]. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43]. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43]. inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ ∗ . The latter is on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than as the value of the 2-norm of the matrix () maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude definite Gram matrix = , the maximum energy growth at time can be calculated as the value of the 2-norm of the matrix () −1 and is equal to the square of its maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43]. To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results. *2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations* can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ , random initial perturbation. The evolution forward is governed by the operator L and is followed by the integration backward governed by the adjoint operator, L ∗ . The latter is derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43].

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations 2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations* To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results. To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results. To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results. To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results. derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding optimal initial vector [43]. a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is followed by the integration backward governed by the adjoint operator, ℒ the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. The second approach is the iterative forward/backward time integration applying a inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than a given from the eigenvalues of the linearized problem. where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive definite Gram matrix = , the maximum energy growth at time can be calculated −1 and is equal to the square of its To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach that involves forward/backward time integrations yields the correct results.

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

that involves forward/backward time integrations yields the correct results.

optimal initial vector [43].

optimal initial vector [43].

followed by the integration backward governed by the adjoint operator, ℒ

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

as the value of the 2-norm of the matrix ()

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

optimal initial vector [43].

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

#### *2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

*Fluids* **2021**, *6*, x FOR PEER REVIEW 6 of 17

We derive the adjoint operator considering the derivation of the adjoint form of the first- and second-order spatial derivatives and the first-order time derivatives. As noted by [44], the direct and adjoint parabolic problems have opposite directions of stable timelike evolution. From the definition (5a,b) and using (9) the adjoint operator L ∗ can be obtained as follows: We derive the adjoint operator considering the derivation of the adjoint form of the first- and second-order spatial derivatives and the first-order time derivatives. As noted by [44], the direct and adjoint parabolic problems have opposite directions of stable timelike evolution. From the definition (5a,b) and using (9) the adjoint operator ℒ ∗ can be obtained as follows:

$$\langle \mathfrak{q}, \mathcal{LM} \mathfrak{q} \rangle = \int\_{\Omega} \begin{bmatrix} -\Delta \left( \left( \mathcal{L}\_{OS} \mathcal{W} \right)^{\*} + \left( \mathcal{L}\_{\eta \mathcal{W}} \eta \right)^{\*} + \left( \mathcal{L}\_{Tw} \theta \right)^{\*} \right) \\\\ \left( \mathcal{L}\_{S\mathfrak{q}} \eta \right)^{\*} \\\\ \frac{Ri}{T\_{x}(\mathbb{O})} \left\{ \left( \mathcal{L}\_{w\theta} \eta \right)^{\*} + \left( \mathcal{L}\_{S\mathfrak{q}\_{\theta}} \theta \right)^{\*} \right\} \end{bmatrix}^{\mathrm{T}} \begin{bmatrix} \eta \\ \eta \\ \theta \end{bmatrix} d\Omega = \langle (\mathcal{LM}\mathfrak{q})^{\*}, \mathfrak{q} \rangle \tag{13}$$

Then, the adjoint system may be obtained by the integration by parts and reads 2 Then, the adjoint system may be obtained by the integration by parts and reads

$$
\Delta \left[ \Delta \frac{\partial}{\partial t} + 2Ll\_z \frac{\partial^2}{\partial x \partial z} + \left( lL \frac{\partial}{\partial x} + \frac{1}{Re} \Delta \right) \Delta \right] w^\* = -lL\_z \frac{\partial}{\partial y} \eta^\* - T^\prime \theta^\* \tag{14a}
$$

$$
\left[\frac{\partial}{\partial t} + \frac{1}{Re}\Delta + \mathcal{U}\frac{\partial}{\partial x}\right]\eta^\* = 0\tag{14b}
$$

1

, ∗ , ∗ ] M [ ] 

(0)≠0

the disturbance which yields this maximal growth value.

as the value of the 2-norm of the matrix ()

a given from the eigenvalues of the linearized problem.

definite Gram matrix =

optimal initial vector [43].

. The latter is

(1<sup>2</sup> <sup>∗</sup> + <sup>1</sup> 

(, , ) = max

〈(0), <sup>∗</sup>

<sup>2</sup> = ∑

,

<sup>∗</sup> = ∑

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

∗ 〈 ∗ , ,

> <sup>∗</sup> ∗

which is the complex-valued Hermitian matrix. Then, using factorization of the positive

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

(0)≠0

() (0)

()()(0)〉 = ‖<sup>∗</sup>

〉

ℳ = ∑

, as follows

,

, the maximum energy growth at time can be calculated

<sup>∗</sup>

−1 and is equal to the square of its

∗

. The latter is

<sup>2</sup> ∗ <sup>+</sup> 1<sup>2</sup>

<sup>∗</sup> +

is the energy amplification at time = , or growth func-

 (0)

−Δ 0 0 0 1 0 0 0 ⁄(0)

1<sup>2</sup> ∗

()()‖ (11)

)] (8)

] (9)

(10)

()()

(12)

()()‖ =

$$\left[\frac{\partial}{\partial t} + \frac{1}{Pe}\Delta + \mathcal{U}\frac{\partial}{\partial \mathbf{x}}\right]\theta^\* = -\Delta\_{2D}Ric^\* \tag{14c}$$

where ∆2= 2 <sup>2</sup> + 2 2 is the 2D Laplacian and is the first derivative of the velocity profile (Equation (2)). Evidently, the adjoint operator ℒ <sup>∗</sup> of the adjoint system <sup>∗</sup>⁄ = −(ℒ) ∗ in terms of a periodical in the stream- and span-wise direction perturbation will be expressed as where ∆2*<sup>D</sup>* = *<sup>∂</sup>* 2 *∂x* <sup>2</sup> + *<sup>∂</sup>* 2 *∂y* 2 is the 2D Laplacian and *U<sup>z</sup>* is the first derivative of the velocity profile (Equation (2)). Evidently, the adjoint operator L ∗ of the adjoint system *∂* where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by 〈 , 〉 = 〈, M〉 = ∫ [ ∗ Ω <sup>Ω</sup> with M = [ ∗ /*∂t* = − L , 〉 = ∫ [ 2(<sup>2</sup> + 2) (1<sup>2</sup> <sup>∗</sup> + <sup>+</sup> 1<sup>2</sup> <sup>∗</sup> + (0) 1<sup>2</sup> )] (8) where \* denotes the complex conjugate. In the following, the inner product of two complex vectors and associates with the energy norm (7), (8) and is given by −Δ 0 0 ∗ in terms of a periodical in the stream- and span-wise direction perturbation will be expressed as

$$
\mathcal{L}^\* = \begin{bmatrix}
\mathcal{L}\_{\text{OS}}^\* & \mathcal{L}\_{\eta w}^\* & \mathcal{L}\_{Tw}^\* \\
0 & \mathcal{L}\_{\text{S}q}^\* & 0 \\
\mathcal{L}\_{w\theta}^\* & 0 & \mathcal{L}\_{\text{S}q\theta}^\*
\end{bmatrix} \tag{15}
$$

where the entries of the matrices are tion, is defined as ([38–40]) where the entries of the matrices are

 ] 

1

∗ , ∗ , ∗ ] M [

2

〈<sup>1</sup>

〉 = 〈, M〉 = ∫ [

So that (, , ) = ‖‖

definite Gram matrix =

optimal initial vector [43].

〈 ,

(, , ) = max

(0)≠0

max 

where <sup>∗</sup>

$$\mathcal{L}\_{\text{CS}}^{\*} = \Delta^{-1} \left( -2i\text{it} \mathcal{L}\_{\text{\tilde{z}}} \frac{\partial}{\partial z} - \left( i\text{t} \mathcal{U} + \text{Re}^{-1} \Delta \right) \Delta \right), \quad \mathcal{L}\_{\text{\tilde{z}}\theta}^{\*} = -\text{Ri} \left( a^2 + \beta^2 \right) \tag{16a}$$

$$\mathcal{L}\_{\eta\mu}^\* = -\Delta^{-1} i\mathcal{J}\mathcal{U}\_{\Sigma} \, \qquad \qquad \mathcal{L}\_{\mathbb{S}\eta}^\* = -\mathrm{Re}^{-1}\Delta - i\mathrm{n}\mathcal{U} \, \tag{16b}$$

$$\mathcal{L}\_{\text{Tur}}^{\*} = -\Delta^{-1} \Gamma\_{\text{D}} \, \qquad \qquad \qquad \mathcal{L}\_{\text{Syl}}^{\*} = -\text{Pe}^{-1} \Delta - \text{ind} \, \text{L} \tag{16c}$$

−1 and is equal to the square of its

∗

, ℒ Finally, in accordance with Equation (11), the time integration forward/backward iterative procedure will result in the optimal initial vector, which maximizes growth at the is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> max = is the largest eigenvalue that corresponds to the maximum possible relative growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is () is the adjoint of the operator () (see Section 2.2). The operator <sup>∗</sup> ()() is self-adjoint and normal, its eigenvalues are real and non-negative, and ‖<sup>∗</sup> ()()‖ = = is the largest eigenvalue that corresponds to the maximum possible relative Finally, in accordance with Equation (11), the time integration forward/backward iterative procedure will result in the optimal initial vector, which maximizes growth at the specified time and a given set of parameters (*α*, *β*, *Re*, *Ri*, *Pe*) [43,45].

#### specified time and a given set of parameters (, , , , ) [43,45]. growth attainable at the time *τ*. The eigenvector corresponding to the eigenvalue is **3. Calculations Techniques**

a given from the eigenvalues of the linearized problem.

followed by the integration backward governed by the adjoint operator, ℒ

that involves forward/backward time integrations yields the correct results.

*2.2. Adjoint form of the Orr-Sommerfeld, Squire, and Energy Equations*

**3. Calculations Techniques** For the non-modal analysis, the linear system (Equation (5)) was discretized via a second-order central finite difference method. The eigenvalues and eigenvectors of matrix ℒ were calculated using the *QR* algorithm. The verification of the method and test calculation for the isothermal mixing layer can be found in [13]. Table 1 presents an example of the convergence of four leading eigenvalues belonging to the discrete part of the spectrum for a non-isothermal case. It is shown that the use of 1000 grid points yields four converged decimal digits for the first mode; however, the convergence slows down for the next modes. Despite this, it is shown that the growth function calculated by the iterative integration To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> ‖‖ <sup>2</sup> = where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, the disturbance which yields this maximal growth value. To examine possible non-modal disturbance growth, we applied two independent approaches. The first one is the Gram matrix factorization, followed by the singular value decomposition (SVD) method [41,42]. The energy norm (7) of an arbitrary perturbation () can be connected with the Euclidean norm of the vector () defined by the eigenvector decomposition of the matrix : ‖‖<sup>2</sup> <sup>2</sup> = ∑ ∗ 〈 ∗ , , 〉 , as follows ‖‖ <sup>2</sup> = <sup>∗</sup> = ∑ <sup>∗</sup> ∗ , ℳ = ∑ <sup>∗</sup> , (12) For the non-modal analysis, the linear system (Equation (5)) was discretized via a second-order central finite difference method. The eigenvalues and eigenvectors of matrix L were calculated using the *QR* algorithm. The verification of the method and test calculation for the isothermal mixing layer can be found in [13]. Table 1 presents an example of the convergence of four leading eigenvalues belonging to the discrete part of the spectrum for a non-isothermal case. It is shown that the use of 1000 grid points yields four converged decimal digits for the first mode; however, the convergence slows down for the next modes. Despite this, it is shown that the growth function calculated by the iterative integration method converges to the value of 21.71 already for 700 grid points.

The second approach is the iterative forward/backward time integration applying a random initial perturbation. The evolution forward is governed by the operator ℒ and is

To validate the results, we apply the full three-dimensional DNS with the optimally perturbed base flow as the initial condition. This validation shows that the second approach

derived in the next section. This approach involves both discrete and continuous parts of the spectrum and results in the growth function value at a given time and the corresponding

maximal singular value. This method allows for the calculation of non-modal growth based on a separated part of the spectrum. Following the arguments of [13], only a discrete part of the spectrum should be included in the analysis. To separate the discrete part of the spectrum, we tried the approach of [13] and extracted the eigenvectors that have finite amplitude inside and are near the mixing zone and vanish far from it. Alternatively, we included in the analysis only those modes whose pseudospectrum is located at a distance smaller than

where is the mass matrix, *p* is the eigenvector, and = 〈,〉 is the Gram matrix, which is the complex-valued Hermitian matrix. Then, using factorization of the positive

method converges to the value of 21.71 already for 700 grid points.


**Table 1.** Convergence of four least stable eigenvalues belonging to the discrete spectrum for stratified flow *α* = 0.5, *β* = 0, *Re* = 1000, *Ri* = 0.8.

The calculated spectrum was examined via analysis of its pseudospectrum [46]. The *ε*-pseudospectrum was computed as the minimal singular value of the matrix (*λI* − L), as was proposed in [47–49]. If *u<sup>s</sup>* is the singular vector of the matrix corresponding to the minimal singular value, then *ε* = ||(*λI* − L)*u<sup>s</sup>* ||2 estimates the accuracy of the calculation of the eigenvalue *λ*. It was found that for the stratified mixing layer configuration, the relatively large part of the eigenvalues corresponds to a relatively large *ε*-pseudospectrum (*ε* = 10−<sup>2</sup> , 10−<sup>3</sup> ), consequently, the eigenvalues calculation accuracy can be even more demanding than in the isothermal case. Results of [13] show that only the eigenmodes with the pseudospectrum *ε* < 10−<sup>6</sup> contributed to the optimal growth in isothermal flow. This is the main limitation of applying the SVD-based method [41] here. A series of tests show that we cannot use the SVD-based method based on a full spectrum for the same reasons as in the case of isothermal flow [10], i.e., we again have to separate part of the spectrum that corresponds to a non-accurate replication of the continuous modes. Figure 1 shows that the growth functions calculated based on extracting the vectors vanishing outside the shear layer (Figure 1a) and based on extracting the vectors corresponding different *ε* values (Figure 1b) do not yield converged consistent results within the SVD-based method. It is emphasized, however, that despite the obvious scatter in the results shown in Figure 1, they do show that significant non-modal growth can be expected even in cases of a strong, stable stratification. It also supports the results of [32], indicating that optimal perturbation can be located outside the shear layer.

We succeeded in obtaining reliable converged quantitative results using the iterative integration method, which, similarly to the isothermal case [13], converges to within three decimal places already on the 700 nodes grid. This method is used for the computation of the results presented below. It should be noted that this method, compared to the SVD-based one, requires much longer computational runs. To validate the non-modal analysis results, the direct numerical simulations taking the optimal vector as an initial condition were carried out. The fully non-linear time-dependent problem was solved using the approach described in [13].

*Re* = 1000, *Ri* = 0.8.

**1st Mode λr, λi**

**N**

**Table 1.** Convergence of four least stable eigenvalues belonging to the discrete spectrum for stratified flow *α* = 0.5, *β* = 0,

 0.0260 ±0.3653 −0.0460 ±0.2653 −0.1379 ±0.2710 −0.1466 ±0.0778 21.70 0.0334 ±0.3986 −0.0460 ±0.2653 −0.1361 ±0.3004 −0.1466 ±0.0684 21.70 0.0324 ±0.4591 −0.0460 ±0.2653 −0.1352 ±0.3017 −0.1466 ±0.0650 21.71 0.0324 ±0.4597 −0.0459 ±0.2654 −0.1334 ±0.3112 −0.1454 ±0.0583 21.71 0.0324 ±0.4597 −0.0459 ±0.2654 −0.1317 ±0.3310 −0.1497 ±0.0287 21.71 0.0324 ±0.4597 −0.0458 ±0.2654 −0.1250 ±0.3302 −0.1497 ±0.0259 21.71 0.0324 ±0.4597 −0.0458 ±0.2654 −0.1250 ±0.3302 −0.1496 ±0.0259 21.71 0.0324 ±0.4597 −0.0458 ±0.2653 −0.1262 ±0.3302 −0.1492 ±0.0256 21.71 0.0324 ±0.4597 −0.0457 ±0.2653 −0.1262 ±0.3302 −0.1489 ±0.0256 21.71 0.0324 ±0.4597 −0.0457 ±0.2652 −0.1262 ±0.3302 −0.1489 ±0.0255 21.71 0.0324 ±0.4597 −0.0457 ±0.2652 −0.1262 ±0.3302 −0.1489 ±0.0255 21.71

**3rd Mode λr, λi**

**4th Mode λr, λi**

The calculated spectrum was examined via analysis of its pseudospectrum [46]. The -pseudospectrum was computed as the minimal singular value of the matrix ( − ℒ),

), consequently, the eigenvalues calculation accuracy can be even more de-

minimal singular value, then = ‖( − ℒ)‖<sup>2</sup> estimates the accuracy of the calculation of the eigenvalue . It was found that for the stratified mixing layer configuration, the relatively large part of the eigenvalues corresponds to a relatively large -pseudospectrum

manding than in the isothermal case. Results of [13] show that only the eigenmodes with

is the main limitation of applying the SVD-based method [41] here. A series of tests show that we cannot use the SVD-based method based on a full spectrum for the same reasons as in the case of isothermal flow [10], i.e., we again have to separate part of the spectrum that corresponds to a non-accurate replication of the continuous modes. Figure 1 shows that the growth functions calculated based on extracting the vectors vanishing outside the shear layer (Figure 1a) and based on extracting the vectors corresponding different values (Figure 1b) do not yield converged consistent results within the SVD-based method. It is emphasized, however, that despite the obvious scatter in the results shown in Figure 1, they do show that significant non-modal growth can be expected even in cases of a strong, stable stratification. It also supports the results of [32], indicating that optimal perturbation can be

is the singular vector of the matrix corresponding to the

contributed to the optimal growth in isothermal flow. This

**Growth Function, G (***t* **= 10)**

**2nd Mode λr, λi**

as was proposed in [47–49]. If

, 10−3

the pseudospectrum < 10−6

located outside the shear layer.

( = 10−2

**Figure 1.** Comparison of growth functions calculated using the singular value decomposition (SVD) method based on extracted vectors (**a**) located in |*z*| < L, i.e., vanishing outside the shear zone and (**b**) corresponding to the values of their pseudospectrum, for *Ri* = 0.9, *Re* = 1000, *α* = 0.5, *β* = 0.

#### **4. Results**

In the following, we describe results on non-modal growth in the stratified mixing layer for *Pr* = 9, *Re* = 1000, and *R<sup>δ</sup>* = √ *Pr* = 3, assuming the fluid properties close to water, where heat diffuses much slower than momentum. The test calculations (Table 2) show that the least unstable eigenmode (or KH mode) is only slightly affected by the increase of the Prandtl number. The results of [48] also show that for a bounded channel flow, heat diffusivity does not affect linear stability. Conversely, the eigenvalues corresponding to the oscillatory Holmboe modes change noticeably with the increase of the Prandtl number, followed by a narrowing of the layer where the temperature varies.


**Table 2.** Least unstable eigenvalues at various Prandtl numbers for *α* = 0.5, *β* = 0, *Re* = 1000.

#### *4.1. Non-Modal Instability as a Function of the Richardson Number*

In this section, we present a characteristic case of non-modal energy growth for varying *Ri* and fixed *α*, *β*, *Pe*, and *R<sup>δ</sup>* . Results for the target time *t* = 10 are summarized in Table 3 (2D case) and Table 4 (3D case). It is seen from Table 2 that the flow is unstable to 2D perturbations (*α* = 0.5, *β* = 0, *Re* = 1000) for *Ri* ≤ 0.8. At small Richardson numbers (*Ri* < 0.1), the monotonic KH mode grows similarly to the isothermal case. With the increase of *Ri*, the two KH modes turn into a pair of Holmboe modes [10]. Further increase of the Richardson number, *Ri* > 0.8, makes the flow linearly stable; however, the non-modal growth remains possible. Although all the eigenvalues are negative, the Holmboe modes remain to be the least stable ones. An examination of the growth function, *G*, at a relatively short time, *t* = 10, shows that the highest amplification occurs at a small Richardson number, and the growth function decreases with the increase of *Ri*. Note that non-modal amplification reaches its value, for example, 76.4, for *Ri* = 0.1, while amplification of the least unstable mode at *t* = 10 is much smaller, approximately 1.4. This result may be important for the choice of the initial vector for non-linear calculation where the least unstable mode is usually

taken as the initial condition. Based on the foregoing results, we argue that at small times the non-modal growth is expected to dominate, i.e., the optimal vector will grow faster than just one least unstable eigenmode. With the steepening stratification, i.e., increase of the Richardson number, a contribution of two-dimensional disturbances into transient non-modal growth diminishes.


**Table 3.** The variations of the leading eigenvalues and values of amplification function, *G*, (at an arbitrary time *t* = 10) with increasing Richardson number for *α* = 0.5, *β* = 0, *Re* = 1000 (2D case).

**Table 4.** The leading eigenvalues, leading monotone and oscillatory (KH-, Holmboe-type) modes, and values of growth function, *G*, at *t* = 10 for the increasing Richardson number at *α* = 0.5, *β* = 1, *Re* = 1000 (3D case). The linearly unstable range is denoted by the gray color.


Based on the results for the isothermal mixing layer [13], where the largest amplification is yielded by 3D optimal disturbance, we examine a particular parameter set, *α* = 0.5, *β* = 1, *Re* = 1000, and vary only the Richardson number. These parameters characterize the perturbation that attains a large non-modal growth in the case of isothermal flow. Table 4 shows that the leading mode belongs to the continuous spectrum, Imag (*λ*) = ±*α*, so that the leading monotone (*λ* = 0, KH) and oscillatory (*λ* 6= 0, Holmboe) modes are presented in additional columns. It is seen that the oscillating Holmboe modes attain larger growth rates than the monotone ones. Computations of the growth function *G* (*t* = 10), using

the adjoint forward/backward time integration technique, show large non-modal growth of linearly stable 3D disturbances for *Ri* < 0.01. The growth function decreases with the increase of the Richardson number for *Ri* ≤ 0.2. Above the value *Ri* = 0.3, the growth function *G* (*t* = 10) starts to increase. We observe that the value of *G* (*t* = 10) continues to increase with increasing *Ri*.

It should be noted that for a relatively strong stratification, in the interval of 0.7 ≤ *Ri* ≤ 1, we found an additional linear instability. A similar result was reported in [19], where it was shown that stratified shear flow could be unstable to 3D disturbances propagating at an angle to the mean flow. Linear instability at large Richardson numbers for 2D perturbation was also reported in [37,50].

We emphasize that, similarly to the 2D case, the 3D non-modal growth at short times can be significantly larger than the modal one. The least unstable modes in the linearly unstable region are characterized by a minimal growth rate, making them almost neutral. For example, for *Ri* = 0.7 and *Ri* = 1, non-modal amplification reaches values 55 and 73.9, while exponential amplification of the least unstable mode at *t* = 10 is about unity in both cases. *Fluids* **2021**, *6*, x FOR PEER REVIEW 10 of 17

> Similar to the isothermal case, the values of the growth function in the 3D case are larger than in the 2D case for stratified flow so that the strongest amplification could be reached via 3D optimal perturbation. For example, for *Ri* = 0.7 and *Ri* = 1, non-modal amplification reaches values 55 and 73.9, while exponential amplification of the least unstable mode at *t* = 10 is about unity in both cases.

Similar to the isothermal case, the values of the growth function in the 3D case are

#### *4.2. Transient Dynamics of Optimal Perturbation* larger than in the 2D case for stratified flow so that the strongest amplification could be

We start exploring the optimal vector evolution from the two-dimensional case, *β* = 0, and consider *Ri* = 1, for which only non-modal transient growth exists. Figure 2 shows the amplitude of the initial optimal disturbances in terms of the temperature and *u*- and *w*- velocity components. We observe that amplitudes of the optimal disturbances are symmetric with respect to the mixing layer midline with the maxima of both of them shifted from the symmetry plane *z* = 0. These shifts indicate that the non-normal growth is governed mainly by the Holmboe modes, whose patterns are discussed in [10]. Contrary to the isothermal case [13], the optimal vector components are not localized only in the shear layer and are noticeably wider. This explains why the selection of vectors lying inside the shear layer is not justified for the stratified model considered here. reached via 3D optimal perturbation. *4.2. Transient Dynamics of Optimal Perturbation* We start exploring the optimal vector evolution from the two-dimensional case, *β* = 0, and consider *Ri* = 1, for which only non-modal transient growth exists. Figure 2 shows the amplitude of the initial optimal disturbances in terms of the temperature and *u*- and *w*- velocity components. We observe that amplitudes of the optimal disturbances are symmetric with respect to the mixing layer midline with the maxima of both of them shifted from the symmetry plane *z* = 0. These shifts indicate that the non-normal growth is governed mainly by the Holmboe modes, whose patterns are discussed in [10]. Contrary to the isothermal case [13], the optimal vector components are not localized only in the shear layer and are noticeably wider. This explains why the selection of vectors lying inside the shear layer is not justified for the stratified model considered here.

**Figure 2.** Amplitudes of the initial optimal disturbance for = 1, = 1000, = 0.5, = 0 yielding the maximal transient growth at = 10. **Figure 2.** Amplitudes of the initial optimal disturbance for *Ri* = 1, *Re* = 1000, *α* = 0.5, *β* = 0 yielding the maximal transient growth at *t* = 10.

To investigate the non-modal growth mechanism, we plot spatial patterns of optimal vector time evolution. Figure 3 shows that the *u*-velocity component's evolution is repre-

To investigate the non-modal growth mechanism, we plot spatial patterns of optimal vector time evolution. Figure 3 shows that the *u*-velocity component's evolution is represented by two waves traveling in the opposite direction to the mean flow. Such a behavior can be a replication of the well-known Holmboe instability mechanism [16]. *Fluids* **2021**, *6*, x FOR PEER REVIEW 11 of 17

**Figure 3.** Evolution of the streamwise velocity component*,* , of a two-dimensional optimal vector calculated for the maximal growth at = 10. = 1, = 1000, = 0.5, = 0. Solid and dashed lines represent positive and negative values of *u*. The plotted values are equally spaced between −2 and 2. **Figure 3.** Evolution of the streamwise velocity component, *u*, of a two-dimensional optimal vector calculated for the maximal growth at *tmax* = 10. *Ri* = 1, *Re* = 1000, *α* = 0.5, *β* = 0. Solid and dashed lines represent positive and negative values of *u*. The plotted values are equally spaced between −2 and 2.

component at four indicative time moments: initial state, growth, maximum at *t* = 10, and further decaying (Figure 4). Note that similar to the isothermal case [13,14], results for the stratified layer show that the initial optimal vector structures are elongated against the

Following [51,52] explanations of the linear stability dynamics using the "buoyancyvorticity wave interaction approach" in an inviscid flow, we presented the *ηy*− vorticity component at four indicative time moments: initial state, growth, maximum at *t* = 10, and further decaying (Figure 4). Note that similar to the isothermal case [13,14], results for the stratified layer show that the initial optimal vector structures are elongated against the shear in the streamwise direction. This allows us to attribute the initial amplification to the Orr mechanism [14,53]. With the increase of time, the perturbation develops into two series of vortices whose origins are located at the lines *z* = ±1, i.e., at the edges of the temperature stratification layer. At the same time, we observe the development of another vortex series with the origins located at the midline. These vortices attain larger amplitude and change their patterns from elongating against the shear to tilting with the shear. Furthermore, the maximal energy growth corresponds to the vertically positioned midline vortices generating vertical velocity. Hence, our calculations confirm the suggestion of [51,52] that vorticity growth occurs due to persistent vertical velocity produced by the propagated waves in the two strong density gradient regions. *Fluids* **2021**, *6*, x FOR PEER REVIEW 12 of 17 shear in the streamwise direction. This allows us to attribute the initial amplification to the Orr mechanism [14,53]. With the increase of time, the perturbation develops into two series of vortices whose origins are located at the lines z = ±1, i.e., at the edges of the temperature stratification layer. At the same time, we observe the development of another vortex series with the origins located at the midline. These vortices attain larger amplitude and change their patterns from elongating against the shear to tilting with the shear. Furthermore, the maximal energy growth corresponds to the vertically positioned midline vortices generating vertical velocity. Hence, our calculations confirm the suggestion of [51,52] that vorticity growth occurs due to persistent vertical velocity produced by the propagated waves in the two strong density gradient regions.

**Figure 4.** Evolution of the spanwise vorticity component, , of a two-dimensional optimal vector calculated for the maximal growth at = 10. = 1, = 1000, = 0.5, = 0. Solid and dashed lines represent positive and negative values of . **Figure 4.** Evolution of the spanwise vorticity component, *ηy*, of a two-dimensional optimal vector calculated for the maximal growth at *tmax* = 10. *Ri* = 1, *Re* = 1000, *α* = 0.5, *β* = 0. Solid and dashed lines represent positive and negative values of *ηy*.

Figure 5 illustrates snapshots of the temperature field where the base temperature profile and the 2D optimal disturbance are superposed. The time evolution of optimal temperature disturbance *θ* exhibits an extension of the initial temperature layer due to two waves propagating in opposite directions along the *x*-axis. Such perturbation leads to elevation of the cold fluid and moving down of the warm fluid, which is driven by the vertical velocity component. During the evolution, the perturbation extracts energy from the mean flow, followed by an increase of the growth function. After passing the maximal growth stage, the cold fluid is naturally descending, which corresponds to the growth function decrease. Owing to the small values of the Holmboe modes decay rates, these up-and-down motions of the fluid persist for a relatively long time. The evolution at larger times results in merged heated structures below and above the midplane, corresponding to an increase in the amplitude of the oscillations and an increase of kinetic and potential energy. At later times, the amplitude of the non-modal disturbance decreases, as expected. Figure 5 illustrates snapshots of the temperature field where the base temperature profile and the 2D optimal disturbance are superposed. The time evolution of optimal temperature disturbance *θ* exhibits an extension of the initial temperature layer due to two waves propagating in opposite directions along the *x*-axis. Such perturbation leads to elevation of the cold fluid and moving down of the warm fluid, which is driven by the vertical velocity component. During the evolution, the perturbation extracts energy from the mean flow, followed by an increase of the growth function. After passing the maximal growth stage, the cold fluid is naturally descending, which corresponds to the growth function decrease. Owing to the small values of the Holmboe modes decay rates, these up-and-down motions of the fluid persist for a relatively long time. The evolution at larger times results in merged heated structures below and above the midplane, corresponding to an increase in the amplitude of the oscillations and an increase of kinetic and potential energy. At later times, the amplitude of the non-modal disturbance decreases, as expected.

**Figure 5.** Snapshots of the base flow superposed with 2D optimal temperature perturbation calculated for the maximal growth at = 10. = 1, = 1000, = 0.5, = 0. The blue color corresponds to cold fluid; the red color corresponds to warm fluid. **Figure 5.** Snapshots of the base flow superposed with 2D optimal temperature perturbation calculated for the maximal growth at *tmax* = 10. *Ri* = 1, *Re* = 1000, *α* = 0.5, *β* = 0. The blue color corresponds to cold fluid; the red color corresponds to warm fluid. **Figure 5.** Snapshots of the base flow superposed with 2D optimal temperature perturbation calculated for the maximal growth at = 10. = 1, = 1000, = 0.5, = 0. The blue color corresponds to cold fluid; the red color corresponds to warm fluid.

Figure 6 illustrates the amplitudes of the optimal vector components for *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1, the parameters characteristic for 3D transient growth. Similar to the 2D case, the optimal vector amplitude has two symmetric maxima, indicating the dominance of the Holmboe-like modes in the initial optimal disturbance. It is shown that the spanwise vorticity component, , has a maximum also at *z* = 0, supporting the suggestions of [52]. The amplitude of all components of the 3D initial optimal vector decays beyond = ±10, making it narrower than the 2D optimal vector. Figure 6 illustrates the amplitudes of the optimal vector components for *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1, the parameters characteristic for 3D transient growth. Similar to the 2D case, the optimal vector amplitude has two symmetric maxima, indicating the dominance of the Holmboe-like modes in the initial optimal disturbance. It is shown that the spanwise vorticity component, *ηy*, has a maximum also at *z* = 0, supporting the suggestions of [52]. The amplitude of all components of the 3D initial optimal vector decays beyond *z* = ±10, making it narrower than the 2D optimal vector. Figure 6 illustrates the amplitudes of the optimal vector components for *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1, the parameters characteristic for 3D transient growth. Similar to the 2D case, the optimal vector amplitude has two symmetric maxima, indicating the dominance of the Holmboe-like modes in the initial optimal disturbance. It is shown that the spanwise vorticity component, , has a maximum also at *z* = 0, supporting the suggestions of [52]. The amplitude of all components of the 3D initial optimal vector decays beyond = ±10, making it narrower than the 2D optimal vector.

**Figure 6.** Amplitudes of the initial optimal disturbance for = 2, = 1000, = 0.5, = 1 yielding the maximal growth at = 10; (**a**) three components of velocity magnitude; (**b**) three vorticity components; (**c**) amplitude of the temperature disturbance. **Figure 6.** Amplitudes of the initial optimal disturbance for = 2, = 1000, = 0.5, = 1 yielding the maximal growth at = 10; (**a**) three components of velocity magnitude; (**b**) three vorticity components; (**c**) amplitude of the temperature disturbance. **Figure 6.** Amplitudes of the initial optimal disturbance for *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1 yielding the maximal growth at *tmax* = 10; (**a**) three components of velocity magnitude; (**b**) three vorticity components; (**c**) amplitude of the temperature disturbance.

Growth function of the three-dimensional optimal temperature perturbation and patterns of its time evolution are illustrated in Figure 7 and Figure 8, respectively, for the case *α* = 0.5, *β* = 1, *Re* = 1000, *Ri* = 2. It is shown that the growth function value at *t* = 10 is similar to those obtained using the forward/backward time integration method, *G* = 122 (Table 4), which justifies the results obtained using two approaches. As already mentioned, the growth rates of the Holmboe modes are low so that the optimal perturbation energy Growth function of the three-dimensional optimal temperature perturbation and patterns of its time evolution are illustrated in Figure 7 and Figure 8, respectively, for the case *α* = 0.5, *β* = 1, *Re* = 1000, *Ri* = 2. It is shown that the growth function value at *t* = 10 is similar to those obtained using the forward/backward time integration method, *G* = 122 (Table 4), which justifies the results obtained using two approaches. As already mentioned, the growth rates of the Holmboe modes are low so that the optimal perturbation energy Growth function of the three-dimensional optimal temperature perturbation and patterns of its time evolution are illustrated in Figures 7 and 8, respectively, for the case *α* = 0.5, *β* = 1, *Re* = 1000, *Ri* = 2. It is shown that the growth function value at *t* = 10 is similar to those obtained using the forward/backward time integration method, *G* = 122 (Table 4), which justifies the results obtained using two approaches. As already mentioned, the growth rates of the Holmboe modes are low so that the optimal perturbation energy slowly

decays at longer times. Similar to the 2D case, the three-dimensional structures exhibit two lines of vortices ordered along with the shear layer and tending to be tilted against the shear slope at the initial stage, constantly changing inclination at later times. At each time moment, corresponding to the maximum of the growth function, the temperature layer's width increases; however, it does not exceed the width observed in the 2D case. Figure 9 shows the evolution of the streamwise velocity at the same times as in Figure 7. The optimal perturbation is composed of two waves traveling in the opposite direction. The wave's amplitude grows and decreases, reflecting in the non-monotonic growth of kinetic energy at longer times. Figure 7. The optimal perturbation is composed of two waves traveling in the opposite direction. The wave's amplitude grows and decreases, reflecting in the non-monotonic growth of kinetic energy at longer times. Contrary to the 2D case, the spanwise vorticity component of initial optimal perturbation appears as three lines of vortices with the centers located in the middle plane and the = ±1 planes. When the growth function increases, e.g., at times *t* = 0, 5, 20, 35, the vortices are always oriented against the shear. The central vortex splits into two smaller vortices, moves from the midline above and below, and scrolls. At times *t* = 10, 25, at which the growth function attains a maximum, the vortices are strongly tilted with the shear, thus supporting the distinctions between 2D and 3D non-modal growths. Figure 7. The optimal perturbation is composed of two waves traveling in the opposite direction. The wave's amplitude grows and decreases, reflecting in the non-monotonic growth of kinetic energy at longer times. Contrary to the 2D case, the spanwise vorticity component of initial optimal perturbation appears as three lines of vortices with the centers located in the middle plane and the = ±1 planes. When the growth function increases, e.g., at times *t* = 0, 5, 20, 35, the vortices are always oriented against the shear. The central vortex splits into two smaller vortices, moves from the midline above and below, and scrolls. At times *t* = 10, 25, at which the growth function attains a maximum, the vortices are strongly tilted with the shear, thus

supporting the distinctions between 2D and 3D non-modal growths.

slowly decays at longer times. Similar to the 2D case, the three-dimensional structures exhibit two lines of vortices ordered along with the shear layer and tending to be tilted against the shear slope at the initial stage, constantly changing inclination at later times. At each time moment, corresponding to the maximum of the growth function, the temperature layer's width increases; however, it does not exceed the width observed in the 2D case. Figure 9 shows the evolution of the streamwise velocity at the same times as in

slowly decays at longer times. Similar to the 2D case, the three-dimensional structures exhibit two lines of vortices ordered along with the shear layer and tending to be tilted against the shear slope at the initial stage, constantly changing inclination at later times. At each time moment, corresponding to the maximum of the growth function, the temperature layer's width increases; however, it does not exceed the width observed in the 2D case. Figure 9 shows the evolution of the streamwise velocity at the same times as in

*Fluids* **2021**, *6*, x FOR PEER REVIEW 14 of 17

*Fluids* **2021**, *6*, x FOR PEER REVIEW 14 of 17

**Figure 7.** Growth function of the 3D optimal temperature perturbation superposed with base flow profile for *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1. The numbers represent different stages of increase and decrease in growth function. The optimal disturbances corresponding to some of these stages are shown in Figure 8. **Figure 7.** Growth function of the 3D optimal temperature perturbation superposed with base flow profile for *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1. The numbers represent different stages of increase and decrease in growth function. The optimal disturbances corresponding to some of these stages are shown in Figure 8. **Figure 7.** Growth function of the 3D optimal temperature perturbation superposed with base flow profile for *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1. The numbers represent different stages of increase and decrease in growth function. The optimal disturbances corresponding to some of these stages are shown in Figure 8.

**Figure 8.** Time evolution of the 3D optimal temperature perturbation superposed with base flow profile for = 2, = 1000, = 0.5, = 1, yielding the maximal values at = 10. The blue color corresponds to cold fluid; the red color **Figure 8.** Time evolution of the 3D optimal temperature perturbation superposed with base flow profile for = 2, = 1000, = 0.5, = 1, yielding the maximal values at = 10. The blue color corresponds to cold fluid; the red color **Figure 8.** Time evolution of the 3D optimal temperature perturbation superposed with base flow profile for *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1, yielding the maximal values at *tmax* = 10. The blue color corresponds to cold fluid; the red color corresponds to warm fluid. The growth function increases at *t* = 0, 5, 20, 35, attains maximal values at *t* = 10, 25, and decays at *t* = 15, 30.

decays at = 15, 30.

**Figure 9.** Time evolution of the streamwise velocity component of a three-dimensional optimal vector calculated for the maximal growth at *t* = 10, = 2, = 1000, = 0.5, = 1. The growth function increases at = 0, 5, 20, 35, attains maximal values at = 10, 25, and decays at = 15, 30. **Figure 9.** Time evolution of the streamwise velocity component of a three-dimensional optimal vector calculated for the maximal growth at *t* = 10, *Ri* = 2, *Re* = 1000, *α* = 0.5, *β* = 1. The growth function increases at *t* = 0, 5, 20, 35, attains maximal values at *t* = 10, 25, and decays at *t* = 15, 30.

**5. Conclusions** The non-modal disturbances growth and transient dynamics of optimal perturbation in stratified viscous mixing layers were investigated. The flow exhibits strong transient growth for streamwise and spanwise wavenumbers, at which the flow is either asymptotically or neutrally stable. Comparing the calculated flow structures with those observed in several previous experimental and numerical studies at the parameters corresponding Contrary to the 2D case, the spanwise vorticity component of initial optimal perturbation appears as three lines of vortices with the centers located in the middle plane and the *z* = ±1 planes. When the growth function increases, e.g., at times *t* = 0, 5, 20, 35, the vortices are always oriented against the shear. The central vortex splits into two smaller vortices, moves from the midline above and below, and scrolls. At times *t* = 10, 25, at which the growth function attains a maximum, the vortices are strongly tilted with the shear, thus supporting the distinctions between 2D and 3D non-modal growths.

to linearly unstable regimes, we conclude that the mixing layer flow at early stages of the

#### linear instability development can be strongly affected by the temporal disturbances **5. Conclusions**

growth when the optimal perturbation is localized inside the shear zone. The main findings of the performed study are as follows. The non-modal disturbances growth and transient dynamics of optimal perturbation in stratified viscous mixing layers were investigated. The flow exhibits strong transient growth for streamwise and spanwise wavenumbers, at which the flow is either asymptotically or neutrally stable. Comparing the calculated flow structures with those observed in several previous experimental and numerical studies at the parameters corresponding to linearly unstable regimes, we conclude that the mixing layer flow at early stages of

the linear instability development can be strongly affected by the temporal disturbances growth when the optimal perturbation is localized inside the shear zone. The main findings of the performed study are as follows.


In the case of stratified flow, we could not find a clear criterion to define which part of the spectrum should be taken into account to obtain a correct and numerically converged growth function. We could not obtain comparable results using SVD decompositionbased approaches for studying the non-modal growth, trying different extractions of the (seemingly) discrete part of the spectrum. This shows that the stratified mixing layer is a significantly more complicated problem, and the question about which part of the spectrum contributes to the non-modal growth is yet to be investigated.

**Author Contributions:** Conceptualization, methodology, validation, formal analysis, investigation, resources, data curation, H.V. and A.G.; writing—original draft preparation, H.V.; writing—review and editing, A.G.; visualization, H.V.; supervision, project administration, funding acquisition, A.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

