The stability of stratified spatially periodic shear flows at low Péclet number

# The stability of stratified spatially periodic shear flows at low Péclet number

## Abstract

This work addresses the question of the stability of stratified, spatially periodic shear flows at low Péclet number but high Reynolds number. This little-studied limit is motivated by astrophysical systems, where the Prandtl number is often very small. Furthermore, it can be studied using a reduced set of “low-Péclet-number equations” proposed by Lignieres [Astronomy & Astrophysics, 348, 933-939, 1999]. Through a linear stability analysis, we first determine the conditions for instability to infinitesimal perturbations. We formally extend Squire’s theorem to the low-Péclet-number equations, which shows that the first unstable mode is always two-dimensional. We then perform an energy stability analysis of the low-Péclet-number equations and prove that for a given value of the Reynolds number, above a critical strength of the stratification, any smooth periodic shear flow is stable to perturbations of arbitrary amplitude. In that parameter regime, the flow can only be laminar and turbulent mixing does not take place. Finding that the conditions for linear and energy stability are different, we thus identify a region in parameter space where finite-amplitude instabilities could exist. Using direct numerical simulations, we indeed find that the system is subject to such finite-amplitude instabilities. We determine numerically how far into the linearly stable region of parameter space turbulence can be sustained.

## I Introduction

The study of the onset of turbulence in stratified shear flows has a long history that dates back to Richardson (1920). He argued that the kinetic energy of turbulent eddies in a stratified shear flow can only decrease if , where is the local buoyancy frequency, and is the local shearing rate of the flow field in the vertical direction . This criterion, derived simply from energetic arguments, is now commonly referred to as Richardson’s criterion, and the local ratio

 J(z)=N2(z)S2(z) (1)

is called the gradient Richardson number. The first linear stability analysis of a stratified shear flow is due to Taylor (1931), who considered both continuously and discretely varying stratification and shear profiles. This work, together with Goldstein (1931), then led to the derivation of the Taylor-Goldstein eigenvalue equation for the complex growth rate of two-dimensional infinitesimal disturbances in stratified shear flows. The solution of this equation for a given shear profile and stratification profile can be obtained either analytically in a few particular cases, or numerically in general. It wasn’t until much later, however, that the first general result on the stability of stratified shear flows was derived by Miles (1961) and Howard (1961): a system is stable to infinitesimal perturbations provided is everywhere larger than 1/4. As discussed by Howard and Maslowe (1973), this theorem should not be viewed as a refinement of Richardson’s argument (i.e. replacing 1 by 1/4), since the latter was specifically interested in determining when turbulence could be sustained, rather than triggered. In this sense, Richardson’s original argument should be viewed more as a nonlinear stability criterion than a linear one.

These results were obtained in the limit of vanishing viscosity and diffusivity. For thermally stratified flows, however, thermal diffusion can have a significant influence on the development of shear instabilities by damping the buoyancy restoring force. This effect was first studied by Townsend (1958) in the context of atmospheric flows. He showed that the thermal adjustment of the fluid parcel to its surroundings, by radiative heating and cooling or by thermal conduction, always acts to destabilize the flow and increases the critical Richardson number for linear stability, , by a factor inversely proportional to the product of the shearing rate with the cooling time (this product is a local Péclet number for the flow), so . Viscosity, meanwhile, has a generally stabilizing influence (Dudis, 1974). Zahn (1974) emphasized the importance of these results for stellar astrophysics: in stellar interiors where the Prandtl number is typically very small () high Reynolds number flows can also have a low Péclet number, or in other words, thermally diffusive shear flows exist when viscosity is nevertheless small enough not to suppress the development of the instability. This combination is ideal for shear instabilities, and is specific to astrophysical systems – it cannot happen for most geophysical flows where the Prandtl number is usually of order unity or larger.

Applying Townsend (1958)’s results to shear-induced turbulence in stellar interiors, Zahn (1974) further argued that the relevant cooling timescale is the radiative timescale based on the size of turbulent eddies, namely where is the thermal diffusivity. He then proposed to take for the smallest length scale for which viscosity is still negligible, that is, one for which the turbulent Reynolds number (where is the kinematic viscosity), where is a constant that he estimates to be around . This would imply , or in other words, . Zahn (1974)’s argument, as in the case of Richardson (1920)’s original argument, should be viewed as a nonlinear stability criterion rather than a linear one, since it relies on the presence of pre-existing turbulent eddies.

Zahn’s work had an enormous impact in the field of stellar evolution. While the standard Richardson criterion is far too stringent to allow for the development of shear instabilities in the absence of thermal diffusion (the typical Richardson number being much larger than one even in the strongest known stellar shear layers), its relaxation allows for the possibility of much-needed mixing in stellar evolution theory. Indeed, models without any form of turbulent mixing in stably stratified regions are not able to account for observations. As reviewed by Pinsonneault (1997) the problem is particularly acute when it comes to explaining the surface chemical abundances of light elements such as lithium and beryllium, as well as products and by-products of nuclear reactions such as helium, carbon, nitrogen and oxygen. More recently, further indication of the need for turbulent mixing was revealed by asteroseisomolgy thanks to the Kepler mission. Measurements of the internal rotation rate of red giant stars (Mosser et al., 2012) are inconsistent with evolution models in which turbulent angular-momentum transport in stably stratified regions is neglected. In both cases, therefore, efficient chemical transport and angular-momentum transport by shear instabilities could be the key to resolving these problems – the question remains, however, of whether these instabilities are indeed triggered, and how efficient mixing is.

In the limit of low-Péclet numbers (i.e., high thermal diffusivity), the temperature fluctuations are slaved to the vertical velocity. The corresponding quasi-static approximation was originally introduced to study low-Prandtl-number thermal convection(Spiegel, 1962; Thual, 1992). In the context of stably-stratified systems, the quasi-static approximation was introduced only recently by Lignières (1999) (see Section II.3 for more detail). He showed that the standard Boussinesq equations can be replaced by a reduced model that is valid in the asymptotic low-Péclet-number limit, and that this model only depends on two parameters: the Reynolds number, and the product of the Richardson number with the Péclet number. As a result, the linear stability properties of the system depend on the product (where with being a characteristic vertical length scale of the laminar flow) rather than on each parameter individually. Since is small by assumption, shear-induced turbulence can be expected even if the Richardson number is much larger than one.

By contrast with linear theory, very little is known to date about the stability of stratified shear flows to finite-amplitude perturbations when viscosity and thermal diffusivity are both taken into account. It is yet a question of crucial importance in stellar astrophysics, since the presence or absence of vertical mixing can strongly affect model predictions. We therefore set out in this work to characterize the domain of instability of low-Péclet-number shear flows to finite-amplitude perturbations. For simplicity, we consider a specific shear profile that is periodic in the vertical direction. We present the model in Section II, and briefly discuss the low-Péclet-number asymptotic equations proposed by Lignières (1999). In Sections III and IV, we study the linear and nonlinear stability of the system respectively, and contrast the results in the low-Péclet-number approximation to those obtained starting from the full set of primitive equations. As we shall demonstrate using an energy stability analysis of the low-Péclet number equations, smooth periodic shear flows are stable to perturbations of arbitrary amplitude for sufficiently large Richardson number, for a given value of the Reynolds number. In Section V, we turn to direct numerical simulations to study the transition to turbulence via linear instabilities and finite-amplitude instabilities. We summarize our results and conclude in Section VI.

## Ii The model

### ii.1 Model setup

Since our intention is to study the energy stability properties of stratified shear flows, it is crucial to start with a model where the mechanism driving the shear is explicit, which guarantees a well-defined energy budget. Two options are available: boundary-forcing, and body-forcing. Having potential applications to stellar astrophysics in mind, we prefer the latter in order to avoid boundary layer dynamics near solid walls, which are rarely present in stars.

A simple and numerically-efficient way of studying body-forced, stratified shear flows is to consider a Boussinesq system (Spiegel and Veronis, 1960), where the forcing and all the perturbations are triply-periodic, and where the background density is linearly stratified1 (see Figure 1). The model equations describing such a system are:

 ∂u∂t+u⋅∇u=−1ρ0∇p+αgTez+ν∇2u+1ρ0F , (2) ∇⋅u=0 , (3) ∂T∂t+u⋅∇T+wT0z=κT∇2T , (4)

where is the triply-periodic velocity field, and are the triply-periodic pressure and temperature perturbations, is the mean density of the region considered, is the coefficient of thermal expansion, is gravity, and are the viscosity and thermal diffusivity (respectively). The quantities , , , , are all assumed to be constant, as in the standard Boussinesq approximation. The use of the latter is justified as long as the vertical height of the domain is much smaller than a density scaleheight. Finally, we assume that there is a constant background temperature gradient2 , and that all thermodynamical and dynamical perturbations have zero mean in the domain.

The applied force should be triply-periodic as well. A natural candidate is a sinusoidal forcing, thus driving a Kolmogorov flow in the laminar regime. In what follows, we assume that is of the form

 F=F0sin(kz)ex , (5)

which defines a typical lengthscale . In the steady laminar regime, this force generates a sinusoidal shear flow along the -direction,

 uL=F0ρ0νk2sin(kz)ex. (6)

Note that while the present paper deals mostly with Kolmogorov forcing, the energy stability of arbitrary smooth velocity profiles is discussed in section IV.2.2.

### ii.2 Non-dimensionalization and model parameters

We non-dimensionalize the equations using the amplitude of the laminar solution, , as a velocity scale. We also use the spatial scale of the laminar solution, , as the unit length scale. This then defines the timescale . With this choice of units, the equations (2)-(4) become

 ∂u∂t+u⋅∇u=−∇p+RiTez+1Re∇2(u−sin(z)ex) , (7) ∇⋅u=0 , (8) ∂T∂t+u⋅∇T+w=1Pe∇2T , (9)

where

 Re=F0ρ0ν2k3,Ri=αgT0zρ20ν2k2F20,Pe=F0ρ0νk3κT. (10)

The laminar solution (6) now takes the dimensionless form . Provided the system remains in the vicinity of this laminar solution, , and are the usual Reynolds, Péclet and Richardson numbers based on the typical velocity of the flow.

### ii.3 Low Péclet number approximation

When a field diffuses on a timescale much shorter than the advective time, it enters a quasi-static regime where the source term and diffusive term instantaneously balance. Such a quasi-static regime has been used for decades in the context of magneto-hydrodynamics of liquid metals: at low magnetic Reynolds number, the induced magnetic field is slaved to the velocity field (see for instance Sommeria and Moreau, 1982). The equivalent approximation for flows of low Prandtl number fluids was originally introduced in the context of thermal convection (Spiegel, 1962; Thual, 1992). Rather surprisingly, this quasi-static approximation appeared only much more recently in the context of stably-stratified flows.

Lignières (1999) proposed that in the low-Péclet-number limit the governing equations (7)-(9) can be approximated by the reduced set of so-called “low-Péclet-number” equations (LPN equations hereafter)

 ∂u∂t+u⋅∇u=−∇p+RiTez+1Re∇2u+F , (11) ∇⋅u=0 , (12) w=1Pe∇2T , (13)

to zeroth and first order in . To derive equation (13), one can assume a regular expansion of in , as in , and further assume that the velocity field is of order unity. At the lowest order, the temperature equation yields which then implies given the applied boundary conditions. The next order then yields the quasi-static balance as required. It is worth mentioning that in Lignières (1999)’s original work the velocity is scaled with its dimensional r.m.s. value , so it is , rather than , that has to be small for the LPN equations to be valid. In Sections III and IV, we shall study the linear and energy stability of a laminar flow for which the r.m.s. velocity is of the same order as the flow amplitude. In that case and we expect the LPN equations to be valid whenever is small. In Section V, however, we shall see that numerical simulations of turbulent shear flows that have large can still be well-described by the LPN equations as long as is small.

It is also worth noting that the LPN equations only allow for temperature fluctuations that have a zero horizontal mean, by contrast with the full equations. As a result, the horizontal mean of the full temperature field (background plus perturbations) is necessarily linear in . To see this, we take the horizontal average of the thermal equation, which, assuming that there is no vertical mean flow (which can be guaranteed by making sure the initial conditions do not have one), results in

 ∂2⟨T⟩h∂z2=0 , (14)

in the LPN equations, where is the horizontal average of . The only solution of this equation which satisfies periodicity is the constant solution; further requiring that the volume-average of be zero then implies that . By contrast, taking the horizontal average of the standard temperature equation under the same assumptions would result in

 ∂⟨T⟩h∂t+∂∂z⟨wT⟩h=1Pe∂2⟨T⟩h∂z2 , (15)

which has solutions with non-zero . These solutions can, for instance, develop into density staircases under the right circumstances(Phillips, 1972; Balmforth, Smith, and Young, 1998). The latter are however prohibited in the LPN equations. Among other effects, this rules out the development of Holmboe modes(Holmboe, 1962), and may explain why it is possible to get simple energy stability results in the LPN limit but not for the full equations.

Finally, combining the momentum and the thermal energy equations yields

 ∂u∂t+u⋅∇u=−∇p+RiPe∇−2wez+1Re∇2u+F , (16)

which formally shows that the Richardson number is no longer the relevant parameter of the system, but that is. The work of Lignières (1999) thus puts the arguments of Townsend (1958) and Zahn (1974) discussed in Section I on a firm theoretical footing. Lignières (1999) and Lignières, Califano, and Mangeney (1999) verified that the LPN equations correctly account for the linear stability properties of various systems in the low-Péclet-number limit. Prat and Lignières (2013) later also verified that they correctly reproduced the low-Péclet dynamics of their 3D nonlinear simulations. In this paper, we continue to verify the validity of the LPN equations through stability analyses and nonlinear simulations.

## Iii Linear stability of a periodic Kolmogorov flow

We first focus on the stability of the laminar solution to infinitesimal perturbations. We solve the linearized versions of equations (7)-(9),

 ∂u′∂t+uL⋅∇u′+u′⋅∇uL=−∇p+RiT′ez+1Re∇2u′ , (17) ∇⋅u′=0 , (18) ∂T′∂t+uL⋅∇T′+w′=1Pe∇2T′ , (19)

where and are infinitesimal perturbations to the linearly stratified background shear flow .

### iii.1 Squire’s transformation

The linear stability of the unstratified Kolmogorov flow was first investigated in detail by Beaumont (1981). Squire’s theorem (Squire, 1933) states that the first unstable mode as the Reynolds number increases is a (-independent) 2D mode. This strong result implies that one can focus on 2D perturbations to determine the stability threshold of the system. Such a 2D analysis is much simpler and computationally less expensive than a 3D one. Beaumont (1981) found that 2D flows are unstable only for .

The linear stability of the stratified Kolmogorov flow to 2D perturbations was studied in detail by Balmforth and Young (2002). The 2D case can be made more generally relevant by noting that Squire’s transformation (Squire, 1933) for the viscous unstratified case can be extended to the stratified case with thermal diffusion to argue that the linear stability of any 3D mode can equivalently be studied by considering that of a 2D mode at lower or equal Reynolds and Péclet numbers, and higher or equal Richardson number. This result, which was summarily discussed by Yih (1955) and clarified by Smyth, Klaassen, and Peltier (1988) and Smyth and Peltier (1990) (see also Drazin and Reid, 2004), states that the growth rate of the 3D normal mode at parameters is related to that of the 2D normal mode at suitably rescaled parameters via

 λ3≡f(l,m;Re,Pe,Ri)=lLλ2≡lLf(L,0;lLRe,lLPe,L2l2Ri) , (20)

where . One can apply the same method to the LPN equations, and the result can readily be deduced from (20). Indeed, the transformation (20) is valid for any Péclet number, so it remains valid for low Péclet numbers. In this limit, we saw that only the product is a relevant parameter: as a consequence, one can replace the last two arguments of the function in (20) by the product of the two. The low-Péclet version of Squire’s transformation therefore gives

 λ3≡f(l,m;Re,RiPe)=lLλ2≡lLf(L,0;lLRe,LlRiPe) . (21)

This relationship between the growth rates of 2D and 3D modes has important implications for the marginal linear stability surface. In order to find the latter in 2D, we first maximize the real part of over all possible values of , yielding the function which returns the growth rate of the fastest growing mode for each parameter set . The marginal linear stability surface is then defined by . Similarly, the marginal linear stability surface for 3D perturbations is obtained by constructing the function and setting . If the functions and are the same, then so are the surfaces and , which implies in turn that the first modes to be destabilized are the 2D modes. This is the case for instance in the limit where stratification is negligible (see above).

In general, the only way to determine whether for a given linear stability problem is to construct these functions by brute force, using their original definition as the growth rates of the fastest growing modes. While this is not too time-consuming in 2D, it can become computationally expensive in 3D. However in this particular problem, since the growth rates of 2D and 3D modes are related, we also have:

 S3(Re,Pe,Ri)=maxχ∈[0,1]χS2(χRe,χPe,Ri/χ2) , (22)

where . A similar relationship applies for the LPN equations:

 S3(Re,RiPe)=maxχ∈[0,1]χS2(χRe,RiPe/χ) . (23)

Note that it is easier to construct from equations (22) or (23) than to do so directly.

Whether or not then simply depends on the properties of . It is quite easy to find sufficient conditions that guarantee . For instance, in the case of the standard equations, if is a strictly increasing function of both and , and a strictly decreasing function of , then the maximum over all possible values of is achieved for , which ensures that . For the LPN equations, it is sufficient to show that is a strictly increasing function of and a strictly decreasing function of . In what follows, we therefore first study the stability of 2D modes, and then use these results to conclude on the stability of the system to 3D modes.

### iii.2 Linear stability analysis using Floquet theory

We use a stream function to describe divergence-free 2D perturbations,

 u=uL+∇×(ψ′ey) , (24)

where is the infinitesimal perturbation. The linearized equations (17)-(19) become

 Missing or unrecognized delimiter for \left ∂T′∂t+sin(z)∂T′∂x+∂ψ′∂x=1Pe∇2T′ . (25)

This set of PDEs for and has coefficients that are independent of and , but periodic in . Normal modes for this system are of the form

 q′(x,z,t)=eiLx+λt^q(z) , (26)

where is either or , and is real. Using Floquet theory, we then seek solutions for given by

 ^q(z)=eiazN∑n=−Nqneinz , (27)

where is real, to satisfy the general periodicity of the system. Substituting this ansatz into the previous equations, we obtain an algebraic system for the and :

 −λ((a+n)2+L2)ψn+L2[(1−(a+n−1)2−L2)ψn−1−(1−(a+n+1)2−L2)ψn+1] =iRiLTn+1Re((a+n)2+L2)2ψn , λTn+L2[Tn−1−Tn+1]+iLψn=−1Pe((a+n)2+L2)Tn , (28)

for . This can be cast as the linear eigenproblem,

 M(L;a;Re,Pe,Ri)x=λx , (29)

where , which can be solved for the complex growth rate . The real part of the latter can then be maximized over all possible values of and for given system parameters to determine the temporal behavior and spatial structure of the most rapidly growing mode of the shear instability, or in other words, to construct (see previous section). In both unstratified and stratified cases studied so far, the first unstable modes at the instability threshold have the same periodicity in as that of the background shear, so that (Beaumont, 1981; Gotoh, Yamada, and Mizushima, 1983; Balmforth and Young, 2002). We verified that this is indeed the case here as well. In what follows, we therefore restrict the presentation of our results to the case .

The equivalent problem for the LPN equations is given by

 −λ((a+n)2+L2)ψn+L2[(1−(a+n−1)2−L2)ψn−1−(1−(a+n+1)2−L2)ψn+1] =L2(a+n)2+L2RiPeψn+1Re((a+n)2+L2)2ψn , (30)

which can be cast as

 M′(L;a;Re,RiPe)y=λy , (31)

where . This time, the fastest growing mode only depends on two system parameters, namely and the product . Again, we restrict the following analysis to the case .

Various aspects of the marginal stability surface for 2D modes are presented in Figure 2. Figure 2a shows the critical Reynolds number as a function of for the standard equations, for various values of the Prandtl number (). The evolution of the shape of these curves as increases is not a priori easy to identify nor explain. However, an obvious result is the existence of unstable modes for reasonably large values of the Richardson number when the Prandtl number is low. This can easily be understood in the light of the work of Townsend (1958) (see also Gage and Miller (1974), Jones (1977),Lignières (1999),Lignières, Califano, and Mangeney (1999) for instance), who showed that stratified shear instabilities can exist beyond the standard Richardson criterion when thermal diffusion is important. Since thermal diffusion increases as decreases, and since , one can naturally expect unstable modes at high Richardson number for fixed and low enough (or vice-versa).

Following Lignières, Califano, and Mangeney (1999), we now show in Figure 2b the same data plotted against , and add the marginal stability curve for the LPN equations. The interpretation of the results is now much clearer. We first see that the LPN equations are indeed a good approximation to the full equations when the Péclet number is small (low ). In the unstratified limit , we find that marginal stability is indeed achieved for , as expected (Beaumont, 1981). We also see that, for the low-Péclet equations, the threshold for linear stability above which the flow is linearly stable becomes independent of the Reynolds number for large enough . The asymptotic value can be estimated numerically, and is roughly . The fact that the critical for linear stability is independent of the Reynolds number for large enough shows that the inviscid limit is a regular limit of this problem. This is, however, in contrast with the findings of Jones (1977) and Lignières, Califano, and Mangeney (1999) for the tanh shear layer. In both cases, they find that for large (albeit using a fairly limited survey of parameter space). Using the LPN equations, Lignières, Califano, and Mangeney (1999) also found that there is no stability threshold in the inviscid limit, that is, unstable modes exist for all values of . The reason for the stark difference between our results and theirs remains to be determined, but could be attributed either to the nature of the boundary conditions used (periodic vs. non-periodic), or to the fact that a sinusoidal velocity profile has shear of both signs while a tanh velocity profile only has shear of one sign.

We now discuss linear stability to 3D perturbations using Squire’s theorem. For the LPN equations, we find that the function is indeed a strictly increasing function of , and a strictly decreasing function of , which implies that the marginal stability of 2D modes is also that of 3D modes (see the previous section). For larger Prandtl number, however, we can immediately see from its null contour that is no longer a monotonic function of which strongly suggests that 3D modes could be the first ones to destabilize the system. Whether this is indeed the case is beyond the scope of this paper, since it belongs to the high-Péclet-number regime. However, this result would be consistent with the work of Smyth and Peltier (1990), who found that 3D modes can be the first ones to be unstable for parallel stratified shear flows which have a tanh profile, albeit in some relatively small region of parameter space.

In preparation for our 3D simulations (see Section V), we are also interested in the spatial structure of the first modes to be destabilized, since the computational domain size must be chosen to be large enough to contain them in order to avoid spurious results. Based on the previous results, we now limit our study entirely to the 2D modes. The range of unstable 2D modes for which marginal stability is achieved is shown in Figure 3, for both the standard equations and for the LPN equations. Again, we see that the results obtained using the LPN equations correctly approximate those obtained using the standard equations at low Péclet number. In all cases, we find that the first mode to be destabilized has , i.e. its horizontal wavelength is larger than the shear lengthscale. For this reason, in the numerical simulations of Section V we use a reasonably long domain size, that can fit at least two wavelengths of the most unstable mode.

Finally, Figure 3 also sheds light on the actual source of the non-monotonicity of the 2D linear stability curves seen in Figure 2 at and above. Indeed, it reveals a new unstable region for low horizontal wavenumber modes, which appears here for ( in this figure). These modes have a growth rate with non-zero imaginary part, by contrast with the standard shearing modes whose growth rates are real3

## Iv Energy stability

### iv.1 Energy stability for stratified shear flows: general ideas

Linear stability only provides information on the stability of a shear flow to infinitesimal perturbations. Energy stability is a much stronger form of stability (e.g. Joseph, 1966, 1976; Drazin and Reid, 2004): when a system is energy stable (also called absolutely stable), perturbations of arbitrarily large amplitudes decay at least exponentially in time, and the laminar flow is the only attractor of the system. Energy stability is thus a sufficient condition for the system to be stable to perturbations of arbitrary amplitude, whereas linear instability is a sufficient condition for the system to be unstable. Often the linear stability limit and the energy stability limit do not coincide in parameter space, an extreme example being the unstratified plane Couette flow, which has a finite threshold Reynolds number for energy stability, but is linearly stable up to infinite Reynolds number: in the region between the two, the system is stable to infinitesimal perturbations, but it may be unstable to perturbations of large amplitude, i.e., it may exhibit finite-amplitude instabilities.

With the goal of further studying the stabilizing effect of background stratification, we now derive an energy stability criterion for forced stratified shear flows. We ask the following question: for a given amplitude of the force, is there a critical strength of the stratification above which the laminar solution is the only attractor of the system?

Again, we insert the decomposition in (7)-(9). However, we do not assume that is small. and then satisfy

 ∂u′∂t+uL⋅∇u′+u′⋅∇uL+u′⋅∇u′=−∇p+RiT′ez+1Re∇2u′ , (32) ∇⋅u′=0 , (33) ∂T′∂t+(uL+u′)⋅∇T′+w′=1Pe∇2T′ , (34)

in the case of the standard equations, and

 ∂u′∂t+uL⋅∇u′+u′⋅∇uL+u′⋅∇u′=−∇p+RiPe∇−2w′ez+1Re∇2u′%, (35) ∇⋅u′=0 , (36)

for the LPN equations.

An energy equation for the perturbations can be obtained by dotting the momentum equation with , adding it to times the temperature equation (where the only constraint on is that it should be a positive scalar), and integrating the result over the domain under consideration. Using the periodicity of the solutions, together with the incompressibility condition greatly simplifies the resulting expressions, which reduce to

 ∂∂t[12⟨u′2⟩+γ2⟨T′2⟩] =(Ri−γ)⟨w′T′⟩−⟨SLw′u′⟩−1Re⟨|∇u′|2⟩−γPe⟨|∇T′|2⟩ , (37) ≡Hγ[u′,T′]

for the standard equations, where denotes a volume integral, and denotes the local vertical shear of the laminar solution. Similarly, for the LPN equations we get

 ∂∂t[12⟨u′2⟩]=RiPe⟨w′∇−2w′⟩−⟨SLw′u′⟩−1Re⟨|∇u′|2⟩≡HLPN[u′] . (38)

This defines the two functionals and , which are both quadratic forms. The task at hand is to determine the region of parameter space or where these quadratic forms are negative definite, i.e., where they are strictly negative for any possible input fields (and ). In this region of parameter space the system is energy stable: the right-hand-side of (37) or (38) is strictly negative and the perturbation decays in time, regardless of its initial amplitude. On the basis of mass and momentum conservation, the only constraints that we place on is that it is divergence-free and has a vanishing average over the whole domain.

Comparing and , we readily see that the task of proving energy stability for stratified shear flows is more involved in the case of the full set of equations than in the case of the LPN equations. We therefore focus on the latter, for which we are able to obtain interesting and generic results on the stability of stratified shear flows.

### iv.2 Energy Stability in the low Péclet number limit

Focussing on the LPN equations, we first compute bounds on the location of the energy stability curve in the () plane. These bounds prove useful, because they correctly describe the scaling behavior of the energy stability limit for large Reynolds number. As we shall see they also validate our numerical results, and provide a simple analytical approximation to the high Reynolds number limit.

#### Lower bound

While the true energy stability boundary can only be obtained by ensuring that for all possible perturbations and , one may also ask the question of when a subset of perturbations is energy stable or unstable. If the subset is unstable, then we know that the system overall is not energy stable either. The critical for energy stability for that subset is therefore a lower bound (called hereafter) on the true energy stability boundary.

We now restrict our attention to perturbations of the following form:

 u′ = Bcos(Ky), (39) v′ = 1Ksin(Ky)sin(z), (40) w′ = cos(z)cos(Ky). (41)

where is the wave number of the perturbation in the direction, and is a free parameter. One can check that such perturbations are divergence-free.

We insert (39)-(41) into the quadratic form (38), recalling that , to obtain

 HLPNLxLyLz = −RiPe4(K2+1)−B4−1Re[K2B22+K2+14(1+1K2)]. (42)

This expression is a quadratic polynomial in . As long as its discriminant is negative, the polynomial is negative and perturbations of the form (39)-(41) decay exponentially. However, when the discriminant is positive, there will be values of corresponding to growing perturbations. The threshold for energy stability of perturbations of the form (39)-(41) is therefore attained when the discriminant vanishes, which gives

 (RiPe)< = K2+18K2Re−(K2+1)2(1+1K2)Re. (43)

For large enough Reynolds number, this value is maximum for the smallest value of that is compatible with the boundary conditions, namely . The corresponding lower bound on the energy stability limit of the system is

 (RiPe)< = ⎡⎢ ⎢ ⎢ ⎢⎣ReL2y32π2−(4π2L2y+1)(1+L2y4π2)Re⎤⎥ ⎥ ⎥ ⎥⎦(4π2L2y+1). (44)

For a given system size, the high- asymptotic behavior of the lower-bound is

 (RiPe)< ≃ L2y32π2(4π2L2y+1)Re, (45)

which shows that the Richardson-Péclet number needs to be at least of the order of to have energy stability. This bound also indicates a strong dependence of the energy stability limit on the size of the domain. Indeed, as the transverse size of the domain increases, expression (45) grows as : larger domains allow for perturbations that are very weakly damped by viscosity.

The lower bound is plotted in Figure 4 for a domain of size , for which

 (RiPe)< = Re4−8Re. (46)

#### Upper bound

Upper bounds on the energy stability limit can be obtained using rigorous estimates of the three terms in . In this subsection, we do not restrict attention to shear flows of the Kolmogorov type. Instead, we consider any smooth shear flow along that is -periodic in . We still use the height of the domain as the characteristic length scale, and consider a domain of size with periodic boundary conditions. We assume that some forcing function with amplitude sustains the laminar flow. The dimensionless profile has an amplitude of order unity. We prove that any such laminar flow is energy stable provided the Richardson-Péclet number is large enough.

To simplify notations and avoid dealing with the inverse Laplacian operator, we introduce , such that . With these notations, the quadratic functional reads

 HLPN[u′]=−RiPe⟨|∇θ|2⟩−⟨duLdzu′∇2θ⟩−1Re⟨|∇u′|2⟩. (47)

Let be the second term of this functional. After integration by parts,

 T=−⟨duLdzu′∇2θ⟩=⟨∇(duLdzu′)⋅∇θ⟩=⟨duLdz∇u′⋅∇θ⟩+⟨d2uLdz2u′∂zθ⟩. (48)

Using classical inequalities, we bound this term according to

 |T| ≤ supz∣∣∣duLdz∣∣∣⟨|∇u′||∇θ|⟩+supz∣∣∣d2uLdz2∣∣∣⟨|u′||∂zθ|⟩ (49) ≤ supz∣∣∣duLdz∣∣∣√⟨|∇u′|2⟩√⟨|∇θ|2⟩+supz∣∣∣d2uLdz2∣∣∣√⟨|u′|2⟩√⟨|∂zθ|2⟩ ≤ RiPe2⟨|∇θ|2⟩+12RiPesupz∣∣∣duLdz∣∣∣2⟨|∇u′|2⟩ +RiPe2⟨|∇θ|2⟩+12RiPesupz∣∣∣d2uLdz2∣∣∣2⟨|u′|2⟩ = RiPe⟨|∇θ|2⟩+12RiPesupz∣∣∣duLdz∣∣∣2⟨|∇u′|2⟩+12RiPesupz∣∣∣d2uLdz2∣∣∣2⟨|u′|2⟩,

where we have used respectively Hölder’s inequality to get the first line, Cauchy-Schwartz inequality to get the second one, and Young’s inequality to get the final expression (see Doering, Spiegel, and Worthing (2000) for another example of the use of these inequalities in a fluid dynamics context). We now wish to express in terms of . To wit, we use Poincaré’s inequality (see Doering, Spiegel, and Worthing, 2000): the divergence-free constraint implies that has vanishing Fourier amplitude on modes with vanishing wave numbers in both the and directions, hence

 ⟨|u′|2⟩≤14π2max{L2y;4π2}⟨|∇u′|2⟩, (50)

where the arguments of the are the maximum allowed values for the squared wavelengths in the and directions. Inserting this inequality into (49), together with , leads to

 HLPN[u′]≤⟨|∇u′|2⟩[12RiPe(supz∣∣∣duLdz∣∣∣2+supz∣∣∣d2uLdz2∣∣∣2max{L2y4π2;1})−1Re]. (51)

As discussed at the beginning of this subsection, the non-dimensionalization is such that

 supz∣∣∣duLdz∣∣∣2 = c1, (52) supz∣∣∣d2uLdz2∣∣∣2 = c2, (53)

where and are constants of order unity that depend on the shape of the laminar profile only (and specifically not on , , , etc). Combining these expressions with (51) leads to

 HLPN[u′]≤⟨|∇u′|2⟩[12RiPe(c1+c2max{L2y4π2;1})−1Re], (54)

hence is a negative quadratic form if the expression inside the square brackets is negative, i.e., if

 RiPe>(RiPe)>=(c1+c2max{L2y4π2;1})Re2. (55)

Because we have used rough but rigorous estimates of the different terms of the quadratic functional, is an upper bound on the actual energy stability limit of the system. It proves that any smooth laminar velocity profile is absolutely stable provided the stratification is strong enough. This has important implications, showing in particular that for fixed Reynolds number and strong enough stratification such a shear flow does not induce any turbulent mixing! Note, however, that since this result is obtained for the LPN equations, it is formally only valid for perturbations that have a low Péclet number. This can be done mathematically by taking the asymptotic limit of the equations for before considering perturbations of arbitrary amplitude. In practice, however, our result does not rule out the possibility of instability for perturbations that have a high Péclet number. A simple example would be perturbations that locally reduce or eliminate the horizontally-averaged vertical stratification the domain: such perturbations are not allowed in the LPN equations, but they are allowed in the full set of equations at very low Péclet number, where they might indeed allow for sustained turbulent solutions localized in the mixed layer.

For domains with large extent in the -direction, the sufficient condition (55) for energy stability becomes approximately . In the particular case of the Kolmogorov velocity profile, we can compare this upper bound to the lower bound (45): both bounds scale as for large Reynolds number, which indicates unambiguously that the actual critical for energy stability obeys the same scaling. This is illustrated in Figure 4, where we plot the upper bound for a Kolmogorov flow in a domain of size . For such a Kolmogorov flow, and the bound becomes

 (RiPe)>=Re. (56)

#### Energy stability boundary using Euler-Lagrange equations

To determine the true energy stability limit of the LPN system, we now consider the variational problem associated with the extremization of the quadratic functional. This gives a set of Euler-Lagrange equations, which can be solved numerically to obtain the critical value of for absolute stability, called hereafter.

Starting from , we separate the viscous dissipation term from the other two, as

 HLPN[u′]=RiPe⟨w′∇−2w′⟩−⟨SLw′u′⟩−D , (57)

where is positive definite for non-trivial flows. We then ask the following question: for fixed viscous dissipation rate , for what values of and is for all incompressible velocity fields ? As we shall see, the selected value of is irrelevant as it merely serves as a general normalization4 of . In order to answer this question, it is sufficient to maximize over all possible incompressible flows , and find out for what values of and this maximum is smaller than . The problem thus reduces to an Euler-Lagrange optimization problem.

Using the notation as in the previous section, we construct the following Lagrangian:

 L[u′]=RiPe⟨w′θ⟩−⟨SLw′u′⟩+⟨π1(x,y,z)∇⋅u′⟩−π2⟨|∇u′|2Re−D0⟩+⟨π3(x,y,z)(w′−∇2θ)⟩ , (58)

where the Lagrange multiplier function enforces incompressibility at every point, enforces equation (13) at every point, and the constant multiplier enforces globally5. Note that we go back here to using the field merely to avoid dealing with inverse Laplacian operators in the variational problem; it is also possible to work through the following derivation without doing it.

The maximizing perturbation field has to solve the Euler-Lagrange equations: three of them (arising from the derivatives of with respect to , and ) simply recover the constraints, and the other four (arising from the derivatives of with respect to , , and ) are

 −SLw′−∂xπ1=−2π2Re∇2u′ , (59) −∂yπ1=−2π2Re∇2v′ , (60) RiPeθ+π3−SLu′−∂zπ1=−2π2Re∇2w′ , (61) RiPew′=∇2π3 . (62)

We see that the multiplier plays a role similar to pressure, a standard result. Comparing the fourth equation with the constraint (13) also reveals that , which implies that we can eliminate both and to get

 2RiPe∇−2w′−S<