## pyqg: Python Quasigeostrophic Model 1
2
3
4
5
6

### 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)
l = fftshift(m.l, **kwargs) * m.radii


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!

1
2
3
4
5
6

This site uses Akismet to reduce spam. Learn how your comment data is processed. From Wikipedia.