Basic Equations

This module solves the non-dimensional surface quasi-geostrophic (SQG) equation for surface buoyancy $b_s = b(x, y, z=0)$, as described by Capet et al. (2008). The buoyancy and the fluid velocity at the surface are related through a streamfunction $\psi$ via:

\[(u_s, v_s, b_s) = (-\partial_y \psi, \partial_x \psi, -\partial_z \psi) .\]

The SQG model evolves the surface buoyancy,

\[\partial_t b_s + \mathsf{J}(\psi, b_s) = \underbrace{-\nu(-1)^{n_\nu} \nabla^{2n_\nu} b_s}_{\textrm{buoyancy diffusion}} + \underbrace{F}_{\textrm{forcing}} .\]

Above, $\mathsf{J}(\psi, b) = (\partial_x \psi)(\partial_y b) - (\partial_y \psi)(\partial_x b)$ is the two-dimensional Jacobian. The evolution of buoyancy is only solved for the surface layer, but $b_s$ is a function of the vertical gradient of $\psi$. In the SQG system, the potential vorticity in the interior of the flow is identically zero. That is, relative vorticity is equal and opposite to the vertical stretching of the buoyancy layers,

\[\underbrace{\left(\partial_x^2 + \partial_y^2 \right) \psi}_{\textrm{relative vorticity}} + \underbrace{\partial_z^2 \psi}_{\textrm{stretching term}} = 0 ,\]

with the boundary conditions $b_s = - \partial_z \psi|_{z=0}$ and $\psi \rightarrow 0$ as $z \rightarrow -\infty$. (We take here the oceanographic convention: $z \le 0$.)

These equations describe a system where the streamfunction (and hence the dynamics) at all depths is prescribed entirely by the surface buoyancy. By taking the Fourier transform in the horizontal ($x$ and $y$), the streamfunction-buoyancy relation is:

\[\widehat{\psi}(k_x, k_y, z, t) = - \frac{\widehat{b_s}}{|𝐤|} \, e^{|𝐤|z} , \]

where $|𝐤| = \sqrt{k_x^2 + k_y^2}$ is the total horizontal wavenumber.


The buoyancy equation is time-stepped forward in Fourier space:

\[\partial_t \widehat{b_s} = - \widehat{\mathsf{J}(\psi, b_s)} - \nu |𝐤|^{2 n_\nu} \widehat{b_s} + \widehat{F} .\]

The surface buoyancy is b. The state variable sol is the Fourier transform of the surface buoyancy, bh.

The Jacobian is computed in the conservative form: $\mathsf{J}(f, g) = \partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]$.

The linear operator is constructed in Equation

Equation(params, grid)

Return the equation for surface QG dynamics with params and grid. The linear opeartor $L$ includes (hyper)-viscosity of order $n_ν$ with coefficient $ν$,

\[L = - ν |𝐤|^{2 n_ν} .\]

Plain old viscocity corresponds to $n_ν=1$.

The nonlinear term is computed via function calcN!().


while the nonlinear terms via

calcN!(N, sol, t, clock, vars, params, grid)

Calculate the nonlinear term, that is the advection term and the forcing,

\[N = - \widehat{𝖩(ψ, b)} + F̂ .\]


which in turn calls calcN_advection! and addforcing!.

Parameters and Variables

All required parameters are included inside Params and all module variables are included inside Vars.

For the decaying case (no forcing, $F = 0$), variables are constructed with Vars. For the forced case ($F \ne 0$) variables are constructed with either ForcedVars or StochasticForcedVars.

Helper functions


Some useful diagnostics are kinetic energy and buoyancy variance.


Return the domain-averaged surface kinetic energy. Since $u² + v² = |{\bf ∇} ψ|²$, we get

\[\int \frac1{2} |{\bf ∇} ψ|² \frac{𝖽x 𝖽y}{L_x L_y} = \sum_{𝐤} \frac1{2} |𝐤|² |ψ̂|² .\]

In SQG, this is identical to half the domain-averaged surface buoyancy variance.


Return the buoyancy variance,

\[\int b² \frac{𝖽x 𝖽y}{L_x L_y} = \sum_{𝐤} |b̂|² .\]

In SQG, this is identical to the velocity variance (i.e., twice the domain-averaged kinetic energy).


Other diagnostic include: buoyancy_dissipation and buoyancy_work.