The Eady model for baroclinic instability

RSES, ANU, 2018

Rotating reference frame

The Coriolis term and the "traditional approximation"

In a rotating frame of reference there is an extra term that appears in the momentum equations: the Coriolis term \begin{align*} 2\boldsymbol{\Omega}\boldsymbol{\times}\boldsymbol{u}. \tag{1} \end{align*} Consider the expression of the Coriolis force in cartesian coordinates where $x$ is along latitude circles, $y$ along longitudes, and $z$ is in the vertical. Then,

\begin{align*} \boldsymbol{\Omega}\boldsymbol{\times}\boldsymbol{u} &= \left|\begin{array}{ccc}\widehat{\boldsymbol{x}}&\widehat{\boldsymbol{y}}&\widehat{\boldsymbol{z}}\\0 &\Omega\cos{\theta}&\Omega\sin{\theta}\\ u & \upsilon& w\end{array}\right| \\&= \big( 2\Omega w \cos\theta - 2\Omega \upsilon\sin\theta\big)\widehat{\boldsymbol{x}} + 2\Omega u\sin\theta\,\widehat{\boldsymbol{y}} - 2\Omega u \cos\theta\,\widehat{\boldsymbol{z}}. \tag{2} \end{align*}

The Coriolis term and the "traditional approximation"

A usual approximation is the so called "traditional approximation", in which terms involving the horizontal component of the Earth’s angular velocity are neglected. With these approximations we have that

\begin{align*} 2\boldsymbol{\Omega}\boldsymbol{\times}\boldsymbol{u} \approx - 2\Omega \upsilon\sin\theta\,\widehat{\boldsymbol{x}} + 2\Omega u\sin\theta\,\widehat{\boldsymbol{y}}, \tag{3} \end{align*}

and the equations of motions simplify to:

\begin{align*} (\partial_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\,u - f \upsilon &= -\partial_x p\big/\rho_m, \tag{4} \\ (\partial_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\,\upsilon + f u &= -\partial_y p\big/\rho_m, \tag{5}\\ (\partial_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\,w &= -\partial_z p\big/\rho_m + b, \tag{6}\\ (\partial_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\,b &= 0, \tag{7} \end{align*}

where we have defined $f\equiv 2\Omega\sin\theta$, where $\theta$ is the latitude.

[There are some intermediate steps still to get the quasi-geostrophic (QG) equations...]

Quasi-geostrophic (QG) equations

In QG approximation the vertical velocity is zero (at leading order) and the horizonal velocities $u,\ \upsilon$ are given in terms of a streamfunction: $u=-\partial_y\psi$ and $\upsilon=\partial_x\psi$. The streamfunction is proportional to the pressure: $\psi=p\big/(\rho_m f_0)$. Furthermore, from hydrostatic balance, we have that buoyancy is also expressed in terms of the streamfunction $b=f_0\partial_z\psi$.

After some fiddling around we get that the equations of motion are:

\begin{align*} (\partial_t + u\partial_x + \upsilon\partial_y)\,q &= 0\quad\text{for }0<z/H<1, \tag{8}\\ (\partial_t + u\partial_x + \upsilon\partial_y)\,b &= 0\quad\text{@ }z=0, H, \tag{9} \end{align*}

where $q = f_0 + \partial_x\upsilon - \partial_y u + \partial_z\big(\frac{f_0}{N^2(z)}b\big)$ is the quasi-geostrophic potential vorticity (QGPV) and $b$ is the buoyancy. The QG approximation is very good for atmospheric flows in the mid-latitude. In that case, $z=0$ corresponds to the Earth's surface and $z=H$ corresponds to the tropopause, which is at about 10 km in the mid-latitudes.

In terms of the streamfunction the QGPV and buoyancy are:

\begin{align*} q = f_0 + \big(\partial_x^2+\partial_y^2 + \partial_z \frac{f_0^2}{N^2} \partial_z\big)\psi \quad\text{and}\quad b=f_0\partial_z\psi. \tag{10} \end{align*}

The domain-averaged energy of the flow is

\begin{align*} E &= \frac1{H L_x L_y}\int \frac1{2}\big( u^2 + \upsilon^2 + \frac{b^2}{N^2}\big)\,\mathrm{d}^3\boldsymbol{x}\\ &= \frac1{H L_x L_y}\int \frac1{2}\big[ (\partial_x\psi)^2 + (\partial_y\psi)^2 + \frac{f_0^2}{N^2} (\partial_z\psi)^2\big]\,\mathrm{d}^3\boldsymbol{x}. \tag{11} \end{align*}

Also, since everything is given in terms of $\psi$ it's also obvious that the following hold:

\begin{align*} f_0\partial_z u = -\partial_y b \quad\text{and}\quad f_0\partial_z \upsilon = \partial_x b. \tag{12} \end{align*}

These are called "thermal wind relations" and they relate the vertical gradients of the flow with horizontal buoyancy gradients.

Eady's problem

Let's consider a constant stratification, $N^2(z)=\textrm{const}$. (This is equivalent with $\mathrm{d}\rho_0/\mathrm{d}z=\textrm{const}$.) In this case we can see that the streamfunction

\begin{align*} \psi^e = -\Lambda y z,\tag{13} \end{align*}

satisfies the QG equations. This streamfunction gives:

\begin{align*} u^e = \Lambda z,\ \upsilon^e=0,\ q^e = 0,\ b^e=-f_0\Lambda y.\tag{14} \end{align*}

Eady's problem

We want to study how perturbations about solution (13) evolve. For simplicity, we will take perturbations that are independent of $y$, i.e., $\psi=\psi^e+\psi'(x,z,t)$. The linearized QG equations about the basic state become:

\begin{align*} \big(\partial_t + u^e_+\partial_x \big)b'_+ + \upsilon'_+\partial_yb^e_+ &= 0\quad\text{@ }z=H,\tag{15}\\ \big(\partial_t + u^e\partial_x \big)q' &= 0\quad\text{for }0\le z/H\le 1,\tag{16}\\ \big(\partial_t + u^e_-\partial_x \big)b'_- + \upsilon'_-\partial_yb^e_- &= 0\quad\text{@ }z=0,\tag{17} \end{align*}

with subscripts $-/+$ denoting the values of the fiends at $z=0$ and $z=H$ respectively. Expressed in terms of streamfuction we have:

\begin{align*} \big(\partial_t + \Lambda H\partial_x \big)\partial_z\psi'_+ - \Lambda\partial_x\psi'_+ &= 0\quad\text{@ }z=H,\tag{18}\\ \big(\partial_t + \Lambda z\partial_x \big)\big(\partial_x^2 + \frac{f_0^2}{N^2}\partial_z^2\big)\psi' &= 0\quad\text{for }0\le z/H\le 1,\tag{19}\\ \partial_t \partial_z\psi'_- - \Lambda\partial_x\psi'_- &= 0\quad\text{@ }z=0,\tag{20} \end{align*}

Non-dimensionalization

It's always very important to try to reduce the number of free parameters by non-dimensionalizing appropriately. Finding the best way to non-dimensionalize a system is a bit of an art. In this case, after fiddling around we get that if we use:

\begin{align*} z = H\widetilde{z},\ x = \frac{NH}{f_0}\widetilde{x},\ t=\frac{N}{\Lambda f_0}\widetilde{t}.\tag{21} \end{align*}

where tildes denote non-dimensional variables. Furthermore, we can scale the streamfunction as

\begin{align*} \psi' = \frac{NH^2\Lambda}{f_0} \widetilde{\psi}'.\tag{22} \end{align*}

However, because our equations are linear in perturbation streamfunction, how we scale $\psi'$ does not really matter...

Putting all these together, the linearized equations of motions become:

\begin{align*} \big(\partial_{\widetilde{t}} + \partial_{\widetilde{x}} \big)\partial_{\widetilde{z}} \widetilde{\psi}'_+ - \partial_{\widetilde{x}} \widetilde{\psi}'_+ &= 0\quad\text{@ }\widetilde{z}=1,\tag{23}\\ \big(\partial_{\widetilde{t}} + \widetilde{z} \partial_{\widetilde{x}} \big)\big(\partial_{\widetilde{x}}^2 + \partial_{\widetilde{z}}^2\big)\widetilde{\psi}' &= 0\quad\text{for }0\le \widetilde{z}\le 1,\tag{24}\\ \partial_{\widetilde{t}} \partial_{\widetilde{z}}\widetilde{\psi}'_- - \partial_{\widetilde{x}} \widetilde{\psi}'_- &= 0\quad\text{@ }\widetilde{z}=0.\tag{25} \end{align*}

Remarkably most free parameters are gone! Only the zonal wavenumber remains. Now we are ready to code-up these equations.

In non-dimensional variables the energy of the perturbations is

\begin{align*} E &= \frac1{H L_x}\int \frac1{2}\big[ (\partial_x\psi)^2 + \frac{f_0^2}{N^2} (\partial_z\psi)^2\big]\,\mathrm{d}^2\boldsymbol{x} \\ &= \frac{N H^3 \Lambda^2}{f_0L_x}\int \frac1{2}\big[ (\partial_{\widetilde{x}}\psi)^2 + (\partial_{\widetilde{z}}\psi)^2\big]\,\mathrm{d}\widetilde{x}\,\mathrm{d}\widetilde{z}. \tag{26} \end{align*}

From here after we will use the nondimensional equations but dropping the tildes.

Homogeneity in $x$

Given that the equations are homogeneous in $x$ we can take

\begin{align*} \psi'(x,z,t)= \widehat{\psi}(z,t)e^{\mathrm{i}kx}.\tag{27} \end{align*}

which implies:

\begin{align*} (\partial_t +\mathrm{i}k)\partial_z\widehat{\psi}_+ - \mathrm{i}k \widehat{\psi}_+ &= 0\quad\text{@ }z=1,\tag{28}\\ (\partial_t + \mathrm{i}k z )(\partial_z^2-k^2)\widehat{\psi} &= 0\quad\text{for }0\le z\le 1,\tag{29}\\ \partial_t \partial_z\widehat{\psi}_- - \mathrm{i}k\widehat{\psi}_- &= 0\quad\text{@ }z=0.\tag{30} \end{align*}

If we compute the energy averaged over a zonal wavelength (i.e., $L_x\to 2\pi/k$) and use our expansion (27) we get:

\begin{align*} E &=\frac{H^2 \Lambda^2}{4} \int \big( |\partial_z\widehat{\psi}|^2 + k^2 |\widehat{\psi} |^2 \big)\,\mathrm{d}\widetilde{z}.\tag{31} \end{align*}

Numerical implementation

Assume a discretization of our $z$-domain as:

\begin{align*} z \to \underbrace{z_1}_{=0}, z_2, \dots, z_{n_z-1}, \underbrace{z_{n_z}}_{=1}, \tag{32} \end{align*}

and similarly

\begin{align*} \widehat{\psi}(z, t) \to \underbrace{{\psi}_1}_{=\widehat{\psi}(0,t)}, {\psi}_2, \dots, {\psi}_{n_z-1}, \underbrace{\widehat{\psi}_{n_z}}_{=\widehat{\psi}(z_{n_z},t)}.\tag{33} \end{align*}

The QGPV equations become

\begin{align*} \partial_{t}\big(\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{\delta^2}-k^2\psi_j\big) &= -\mathrm{i} k z_j\big(\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{\delta^2}-k^2\psi_j\big)\quad\text{for }j=1,\dots,n_z.\tag{34} \end{align*}

But this means the QGPV equations at the boundaries require $\psi_{-1}$ and $\psi_{n_z+1}$ which are outside of our domain! However we can get around that using the extra two buoyancy equations at the boundaries. For example, the buoyancy equation at $z=0$ implies

\begin{align*} \partial_{t} \frac{\psi_2-\psi_{0}}{2\delta}& = \mathrm{i}k \psi_1.\tag{35} \end{align*}

We can eliminate $\psi_{0}$ by multiplying (35) with $2/\delta$ and add to it equation (34) for $j=1$. We then have:

\begin{align*} \partial_{t}\big(\frac{2\psi_{2}-2\psi_1}{\delta^2}-k^2\psi_1\big) &= \frac{2\mathrm{i}k}{\delta}\psi_1. \tag{36} \end{align*}

Similarly, we can eliminate $\psi_{n_z+1}$ by using the buoyancy equation at $z=1$. That gives:

\begin{align*} \partial_{t}\big(\frac{-2\psi_{n_z}+2\psi_{n_z-1}}{\delta^2}-k^2\psi_{n_z}\big) &= -\mathrm{i}k \big(\frac{-2\psi_{n_z}+2\psi_{n_z-1}}{\delta^2}-k^2\psi_{n_z}\big) - \frac{2\mathrm{i}k}{\delta}\psi_{n_z}. \tag{37} \end{align*}

Thus our equations are (36), (37) and

\begin{align*} \partial_{t}\big(\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{\delta^2}-k^2\psi_j\big) &= -\mathrm{i} k z_j\big(\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{\delta^2}-k^2\psi_j\big)\quad\text{for }j=2,\dots,n_z-1. \tag{38} \end{align*}

We have $n_z+1$ equations for $n_z+1$ variables. (It always feels good when one realises that.)

Next we want to express the energy integral in discrete variables. For the potential energy term we have:

\begin{align*} \int |\partial_z\widehat{\psi}|^2 \,\mathrm{d}\widetilde{z} &\to \frac1{2} \big[\partial_z\widehat{\psi}^*\big]_1\big[\partial_z\widehat{\psi}\big]_1 \delta + \big[\partial_z\widehat{\psi}^*\big]_2\big[\partial_z\widehat{\psi}\big]_2 \delta + \dotsb +\big[\partial_z\widehat{\psi}^*\big]_{n_z-1}\big[\partial_z\widehat{\psi}\big]_{n_z-1} \delta + \frac1{2}\big[\partial_z\widehat{\psi}^*\big]_{n_z}\big[\partial_z\widehat{\psi}\big]_{n_z} \delta.\tag{39} \end{align*}

We'll approximate the derivatives at the boundaries as, e.g.,

\begin{align*} \big[\partial_z\widehat{\psi}^*\big]_1\big[\partial_z\widehat{\psi}\big]_1 = \frac{\psi_2^*-\psi_1^*}{\delta}\frac{\psi_2-\psi_1}{\delta}\quad\text{and}\quad\big[\partial_z\widehat{\psi}^*\big]_{n_z}\big[\partial_z\widehat{\psi}\big]_{n_z} = \frac{\psi^*_{n_z}-\psi^*_{n_z-1}}{\delta}\frac{\psi_{n_z}-\psi_{n_z-1}}{\delta},\tag{40} \end{align*}

while for the interior terms we will use the average of the forward and backwards finite difference, e.g.,

\begin{align*} \big[\partial_z\widehat{\psi}^*\big]_j\big[\partial_z\widehat{\psi}\big]_j = \frac1{2}\Big(\frac{\psi_{j+1}^*-\psi_j^*}{\delta}\frac{\psi_{j+1}-\psi_j}{\delta} + \frac{\psi_{j}^*-\psi_{j-1}^*}{\delta}\frac{\psi_{j}-\psi_{j-1}}{\delta}\Big).\tag{41} \end{align*}

Therefore we have that

\begin{align*} \frac{E}{2\pi/k} = \int |\partial_z\widehat{\psi}|^2 + k^2|\psi|^2 \,\mathrm{d}\widetilde{z} &\to \big[ \psi_1^*,\dots,\psi^*_{n_z}\big]\,\mathbb{M}\begin{bmatrix} \psi_1\\\vdots\\\psi_{n_z}\end{bmatrix},\tag{42} \end{align*}

with

\begin{align*} \mathbb{M} \equiv \frac1{\delta} \begin{pmatrix} 1+\frac1{2}k^2\delta^2 & -1 & 0 & \dotsb & 0 \\ -1 & 2+k^2\delta^2 & -1 & & \vdots \\ 0 & \ddots & \ddots & \ddots & 0\\ \vdots & & -1 & 2+k^2\delta^2 & -1\\ 0 & \dotsb & 0 & -1 & 1+\frac1{2}k^2\delta^2 \end{pmatrix}.\tag{43} \end{align*}

Note that $\mathbb{M}$ is self-adjoint.

Exercise: Would that be the case if we hadn't taken the average of forward and backwards finite difference approximation of the derivatives for interior points?

Exercise: Confirm (numerically) that $\mathbb{M}$ is positive definite.

Let's code all these up!

In [1]:
using LinearAlgebra, PyPlot
In [2]:
nz = 101
z = range(0, stop=1, length=nz)
dz = z[2]-z[1]

ue = Diagonal(z)
be = Diagonal([2/dz; zeros(nz-2); -2/dz]) # for j=1,nz equations
Inz = Matrix(Diagonal(ones(nz)))
Inz2 = Matrix(Diagonal([1/2; ones(nz-2); 1/2]))
diag, sidediag = -2ones(nz), ones(nz-1)
D2z = Matrix(Tridiagonal(sidediag, diag, sidediag))
D2z = D2z/dz^2

D2zmod = D2z
D2zmod[1, 1:2] = [-2, 2]/dz^2       # for j=1  equation
D2zmod[nz, nz-1:nz] = [2, -2]/dz^2  # for j=nz equation

function create_operators(k)
    Laplacian = D2zmod - k^2*Inz
    invLaplacian = inv(Laplacian)
    # the operator A
    A =invLaplacian*(-im*k*ue*Laplacian + im*k*be)
    # the metric M
    M = Matrix(Tridiagonal(-sidediag, -diag, -sidediag))    
    M[1, 1:2], M[nz, nz-1:nz] = [1, -1], [-1, 1]
    M = M/dz + k^2*dz*Inz2
    Msqrt = sqrt(M)
    invMsqrt = inv(Msqrt)
    # the operator A_M
    AM = Msqrt*A*invMsqrt
    return A, M, AM, Msqrt, invMsqrt
end;

First of all let's compute the eigenvalues of $\mathbb{A}$ for various $k$ and compare with the analytic expressions for the growth rates.

In [3]:
k = range(0, stop=3, length=51)

# the analytic growth rates
mu = k;
growthrate_analytic = real.( k ./ mu .* sqrt.(Complex.(1 ./ tanh.(mu/2)-mu/2).*(mu/2-tanh.(mu/2))));
growthrate_analytic[1] = 0

# the numerical growth rates
growthrate = zeros(length(k))
for j in 1:length(k)
    if k[j] == 0
        growthrate[j] = 0
    else
        A, M, AM = create_operators(k[j])
        growthrate[j] = maximum(real(eigvals(A)))        
    end
end
In [4]:
plot(k, growthrate_analytic, label="analytic")
plot(k, growthrate, "*", label="numerical")
xlabel(L"$k$", fontsize=18)
ylabel(L"$\mathrm{Re}(\sigma)$", fontsize=18)
title("Eady growth rates", fontsize=18)
legend();

Now let's find the eigenvectors of $\mathbb{A}_{\mathbb{M}}$, $\mathbb{A}_{\mathbb{M}}^\dagger$, $\mathbb{A}_{\mathbb{M}}+\mathbb{A}_{\mathbb{M}}^\dagger$, and $e^{\mathbb{A}^\dagger t_*}e^{\mathbb{A} t_*}$ for some finite time $t_*$. Remember that to plot we need to return back to streamfunction variables by acting with $\mathbb{M}^{-1/2}$.

In [5]:
function eigenanalysis(A)
    S, U = eigen(A);
    mostunstable = findall(real(S) .== maximum(real(S)))
    return S[mostunstable], U[:, mostunstable]
end

function find_phi0s(k, tstar)
    # finds the eigenvector of A, A^\dagger, A+A^\dagger,
    # and e^{A^dagger t}e^{A t} for t=tstar
    
    A, M, AM, Msqrt, invMsqrt = create_operators(k)
    maxgrowthrate = maximum(real(eigvals(A)))

    s1, u1 = eigenanalysis(AM);
    s1adj, v1 = eigenanalysis(Matrix(adjoint(AM)));
    s1inst, w1 = eigenanalysis(AM + Matrix(adjoint(AM)));
    s1opt, y1 = eigenanalysis(Matrix(adjoint(exp(AM*tstar))*exp(AM*tstar)))
    return s1, u1, s1adj, v1, s1inst, w1, s1opt, y1
end;
In [6]:
function plot_phi0s(;k=1.6, tstar=1)
    # plots the eigenvector of A, A^\dagger, A+A^\dagger,
    # and e^{A^dagger t}e^{A t} for t=tstar
    
    A, M, AM, Msqrt, invMsqrt = create_operators(k)
    s1, u1, s1adj, v1, s1inst, w1, s1opt, y1 = find_phi0s(k, tstar)
    
    nx = 100
    x = range(0, stop=2π, length=nx)
    eikx = Matrix(transpose(exp.(im*k*x)))
    
    # go back to 'normal' coordinates by multiplying with invMsqrt
    u1real = real.((invMsqrt*u1)*eikx)
    v1real = real.((invMsqrt*v1)*eikx)
    w1real = real.((invMsqrt*w1)*eikx)
    y1real = real.((invMsqrt*y1)*eikx)

    fig, axs = subplots(ncols=2, nrows=2, figsize=(6.5, 6.5))
    sca(axs[1]); contour(x, z, u1real)
    title(L"eigvector of $A$", fontsize=14)
    ylabel(L"$z$", fontsize=14)
    sca(axs[2]); contour(x, z, v1real)
    title(L"eigvector of $A^\dagger$", fontsize=14)
    xlabel(L"$x$", fontsize=14); ylabel(L"$z$", fontsize=14)
    sca(axs[3]); contour(x, z, w1real)
    title(L"eigvector of $A+A^\dagger$", fontsize=14)
    sca(axs[4]); contour(x, z, y1real)
    title("eigvector of \$e^{A^\\dagger t_*}e^{A t_*}\$, \$t_*= $tstar\$", fontsize=14)
    xlabel(L"$x$", fontsize=14)
end;
In [7]:
plot_phi0s(k=1.6, tstar=5);

Now let's see how the perturbation energy evolves if we use each of the previous eigenvectors as an initial condition.

In [8]:
function compute_energy(phi0, tfinal, AM)
    t = range(0, stop=tfinal, length=101)
    energy = zeros(length(t))
    dt = t[2]-t[1]
    expAMdt = exp(AM*dt)
    phit = phi0/sqrt(dot(phi0, phi0))
    for j in 1:length(t)
        energy[j] = dot(phit, phit)
        phit .= expAMdt*phit
    end
    return t, energy
end;
In [9]:
k=1.6
A, M, AM, Msqrt, invMsqrt = create_operators(k)

tstar, tfinal = 5, 10
s1, u1, s1adj, v1, s1inst, w1, s1opt, y1 = find_phi0s(k, tstar)

figure(figsize=(6, 2.8))
t, energy = compute_energy(u1, tfinal, AM)
semilogy(t, energy , label="eigenvector")
t, energy = compute_energy(v1, tfinal, AM)
semilogy(t, energy , label="adjoint")
t, energy = compute_energy(w1, tfinal, AM)
semilogy(t, energy , label="instantaneous")
t, energy = compute_energy(y1, tfinal, AM)
semilogy(t, energy , label="optimal \$t_*= $tstar\$")
semilogy(t, exp.(2*real(s1[1])*t), "--k")
xlabel(L"$t$", fontsize=16); ylabel(L"$E(t)$", fontsize=16)
legend();

We see that the eigenmode grows the least of all!

What's even more interesting is to repeat all these but for a zonal wavenumber $k$ for which there is not any unstable eigenmode of $\mathbb{A}$. For example, let's do it for $k=2.6$.

In [10]:
plot_phi0s(k=2.6, tstar=5);
In [11]:
k=2.6
A, M, AM, Msqrt, invMsqrt = create_operators(k)

tstar, tfinal = 5, 10
s1, u1, s1adj, v1, s1inst, w1, s1opt, y1 = find_phi0s(k, tstar)

figure(figsize=(6, 2.8))
t, energy = compute_energy(u1, tfinal, AM); semilogy(t, energy , label="eigenvector")
t, energy = compute_energy(v1, tfinal, AM); semilogy(t, energy , label="adjoint")
t, energy = compute_energy(w1, tfinal, AM); semilogy(t, energy , label="instantaneous")
t, energy = compute_energy(y1, tfinal, AM); semilogy(t, energy , label="optimal \$t_*=10\$")
xlabel(L"$t$", fontsize=16); ylabel(L"$E(t)$", fontsize=16)
legend();

What if we plot for larger times?

In [12]:
k=2.6
A, M, AM, Msqrt, invMsqrt = create_operators(k)

tstar, tfinal = 5, 100
s1, u1, s1adj, v1, s1inst, w1, s1opt, y1 = find_phi0s(k, tstar)

figure(figsize=(6, 2.8))
t, energy = compute_energy(u1, tfinal, AM); semilogy(t, energy , label="eigenvector")
t, energy = compute_energy(v1, tfinal, AM); semilogy(t, energy , label="adjoint")
t, energy = compute_energy(w1, tfinal, AM); semilogy(t, energy , label="instantaneous")
t, energy = compute_energy(y1, tfinal, AM); semilogy(t, energy , label="optimal \$t_*=10\$")
xlabel(L"$t$", fontsize=16); ylabel(L"$E(t)$", fontsize=16)
legend();

The adjoint eigenmode grows in energy up to $10^5$ times...

Exercise: Try looking for what happens for even longer times.

So it's apparent that for non-normal operators eigenanalysis does not tell us much for how the perturbations will actually evolve.

As a last example let's consider the following. From the eigenanalysis of the problem we see that there is a high-wavenumber cutoff of the Eady instability. That is for $k_1<k_{\mathrm{crit}}\approx2.4$ there is instability while for $k_2>k_{\mathrm{crit}}$ all eigenvalues of $\mathbb{A}$ have negative real part.

Let's pick two wavenumbers close to the $k_{\mathrm{crit}}$ and investigate what's the maximum perturbation growth that we can get.

In [13]:
tstar, tfinal = 10, 50
k1 = 2.38 # this is <k_crit
A1, M1, AM1, M1sqrt, invM1sqrt = create_operators(k1)
s1, u1, s1adj, v1, s1inst, w1, s1opt, y1 = find_phi0s(k1, tstar)
t1, energy1 = compute_energy(y1, tfinal, AM1)
k2 = 2.42 # this is >k_crit
A2, M2, AM2, M2sqrt, invM2sqrt = create_operators(k2)
s2, u2, s2adj, v2, s2inst, w2, s2opt, y2 = find_phi0s(k2, tstar);
t2, energy2 = compute_energy(y2, tfinal, AM2)
gr1, gr2 = round(real(s1[1]), digits=3), round(real(s2[1]), digits=3)

figure(figsize=(6, 2.8))
semilogy(t1, energy1, label="\$k_1= $k1 \$, \$\\sigma_1 = $gr1\$")
semilogy(t2, energy2, label="\$k_2= $k2 \$, \$\\sigma_2 = $gr2\$")
semilogy(t1, 2e3*exp.(2*real(s1[1])*t1), "--k", linewidth=0.8, label="\$\\propto e^{2\\mathrm{Re}(\\sigma_1)t}\$")
semilogy(t1, 1e3*exp.(2*real(s2[1])*t1), "-.k", linewidth=0.8, label="\$\\propto e^{2\\mathrm{Re}(\\sigma_2)t}\$")
xlabel(L"$t$", fontsize=16); ylabel(L"$E(t)/E(0)$", fontsize=16); legend();

We see that for both wavenumbers if we start with the eigenvector of $\mathbb{A}_{\mathbb{M}}^\dagger$ as initial condition the energy evolution is same for up to $t\approx 10$. Beyond that point the energy evolution goes as $e^{2\mathrm{Re}(\sigma_1)t}$ for $k_1<k_{\mathrm{crit}}$, while it oscillates about a constant value for $k_2$.

But regardless of whether there was an instability or not, both wavenumbers showed perturbation energy growth of about $10^3$ before they diverge.

This brings us back to what we were discussing in the first lectures: Eigenanalysis or modal stability only tells you what will happen at $t\to\infty$. Indeed, at $t\to\infty$ the energy for $k_1$ and $k_2$ is very different. But that's not the case at all for intermediate times (in this example for up to $t\approx 10$).

The following code evoles an initial condition and plots snapshots of the streamfunction and the energy evolution.

In [ ]:
k, tstar = 2.6, 5
A, M, AM, Msqrt, invMsqrt = create_operators(k)
s1, u1, s1adj, v1, s1inst, w1, s1opt, y1 = find_phi0s(k, tstar)

phi0 = y1
phit = phi0 / sqrt(dot(phi0, phi0))

tfinal = 10
t = range(0, stop=tfinal, length = 51)
dt = t[2]-t[1]
energy = zeros(length(t))
expAMdt = exp(AM*dt)

nx = 100
x = range(0, stop=2π, length=nx)
eikx = Matrix(transpose(exp.(im*k*x)))

fig, axs = subplots(ncols=2, nrows=1, figsize=(12, 2.8))

for j in 1:length(t)
    energy[j] = dot(phit, phit)
    phit = expAMdt*phit
    
    phitreal = real.((invMsqrt*phit)*eikx)
    
    sca(axs[1]); cla();
    contour(x, z, phitreal)
    xlabel(L"$x$", fontsize=14)
    ylabel(L"$z$", fontsize=14)
    sca(axs[2]); cla();
    semilogy(t[1:j], energy[1:j])
    xlabel(L"$t$", fontsize=14)
    ylabel(L"$E(t)$", fontsize=14)
    sleep(0.01)
    IJulia.clear_output(true)
    display(fig)
end

References

  • Vallis (2017). Atmospheric and oceanic fluid dynamics. Cambridge University Press.
  • Farrell (1984). Modal and non-modal baroclinic waves, J. Atmos. Sci., 41 (4), 668-673.