Kelvin-Helmholtz billows above Richardson number $1/4$

TL;DR

This study investigates the existence of finite amplitude Kelvin-Helmholtz billows above the classical Richardson number threshold, revealing complex bifurcation structures and new instabilities in stratified shear flows at finite Reynolds numbers.

Contribution

It demonstrates that Kelvin-Helmholtz billows can exist at Richardson numbers greater than 1/4, challenging classical linear stability results, and explores their bifurcation behavior at finite Reynolds numbers.

Findings

Kelvin-Helmholtz billows exist above Ri=1/4 at finite Re.

Bifurcation diagrams reveal complex flow structures and bistability.

A new slow-growing linear instability is identified in hyperbolic tangent stratification.

Abstract

We study the dynamical system of a forced stratified mixing layer at finite Reynolds number $Re$, and Prandtl number $Pr=1$. We consider a hyperbolic tangent background velocity profile in the two cases of hyperbolic tangent and uniform background buoyancy stratifications. The system is forced in such a way that these background profiles are a steady solution of the governing equations. As is well-known, if the minimum gradient Richardson number of the flow, $Ri_m$, is less than a certain critical value $Ri_c$, the flow is linearly unstable to Kelvin-Helmholtz instability in both cases. Using Newton-Krylov iteration, we find steady, two-dimensional, finite amplitude elliptical vortex structures, i.e. `Kelvin-Helmholtz billows', existing above $Ri_c$. Bifurcation diagrams are produced using branch continuation, and we explore how these diagrams change with varying $Re$. In particular, when $Re$ is sufficiently high we find that finite amplitude Kelvin-Helmholtz billows exist at $Ri_m>1/4$, where the flow is linearly stable by the Miles-Howard theorem. For the uniform background stratification, we give a simple explanation of the dynamical system, showing the dynamics can be understood on a two-dimensional manifold embedded in state space, and demonstrate the cases in which the system is bistable. In the case of a hyperbolic tangent stratification, we also describe a new, slow-growing, linear instability of the background profiles at finite $Re$, which complicates the dynamics.

Equations

$${\partial }_{t}u+\left(U+u\right){\partial }_{x}u+w{\partial }_z\left(U+u\right)$$

$${\partial }_{t}w+\left(U+u\right){\partial }_{x}w+w{\partial }_{z}w$$

$${\partial }_{t}b+\left(U+u\right){\partial }_{x}b+w{\partial }_z\left(B+b\right)$$

$${\partial }_{x}u+{\partial }_{z}w$$

$$X(t_0+t)=F\left(X(t_0),t;R{i}_b,Re\right),$$

$$F(X,T;R{i}_b,Re)-X=0$$

$$F(X,T;R{i}_b,Re)-X$$

$$F_X(X,Y,T;R{i}_b,Re)-Y$$

$$Y\cdot A-1$$

$$F(X,T;R{i}_b,Re)-X$$

$$F_X(X,Y_1,T;R{i}_b,Re)-\cos \theta Y_1+\sin \theta Y_2$$

$$F_X(X,Y_2,T;R{i}_b,Re)-\sin \theta Y_1-\cos \theta Y_2$$

$$Y_1\cdot A-1$$

$$Y_2\cdot A$$

$$E=\frac{1}{2L_x}{\int }_0^{L_x}\mathrm{d}x{\int }_{-L_z}^{L_z}\mathrm{d}z\left(u^2+w^2+R{i}_b{b}^2\right).$$

Full text

1 Introduction

The Miles-Howard theorem (Miles, 1961; Howard, 1961) tells us that for inviscid, infinitesimal perturbations to steady, one-dimensional, parallel shear flows, the minimum gradient Richardson number $Ri_{m}$ of the flow must be less than $1/4$ for such ‘linear’ perturbations to grow exponentially. From this, it is often argued that oceanic measurements will always find a Richardson number greater than or equal to $1/4$, otherwise turbulence will ensue (see Smyth et al., 2019, and references therein), despite the very specific restrictions on the applicability of the theorem. In this paper we will examine two aspects of these restrictions, namely that perturbations are infinitesimal, and that $Re=\infty$.

With finite amplitude perturbations, nonlinear effects can no longer be neglected. There is various evidence that for flows susceptible to Kelvin-Helmholtz instability (KHI), complex nonlinear behaviour exists when $Ri_{m}>1/4$. Kaminski et al. (2017) showed that perturbations which grow transiently before decaying in the linearised setting can lead to turbulent-like irreversible mixing with $Ri_{m}>1/4$ when nonlinearity is included. Howland et al. (2018) showed that as $Ri_{m}\to 1/4$ from below, the maximum amplitude of a saturated Kelvin-Helmoltz billow does not tend to zero, but to some finite value. One possible cause of these observations is that the pitchfork bifurcation, generically expected to occur at the critical Richardson number, $Ri_{c}$, is subcritical, so that finite amplitude states exist above $Ri_{c}$. This could mean that the system is bistable in a certain range of $Ri_{m}$ with $Ri_{m}>Ri_{c}$. (Note that by ‘subcritical’ here we mean those regions where the laminar flow is linearly stable, above $Ri_{c}$, consistent with normal dynamical systems terminology, as opposed to the occasional oceanographic usage meaning below $Ri_{c}$.)

Historically, the best way to determine the nature of the bifurcation has been to consider the next order nonlinear effects, a so-called weakly nonlinear analysis. Such analysis was performed by Maslowe (1977) and Brown et al. (1981) for stratified shear layers, finding subcriticality (in the above sense) for $Pr<1$. However, our results suggest that the weakly nonlinear analysis can potentially be misleading, as discussed in section 4, since higher order effects can quickly dominate.

More recently, as it has become possible computationally to solve the Navier-Stokes equations directly, finding the finite amplitude states which arise from bifurcations has emerged as an alternative. Newton’s method can be used to find solutions of nonlinear problems, such as steady states, iteratively. The introduction of Newton-Krylov methods (Edwards et al., 1994), where a Krylov-subspace method such as generalised minimal residuals (GMRES) (Saad & Schultz, 1986) is used to solve the linear system inexactly at each Newton step, has allowed this to be applied to very high dimensional systems for which it is prohibitively expensive to work with matrices directly (for a comprehensive review, see Dijkstra, 2014). It is also possible to use Newton’s method to find and track bifurcation points of high dimensional dynamical systems (Salinger et al., 2002; Haines et al., 2011). Net & Sánchez (2015) used a matrix-free bifurcation tracking technique with a Newton-Krylov method, as employed in this paper, and further extended this to find bifurcations of periodic orbits.

In this paper, we find the exact coherent states that bifurcate from the laminar flow at $Ri_{c}$, and track these as both $Ri_{m}$ and $Re$ vary, to build a picture of the dynamical system near $Ri_{m}=1/4$, and, crucially, answer the question of whether the system can be bistable above $Ri_{c}$. Two different models susceptible to KHI are considered. The first, the ‘Holmboe’ model (Holmboe, 1962), with a hyperbolic tangent buoyancy profile, is the standard model in this field (Hazel, 1972; Klaassen & Peltier, 1985; Smyth & Peltier, 1991; Mallier, 2003), but we demonstrate that complex behaviour arises–associated with what we believe to be a previously unreported linear instability–and dominates at long times when this model is forced onto the system at finite $Re$, obscuring the KHI. We then examine an alternative ‘Drazin’ model (Drazin, 1958), with a uniform stratification, which shares many of the features of the Holmboe model but does not exhibit this complex behaviour. Note that, with the parameters studied, both models are only known to be susceptible to stationary KHI, and not the propagating Holmboe wave instability. The paper proceeds as follows: in section 2, we describe the methodology and code used. In section 3.1 a bifurcation diagram is presented for the Holmboe model, as well as a description of the newly discovered linear instability. In section 3.2, a bifurcation diagram and a full description of the dynamics is given for the Drazin model. Section 4 gives a brief discussion of these results.

2 Methodology

We consider the Boussinesq equations in two dimensions, and study the nonlinear evolution of perturbations away from a steady parallel velocity profile $U(z)$ and buoyancy stratification $B(z)$. Solving for the perturbation away from these constant-in-time profiles is equivalent to solving for the full system, with an artificial body force to counteract diffusion. In non-dimensional form, the equations are:

[TABLE]

Here $u$ is the fluid velocity in the horizontal ($x$) direction, and $w$ is the velocity in the vertical ($z$) direction. Buoyancy acts in the positive $z$ direction. We impose periodic boundary conditions at $x=0$ and $x=L_{x}$, and at $z=\pm L_{z}$ we enforce no-penetration ($w=0$), stress-free ($\partial u/\partial z=0$), and insulating ($\partial b/\partial z=0$) boundary conditions. Given the dimensional shear layer depth, $2L$, velocity difference $2\Delta U$, density difference $2\Delta\rho$, typical density $\rho^{*}$, and diffusivities of momentum $\nu$ and density $\kappa$, the Reynolds number is defined as $Re=\frac{\Delta UL}{\nu}$, the Prandtl number $Pr=\frac{\nu}{\kappa}$, and the bulk Richardson number $Ri_{b}=\frac{g}{\rho^{*}}\frac{L\Delta\rho}{\Delta U^{2}}$. Throughout, we take $Pr=1$ for simplicity. Two different choices of $U$ and $B$ are considered in sections 3.1 and 3.2 respectively. For both background flows studied, the minimum gradient Richardson number $Ri_{m}$, as relevant to the Miles-Howard theorem, is equal to the bulk Richardson number $Ri_{b}$.

2.1 Time-stepping

A new solver was developed to solve the Boussinesq equations around arbitrary background flows. Time integration uses a third order Runge-Kutta-Wray scheme, and spatial derivatives are handled pseudo-spectrally in the periodic horizontal direction, and with explicitly conservative quasi-second order finite differences in the vertical, on a non-uniform staggered grid with more points closer to the central shear layer. The code was validated against DIABLO (Taylor, 2008). Further, a linearised version of the same timestepper was produced, and validated against very low amplitude states in the full nonlinear solver. For the system studied in section 3.1, a grid is used with 256 equispaced points in the streamwise direction, and 512 points in the vertical direction, with a greater density of points in the middle of the domain. For the system studied in section 3.2, 128 points are used in the streamwise direction, covering a shorter domain, and 768 vertically, in order to accurately capture behaviour at higher $Re$. The results are validated by reconverging certain solutions at a higher resolution of $384\times 768$ in section 3.1 and $256\times 1024$ in section 3.2.

2.2 Steady states and bifurcation points

Formally, we may describe our dynamical system as the evolution of a state $X$ by a time $t$ through

[TABLE]

where $Ri_{b}$ and $Re$ are the constant parameters at which we are considering the evolution. Finding steady states of the flow is then equivalent to finding solutions to

[TABLE]

for some arbitrary fixed $T$. A larger $T$ acts to precondition the equations, but if it is too large, computation will be prohibitively expensive. For our system, we found $T=11$ to be a good compromise. It is possible, though extremely unlikely, that this will also find a periodic orbit of period $T$.

Solving (6) is done by using Newton-GMRES (generalised minimum residual) iteration on an initial guess. Our implementation closely matches that employed by Chandler & Kerswell (2013), including the use of a trust region to make the algorithm globally convergent. The GMRES iteration at each Newton step is continued until the residual is less than $10^{-2}$, and the Newton iteration is continued until its residual, the norm of the left-hand side of (6), is less than $10^{-8}$. Through trial and error, we converge a steady state solution, the result of a very long time integration of equations (1-4), at $Re=1000$ and $Ri_{b}=0.2$, in both the flows studied in this paper. Once one state is found at these particular $Ri_{b}$ and $Re$, we converge another very close by at a different $Ri_{b}$ but the same $Re$. We then follow the solution branch at this $Re$ over a range of $Ri_{b}$ using pseudo-arclength continuation (Keller, 1977). We examine the stability of the branch with Arnoldi iteration, using a linearised version of the same timestepping code.

The stability analysis reveals the existence of bifurcation points, where eigenvalues of the state cross a stability boundary. To continue these bifurcation points to different $Re$, we use the states found by stability analysis as an initial guess in a different iterative solver. The system we solve is similar to that implemented in LOCA (Salinger et al., 2002), but we use a matrix free method, as discussed in detail in Sánchez & Net (2016). We look for solutions to

[TABLE]

with Newton-GMRES. In this case we allow $X$, $Y$ and $Ri_{b}$ to be found by the iteration, but hold $Re$ fixed. Here $F_{X}(X,Y,t;Ri_{b},Re)$ is the linearised time evolution of a state $Y$ about a nonlinear state $X$, computed using the linearised timestepper. Equation (7b) enforces that $Y$ is a neutral eigenmode of the Jacobian at $X$. We normalise $Y$ using (7c), with some fixed arbitrary state $A$. Once bifurcation points are found at a particular $Re$, they are reconverged at higher $Re$. We are particularly interested in how the $Ri_{b}$ value of the bifurcation point varies with $Re$.

Equations (7) find bifurcation points with purely real neutral eigenmodes, i.e. pitchfork and saddle-node bifurcations. For Hopf bifurcations, a set of five equations is needed, including two different linearised time evolutions. These arise from the real and imaginary parts of the eigenvalue $e^{i\theta}$. The following are solved for the unknowns $X$, $Y_{1}$, $Y_{2}$, $Ri_{b}$ and $\theta$:

[TABLE]

The additional computational requirements of (8) mean that we are unable to track Hopf bifurcations to as high Reynolds numbers as pitchfork and saddle-node bifurcations.

3 Results

3.1 Hyperbolic tangent stratification: the Holmboe model

First we consider a background profile of $U=\tanh{z}$, $B=\tanh{z}$. This is a commonly used model of a mixing layer, introduced by Holmboe (1962). It has the useful property that, at infinite $Re$, the linear stability analysis can be performed analytically by hand (Miles, 1963). With this choice, we find that the minimum gradient Richardson number $Ri_{m}$ is equal to $Ri_{b}$, and so the Miles-Howard theorem tells us that the flow is certainly stable for $Ri_{b}>1/4$. We choose $L_{x}=4\pi$, which is one wavelength of the most unstable mode at $Ri_{b}=1/4$ as $Re\to\infty$, assuming a domain of infinite height vertically. We take $L_{z}=10$. This is a compromise between being large enough that the boundaries do not significantly affect the flow but small enough to keep computation costs down. The finite value of $L_{z}$ means that $Ri_{c}$ tends to a value slightly less than $1/4$ as $Re\to\infty$.

Following Howland et al. (2018) we define the energy of perturbations to be

[TABLE]

State space is taken as the space of all possible incompressible perturbation flows $X=\left({u},{w},{b}\right)$, with norm $\left\|X\right\|:=\sqrt{2E}$. Note that $p$ is not a dynamical variable as it can be calculated from a Poisson equation forced by the velocity field.

Figure 1(a) shows a bifurcation diagram at $Re=4000$. Where the background state becomes unstable to KHI at $Ri_{c}\approx 0.2494$, a pitchfork bifurcation occurs (the green dot on figure 1(a)), giving rise to a branch of finite amplitude, billow-like states. This branch is initially stable–except to the unrelated instability discussed below–and decreasing in $Ri_{b}$, but there is soon a saddle-node bifurcation (see inset in figure 1(a)) and it then increases in $Ri_{b}$. As the unstable branch increases in amplitude, $Ri_{b}$ increases, and we find steady, though unstable, states above $Ri_{b}=1/4$. There is another saddle-node bifurcation at $Ri_{s}$ (blue dot), adding a second unstable direction to the branch. The Hopf bifurcation, at $Ri_{h}\approx 0.244$ (red dot), has two neutral eigenmodes which are connected with the eigenmodes of the two saddle-node bifurcations via a periodic orbit.

A very weak linear instability, apparently hitherto unreported, is present in all states on the bifurcation diagram. Figure 1(b) shows the maximum growth rate of linear instability of the background state, as $Ri_{b}$ and $Re$ vary. For $Ri_{b}>1/4$, the new instability is the dominant one. This has a phase speed of less than one, and manifests as convective rolls, advected through the domain, above and below the interface at some critical layer, as shown in figure 2. As $Re\to\infty$, the growth rate tends to zero, as required by the Miles-Howard theorem. Close agreement of growth rates, to one part in $10^{3}$, was found for this instability between the Arnoldi stability algorithm of our code, and a direct solution of the Orr-Sommerfeld equations, using a MATLAB code by W. D. Smyth. Despite the small growth rate, at long times this instability leads to significant nonlinear behaviour in the forced problem, which eventually dominates and obscures all signature of KHI below $Ri_{b}=1/4$. This means that the dynamical behaviour of KHI is very difficult to describe in this model.

3.2 Uniform stratification: the Drazin model

We now consider the case with a uniform background stratification, so that $U=\tanh{z}$ but $B=z$. This is also a commonly studied problem (Drazin, 1958; Churilov & Shukhman, 1987; Kaminski et al., 2014) as again, linear stability analysis can be performed analytically. As before, $Ri_{m}=Ri_{b}$ for this flow. Linear stability analysis on a domain of infinite height tells us we should now take $L_{x}=2\sqrt{2}\pi$ to achieve $Ri_{c}\to 1/4$ as $Re\to\infty$. As before, $L_{z}=10$. We use the same definition of energy $E$ as in the hyperbolic tangent case, though we note that only for uniform stratification does this form of energy correspond precisely to the sum of perturbation kinetic and potential energies.

Qualitatively, the bifurcation diagram is very similar to the $\tanh$ stratification case. Figure 3(a) shows the diagram for $Re=4000$. The main difference from figure 1(a) is the lack of the first saddle-node bifurcation near the pitchfork. The values of the various bifurcation Richardson numbers are different, for example the Hopf bifurcation at $Ri_{h}$ (shown in red) occurs at somewhat lower $Ri_{b}$ than before. Also crucially, the propagating linear instability described in section 3.1 is no longer present, and consequently we can study the long-time behaviour of KHI.

The period of the Hopf bifurcation at $Re=4000$ is about $1690$ advective time units, which is much too high to allow us to converge the resulting periodic orbit directly, but long time integrations at a range of $Ri_{b}$ give us an idea of the behaviour, since it appears to be stable in this case. Even this simple method becomes useless as we approach $Ri_{c}$, since the period increases towards infinity. This is the generic behaviour near a homoclinic bifurcation (Strogatz, 2014), which we believe occurs somewhere between $Ri_{c}$ and $Ri_{s}$: the periodic orbit collides with the lower branch state.

The behaviour of the system, which is generic for sufficiently high $Re$, can be completely understood on a two-dimensional manifold described by the two most unstable eigenmodes, as shown schematically in figure 4. In region 1, where $Ri_{b}<Ri_{h}$, the laminar state is unstable, and the instability saturates and eventually leads to the upper branch state, which is stable. For $Ri_{h}<Ri_{b}<Ri_{c}$, region 2, the laminar state and upper branch are both unstable, and perturbations lead to a stable periodic orbit. Immediately above the pitchfork bifurcation $Ri_{c}$ in the region 3, the laminar state is stable and there exists a lower branch edge state, which is unstable. If finite amplitude perturbations to the laminar state are past this edge, they are attracted to the periodic orbit, and we have subcritical ‘transition’. Region 4 is between the homoclinic bifurcation of the periodic orbit with the lower branch state, and the saddle-node bifurcation of the lower and upper branches. Here, there are unstable finite amplitude states and large transient trajectories, but the laminar state is the only attractor. In region 5, past the saddle-node bifurcation, $Ri_{b}>Ri_{s}$, the laminar state is the only known exact coherent structure. Of course, in reality the finite amplitude states break the translational symmetry of the laminar state, and there are in fact a continuum of upper branch states, periodic orbits and so on, with a shift of origin. Which of these the system is attracted to depends on the phase of the initial perturbation.

Figure 5 shows the vorticity structure of the steady states at two different values of $Ri_{b}$. In the case of the Hopf bifurcation, billow-like structures are clearly seen, bearing a strong resemblance to the saturated, unsteady billows found by Howland et al. (2018). Increasing $Ri_{b}$ along the upper branch to the saddle-node bifurcation, these structures remain but become significantly less pronounced. Baroclinic effects mean that the height of the billows decreases with increasing $Ri_{b}$.

We track the values of $Ri_{c}$ and $Ri_{s}$ for $Re$ from $1000$ to $10000$ using the method described in section 2.2, and the results are shown on figure 3(b). As $Re\to\infty$, extrapolation, assuming linearity in $1/Re$, suggests $Ri_{c}\to 0.25-1.4\times 10^{-5}$, slightly less than $1/4$ because of the finite height of the domain. Nevertheless, for $Re\gtrsim 9000$, we find that $Ri_{s}>1/4$. It has been difficult to extrapolate $Ri_{s}$ convincingly, though at $Re=\infty$, it seems likely that $Ri_{s}\approx 0.251$. Since we have been unable to find the location of the conjectured homoclinic bifurcation, we are unable to say whether region 3, with a stable periodic orbit, extends above $Ri_{b}=1/4$, and hence whether the system is bistable here. Nevertheless, region 4 certainly exists above $Ri_{b}=1/4$, so there will be nonlinear transient behaviour, with the development of Kelvin-Helmholtz style billows as shown in figure 5. We have also tracked the Hopf bifurcation (omitted from figure 3(b) for scale reasons) and this shows a similar trend to the saddle-node bifurcation.

4 Discussion and Conclusions

The Miles-Howard theorem is an important result in the theory of linear stability of inviscid flows. However, the fact it seems to work in more general conditions than those for which it is proven means it has been informally applied as a ‘rule of thumb’ at high $Re$. We have shown that subcritical instability can exist in such flows, so that complex nonlinear behaviour can occur even when the flow is linearly stable. This is not a new result; Maslowe (1977) found subcritical instability in the Holmboe model with $Pr=0.72$ and $Re=100$ using a weakly nonlinear analysis. We note however, that this technique of finding the first order correction to the linear theory would have given misleading results applied to the parameters we study, since in the Holmboe model, we find a saddle-node bifurcation very close to the pitchfork, leading to subcritical instability instead of the apparent supercriticality. Furthermore, the technique presented in this paper allows us to precisely find the location of the saddle-node bifurcations, and demonstrate explicitly that finite ampltiude states exist at $Ri_{b}>1/4$, which has only been inferred previously.

We have been able to give a simple description of the dynamics in the Drazin model. It is not immediately clear that the dynamics of the forced system studied here will be relevant to those of an unforced system, which has traditionally been used as a model for geophysical flows. The incredibly long periods of the orbits born from the Hopf bifurcation discussed earlier, for example, mean that in an unforced problem, the background flow would have diffused almost entirely away before one complete cycle. Nevertheless, the instability of the unforced flows still leads to saturated states very similar to the steady solutions we have found, and the subcriticality we have demonstrated would certainly lead to nontrivial transient behaviour.

Despite our results, it seems that $Ri_{b}=1/4$ is indeed a useful ‘rule of thumb’ for stability in physical flows. The subcriticality we have found extends only very slightly about $1/4$ in both cases studied. However, this is in apparent disagreement with the results of Maslowe (1977), who found large subcriticality. This is likely because of the different value of $Pr$ studied. Indeed, in a follow-up work Brown et al. (1981) found subcriticality when $Pr<1$ but supercriticality when $Pr>1$, and showed that higher order terms must be considered at our choice of $Pr=1$, so we would like to extend this work to the more oceanographically relevant range $Pr\sim O(10)$.

In addition to these finite amplitude nonlinear states, we have found linear instability with $Ri_{b}>1/4$ in the Holmboe model (see figure 2), which disappears as $Re\to\infty$, as required by the Miles-Howard theorem. A similar phenomenon was found by Miller & Lindzen (1988). However, their instability had large growth rates and required a carefully constructed flow. We have found an instability in a widely used model, hitherto unreported to the best of our knowledge. The new instability has a tiny growth rate at physically realistic $Re$. This suggests it can be ignored in oceanic problems, but does not entirely explain why it has not been discussed before. It is commonly assumed that finite $Re$ effects are always stabilising compared to inviscid behaviour. This instability demonstrates that such assumptions should be checked carefully. While it is not appropriate to classify this instability as ‘classic’ Holmboe wave instability, since it lacks the characteristic ‘wave-interaction’ resonance between an interfacial gravity wave and vorticity waves localised at the edge of the shear layer, we conjecture that it may be homotopically connected to Holmboe instability as parameters are varied, since it has a similar phase speed and occurs at similar values of $Ri_{b}$. This is an area for future research.