That's Geophysical Fluid Dynamics.

Surface Quasi-Geostrophic (SQG) Model

pyqg will be used to reproduce the results of this seminal paper by I. M. Held, R. T. Pierrehumbert, S. T. Garner and K. L. Swanson (1995).

SQG describes surface-intensified flows due to buoyancy. Assume we’re on a semi-infinite domain from \(-\infty\) to 0 in \(z\). We have buoyancy conservation at \(z=0\):

$$
\frac{\partial b}{\partial t} + J(\psi, b) = 0
$$

And \(b = f_0 \frac{\partial \psi}{\partial z}\) there. The potential vorticity in the interior is zero. If we assume we’re on an f-plane and in a Boussinesq framework, PV is given by:

$$q = \nabla_h^2 \psi + \frac{\partial}{\partial z}\left(\frac{f_0^2}{N^2}\frac{\partial \psi}{\partial z}\right) = 0$$

The domain is periodic in x and y, so we can use replace differentiation with multiplication by wavenumbers in Fourier space. Likewise, we will assume that the streamfunction decays as \(z \to -\infty\), so the solution will have an exponential factor decaying in \(z\). We can assume plane wave solutions of the form:

$$\psi = \sum\limits_{k,l}\Re\left\{\hat{\psi}\exp(i[kx + ly])\right\}\exp\left(\frac{N}{f_0}Kz\right)$$

where \(K = \sqrt{k^2+l^2}\).

Note that this satisfies the zero PV condition in the interior:

$$\hat{q}_{k,l} = \left((ik)^2+(il)^2 + \frac{f_0^2}{N^2}\left(\frac{N}{f_0}\sqrt{k^2+l^2}\right)^2\right)\hat{\psi}_{k,l} = 0$$

From this relation we get that:

$$\hat{b}_{k,l} = N K \hat{\psi}_{k,l}$$

Now it’s time to start coding! Let’s create the model object:

year = 1.
kwargs = {'L': 2.*np.pi, 'nx': 512, 'tmax': 26.005, 'beta': 0., 'Nb': 1., 'H': 1., 'rek': 0., 'rd': None,
'dt': 0.005, 'taveint': 1., 'twrite': 400., 'ntd': 4}
m = sqg_model.SQGModel(**kwargs)
# in this example we used ntd=4, four threads
# if your machine has more (or fewer) cores available, you could try changing it

The initial condition will be an elliptical vortex:

$$b = \exp\left( -(x^2+(4y)^2) \right)$$

The eccentricity of this ellipse is 4, and this results in an instability. If the eccentricity is lowered (say the initial condition is a circle rather than an ellipse), such an instability will not appear.

Let’s set the initial condition (note that qi represents the surface buoyancy here and NOT potential vorticity):

# Choose ICs from Held et al. (1995)
# case i) Elliptical vortex
x = np.linspace(m.dx/2,2*np.pi,m.nx) - np.pi
y = np.linspace(m.dy/2,2*np.pi,m.ny) - np.pi
x,y = np.meshgrid(x,y)

qi = np.exp(-(x**2 + (4.0*y)**2))

# initialize the model with that initial condition
m.set_q(qi[np.newaxis,:,:])

We will run the model and make a movie using the snapshots feature, like we did for the barotropic model (see the previous page):

extent = [-np.pi, np.pi, np.pi, -np.pi]
fig, ax = plt.subplots()
l = ax.imshow(m.q.squeeze(), cmap='RdBu_r', extent=extent)
l.set_clim([0, 1])
ax.set_aspect(1)
plt.colorbar(l)
plt.title('b(x,y,t={})'.format(m.t))

import matplotlib.animation as manimation
FFMpegWriter = manimation.writers['ffmpeg']
writer = FFMpegWriter(fps=5)
vid_title = 'sqg_movie'
with writer.saving(fig, "{}.mp4".format(vid_title), 100): 
    writer.grab_frame() # initial data 
    for _ in m.run_with_snapshots(tsnapstart=0, tsnapint=1): 
        l.set_data(m.q.squeeze()) 
        plt.title('b(x,y,t={})'.format(m.t))
        writer.grab_frame()

By the time our GIF ends, we can see the development of secondary instabilities within the filaments. This is part of what sets surface buoyancy evolution in SQG apart from the 2D vorticity equations. The SQG model can form singular vorticity gradients in finite time, whereas 2D vorticity cannot. For the 2D vorticity equation, the growth rates of filamentary instabilities remains constant due to vorticity conservation, and . In SQG, it is buoyancy that is conserved, not vorticity, so the growth rate of the instability may continue to increase until there is breakdown.

The last thing we’ll do is talk about the built-in linear stability analysis feature.