That's Geophysical Fluid Dynamics.

Built-in Linear Stability Analysis

There is a built-in tool for computing the baroclinic linear stability analysis from a model’s initial conditions. It looks for solutions for \(\vec{\Phi}\) and \(\omega\) of the following equation:

$$\mathbf{A}\vec{\Phi} = \omega \mathbf{B}\vec{\Phi}$$

where:

$$
\begin{align*}
\mathbf{B} &= \mathbf{S} – \kappa^2\mathbf{I} \\
\mathbf{A} &= (Uk + Vl)\mathbf{B} + \mathbf{I}(kQ_y – lQ_x) + i(\delta_{N N}r_{ek} \kappa^2)\mathbf{I}
\end{align*}
$$

Recall that \(\mathbf{S}\) is the “stretching matrix” in a layered model. In the 3-layer case:

$$\mathbf{S} = \begin{pmatrix}
-\frac{f_0^2}{g’_1 H_1} & \frac{f_0^2}{g’_1 H_1} & 0 \\
\frac{f_0^2}{g’_1 H_2} & – (\frac{f_0^2}{g’_1 H_2} + \frac{f_0^2}{g’_2 H_2}) & \frac{f_0^2}{g’_2 H_2} \\
0 & \frac{f_0^2}{g’_2 H_3} & -\frac{f_0^2}{g’_2 H_3}
\end{pmatrix}$$

Growth rates are given by \(\Im(\omega)\). We can use this to find the most unstable mode. We will consider an example from the 2-layer Phillips model of baroclinic instability, without bottom friction. The initial conditions are as follows:

kwargs = {'nx': 256, 'nz': 2, 'U': [0.01, -0.01], 'V': [0.,0.], 'H': [1.,1.], 'L': 2*np.pi, 'beta': 1.5, 'rd': 1./20., 'rek': 0.05, 'f': 1., 'delta': 1.}
m = pyqg.LayeredModel(**kwargs)

We call the stability_analysis function which returns the eigenvectors and eigenvalues from the linear algebra problem posed above:

evals, evecs = m.stability_analysis()

For plotting purposes, use fftshift to reorder the entries. The wavenumbers are also multiplied by the first baroclinic deformation radius to make them dimensionless:

from numpy.fft import fftshift
kwargs = {'axes':(0,)}
evals = fftshift(evals.imag, **kwargs)
k = m.k*m.radii[1]
l = fftshift(m.l, **kwargs) * m.radii[1]

We can analyze the fastest growing mode:

# we want the l=0 slice
kwargs = {'axes':(1,)}
l0ind = int(m.Ny/2)
argmax = evals[l0ind, :].argmax()
evecs_sh = fftshift(evecs, **kwargs)
evec = evecs_sh[:, l0ind, argmax]
kmax = k[l0ind, argmax]

x = np.linspace(0, 4 * np.pi / kmax, 100)
mag = np.abs(evec)
phase = np.arctan2(evec.imag, evec.real)

And now let’s plot:

plt.contour(k, l, evals, colors='k')
plt.pcolormesh(k, l, evals, cmap='Blues')
plt.colorbar()
plt.xlim(0,2.); plt.ylim(-2., 2.)
plt.clim([0., .1])
plt.xlabel(r'$k \, L_d$'); plt.ylabel(r'$l \, L_d$')
plt.title('without bottom friction')
plt.show()

Consider the \(l=0\) slice of the above plot. We clearly must have \(k < k_d\), but there is also a low wavenumber cutoff. This example model has \(U = 0.01\),  \(\beta = 1.5 \) and \(L_d = \frac{1}{20}\). Our cutoff is set by the Rhines scale:

$$k L_d >L_d \sqrt{\frac{\beta}{2U}} \approx 0.433$$

Which confirms what we see on the plot. If there were no \(\beta\) effect, there would be no low wavenumber cutoff:

That’s about everything I have to say about pyqg. You’ve got plenty of models to choose from and lots of neat experiments to come up with. Have fun!