Hydrodynamic stability: Parrallel stratified shear flows

RSES, ANU, 2018

The equations of motion

Fluids are governed by the Navier-Stokes equations. These describe the evolution of the fluid variables $\boldsymbol{u}(\boldsymbol{x}, t)$, $\rho(\boldsymbol{x}, t)$, and $T(\boldsymbol{x}, t)$. Under some simplifications these read:

\begin{align} \partial_t\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} = -\frac{\boldsymbol{\nabla}p}{\rho} + \nu \nabla^2\boldsymbol{u} + \frac{\boldsymbol{F}}{\rho} &\quad\text{N-S momentum equation (Newton's second law)}\label{1}\tag{1}\\ \partial_t \rho + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}\rho) = 0&\quad\text{conservation of mass}\label{2}\tag{2}\\ \partial_t T + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} T = \kappa\nabla^2 T &\quad\text{thermodynamic or "heat" equation}\label{3}\tag{3}\\ p = \mathrm{function}(\rho, T) &\quad\text{equation of state}\label{4}\tag{4} \end{align}

Hereafter, $\boldsymbol{u}\equiv (u,v,w)$ and $\boldsymbol{x} \equiv (x, y, z)$. Sometimes we may write $D/Dt$ as a shorthand for the advective derivative, i.e., $$ \frac{D\phi}{Dt} \equiv \big(\partial_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\big)\phi. $$

Above in the thermodynamic equation we have neglected on the right-hand-side the forcing term that comes from conversion of kinetic energy into heating. This is very good assumption for (most of) the ocean or for fluids in the lab. See the discussion in the book by Salmon (1998) or the paper by Spiegel & Veronis (1959).

The term $\boldsymbol{F}$ on the right-hand-side of the momentum equation denotes external forces per unit volume. For example, gravitational force appears as $\boldsymbol{F} = -\rho g \widehat{\boldsymbol{z}}$.

Boussinesq approximation

A very usefull and often applicable approximation to the full Navier-Stokes equations under the influence of gravity is the, so-called, Boussinesq approximation. The approximaton exploit the smallness of density variations in liquids. That is, the approximation is valid if we can decompose the density field as

$$ \rho(\boldsymbol{x}, t) = \rho_m + \rho_0(z) + \rho'(\boldsymbol{x}, t)\label{5}\tag{5} $$

with \begin{equation} |\rho'|\ll|\rho_0|\ll|\rho_m|. \label{6}\tag{6} \end{equation}

If variations of density are very small, then it is also reasonable to assume that the equation of state is linear in all arguments. Thus: $$ \rho = \rho_m\left( 1-\frac{T-T_m}{T_m} +\frac{p-p_m}{p_m} \right). \label{7}\tag{7} $$

Above, $\rho_m$ is the spatial and time average of $\rho$: $\rho_m = \langle \rho \rangle^{x,y,z,t}$, where with superscript on angle brackets we denote averages. The vertical profile $\rho_0(z)$ is the $x$, $y$ and time-average of $\rho-\rho_m$: $\rho_0(z) = \langle \rho - \rho_m \rangle^{x,y,t}$. Last, $\rho'$ is what's left.

For the Earth's ocean (6) is a very good since the density components are about two orders of magnitude appart.

Let's go here through a quick derivation of the Boussinesq equations. For more detailed and carefull derivation of the Boussinesq approximation refer to the books by Salmon (1998) and Vallis (2017) or in the short and elegant paper by Spiegel and Veronis (1959). Note that this was Ed Spiegel's summer project at GFD summer school 1959 --- the first GFD summer school held at WHOI.

To the background density field $\rho_m+\rho_0(z)$ there corresponds a background pressure field such that pressure balances buoyance when the fluid is at rest ($\boldsymbol{u}=0$). From (1), setting $\boldsymbol{u}=0$ we get that the background pressure $p_0$ satisfies: $$ \partial_z p_0(z) = - g \big[\rho_m+\rho_0(z)\big]. $$

Only departures from $\rho_m+\rho_0(z)$ and $p_0(z)$ will induce fluid motion. (Actually $p_0(z)$ is determined up to a constant but a constant does not matter for pressure gradients.)

Consider now the full density and pressure fields. From (1) if we multiply by $\rho$ we get:

$$ (\rho_b+\rho') \left(\partial_t\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}\right) = -\boldsymbol{\nabla}p_0 - \boldsymbol{\nabla}p' + (\rho_m+\rho_0+\rho')\nu \nabla^2\boldsymbol{u} - g \rho_b\widehat{\boldsymbol{z}} - g\rho'\widehat{\boldsymbol{z}} $$

The first and fourth terms on the right-hand-side cancel. We factor out $\rho_m$ and get: $$ \left(1+\frac{\rho_0}{\rho_m}+\frac{\rho'}{\rho_m}\right) \left(\partial_t\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}\right) = - \frac{\boldsymbol{\nabla}p'}{\rho_m} + \left(1+\frac{\rho_0}{\rho_m}+\frac{\rho'}{\rho_m}\right)\nu \nabla^2\boldsymbol{u} - g\frac{\rho'}{\rho_m}\widehat{\boldsymbol{z}} $$

Using (5) we simplify terms $1+\rho_0/\rho_m+\rho'/\rho_m$ to unity: $$ \partial_t\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} = - \frac{\boldsymbol{\nabla}p'}{\rho_m} + \nu \nabla^2\boldsymbol{u} - g\frac{\rho'}{\rho_m}\widehat{\boldsymbol{z}}. \label{8}\tag{8} $$

The mass conservation equation reads:

$$ \partial_t \rho' + \rho_m \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0\boldsymbol{u}) + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho'\boldsymbol{u}) = 0, $$

but again using (5) we can deduce that term $\rho_m \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}$ is much larger than all other terms in the above. Thus, to a good approximation the mass conservation becomes simply:

$$ \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0. \label{9}\tag{9} $$

The thermodynamic equation derivation is a bit more elaborate and outside of our scope here. It turns out that from (3) using decomposition (5) and the linear equation of state (7) we have

$$ \partial_t \rho' + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \rho' + w \frac{\mathrm{d}\rho_0}{\mathrm{d}z}= \kappa\nabla^2 \rho'. \label{10}\tag{10} $$

The term $\mathrm{d}\rho_0/\mathrm{d}z$ is often written in terms of the Brunt-Väisälä frequency, which up to the approximations already made in (3) is defined as:

$$N^2(z)\equiv -\frac{g}{\rho_m}\frac{\mathrm{d}\rho_0}{\mathrm{d}z}.$$

Gathering everything together, the Boussinesq equations are:

\begin{align} \partial_t\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} = -\boldsymbol{\nabla}\left(\frac{p'}{\rho_m}\right) + \nu \nabla^2\boldsymbol{u} - g\frac{\rho'}{\rho_m}\widehat{\boldsymbol{z}} &\quad\text{N-S momentum equation (Newton's second law)} \tag{11}\\ \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0&\quad\text{conservation of mass}\tag{12}\\ \partial_t \rho' + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \rho' + w \frac{\mathrm{d}\rho_0}{\mathrm{d}z}= \kappa\nabla^2 \rho' &\quad\text{thermodynamic equation}\tag{13}\\ \rho = \rho_m\left( 1-\frac{T-T_m}{T_m} +\frac{p-p_m}{p_m} \right) &\quad\text{equation of state}\tag{14} \end{align}

or, if we write them in terms of the total density and pressure fields:

\begin{align} \partial_t\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} = -\boldsymbol{\nabla}\left(\frac{p}{\rho_m}\right) + \nu \nabla^2\boldsymbol{u} - g\frac{\rho}{\rho_m}\widehat{\boldsymbol{z}} &\quad\text{N-S momentum equation (Newton's second law)} \tag{15}\\ \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0&\quad\text{conservation of mass}\tag{16}\\ \partial_t \rho + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \rho = \kappa\nabla^2 \rho &\quad\text{thermodynamic equation}\tag{17}\\ \rho = \rho_m\left( 1-\frac{T-T_m}{T_m} +\frac{p-p_m}{p_m} \right) &\quad\text{equation of state}\tag{18} \end{align}

Note #1

We will denote $p/\rho_m$ simply as $p$ (or similarly $p'/\rho_m$ simply as $p'$). This is often done in the literature. This way "pressure" has dimensions of $(\textrm{length})^2(\textrm{time})^{-2}$.

Note #2

You might have seen a different version of the Boussinesq equations. It is customary to write down the buoyancy term in a different way, usually depending ones background. For example, \begin{align} \textrm{astrophysics} & & -g \rho/\rho_0 \\ \textrm{engineering} & & -g \alpha T \\ \textrm{oceanography} & & b \end{align}

All definitions are equivalent. You can read off the conversions between the different definitions from the table. In each field they write the thermodynamic equation in the corresponding variable ($\rho$, $T$ or $b$). Once more we stress that in all these definitions what comes into play in the dynamical equations are the deviations from the hydrostatic balance.

Equilibrium solutions of the inviscid Boussinesq equations

Let's start with the inviscid Boussinesq equations (drop viscosity $\nu$ and diffusivity $\kappa$):

\begin{gather} \partial_t\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} = -\boldsymbol{\nabla} p - g\frac{\rho}{\rho_m}\widehat{\boldsymbol{z}} \tag{19}\\ \partial_t \rho + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \rho = 0 \tag{20}\\ \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0 \tag{21} \end{gather}

It's easy to find non-trivial equilibrium solutions of the above. In particular, any flow fields of the form:

$$ \boldsymbol{u} = U(z)\,\widehat{\boldsymbol{x}},\quad\rho = \rho(z), \tag{22} $$

accompanied with an appropriate pressure field.

Exercise: Confirm that indeed (22) satisfies (19)-(21). What's the pressure field needed so that (22) satisfies (19)-(21)?

Exercise: Can you generalize (22) so that it still remains a solution of (19)-(21)? For example, would any flow $\boldsymbol{u} = U(x, y, z)\,\widehat{\boldsymbol{x}}$ also work? How about flow with more than one non-zero components?

The fact that any basic state in the form of (22) is an equilibrium solution of the inviscid Boussinesq equations has triggered a lot of instability studies. Some of the simplest example of basic states in the form of (22) are shown below:

Each of these basic states are so well studied that they have a name; usually the name of the person who studied them first.

In the next lecture we will study the stability of the unstratified Kelvin-Helmholtz profile. We will study the stability of this profile both analytically and numerically. The rest of the basic states will be some of the proposed project.

There are theorems that help you identify when a profile like those above is a good candidate for instability (Rayleigh's and Fjortoft's criterion or Miles & Howard's criterion). These criteria would be the subject of another project.


  • Salmon, R., 1998. Lectures on geophysical fluid dynamics. Oxford University Press.
  • Spiegel, E. A., & Veronis, G. (1960). On the Boussinesq approximation for a compressible fluid. Astrophys. J., 131, 442-447.
  • Vallis, G. K. (2017). Atmospheric and oceanic fluid dynamics. Cambridge University Press.