QG Equations

This page details the quasi-geostrophic (QG) equations implemented in QGYBJ+.jl.

Potential Vorticity Evolution

The core QG equation is the conservation of potential vorticity:

\[\frac{\partial q}{\partial t} + J(\psi, q) = \mathcal{D}_q + \mathcal{F}_q\]

where $\mathcal{D}_q$ is dissipation and $\mathcal{F}_q$ is forcing.

Potential Vorticity Definition

\[q = \underbrace{\nabla^2\psi}_{\text{relative vorticity}} + \underbrace{\frac{\partial}{\partial z}\left(\frac{f_0^2}{N^2}\frac{\partial\psi}{\partial z}\right)}_{\text{stretching term}}\]

For uniform stratification ($N^2 = \text{const}$):

\[q = \nabla^2\psi + \frac{f_0^2}{N^2}\frac{\partial^2\psi}{\partial z^2}\]

Physical Interpretation

Term	Physical Meaning
$\nabla^2\psi$	Relative vorticity $\zeta$
$(f_0^2/N^2)\partial_z^2\psi$	Vortex stretching due to vertical motion
$J(\psi, q)$	Advection of PV by geostrophic flow

Streamfunction Inversion

Given $q$, we solve for $\psi$ via the elliptic equation:

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

Spectral Representation

In spectral space (horizontal) with vertical finite differences:

\[-k_h^2 \hat{\psi} + \frac{\partial}{\partial z}\left(a(z)\frac{\partial\hat{\psi}}{\partial z}\right) = \hat{q}\]

where:

$k_h^2 = k_x^2 + k_y^2$ (horizontal wavenumber squared)
$a(z) = f_0^2/N^2(z)$ (stretching coefficient)

Tridiagonal System

For each horizontal wavenumber $(k_x, k_y)$, the vertical discretization gives:

\[a_k \hat{\psi}_{k-1} + b_k \hat{\psi}_k + c_k \hat{\psi}_{k+1} = \hat{q}_k\]

This is solved efficiently with the Thomas algorithm in O(nz) operations.

Boundary Conditions

Top and Bottom: $\frac{\partial\psi}{\partial z} = 0$ (no buoyancy flux)

In the code:

# Called for each time step
invert_q_to_psi!(state, grid, params, a_ell)

Velocity Fields

Geostrophic Velocities

From geostrophic balance:

\[u = -\frac{\partial\psi}{\partial y}, \quad v = \frac{\partial\psi}{\partial x}\]

In spectral space:

\[\hat{u} = -ik_y\hat{\psi}, \quad \hat{v} = ik_x\hat{\psi}\]

Vertical Velocity

The QG omega equation gives the ageostrophic vertical velocity:

\[N^2 \nabla^2 w + f_0^2\frac{\partial^2 w}{\partial z^2} = 2f_0 J\left(\frac{\partial\psi}{\partial z}, \nabla^2\psi\right)\]

or equivalently (dividing by $N^2$):

\[\nabla^2 w + \frac{f_0^2}{N^2}\frac{\partial^2 w}{\partial z^2} = \frac{2f_0}{N^2}J\left(\frac{\partial\psi}{\partial z}, \nabla^2\psi\right)\]

The coefficient $f_0^2/N^2 \ll 1$ in typical oceanic conditions, reflecting that stratification strongly suppresses vertical motion relative to horizontal. The RHS represents frontogenesis/frontolysis forcing through the interaction of vertical shear (thermal wind) with relative vorticity gradients.

Jacobian Operator

The Jacobian $J(a, b)$ is computed pseudo-spectrally:

\[J(a, b) = \frac{\partial a}{\partial x}\frac{\partial b}{\partial y} - \frac{\partial a}{\partial y}\frac{\partial b}{\partial x}\]

Algorithm

Compute $\partial a/\partial x$, $\partial a/\partial y$ in spectral space
Transform to physical space
Multiply in physical space
Transform back to spectral space
Apply dealiasing (2/3 rule)

Conservation Properties

The Jacobian satisfies:

$\int J(a, b) \, dA = 0$ (integral vanishes)
$J(a, a) = 0$ (anti-symmetry)

These ensure energy and enstrophy conservation in the inviscid limit.

Dissipation

Hyperdiffusion

The model uses scale-selective hyperdiffusion:

\[\mathcal{D}_q = -\nu_{h1}(-\nabla^2)^{p_1} q - \nu_{h2}(-\nabla^2)^{p_2} q - \nu_z\frac{\partial^2 q}{\partial z^2}\]

where:

$\nu_{h1}, p_1$: Large-scale dissipation (drag)
$\nu_{h2}, p_2$: Small-scale dissipation (hyperviscosity)
$\nu_z$: Vertical diffusion

Integrating Factor Method

To handle stiff diffusion, we use integrating factors:

\[\tilde{q} = q \cdot e^{\nu k^{2p} \Delta t}\]

This allows larger time steps while maintaining stability.

Wave Feedback

When waves are present, the QG equation includes a feedback term through a modified effective PV.

The Wave Feedback Mechanism

The wave-induced PV $q^w$ is computed from the wave envelope $B$:

\[q^w = \frac{i}{2f_0} J(B^*, B) + \frac{1}{4f_0} \nabla_h^2 |B|^2\]

where $B = B_R + i B_I$ is the complex wave envelope with units of velocity (m/s) and $f_0$ is the Coriolis parameter.

Dimensional Equations

The model solves dimensional equations where $B$ has actual velocity amplitude. The $1/f_0$ factor ensures correct dimensionality of the wave feedback term.

Effective PV for Inversion

The streamfunction is obtained by inverting the effective PV:

\[q^* = q - q^w\]

\[\nabla^2\psi + \frac{\partial}{\partial z}\left(\frac{f_0^2}{N^2}\frac{\partial\psi}{\partial z}\right) = q^*\]

This means the wave feedback modifies the inversion rather than appearing as an explicit advection term.

See Wave-Mean Interaction for detailed formulas and implementation.

Implementation

Key Functions

The QG equation implementation uses these core functions:

invert_q_to_psi! - Solve elliptic equation for streamfunction
jacobian_spectral! - Compute Jacobian pseudo-spectrally
compute_velocities! - Get (u, v) from streamfunction

See the Physics API Reference for detailed documentation.

Code Location

elliptic.jl: Streamfunction inversion (invert_q_to_psi!)
nonlinear.jl: Jacobian computation (jacobian_spectral!)
operators.jl: Velocity computation (compute_velocities!)

References

Vallis, G. K. (2017). Atmospheric and Oceanic Fluid Dynamics. Cambridge University Press.
Pedlosky, J. (1987). Geophysical Fluid Dynamics. Springer.