1

2

3

4

5

6

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!

1

2

3

4

5

6

- pyqg: Python Quasigeostrophic Model
- Vorticity and Quasi-Geostrophy
- GFD = Rotation + Stratification
- My first blog post!

February 2018 (1)

August 2017 (1)

Modeling (1)

Pedagogical (2)