2.1 Problem set-up

The fluid is assumed to be contained in a rectangular domain of aspect ratio $A\equiv L_{x}/H$ , where $H$ is the height and $L_{x}$ the length ( $L_{y}$ will denote the width in the three-dimensional case). Using $H$ as a unit of length, the non-dimensional two-dimensional domain is $[-A/2,A/2]\times [0,1]$ . Here, similar to previous studies (Paparella & Young Reference Paparella and Young2002; Scotti & White Reference Scotti and White2011) $A=4$ will be used as a baseline geometry.

The flow is described by the incompressible Navier–Stokes equations in the Boussinesq limit with a linear equation of state. Under these assumptions, the density $\unicode[STIX]{x1D70C}$ is given by $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}(1+\unicode[STIX]{x1D6FC}T)$ , where $T$ is the quantity that describes the stratifying agent and $\unicode[STIX]{x1D70C}_{0}$ is a reference density. In experiments involving water, stratification is normally controlled either via changes in temperature, in which case $\unicode[STIX]{x1D6FC}$ is the (negative) coefficient of thermal expansion (Niemela et al. Reference Niemela, Skrbek, Sreenivasan and Donnelly2000; Roche et al. Reference Roche, Castaing, Chabaud and Hébral2002; Mullarney, Griffiths & Hughes Reference Mullarney, Griffiths and Hughes2004), or via varying the concentration of a salt (Stewart et al. Reference Stewart, Graham and Griffiths2012; Passaggia et al. Reference Passaggia, Hurley, Scotti and White2017), and $\unicode[STIX]{x1D6FC}$ is the positive haline contraction coefficient. Since, from a theoretical point of view, the precise nature of the stratifying agent is irrelevant, we describe stratification in terms of the buoyancy $b\equiv g(\unicode[STIX]{x1D70C}_{0}-\unicode[STIX]{x1D70C})/\unicode[STIX]{x1D70C}_{0}$ , where $g$ is the acceleration of gravity. $b$ satisfies the same advection–diffusion equation of the stratifying agent. Let $\unicode[STIX]{x0394}b$ be the buoyancy contrast imposed by the boundary conditions. We can define a velocity scale $\sqrt{\unicode[STIX]{x0394}bH}$ , which allows to write non-dimensionalize the governing equations as follows

where $\boldsymbol{u}=(u,v,w)^{\text{T}}$ and $p$ are the non-dimensional velocity vector and pressure. The Rayleigh number is defined by $Ra=\unicode[STIX]{x0394}bL^{3}/(\unicode[STIX]{x1D708}\unicode[STIX]{x1D705})$ and the Prandtl number is given by $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ where $\unicode[STIX]{x1D708}$ is the kinematic viscosity and $\unicode[STIX]{x1D705}$ is the diffusivity of the stratifying agent (calling $\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ the Prandtl number should not be construed as implying that we are thinking of heat as the stratifying agent). A value of $Pr=1$ was used throughout the present study. While it would be more realistic to use higher values, it would limit the range of Rayleigh numbers attainable. The unit vector $\boldsymbol{e}_{\boldsymbol{z}}$ points opposite the local gravity.

Figure 1.

Sketch of the flow set-up.

The base flow is assumed to be two-dimensional whereas the perturbation is sought to be periodic in the transverse direction $y$ . The bottom and the side walls are insulating whereas a piecewise constant buoyancy profile $B(x<0,1)=0$ and $B(x>0,1)=1$ is imposed at the upper boundary ( $z=1$ ), as shown in figure 1. The choice of such idealized boundary conditions for the buoyancy was essentially motivated by previous studies (Gayen et al. Reference Gayen, Griffiths, Hughes and Saenz2013, Reference Gayen, Griffiths and Hughes2014). The step profile here compared for instance to the linear profile provides more circulation as the amount of buoyancy provided for the circulation is larger for the step case than the linear case. More importantly, a step profile does not add an extra length scale to the problem, simplifying the analysis.

We consider both free- and no-slip boundary conditions for the velocity along the boundaries. The latter is required to model actual laboratory experiments, whereas the former is more appropriate if the aim is to extrapolate the results to geophysical situations, where no-slip boundary layers should play a negligible role. In the literature we find examples of both (Paparella & Young Reference Paparella and Young2002; Scotti & White Reference Scotti and White2011; Gayen et al. Reference Gayen, Griffiths and Hughes2014). The flow domain was discretized with $400$ points using fourth-order finite differences for the Laplacian and eighth-order differences for the gradient in the longitudinal direction and $60$ Chebyshev collocation points in the vertical direction.

2.2 Steady states method

Steady solutions and the associated global stability results were obtained across the range $6.4\times 10^{6}\leqslant Ra\leqslant 1.28\times 10^{10}$ for free-slip boundary conditions and $6.4\times 10^{6}\leqslant Ra\leqslant 1.28\times 10^{9}$ for the no-slip case. The choice was made to use a time stepper approach to compute both the steady states and the eigenvalue problem associated with the stability. The same procedure was used in Marquillie & Ehrenstein (Reference Marquillie and Ehrenstein2003) and Passaggia, Meunier & Le Dizès (Reference Passaggia, Meunier and Le Dizès2014) and details about the numerical methods can be found in Marquillie & Ehrenstein (Reference Marquillie and Ehrenstein2002). Equilibrium states, steady solutions of the system (2.1), are computed using the selective frequency damping method (Akervik et al. Reference Akervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006)

where $\boldsymbol{q}=(\boldsymbol{u},p,b)^{\text{T}}$ . For the stabilization coefficients, the optimal cutoff frequency $\unicode[STIX]{x1D6EC}=12.7423$ and damping factor $\unicode[STIX]{x1D701}=0.0534$ computed at $Ra=1.28\times 10^{9}$ for the no-slip case were used to converge to a steady state (Cunha, Passaggia & Lazareff Reference Cunha, Passaggia and Lazareff2015). Simulations of system (2.1) were considered to be converged to a steady state when the $l_{2}$ -norm of the difference of the velocity and buoyancy fields renormalized with the time step, was smaller than $10^{-8}$ . Here the time step $\unicode[STIX]{x0394}t=2\times 10^{-3}$ was used for all the simulations and the forthcoming stability analysis.

2.3 Global instability modes

The flow solution $\boldsymbol{q}$ is decomposed into a steady state $\boldsymbol{Q}=(\boldsymbol{U},P,B)$ , and a perturbation with the special form $\boldsymbol{q}^{\prime }=(\boldsymbol{u}^{\prime }(x,z,t),p^{\prime }(x,z,t),b^{\prime }(x,z,t))\exp (\text{i}\unicode[STIX]{x1D6FD}y)$ . The linearized perturbation dynamics is governed by

with $\tilde{\unicode[STIX]{x1D735}}=(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x,\text{i}\unicode[STIX]{x1D6FD},\unicode[STIX]{x2202}/\unicode[STIX]{x2202}z)^{\text{T}}$ and $\tilde{\unicode[STIX]{x1D6FB}}^{2}=(\unicode[STIX]{x2202}^{2}/\unicode[STIX]{x2202}x^{2}-\unicode[STIX]{x1D6FD}^{2}+\unicode[STIX]{x2202}^{2}/\unicode[STIX]{x2202}z^{2})$ . The free-slip and no-slip boundary conditions for the velocity of the perturbation were considered. A Dirichlet condition $b^{\prime }=0$ was used at the upper boundary whereas the buoyancy perturbation is insulated in the rest of the domain. Note that the numerical methods and discretization are the same as described in § 2.2. The linear system (2.3) may be equivalently written

where $\unicode[STIX]{x1D63C}$ depends on $(\boldsymbol{U},\boldsymbol{B},Ra,Pr,\unicode[STIX]{x1D6FD})$ and $\boldsymbol{E}$ is the projection operator onto the velocity and buoyancy field such that $\boldsymbol{E}\boldsymbol{q}=(\boldsymbol{u},0,b)^{\text{T}}$ . The stability is computed by considering the three-dimensional temporal modes $\hat{\boldsymbol{q}}(x,z)\exp (-\text{i}(\unicode[STIX]{x1D714}t-\unicode[STIX]{x1D6FD}y))$ , where $\hat{\boldsymbol{q}}$ and $\unicode[STIX]{x1D714}$ satisfy the generalized eigenvalue problem

Here, the matrix $\unicode[STIX]{x1D63C}$ is not computed explicitly and it is advantageous for large eigenvalue problems to only consider the operator that takes a state $\boldsymbol{q}^{\prime }(t)$ as initial condition for (2.4) and returns to the state $\boldsymbol{q}^{\prime }(t+\unicode[STIX]{x0394}t)$ after a given time $\unicode[STIX]{x0394}t$ (Edwards et al. Reference Edwards, Tuckerman, Friesner and Sorensen1994). Since Krylov-based eigenvalue solvers only require the storage of several flow states $\boldsymbol{q}^{\prime }(t)$ , the standard time marching algorithm is used to propagate the solution. However, such an approach will only identify a subset of eigenmodes with the largest growth rates $\unicode[STIX]{x1D714}_{i}$ . To circumvent this limitation, the generalized eigenvalue problem (2.5) was solved using the shift-and-relax method of Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Chomaz2012). The linear perturbation dynamics (2.4) was coupled to a lowpass filter $\bar{\boldsymbol{q}}$ such that

This method allows for computing eigenmodes with different frequency ranges, such that real frequencies $\unicode[STIX]{x1D714}_{r}$ , far from the target frequency $\unicode[STIX]{x1D714}_{0,r}$ are attenuated (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Chomaz2012). The parameters for the shift-and-relax method were $\unicode[STIX]{x1D712}=\unicode[STIX]{x1D70F}=0.1$ and all the computed eigenvalues were checked to converge with a relative precision of $10^{-3}$ , which is the bound given by the accuracy of the temporal scheme. Note that the system (2.6) is a generalization of (2.2) where the linear perturbation dynamics is forced at a particular frequency $\unicode[STIX]{x1D714}_{0}$ and damped everywhere else in the spectrum. The eigenvalues $\unicode[STIX]{x1D6FA}$ , associated with (2.6) can be transformed back to their original counterparts $\unicode[STIX]{x1D714}$ in (2.5) using the transformation (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Chomaz2012)

The eigenvalue problem is solved using a Krylov subspace of $80$ snapshots. Another simulation is also performed with a shift in the complex plane at $\unicode[STIX]{x1D714}_{0}=1$ to ensure the convergence of all the eigenvalues. All eigenvalues were computed with a relative error below $10^{-3}$ for each eigenvalue.

Figure 2.

Comparison between the local vertical profiles for the three-dimensional DNS code (——) (Scotti Reference Scotti2008) computed in a narrow domain (i.e. $W=0.125$ ) and the present two-dimensional DNS ( $\Box$ ) at $x=1$ and $Ra=2\times 10^{6}$ and $A=4$ for the horizontal velocity $U(z)$ (a) and the buoyancy $B(z)$ (b), in the free-slip case.

2.4 Validation of the numerical procedure

Results have been compared with a separate, well-validated, three-dimensional DNS code Scotti (Reference Scotti2008), and profiles of buoyancy and velocity were found to differ by less than 1 % between the two codes (e.g. see figure 2 a,b). The fact that the local and global stability analyses provide comparable results offers a validation of the stability analysis procedure. Also, verification of the bi-orthogonality condition (see (6.1) introduced in § 6.2) proves that the stability codes were implemented correctly. Furthermore, the local stability results offer a cross-validation for the stability methods used in the present study. We also simulated the flow using a three-dimensional DNS code for a transverse box width $W=1$ at $Ra=2\times 10^{6}$ in the free-slip case, and recovered the three-dimensional instability predicted by the stability analyses. The same code was also used and validated in previous studies (Marquillie & Ehrenstein Reference Marquillie and Ehrenstein2002, Reference Marquillie and Ehrenstein2003; Passaggia, Leweke & Ehrenstein Reference Passaggia, Leweke and Ehrenstein2012; Passaggia et al. Reference Passaggia, Meunier and Le Dizès2014; Cunha et al. Reference Cunha, Passaggia and Lazareff2015). Further detail on validation of the codes used here may be found therein.

Figure 3.

Streamlines (a–c) and buoyancy distribution (d–f) of the steady base flows at $Ra=1.92\times 10^{7}$ (a,d), $Ra=1.44\times 10^{8}$ (b,e) and $Ra=6.4\times 10^{8}$ (c,f) at $Pr=1$ , $A=4$ using free-slip boundary conditions.

4.1 Three-dimensional global modes

The DNS results of Gayen et al. (Reference Gayen, Griffiths and Hughes2014) suggest that the most amplified global modes are three-dimensional and are characterized by steady counter-rotating vortices with axes aligned along the $x$ direction. These modes have been observed for a large range of Rayleigh numbers and mark the onset of instability in horizontal convection. To study these modes, a global stability analysis was performed over a range of Rayleigh numbers $Ra=[2.02,2.56]\times 10^{7}$ and transverse wavenumbers $\unicode[STIX]{x1D6FD}=[4,8]$ where a family of unstable modes was found. These modes are characterized by zero frequency, meaning that they are steady with growth rates in the range $\unicode[STIX]{x1D714}_{i}=[-8.6\times 10^{-3},1.194\times 10^{-2}]$ . These growth rates are shown in figure 4(a) as a function of the transverse wavenumber $\unicode[STIX]{x1D6FD}$ . Increasing the Rayleigh number appears to increase both the growth rate and the most unstable wavenumber, as well as the unstable wavenumber range. The critical Rayleigh number for this mode was found at $Ra_{c,3D}\approx 1.24\times 10^{8}$ and gives the onset of instability for an aspect ratio $A=4$ (cf. figure 4 b).

Figure 4.

Growth rate $\unicode[STIX]{x1D714}_{i}$ (a) of the most unstable three-dimensional global mode as a function of the transverse wavenumber $\unicode[STIX]{x1D6FD}$ at $Ra=2\times 10^{7}$ (-▪-), $Ra=2.08\times 10^{7}$ (-●-), $Ra=2.24\times 10^{7}$ (-▴-), $Ra=2.56\times 10^{7}$ (-▾-). Most unstable growth rate $\unicode[STIX]{x1D714}_{i}$ (b) as a function of the Rayleigh number $Ra$ , depicting also the critical Rayleigh number $Ra_{c,3D}=2.02\times 10^{7}$ at $A=4$ , for free-slip boundary conditions.

Figure 5.

Visualization of the most unstable three-dimensional eigenmode by means of iso-buoyancy surfaces $\hat{b}=\pm 44\,\%$ normalized by the maximum of $\hat{b}$ (red/blue) at $Ra=1.44\times 10^{8}$ , $Pr=1$ , $\unicode[STIX]{x1D6FD}=6$ and $\unicode[STIX]{x1D714}_{i}=0.0026259$ , for free-slip boundary conditions.

Figure 6.

Velocity distribution $\hat{u} (x,z)$ (a), $\hat{v}(x,z)$ (b), ${\hat{w}}(x,z)$ (c), $\hat{b}(x,z)$ (d), of the most unstable three-dimensional mode ( $\unicode[STIX]{x1D6FD}=6$ ) at $Ra=1.44\times 10^{8}$ and $Pr=1$ , using free-slip boundary conditions (the amplitude was normalized using the maximum of transverse velocity $\hat{v}$ ).

The structure of the unstable modes, shown in figure 5, is characterized by convection rolls, originating in the vicinity of the statically unstable fluid in the upper left boundary. In the near-surface region, the global instability mode, shown in figure 5, is essentially characterized by an increase of buoyancy and longitudinal vorticity. Within the descending plume, all quantities appear to switch sign (cf. figure 6). The same observation can be made near the bottom of the domain where the counter-rotating vortices span most of the lower boundary. Finally, its structure vanishes in the upwelling region.

The growth rates and the most unstable wavenumber for larger values of the Rayleigh number are shown in figure 7(a,b). The growth rate follows the inverse scaling of the e-folding time scale given in (4.1), suggesting that the instability is of the Rayleigh–Taylor type. This scaling however changes at $Ra\approx 6.4\times 10^{8}$ and this interesting feature will be discussed in a forthcoming manuscript. The scaling for the most unstable wavelength follows a $Ra^{1/3}$ for all values of $Ra$ , confirming the hypothesis of a Rayleigh–Taylor instability.

Figure 7.

Most unstable growth rate $\unicode[STIX]{x1D714}_{i}$ (a) as a function of $Ra$ (○) compared with the $Ra^{2/3}$ scaling (——). Most unstable wavenumber $\unicode[STIX]{x1D6FD}$ (b) as a function of $Ra$ (○) compared with the $Ra^{1/3}$ scaling (——), using free-slip boundary conditions.

6.1 Global stability analysis

The growth rate of the pure two-dimensional mode ( $\unicode[STIX]{x1D6FD}=0$ ) is shown in figure 12(a) as a function of the Rayleigh number. The critical value for the onset of two-dimensional instability is found to be $Ra_{c,2D}\approx 1.65\times 10^{8}$ . The eigenfunction for the unstable two-dimensional mode ( $\unicode[STIX]{x1D6FD}=0$ ) is depicted in figure 13. The structure of the mode mostly spans the left boundary of the flow and takes the form of regular transverse vortices shed inside the domain, which effectively transport buoyancy as well as kinetic and potential energy. While the eigenfunction shown in figure 13 is for $\unicode[STIX]{x1D6FD}=0$ , pure Rayleigh–Taylor modes are found in the global spectrum for $\unicode[STIX]{x1D6FD}=9$ corresponding to transverse wavelengths smaller than the depth of the domain. Three-dimensional modes (cf. figure 14 b,c) are also found with a frequency and have a maximum growth rate $\unicode[STIX]{x1D714}_{i}$ , found for $\unicode[STIX]{x1D6FD}\sim 12$ . This is depicted in figure 15 where the global spectrum at $Ra=6.4\times 10^{8}$ is shown for several values of the transverse wavenumber. These modes also exhibit a lower frequency $\unicode[STIX]{x1D714}_{r}$ than the two-dimensional instability, and the latter seems to be independent of the Rayleigh number (cf. figure 12(b) showing the most unstable quasi two-dimensional mode at $Ra=6.4\times 10^{8}$ ). Tsai et al. (Reference Tsai, Hussam, Fouras and Sheard2016) found similar results in their local stability analyses. They showed that both two and three-dimensional instabilities had comparable temporal growth rates using a local stability analysis of the parallel base flow near the plume.

They also computed the marginal stability and showed that the longitudinal and transverse rolls were marginally unstable at $Ra\approx 1.23\times 10^{7}$ but for $Pr=6.14$ and $A=6.25$ . The most unstable wavelength are also in qualitative agreements with the present study. Tsai et al. (Reference Tsai, Hussam, Fouras and Sheard2016) found that the onset of three-dimensional local instability was found near the plume, for $\unicode[STIX]{x1D6FD}\approx 19$ (in our units) compared to $\unicode[STIX]{x1D6FD}\approx 12$ in our global stability analysis. Keeping in mind that the local temporal stability analysis defines regions where the flow is convectively unstable, we perform a spatio-temporal stability analysis with the support of a sensitivity analysis to link the local stability properties of the flows to the global stability analysis.

Gayen et al. (Reference Gayen, Griffiths and Hughes2014) and Tsai et al. (Reference Tsai, Hussam, Fouras and Sheard2016) reported that for high Rayleigh numbers, convective instabilities dominate the unstable dynamics as well as the transfer from available potential energy to kinetic energy, responsible for mixing (Gayen et al. Reference Gayen, Griffiths, Hughes and Saenz2013). In the next part, the convective/absolute nature of the flow is examined to determine the mechanisms responsible for this instability.

6.2 Linking the global to the local instability properties

In this non-parallel flow set-up, it is important to define the regions where the flow is most sensitive to a modification of the flow. The structural sensitivity as developed by Hill (Reference Hill1992) and Giannetti & Luchini (Reference Giannetti and Luchini2007) is given by the dyadic product of the direct and adjoint momentum vectors $\unicode[STIX]{x1D61A}_{i,j}$ (i.e. $\unicode[STIX]{x1D61A}_{i,j}=\hat{\boldsymbol{u}_{i}}\hat{\boldsymbol{u}}_{j}^{+,H}$ ) where the superscript $H$ denotes the complex conjugate transpose, the superscript $\dagger$ defines the adjoint eigenmode and the subscripts $i$ and $j$ are the indices $[1,2,3]$ associated with each component of the velocity vector $\boldsymbol{u}$ . The components of this complex valued $3\times 3$ tensor shows how feedback between the components of the linearized momentum vector affects the growth rate and frequency of the global mode (Qadri, Mistry & Juniper Reference Qadri, Mistry and Juniper2013). Here the adjoint global modes $\hat{\boldsymbol{q}}^{\dagger }(x,0,z)$ can be defined using the bi-orthogonal property

where $\unicode[STIX]{x1D6FF}$ is the Kronecker symbol and $\hat{\boldsymbol{q}}_{j}^{\dagger }$ is the adjoint eigenvector associated with the $j$ th direct eigenvector. Using equation (2.5), the inner product $\langle \boldsymbol{a},\boldsymbol{c}\rangle =\int _{\unicode[STIX]{x1D6E4}}\hspace{2.22198pt}\boldsymbol{a}^{H}\boldsymbol{\cdot }\boldsymbol{c}\,\text{d}\unicode[STIX]{x1D6E4}$ , where the superscript $H$ denotes the complex conjugate transpose and $\unicode[STIX]{x1D6E4}$ is the computational domain. We may then express the adjoint identity

where the adjoint operator $\unicode[STIX]{x1D63C}^{H}$ is defined by

and where the boundary terms

are cancelled by considering the same boundary conditions between the linear perturbation dynamics and the adjoint solver (Kaminski, Caulfield & Taylor Reference Kaminski, Caulfield and Taylor2014). Finally, the adjoint global modes are computed using the shift-and-relax method and the difference between the direct and the adjoint eigenvalues were obtained with a relative error of $5\times 10^{-3}$ .

The expression for $\unicode[STIX]{x1D61A}_{i,j}$ can be simplified to a scalar quantity (Giannetti & Luchini Reference Giannetti and Luchini2007), namely the Frobenius norm of the sensitivity tensor $\unicode[STIX]{x1D61A}_{i,j}$ which identifies where a modification in the linearized equations produces the largest drift of the eigenvalue and thereby, reveals the region of the flow that acts as a wavemaker. This expression writes

and is maximum in regions that drives the dynamics associated with the global mode (Juniper & Pier Reference Juniper and Pier2015). Here $\Vert \hat{\boldsymbol{u}}^{\dagger }(x,0,z)\Vert$ and $\Vert \hat{\boldsymbol{u}}(x,0,z)\Vert$ represent the $l^{2}$ -norm of the adjoint and direct global momentum vectors at each point in space (Chandler et al. Reference Chandler, Juniper, Nichols and Schmid2012).

6.3 Local stability analysis

The question naturally arises whether the unstable global mode is connected to local linear stability properties, that is whether specific regions of the flow are responsible for triggering the onset of global instability. We consider the unstable state at $Ra=6.4\times 10^{8}$ , which is a globally unstable two-dimensional steady state computed in the § 2.2. In light of the global mode analysis, the structural sensitivity $\unicode[STIX]{x1D709}$ is used to locate regions likely to be absolutely unstable. The wavemaker $\unicode[STIX]{x1D709}(x,0,z)$ is shown in figure 16 where region near the unstable boundary seems to play an important role. These contours have their greatest magnitude bellow the boundary $b=0$ , near to the top-left corner, where the descending plume originates. These results are in agreements with the local stability analysis of Tsai et al. (Reference Tsai, Hussam, Fouras and Sheard2016) where this region was found to be the most convectively unstable.

Figure 16.

Structural sensitivity vector $\unicode[STIX]{x1D709}(x,0,z)$ associated with the most unstable two-dimensional mode ( $\unicode[STIX]{x1D6FD}=0$ ), at $Ra=6.4\times 10^{7}$ and $Pr=1$ for no-slip boundary conditions.

In the following, we define the streamfunction for the perturbation such that

Different velocity profiles $U(z)$ and buoyancy profiles $B(z)$ for the streamwise positions between $x_{p}=-1.375$ and $x_{p}=-1.525$ have been considered where the wavemaker dominates. A local linear stability analysis is performed for parallel flow profiles $(U(x_{p},z),0,B(x_{p},z))$ using normal modes of the form

Linearizing the Navier–Stokes system at $(U(z),0,B(z))$ and taking the curl of (2.1) in the $(x{-}z)$ plane, the Orr–Sommerfeld equation is coupled to the buoyancy such that

Once the modes are discretized using Chebyshev collocation, a matrix eigenvalue problem is recovered and solved numerically for the dispersion relation $D(\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FC})=0$ . A spatio-temporal linear stability analysis has been performed, solving $\unicode[STIX]{x1D714}(\unicode[STIX]{x1D6FC})$ for both complex wavenumbers $\unicode[STIX]{x1D6FC}$ and frequencies $\unicode[STIX]{x1D714}$ . Solving $\unicode[STIX]{x1D714}(\unicode[STIX]{x1D6FC})$ , the complex local absolute frequency is given in classical fashion by

The locally parallel flow defined by the local profile $(U(x_{p},z),0,B(x_{p},z))$ is absolutely unstable if the imaginary part of the complex absolute frequency is positive (Schmid & Henningson Reference Schmid and Henningson2001; Chomaz Reference Chomaz2005). Note that in the present vertical profiles, the fluid near the top boundary is unstable (cf. figure 17 a). This was confirmed by the local stability analysis where the unstable eigenfunctions were only found near the upper boundary. Figure 17(b) depicts the analysis for the local profiles from $x=-1.375$ (dotted line) until $x=-1.525$ (continuous line). The $\unicode[STIX]{x1D6FC}_{i}=\text{Const.}$ contours (for $0.45<\unicode[STIX]{x1D6FC}_{r}<0.7$ ) in the complex frequency plane are shown and the cusp (the point $(\unicode[STIX]{x2202}\unicode[STIX]{x1D714})/(\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC})(\unicode[STIX]{x1D6FC}_{0})=0$ ) appears for $\unicode[STIX]{x1D6FC}_{r}\approx 5.98$ , $\unicode[STIX]{x1D6FC}_{i}\approx 5.25$ , under and above the real axis for $x=-1.375$ and $x=-1.425$ respectively. The saddle point described in (6.9) is shown in figure 17(c) where the pinching point between two distinct branches $\unicode[STIX]{x1D6FC}+$ and $\unicode[STIX]{x1D6FC}-$ ( $\unicode[STIX]{x1D6FC}$ complex and $\unicode[STIX]{x1D714}_{i}$ constant) originate from distinct halves of the $\unicode[STIX]{x1D6FC}-$ plane (Schmid & Henningson Reference Schmid and Henningson2001). These results demonstrate that the flow becomes absolutely unstable near the upper boundary, in the vicinity of the plume. Prior to the discussion in § 4, we investigated the effect of both shear and viscosity on the local stability properties of the flow. In the present local stability analysis at $x_{p}=-1.525$ , decreasing viscosity only has a weak effect on the local stability, but increasing the magnitude of the shear stabilizes the flow.

Figure 17.

Longitudinal velocity $U(x_{p},z)$ (- - -) and Brunt–Väisälä frequency profiles $\unicode[STIX]{x2202}B(x_{p},z)/\unicode[STIX]{x2202}z$ (——) (a) located near the plume at $x_{p}=-1.525$ , using no-slip boundary conditions. Contours in the plane ( $\unicode[STIX]{x1D714}_{r},\unicode[STIX]{x1D714}_{i}$ ) at $x_{p}=-1.525$ (——) and $x_{p}=-1.475$ (- - -) at $x_{p}=-1.425$ (- - -) and $x_{p}=-1.375$ ( $\ldots$ ) at $\unicode[STIX]{x1D6FC}_{i}=5.25$ , (cusp map $(\unicode[STIX]{x2202}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D6FC}_{0}))/(\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC})=0$ ) (b). Spatial branches in the complex $\unicode[STIX]{x1D6FC}$ -plane for $Ra=6.4\times 10^{7}$ showing the saddle point for iso-contours of $\unicode[STIX]{x1D714}_{r}$ (c) and where the white lines show the saddle point associated with the growth rate $\unicode[STIX]{x1D714}_{i}$ .

Comparing these local results with the global stability analysis, both the growth rate and the longitudinal wavenumber of the instability compare favourably ( $\unicode[STIX]{x1D714}_{0,i}=0.0291$ versus $\unicode[STIX]{x1D714}_{global,i}=0.0307$ (cf. figure 12 a) and $\unicode[STIX]{x1D6FC}_{0}\approx 5.42$ versus $\unicode[STIX]{x1D6FC}_{global}\approx 5.2$ ). Despite this good agreement in the wavenumber and growth rate, the global frequency $\unicode[STIX]{x1D714}_{0,r}=0.498$ is smaller than the global frequency $\unicode[STIX]{x1D714}_{global,r}=0.574$ The spatial growth rate $\unicode[STIX]{x1D6FC}_{i}\approx 5.25$ is large, meaning that the mode decays very fast towards the centre of the domain, in agreement with the modal structure of the global instability, but the global mode is influenced by the non-parallel nature of the flow as $x\rightarrow -2$ . Thus it seems unlikely that non-parallel effects can be neglected. A fundamental assumption of the local stability analysis is the existence of two separate length scales in the flow: the evolution length scale $\unicode[STIX]{x1D706}$ , characterizing the spatial development of the base flow, and the instability wavelength, $2\unicode[STIX]{x03C0}/(\unicode[STIX]{x1D6FC}_{0})\approx 1$ . The length scale $\unicode[STIX]{x1D706}$ can be defined by means of the momentum or vorticity thickness of the base flow profiles $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D6E9}(x)$ and writes

The small parameter of interest is then

In the present investigations, both the velocity $U(z)$ and the buoyancy profiles $B(z)$ of the base flow evolve rapidly near the upper boundary. In this region, the parameter $\unicode[STIX]{x1D716}\approx 2$ shows that a weakly non-parallel asymptotic analysis is impossible. The same issue was reported for instance in the case of a separated flow behind a bump where the absolute frequency $\unicode[STIX]{x1D714}_{0,r}$ (Marquillie & Ehrenstein Reference Marquillie and Ehrenstein2003) did not match the global instability $\unicode[STIX]{x1D714}_{global,r}$ (Ehrenstein & Gallaire Reference Ehrenstein and Gallaire2008) due to the rapid evolution of the base flow near the reattachment region. We conclude that the two-dimensional instability is associated with a local absolutely unstable sheared Rayleigh–Taylor instability controlled by non-parallel effects near the plume region.

Figure 18.

Growth rate $\unicode[STIX]{x1D714}_{i}$ extracted from the global stability as a function of the domain’s aspect ratio $A=L/H$ , for the most unstable three-dimensional mode (-▪-) and the two-dimensional mode (-●-) at $Ra=A^{3}2.26\times 10^{6}$ for free-slip boundary conditions (a) and $Ra=A^{3}3\times 10^{6}$ for no-slip boundary conditions (b) at $Pr=1$ .