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 in 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 decaying case (no forcing, $F = 0$), vars can be constructed with DecayingVars. For the forced case ($F \ne 0$) the vars struct is with ForcedVars or StochasticForcedVars.

Helper functions



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.


  • examples/surfaceqg_decaying.jl: A script that simulates decaying surface quasi-geostrophic flow with a prescribed initial buoyancy field, producing an animation of the evolution of the surface buoyancy.

    Capet, X. et al., (2008). Surface kinetic energy transfer in surface quasi-geostrophic flows. J. Fluid Mech., 604, 165-174.