\documentclass[12pt]{article}

\renewcommand\thesubsection{\alph{subsection}}

\usepackage{float}
\usepackage[english]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage[colorinlistoftodos]{todonotes}
\usepackage{fancyhdr}
\usepackage{enumerate}
\usepackage{amsmath}
\usepackage{bm}
\usepackage{cancel}
\usepackage[justification=centering]{caption}
\usepackage{parskip}
\usepackage{array}
\usepackage{tabu}
\usepackage{subcaption}
\usepackage[style=ieee]{biblatex}
\usepackage{setspace}
\usepackage[section]{placeins}
\pagestyle{fancy}
\fancyhf{}
\rhead{CEE 262B Supplementary Section}
\lhead{Section Notes}
\rfoot{\thepage}
\usepackage[nomarginpar, margin=1in, foot=.25in]{geometry}


\begin{document}

\section{Section Agenda}

\begin{enumerate}
    \item KH instability background
    \item Taylor-Goldstein equation derivation
    \item Stability analysis: problem statement and solution
\end{enumerate}


\section{Kelvin-Helmholtz Instability}

Hydrodynamic instabilities occur in a wide range of geophysical conditions, and are generally characterized by unrestrained growth and a transition to turbulence in the flow field. The Kelvin-Helmholtz (KH) instability specifically occurs in stratified shear flows. It can be seen in the atmosphere and ocean along density interfaces (pycnoclines). It often results in billow structures as seen in the figures below.

\begin{center}
\includegraphics[width=5 in, height=2.5 in]{kh1.jpg}
\end{center}

\begin{center}
\includegraphics[width=5 in, height=2.5 in]{kh2.jpg}
\end{center}

\section{Taylor-Goldstein equation derivation}

KH instabilities can be examined through the Taylor-Goldstein equation. We will derive this equation as follows:

Assume a two-dimmensional shear flow with continuous density stratification $\overline{\rho}(z)$ and buoyancy frequency $N^2 =  -\frac{g}{\rho_0} \frac{\partial \overline{\rho}}{\partial z}$. Ignoring viscous effects and invoking the Boussinesq approximation, our equations of motion are:

\begin{equation}\label{eq:umom}
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + w \frac{\partial u}{\partial z} = -\frac{1}{\rho_0} \frac{\partial p}{\partial x}
\end{equation}

\begin{equation}\label{eq:wmom}
\frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + w \frac{\partial w}{\partial z} = -\frac{1}{\rho_0} \frac{\partial p}{\partial z} - \frac{\rho g}{\rho_0}
\end{equation}

\begin{equation}\label{eq:cont}
\frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0
\end{equation}

\begin{equation}\label{eq:incomp}
\frac{\partial \rho}{\partial t} + u \frac{\partial \rho}{\partial x} + w\frac{\partial \rho}{\partial z} = 0 
\end{equation}

We can decompose the flow variables into mean and perturbation components. Assuming small mean vertical velocity:

\begin{equation}\label{eq:udec}
    u = \overline{u}(z) + u'(x,z,t)
\end{equation}

\begin{equation}\label{eq:wdec}
    w = w'(x,z,t)
\end{equation}

We can further decompose the pressure and density fields into background and perturbation components:

\begin{equation}\label{eq:pdec}
    p = p_0(z) + p'(x,z,t)
\end{equation}

\begin{equation}\label{eq:rhodec}
    \rho = \rho_0 + \overline{\rho}(z) + \rho'(x,z,t)
 \end{equation}

 Assuming that the background pressure is hydrostatic, i.e. 
 \begin{equation}\label{eq:hydrostatic}
     \frac{\partial p_0}{\partial z} = -(\rho_0 + \overline{\rho})g
 \end{equation}

Plugging equations \ref{eq:udec} - \ref{eq:hydrostatic} into equations \ref{eq:umom} - \ref{eq:incomp} and ignoring small terms, we obtain the following equations for small amplitude motions:

\begin{equation}\label{eq:newumom}
    \frac{\partial u'}{\partial t} + \overline{u} \frac{\partial u'}{\partial x} + w' \frac{d \overline{u}}{d z} = -\frac{1}{\rho_0}\frac{\partial p'}{\partial x}
\end{equation}

\begin{equation}\label{eq:newwmom}
    \frac{\partial w'}{\partial t} + \overline{u} \frac{\partial w'}{\partial x}  = -\frac{1}{\rho_0}\frac{\partial p'}{\partial z} - \frac{\rho' g}{\rho_0}
\end{equation}

\begin{equation}\label{eq:newcont}
    \frac{\partial u'}{\partial x} + \frac{\partial w'}{\partial z } = 0 
\end{equation}

\begin{equation}\label{eq:newdens}
    \frac{\partial \rho'}{\partial t} + \overline{u} \frac{\partial \rho'}{\partial x} + w' \frac{d \overline{\rho}}{d z} = 0
\end{equation}

Because this flow is 2D and divergence free, we can define a streamfunction $\psi$ such that:

\begin{equation}\label{eq:stream}
    u' = -\frac{\partial \psi}{\partial z}, \quad w' = \frac{\partial \psi}{\partial x}
\end{equation}

Using the definitions in equation \ref{eq:stream}, equation \ref{eq:newcont} is automatically satisfied. Plugging equation \ref{eq:stream} into equations \ref{eq:newumom} and \ref{eq:newwmom}, cross-differentiating, and then adding together to eliminate the pressure term, we obtain:

\begin{equation}\label{eq:streamlaplace}
    \frac{\partial}{\partial t} \nabla^2 \psi + \overline{u} \frac{\partial}{\partial x} \nabla^2 \psi - \frac{d^2 \overline{u}}{d z^2} \frac{\partial \psi}{\partial x} = -\frac{g}{\rho_0} \frac{\partial \rho'}{\partial x}
\end{equation}

The density equation becomes:

\begin{equation}\label{eq:density}
    \frac{\partial \rho'}{\partial t} + \overline{u} \frac{\partial \rho'}{\partial x} + \frac{\partial \psi}{\partial x} \frac{d \overline{\rho}}{d z} = 0
\end{equation}

Now assuming plane-wave solutions of the form:

\begin{equation}\label{eq:planewave}
\psi = \psi_0 (z) e^{i(kx - \omega t)}, \quad \rho' = R_0(z) e^{i(kx-\omega t)}
\end{equation}

We can plug equation \ref{eq:planewave} into equations \ref{eq:streamlaplace} and \ref{eq:density}, which reduces them to: 

\begin{equation}
(\overline{u} - c)\left(\frac{d^2 \psi_0}{d z^2} - k^2 \psi_0 \right) - \frac{d^2 \overline{u}}{d z^2} \psi_0 = -\frac{g}{\rho_0} R_0
\end{equation}

\begin{equation}
(\overline{u} - c)R_0 + \frac{d \overline{\rho}}{d z} \psi_0 = 0
\end{equation}

We can combine the previous two equations and eliminate $R_0$, resulting in the Taylor-Goldstein Equation: 

\begin{equation}\label{eq:TG}
(\overline{u} - c)\left(\frac{d^2 \psi_0}{d z^2} - k^2 \psi_0 \right) + \left( \frac{N^2}{\overline{u} - c} - \frac{d^2 \overline{u}}{d z^2}\right) \psi_0 = 0
\end{equation}

What you'll be doing on your problem set is a stability analysis, i.e. finding a condition for which the solution experiences unstable growth. This occurs when the wave speed has an imaginary part. Let's see why, allowing $c = c_r + i c_i$:

\begin{gather*}
\omega = c k \\
\omega  = (c_r + i c_i)k = \underbrace{c_rk}_{\omega_r} + i \underbrace{c_i k}_{\omega_i}
\end{gather*}

Next, assume we have a wavelike solution $u = u_0 e^{i(kx -\omega t)}$:

\begin{gather*}
u = u_0 e^{i(kx -\omega t)}\\
u = u_0 e^{i(kx - (\omega_r + i \omega_i) t)}\\
u = u_0 e^{(ikx - i \omega_r + \omega_i t)}
\end{gather*}

As we can see, there is a $u_0 e^{\omega_i t}$ term here due to the imaginary part of the wavespeed, which will cause exponential growth in time, i.e. an instability. 

\newpage

\section{Stability Analysis}

\subsection{Problem statement}
Using equation \ref{eq:TG} above and what you know about instabilities, show that instabilities may occur in a continuously stratified shear flow for Ri$_g = N^2/S^2  < \frac{1}{4}$, where $N^2 = -\frac{g}{\rho_0}\frac{\partial \overline{\rho}}{\partial z}$ and $S^2  = (\partial u/\partial z)^2$

Hint: introduce a new variable $\phi$:

\begin{equation}
\phi = \frac{\psi_0}{\sqrt{\overline{u}-c}}
\end{equation}

And write equation \ref{eq:TG} in terms of $\phi$. Then, multiply your new equation by the complex conjugate $\phi^*$ and integrate each term over the vertical. 

\subsection{Solution}

We can write equation \ref{eq:TG} in terms of $\phi$. This is slightly tricky, so let's go through it term by term. 

First we have $(\overline{u} - c) \frac{d^2 \psi_0}{d z^2}$. Plugging in our new variable yields: 

\begin{gather*}
(\overline{u}-c)\frac{d^2}{d z^2}(\phi (\overline{u} - c)^{1/2})
\end{gather*}

Taking this derivative, we obtain:

\begin{gather*}
(\overline{u}-c)\left[(\overline{u}-c)^{1/2} \frac{d^2 \phi}{d z^2} + \frac{d \phi}{d z}\frac{d \overline{u}}{d z}(\overline{u}-c)^{-1/2} + \phi\left[ \frac{1}{2} \frac{d^2 u}{d z^2} (\overline{u}-c)^{-1/2} - \frac{1}{4} \left(\frac{d u}{d z}\right)^2 (\overline{u}-c)^{-3/2}\right] \right]
\end{gather*}

Next we have $-(\overline{u}-c)k^2 \psi_0$. In terms of $\phi$:

\begin{gather*}
-(\overline{u}-c)^{1/2}k^2 \phi
\end{gather*}

Next we have $\frac{N^2 \psi_0}{(\overline{u}-c)}$. In terms of $\phi$:

\begin{gather*}
N^2 \phi(\overline{u}-c)^{-3/2}
\end{gather*}

Finally we have $-\frac{d^2 \overline{u}}{d z^2}\psi_0$. In terms of $\phi$:

\begin{gather*}
-\frac{d^2 \overline{u}}{d z^2}\phi (\overline{u}-c)^{-1/2}
\end{gather*}

Combining all of these terms and grouping by $\phi$, we obtain:

\begin{gather*}
\frac{d^2 \phi}{d z^2}(\overline{u}-c)^{3/2} + \frac{d \phi}{dz} \frac{d \overline{u}}{d z}(\overline{u}-c)^{1/2} + \dots\\
\phi\left[ \frac{1}{2} \frac{d^2 \overline{u}}{d z^2} (\overline{u}-c)^{1/2} - \frac{1}{4}\left(\frac{d^2 u}{d z^2}\right)(\overline{u}-c)^{-1/2} -(\overline{u}-c)^{3/2}k^2 + N^2 (\overline{u}-c)^{-1/2} - \frac{d^2 \overline{u}}{d z^2}(\overline{u}-c)^{1/2} \right] = 0
\end{gather*}

Next, we divide through by a factor of $(\overline{u}-c)^{1/2}$, and rearrange terms to obtain:

\begin{equation}
\frac{d}{dz}\left[ (\overline{u}-c) \frac{d \phi}{d z}\right] - \phi \left[(\overline{u}-c) k^2 + \frac{1}{2} \frac{d^2 \overline{u}}{d z^2} + \frac{\left(\frac{d \overline{u}}{dz}\right)^2 - 4N^2}{4 (\overline{u}-c)} \right] = 0
\end{equation}

Next, we multiply by $\phi^*$, the complex conjugate of $\phi$. Recalling that $\phi \phi^* = |\phi|^2$, we obtain:

\begin{equation}\label{eq:compconj}
-\frac{\partial }{\partial z}\left(\phi^* (\overline{u}-c) \frac{\partial \phi}{\partial z}\right) + (\overline{u}-c)\frac{\partial \phi}{\partial z}\frac{\partial \phi^*}{\partial z}+ (\overline{u}-c)k^2 |\phi|^2+ \frac{1}{2} \frac{d^2 \overline{u}}{d z^2} |\phi|^2 + \frac{1}{4} \frac{\left(\frac{d \overline{u}}{dz}\right)^2 - 4N^2}{ (\overline{u}-c)}|\phi|^2 = 0
\end{equation}

Next we can integrate over the vertical, from $z_1$ to $z_2$. Given that $\psi_0(z=z_1) = \psi_0(z=z_2) = 0$, this eliminates the first term of equation \ref{eq:compconj}: 
\begin{equation}
\int_{z_1}^{z_2}(\overline{u}-c)\left(\left|\frac{d \phi}{d z}\right|^2 + k^2 |\phi|^2 \right) dz + \frac{1}{2} \int_{z_1}^{z_2} \frac{d^2 \overline{u}}{d z^2} |\phi|^2 dz + \frac{1}{4} \int_{z_1}^{z_2} \frac{\left(\frac{d \overline{u}}{dz}\right)^2 - 4N^2}{ (\overline{u}-c)}|\phi|^2 dz = 0
\end{equation}

An instability can only occur if we have $c_i > 0$, so let's assume we have $c_i>0$ and look at the resulting imaginary part of the solution. The second integral is always real, leaving us with the first and third integrals:

\begin{equation}
\int_{z_1}^{z_2}\left(\left|\frac{d \phi}{d z}\right|^2 + k^2 |\phi|^2 \right) dz = \frac{1}{4} \int_{z_1}^{z_2} \frac{\left(\frac{d \overline{u}}{dz}\right)^2 - 4N^2}{ |\overline{u}-c|^2}|\phi|^2 dz 
\end{equation}

The left hand side is always positive, meaning the only way this relationship holds is if $\left(\frac{d \overline{u}}{dz}\right)^2 - 4 N^2 > 0 $ somewhere. Put another way, a \textit{necessary} condition for \textit{instability} is that $\left(\frac{d \overline{u}}{dz}\right)^2 - 4N^2 > 0$.  We can state this in terms of the gradient Richardson number:

\begin{gather*}
\left(\frac{d \overline{u}}{dz}\right)^2 - 4N^2 > 0\\
\left(\frac{d \overline{u}}{dz}\right)^2 > 4N^2 \\
\frac{N^2}{\left(\frac{d \overline{u}}{dz}\right)^2} < \frac{1}{4}\\
\textrm{Ri}_g < \frac{1}{4}
 \end{gather*}



\end{document}