That's Geophysical Fluid Dynamics.

The Python Quasigeostrophic Model, or pyqg, is a Python library that allows one to model quasigeostrophic systems. The QG equations are an approximation to the full fluid equations of motion in the limit of strong rotation and stratification. pyqg enables one to solve the QG equations in one or multiple layers. They use a pseudo-spectral doubly-periodic model. The QG equations are transformed to Fourier space, where the potential vorticity Fourier coefficients are advanced via a third-order Adams Bashforth scheme. The streamfunction Fourier coefficients can be calculated from these. Then the model moves back to physical space, where the process repeats.

We will describe the model in more detail below.

N-Layer QG Model

From Vallis

In each layer of the QG model other than the bottom layer, the evolution equation is \(\frac{Dq_i}{Dt} = ssd\), where \(ssd\) stands for “small-scale dissipation”. In the bottom layer, there is an additional linear bottom drag term on the right hand side of the form \(-r_{ek}\nabla^2\psi_N\). The form of \(q\) in each layer is:

$$
\begin{align*}
q_1 &= \nabla^2 \psi_1 + \frac{f_0^2}{H_1}\left(\frac{\psi_2-\psi_1}{g’_1} \right) \\
q_k &= \nabla^2 \psi_k + \frac{f_0^2}{H_k}\left(\frac{\psi_{k-1} – \psi_k}{g’_{k-1}} – \frac{\psi_k – \psi_{k+1}}{g_k’} \right) \\
q_N &= \nabla^2 \psi_N + \frac{f_0^2}{H_N}\left(\frac{\psi_{N-1} – \psi_{N}}{g_{N-1}’}\right) + \frac{f_0}{H_N}h_b(x,y)
\end{align*}
$$

where \(g_k’ = g\frac{\rho_{k+1} – \rho_{k}}{\rho_k}\). Assume we have a background velocity (\(U_k,V_k\)) and a background PV \(Q_k\) in each layer. The total velocities and PV in each layer are then:

$$
\begin{align*}
u_k^{tot} &= U_n – \psi_{ky} \\
v_k^{tot} &= V_n + \psi_{kx} \\
q_k^{tot} &= Q_k + \delta_{kN} \frac{f_0}{H_N}h_b + q_k
\end{align*}
$$

You can substitute this into the evolution equations to get a more complete expression.

Spectral Space

In spectral space, operators such as the Laplacian reduce to products of wavenumbers, allowing for simpler computations. Nonlinear products are computed and then converted to spectral space, making this a pseudospectral solution. One can arrive at a relationship between the Fourier coefficients for the streamfunction and potential vorticity:

$$\hat{q}_i = (\mathbf{S}-\kappa^2\mathbf{I})\hat{\psi}_i$$

where \(\kappa^2 = k^2+l^2\) and \(\mathbf{S}\) is copied below from the pyqg documentation page:

Energy and Enstrophy

The energy spectrum looks like:

$$E(k,l) = \frac{1}{2H}\sum\limits_{i=1}^N H_i \kappa^2 |\hat{\psi}_i|^2 + \frac{1}{2H}\sum\limits_{i=1}^{N-1} \frac{f_0^2}{g’_i} | \hat{\psi}_i – \hat{\psi}_{i+1} |^2$$

In the special case of the two-layer model, the energy in physical space is:

$$E = \frac{1}{HS}\int\frac{1}{2}H_1\|\nabla\psi_1\|^2 + \frac{1}{2}H_2\|\nabla\psi_2\|^2 dS$$

The enstrophy spectrum (N-layers) is:

$$
Z(k,l) = \frac{1}{2H}\sum\limits_{i=1}^N H_i | \hat{q}_i |^2
$$

And in the two-layer case, the enstrophy in physical space is:

$$Z = \frac{1}{HS}\int \frac{1}{2}H_1 q_1^2 + \frac{1}{2}H_2 q_2^2 dS$$

The eddy turnover time can then be estimated as \( \frac{2\pi}{\sqrt{Z}}\).

Surface Quasigeostrophy

Another use of pyqg is to model SQG, in which the evolution of buoyancy at the surface induce a fully three-dimensional flow. Buoyancy conservation describes the flow at \(z=0\):

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

The interior is continuously stratified, so the PV is:

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

The potential vorticity and buoyancy are related by \(b = f_0 \frac{\partial \psi}{\partial z}\) at the surface, and \(\psi\) goes to zero as \(z \to -\infty\).

Let’s start looking at some code examples to demonstrate the power of the model.