III Hydrodynamic Stability - Linear stability analysis

1Linear stability analysis

III Hydrodynamic Stability

1.5 Stratified flows

In the Kelvin–Helmholtz scenario, we had a varying density. Let’s now try to

model the situation in complete generality.

Recall that the set of (inviscid) equations for Boussinesq stratified fluids is



∂u

∂t

+ u · ∇u



= −∇p

− gρ

with incompressibility equations

∇ · u = 0,

Dρ

= 0.

We consider a mean base flow and a 2D perturbation as a normal mode:

u =

U(z)

x + u

(x, z, t)

p = ¯p(z) + p

(x, z, t)

ρ = ¯ρ(z) + ρ

(x, z, t),

with

, p

, ρ

) = [

u(z), ˆp(z), ˆρ(z)]e

i(kx−ωt)

= [

u(z), ˆp(z), ˆρ(z)]e

ik(x−ct)

We wish to obtain an equation that involves only

ˆw

. We first linearize our

equations to obtain

¯ρ



∂u

∂t

∂u

∂x

+ w



= −

∂p

∂x

¯ρ



∂w

∂t

∂w

∂x



= −

∂p

∂z

− gρ

∂ρ

∂t

∂ρ

∂x

+ w

d¯ρ

= 0

∂u

∂x

∂w

∂z

= 0.

Plugging in our normal mode solution into the four equations, we obtain

ik¯ρ(

U − c)ˆu + ¯ρ ˆw

U = −ikˆp

ik¯ρ(

U − c) ˆw = −

dˆp

− g ˆρ. (∗)

ik(

U − c)ˆρ + w

d¯ρ

= 0 (†)

ikˆu +

∂w

∂z

= 0.

The last (incompressibility) equation helps us eliminate

ˆu

from the first equation

to obtain

−¯ρ(

U − c)

d ˆw

+ ¯ρ ˆw

= −ikˆp.

To eliminate

ˆp

as well, we differentiate with respect to

and apply (

∗

), and then

further apply (†) to get rid of the ˆρ term, and end up with

−¯ρ(

U − c)

ˆw

+ ¯ρ ˆw

+ k

¯ρ(

U − c) ˆw = −

g ˆw

U − c

d¯ρ

This allows us to write our equation as the Taylor–Goldstein equation



− k



ˆw −

ˆw

(

U − c)

ˆw

(

U − c)

= 0,

where

= −

¯ρ

d¯ρ

This

is known as the buoyancy frequency, and has dimensions

−2

. This is

the frequency at which a slab of stratified fluid would oscillate vertically.

Indeed, if we have a slab of volume

, and we displace it by an infinitesimal

amount ζ in the vertical direction, then the buoyancy force is

F = V gζ

∂ρ

∂z

Since the mass of the fluid is

, by Newton’s second law, the acceleration of

the parcel satisfies



−

∂ρ

∂z



ζ = 0.

This is just simple harmonic oscillation of frequency N .

Crucially, what we see from this computation is that a density stratification

of the fluid can lead to internal waves. We will later see that the presence of

these internal waves can destabilize our system.

Miles–Howard theorem

In 1961, John Miles published a series of papers establishing a sufficient condition

for infinitesimal stability in the case of a stratified fluid. When he submitted

this to review, Louis Howard read the paper, and replied with a 3-page proof of

a more general result. Howard’s proof was published as a follow-up paper titled

“Note on a paper of John W. Miles”.

Recall that we just derived the Taylor–Goldstein equation



− k



ˆw −

ˆw

(

U − c)

ˆw

(

U − c)

= 0,

The magical insight of Howard was to introduce the new variable

H =

(

U − c)

1/2

We can regretfully compute the derivatives

ˆw = (

U − c)

1/2

ˆw =

(

U − c)

−1/2

+ (

U − c)

1/2

ˆw = −

(

U − c)

−3/2





(

U − c)

−1/2

+ (

U − c)

−1/2

+ (

U − c)

1/2

We can substitute this into the Taylor–Goldstein equation, and after some

algebraic mess, we obtain the decent-looking equation



(

U − c)



− H







(

U − c) +





− N

U − c







= 0.

This is now self-adjoint. Imposing boundary conditions

ˆw,

d ˆw

→

0 in the far

field, we multiply by H

∗

and integrate over all space. The first term gives

∗



(

U − c)



dz = −

(

U − c)



dz,

while the second term is given by







−k

|H|

(

U − c) −

|H|

− |H|







− N



(

U − c

∗

)

U − c|







dz.

Both the real and imaginary parts of the sum of these two terms must be zero.

In particular, the imaginary part reads









+ k

|H|

+ |H|

−





U − c|







dz = 0.

So a necessary condition for instability is that

−





< 0

somewhere. Defining the Richardson number to be

Ri(z) =

U/dz)

the necessary condition is

Ri(z) <

Equivalently, a sufficient condition for stability is that Ri(z) ≥

everywhere.

How can we think about this? When we move a parcel, there is the buoyancy

force that attempts to move it back to its original position. However, if we move

the parcel to the high velocity region, we can gain kinetic energy from the mean

flow. Thus, if the velocity gradient is sufficiently high, it becomes “advantageous”

for parcels to move around.

Piecewise-linear profiles

If we want to make our lives simpler, we can restrict our profile with piecewise

linear velocity and layered density. In this case, to solve the Taylor–Goldstein

equation, we observe that away from the interfaces, we have

= 0. So

the Taylor–Goldstein equation reduces to



− k



ˆw = 0.

This has trivial exponential solutions in each of the individual regions, and to

find the correct solutions, we need to impose some matching conditions at the

interface.

So assume that the pressure is continuous across all interfaces. Then using

the equation

−¯ρ(

U − c)

d ˆw

+ ρ ˆw

= −ikˆp,

we see that the left hand seems like the derivative of a product, except the sign

is wrong. So we divide by

(

U−c)

, and then we can write this as

¯ρ

ˆp

(

U − c)



ˆw

U − c



For a general

, we integrate over an infinitesimal distance at the interface. We

assume

ˆp

is continuous, and that

¯ρ

and (

U − c

) have bounded discontinuity.

Then the integral of the LHS vanishes in the limit, and so the integral of the

right-hand side must be zero. This gives the matching condition



ˆw

U − c



−

= 0.

To obtain a second matching condition, we rewrite the Taylor–Goldstein equa-

tion as a derivative, which allows for the determination of the other matching

condition:



(

U − c)

d ˆw

− ˆw

−

g ¯ρ



ˆw

U − c



= k

(

U − c) ˆw −

g ¯ρ



ˆw

U − c



Again integrating over an infinitesimal distance over at the interface, we see

that the LHS must be continuous along the interface. So we obtain the second

matching condition.



(

U − c)

d ˆw

− ˆw

−

g ¯ρ



ˆw

U − c



−

= 0.

We begin by applying this to a relatively simple profile with constant density,

and whose velocity profile looks like

∆U

For convenience, we scale distances by

and speeds by

∆U

, and define

˜c

and α by

c =

∆U

˜c, α =

These quantities

and

˜c

(which we will shortly start writing as

instead) are

nice dimensionless quantities to work with. After the scaling, the profile can be

described by

z = 1

z = −1

U = −1

U = z

U = 1

III

ˆw = Ae

−α(z−1)

ˆw = Be

αz

+ Ce

−αz

ˆw = De

α(z+1)

We have also our exponential solutions in each region between the interfaces,

noting that we require our solution to vanish at the far field.

We now apply the matching conditions. Since

U − c

is continuous, the first

matching condition just says ˆw has to be continuous. So we get

A = Be

+ Ce

−α

D = Be

−α

+ Ce

The other matching condition is slightly messier. It says



(

U − c)

d ˆw

− ˆw

−

g ¯ρ



ˆw

U − c



−

= 0,

which gives us two equations

(Be

+ Ce

−α

)(1 − c)(−α) = (Be

− Ce

−α

)(1 − c)(α) − (Be

+ Ce

−α

)

(Be

−α

+ Ce

)(−1 − c)(α) = (Be

−α

− Ce

)(−1 − c)(α) − (Be

−α

+ Ce

Simplifying these gives us

(2α(1 − c) − 1)Be

= Ce

−α

(2α(1 + c) − 1)Ce

= Be

−α

Thus, we find that

(2α − 1)

− 4α

= e

−4α

and hence we have the dispersion relation

(2α − 1)

− e

−4α

4α

We see that this has the possibility of instability. Indeed, we can expand the

numerator to get

(1 − 4α + 4α

) − (1 − 4α + 8α

+ O(α

))

4α

= −1 + O(α).

So for small

, we have instability. This is known as Rayleigh instability. On

the other hand, as

grows very large, this is stable. We can plot these out in a

graph

0.64

We see that the critical value of

is 0

64, and the maximally instability point

is α ≈ 0.4.

Let’s try to understand where this picture came from. For large

, the

wavelength of the oscillations is small, and so it seems reasonable that we can

approximate the model by one where the two interfaces don’t interact. So

consider the case where we only have one interface.

z = 1

U = z

U = 1

−α(z−1)

α(z−1)

We can perform the same procedure as before, solving



ˆw

U − c



−

= 0.

This gives the condition A = B. The other condition then tell us

(1 − c)(−α)A = (1 − c)αA − Ae

α(z−1)

We can cancel the A, and rearrange this to say

c = 1 −

2α

So we see that this interface supports a wave at speed

= 1

−

2α

at a speed

lower than

. Similarly, if we treated the other interface at isolation, it would

support a wave at c

−

= −1 +

2α

Crucially these two waves are individually stable, and when

is large, so that

the two interfaces are effectively in isolation, we simply expect two independent

and stable waves. Indeed, in the dispersion relation above, we can drop the

−4α

term when α is large, and approximate

c = ±

2α − 1

2α

= ±



1 −

2α



When

is small, then we expect the two modes to interact with each other,

which manifests itself as the

−e

4α

term. Moreover, the resonance should be the

greatest when the two modes have equal speed, i.e. at

α ≈

, which is quite

close to actual maximally unstable mode.

Density stratification

We next consider the case where we have density stratification. After scaling,

we can draw our region and solutions as

z = −1

z = 0

z = 1

U = −1

U = z

U = 1

III

−α(z−1)

α(z−1)

+ Ce

−α(z−1)

α(z+1)

+ Ee

−α(z+1)

F e

α(z+1)

¯ρ = −1

¯ρ = +1

Note that the

¯ρ

is the relative density, since it is the density difference that

matters (fluids cannot have negative density).

Let’s first try to understand how this would behave heuristically. As before,

at the I-II and III-IV interfaces, we have waves with



1 −

2α



. At the

density interface, we previously saw that we can have internal gravity waves

with c

igw

= ±

, where

is the bulk Richardson number (before scaling, it is Ri

g∆ρh

∆U

We expect instability to occur when the frequency of this internal gravity

wave aligns with the frequency of the waves at the velocity interfaces, and this

is given by



1 −

2α



It is easy to solve this to get the condition

' α − 1.

This is in fact what we find if we solve it exactly.

So let’s try to understand the situation more responsibly now. The procedure

is roughly the same. The requirement that ˆw is continuous gives

A = B + C

−α

+ Ce

= De

+ Ee

−α

D + E = F.

If we work things out, then the other matching condition



(

U − c)

d ˆw

− ˆw

− Ri

¯ρ



ˆw

U − c



−

= 0

gives

(1 − c)(−α)A = (1 − c)(α)(B − C) − (B + C)

(−1 − c)(α)F = (−1 − c)(α)(D − E) − (D + E)

(−c)(α)(Be

−α

− Ce

) +

(−c)

(Be

−α

+ Ce

)

= (−c)(α)(De

− Ee

−α

) −

(−c)

(De

+ Ee

−α

This defines a 6

6 matrix problem, 4th order in

. Writing it as

= 0 with

= (

A, B, C, D, E, F

)

, the existence of non-trivial solutions is given by the

requirement that det S = 0, which one can check (!) is given by

+ c



−4α

− (2α − 1)

4α

−





−2α

+ (2α − 1)

2α



= 0.

This is a biquadratic equation, which we can solve.

We can inspect the case when

is very small. In this case, the dispersion

relation becomes

+ c



−4α

− (2α − 1)

4α



Two of these solutions are simply given by

= 0, and the other two are just those

we saw from Rayleigh instability. This is saying that when the density difference

is very small, then the internal gravitational waves are roughly non-existent, and

we can happily ignore them.

In general, we have instability for all

! This is Holmboe instability. The

scenario is particularly complicated if we look at small

, where we expect the

effects of Rayleigh instability to kick in as well. So let’s fix some 0

< α <

64,

and see what happens when we increase the Richardson number.

As before, we use dashed lines to denote the imaginary parts and solid lines

to denote the real part. For each

, any combination of imaginary part and

real part gives a solution for

, and, except when there are double roots, there

should be four such possible combinations.

We can undersatnd this as follows. When

= 0, all we have is Rayleigh

instability, which are the top and bottom curves. There are then two

= 0

solutions, as we have seen. As we increase

to a small, non-sero amount, this

mode becomes non-zero, and turns into a genuine unstable mode. As the value

increases, the two imaginary curves meet and merge to form a new curve.

The solutions then start to have a non-zero real part, which gives a non-zero

phase speed to our pertrubation.

Note that when

becomes very large, then our modes become unstable

again, though this is not clear from our graph.

Taylor instability

While Holmboe instability is quite interesting, our system was unstable even

without the density stratification. Since we saw that instability is in general trig-

gered by the interaction between two interfaces, it should not be surprising that

even in a Rayleigh stable scenario, if we have two layers of density stratification,

then we can still get instability.

Consider the following flow profile:

z = 1

z = −1

U = z

III

−α(z−1)

αz

+ Ce

−αz

α(z+1)

¯ρ = R − 1

¯ρ = R + 1

As before, continuity in ˆw requires

A = Be

+ Ce

−α

D = Be

−α

+ Ce

Now there is no discontinuity in vorticity, but the density field has jumps. So

we need to solve the second matching condition. They are

(1 − c)(−α)(Be

+ Ce

−α

) +

1 − c

(Be

+ Ce

−α

) = (1 − c)(α)(Be

− Ce

−α

)

(1 − c)(α)(Be

−α

+ Ce

) +

1 + c

(Be

−α

+ Ce

) = (−1 − c)(α)(Be

− Ce

)

These give us, respectively,

(2α(1 − c)

− Ri

)B = Ri

−2α

(2α(1 + c)

− Ri

)C = Ri

−2α

So we get the biquadratic equation

− c



2 +





1 −

2α



−

−4α

4α

= 0.

We then apply the quadratic formula to say

= 1 +

2α

2Ri

−4α

4α

So it is possible to have instability with no inflection! Indeed, instability occurs

when c

< 0, which is equivalent to requiring

2α

1 + e

−2α

< Ri

2α

1 − e

−2α

Thus, for any

0, we can have instability. This is known as Taylor

instability.

Heuristically, the density jumps at the interface have speed

igw

= ±1 ∓



2α



1/2

We get equality when

= 0, so that

= 2

. This is in very good agreement

with what our above analysis gives us.