Basic Equations

This module solves the layered quasi-geostrophic equations on a beta plane of variable fluid depth $H - h(x, y)$. The flow in each layer is obtained through a streamfunction $\psi_j$ as $(u_j, v_j) = (-\partial_y \psi_j, \partial_x \psi_j)$, $j = 1, \dots, n$, where $n$ is the number of fluid layers.

The QGPV in each layer is

\[\mathrm{QGPV}_j = q_j + \underbrace{f_0 + \beta y}_{\textrm{planetary PV}} + \delta_{j, n} \underbrace{\frac{f_0 h}{H_n}}_{\textrm{topographic PV}}, \quad j = 1, \dots, n .\]

where $q_j$ incorporates the relative vorticity in each layer $\nabla^2 \psi_j$ and the vortex stretching terms:

\[\begin{aligned} q_1 &= \nabla^2 \psi_1 + F_{3/2, 1} (\psi_2 - \psi_1) ,\\ q_j &= \nabla^2 \psi_j + F_{j-1/2, j} (\psi_{j-1} - \psi_j) + F_{j+1/2, j} (\psi_{j+1} - \psi_j) , \quad j = 2, \dots, n-1 ,\\ q_n &= \nabla^2 \psi_n + F_{n-1/2, n} (\psi_{n-1} - \psi_n) . \end{aligned}\]


\[F_{j+1/2, k} = \frac{f_0^2}{g'_{j+1/2} H_k} \quad \text{and} \quad g'_{j+1/2} = g \frac{\rho_{j+1} - \rho_j}{\rho_{j+1}} .\]

In view of the relationships above, when we convert to Fourier space $q$'s and $\psi$'s are related via the matrix equation

\[\begin{pmatrix} \widehat{q}_{๐ค, 1}\\\vdots\\\widehat{q}_{๐ค, n} \end{pmatrix} = \underbrace{\left(-|๐ค|^2 \mathbb{1} + \mathbb{F} \right)}_{\equiv \mathbb{S}_{๐ค}} \begin{pmatrix} \widehat{\psi}_{๐ค, 1}\\\vdots\\\widehat{\psi}_{๐ค, n} \end{pmatrix}\]


\[\mathbb{F} \equiv \begin{pmatrix} -F_{3/2, 1} & F_{3/2, 1} & 0 & \cdots & 0\\ F_{3/2, 2} & -(F_{3/2, 2}+F_{5/2, 2}) & F_{5/2, 2} & & \vdots\\ 0 & \ddots & \ddots & \ddots & \\ \vdots & & & & 0 \\ 0 & \cdots & 0 & F_{n-1/2, n} & -F_{n-1/2, n} \end{pmatrix} .\]

Including an imposed zonal flow $U_j(y)$ in each layer, the equations of motion are:

\[\partial_t q_j + \mathsf{J}(\psi_j, q_j ) + (U_j - \partial_y\psi_j) \partial_x Q_j + U_j \partial_x q_j + (\partial_y Q_j)(\partial_x \psi_j) = -\delta_{j, n} \mu \nabla^2 \psi_n - \nu (-1)^{n_\nu} \nabla^{2 n_\nu} q_j ,\]


\[\begin{aligned} \partial_y Q_j &\equiv \beta - \partial_y^2 U_j - (1-\delta_{j,1}) F_{j-1/2, j} (U_{j-1} - U_j) - (1 - \delta_{j,n}) F_{j+1/2, j} (U_{j+1} - U_j) + \delta_{j,n} \partial_y \eta , \\ \partial_x Q_j &\equiv \delta_{j, n} \partial_x \eta . \end{aligned}\]


Matrices $\mathbb{S}_{๐ค}$ as well as $\mathbb{S}^{-1}_{๐ค}$ are included in params as params.S and params.Sโปยน respectively. Additionally, the background PV gradients $\partial_x Q$ and $\partial_y Q$ are also included in the params as params.Qx and params.Qy.

One can get the $\widehat{\psi}_j$ from $\widehat{q}_j$ via streamfunctionfrompv!(psih, qh, params, grid), while the inverse, i.e. obtain $\widehat{q}_j$ from $\widehat{\psi}_j$, is done via pvfromstreamfunction!(qh, psih, params, grid).

The equations of motion are time-stepped forward in Fourier space:

\[\partial_t \widehat{q}_j = - \widehat{\mathsf{J}(\psi_j, q_j)} - \widehat{U_j \partial_x Q_j} - \widehat{U_j \partial_x q_j} + \widehat{(\partial_y \psi_j) \partial_x Q_j} - \widehat{(\partial_x \psi_j)(\partial_y Q_j)} + \delta_{j, n} \mu |๐ค|^{2} \widehat{\psi}_n - \nu |๐ค|^{2n_\nu} \widehat{q}_j .\]

In doing so 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 state variable sol consists of the Fourier transforms of $q_j$ at each layer, i.e., qh.

The linear operator is constructed in Equation

GeophysicalFlows.MultiLayerQG.Equation โ€” Function
Equation(params, grid)

Return the equation for a multi-layer quasi-geostrophic problem with params and grid. The linear opeartor $L$ includes only (hyper)-viscosity and is computed via hyperviscosity(params, grid).

The nonlinear term is computed via function calcN!.

GeophysicalFlows.MultiLayerQG.hyperviscosity โ€” Function
hyperviscosity(params, grid)

Return the linear operator L that corresponds to (hyper)-viscosity of order $n_ฮฝ$ with coefficient $ฮฝ$ for $n$ fluid layers.

\[L_j = - ฮฝ |๐ค|^{2 n_ฮฝ}, \ j = 1, ...,n .\]


The nonlinear terms are computed via

GeophysicalFlows.MultiLayerQG.calcN! โ€” Function
calcN!(N, sol, t, clock, vars, params, grid)

Compute the nonlinear term, that is the advection term, the bottom drag, and the forcing:

\[N_j = - \widehat{๐–ฉ(ฯˆ_j, q_j)} - \widehat{U_j โˆ‚_x Q_j} - \widehat{U_j โˆ‚_x q_j} + \widehat{(โˆ‚_y ฯˆ_j)(โˆ‚_x Q_j)} - \widehat{(โˆ‚_x ฯˆ_j)(โˆ‚_y Q_j)} + ฮด_{j, n} ฮผ |๐ค|^2 ฯˆฬ‚_n + Fฬ‚_j .\]


which in turn calls calcN_advection! and addforcing!.

Linearized MultiLayerQG dynamics

The MultiLayerQG module includes also a linearized version of the dynamics about a base flow $U_j(y)$, $j = 1, \dots, n$; see LinearEquation, calcNlinear!, and calcN_linearadvection!.

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

Helper functions

GeophysicalFlows.MultiLayerQG.set_ฯˆ! โ€” Function
set_ฯˆ!(params, vars, grid, sol, ฯˆ)
set_ฯˆ!(prob, ฯˆ)

Set the solution prob.sol to the transform qh that corresponds to streamfunction ฯˆ and update variables.



The eddy kinetic energy in each layer and the eddy potential energy that corresponds to each fluid interface is computed via energies:

GeophysicalFlows.MultiLayerQG.energies โ€” Function
energies(vars, params, grid, sol)

Return the kinetic energy of each fluid layer KE$_1, ...,$ KE$_{n}$, and the potential energy of each fluid interface PE$_{3/2}, ...,$ PE$_{n-1/2}$, where $n$ is the number of layers in the fluid. (When $n=1$, only the kinetic energy is returned.)

The kinetic energy at the $j$-th fluid layer is

\[๐–ช๐–ค_j = \frac{H_j}{H} \int \frac1{2} |{\bf โˆ‡} ฯˆ_j|^2 \frac{๐–ฝx ๐–ฝy}{L_x L_y} = \frac1{2} \frac{H_j}{H} \sum_{๐ค} |๐ค|ยฒ |ฯˆฬ‚_j|ยฒ, \ j = 1, ..., n ,\]

while the potential energy that corresponds to the interface $j+1/2$ (i.e., the interface between the $j$-th and $(j+1)$-th fluid layer) is

\[๐–ฏ๐–ค_{j+1/2} = \int \frac1{2} \frac{fโ‚€^2}{g'_{j+1/2} H} (ฯˆ_j - ฯˆ_{j+1})^2 \frac{๐–ฝx ๐–ฝy}{L_x L_y} = \frac1{2} \frac{fโ‚€^2}{g'_{j+1/2} H} \sum_{๐ค} |ฯˆฬ‚_j - ฯˆฬ‚_{j+1}|ยฒ, \ j = 1, ..., n-1 .\]


The lateral eddy fluxes in each layer and the vertical fluxes across fluid interfaces are computed via fluxes:

GeophysicalFlows.MultiLayerQG.fluxes โ€” Function
fluxes(vars, params, grid, sol)

Return the lateral eddy fluxes within each fluid layer, lateralfluxes$_1,...,$lateralfluxes$_n$ and also the vertical eddy fluxes at each fluid interface, verticalfluxes$_{3/2},...,$verticalfluxes$_{n-1/2}$, where $n$ is the total number of layers in the fluid. (When $n=1$, only the lateral fluxes are returned.)

The lateral eddy fluxes within the $j$-th fluid layer are

\[\textrm{lateralfluxes}_j = \frac{H_j}{H} \int U_j v_j โˆ‚_y u_j \frac{๐–ฝx ๐–ฝy}{L_x L_y} , \ j = 1, ..., n ,\]

while the vertical eddy fluxes at the $j+1/2$-th fluid interface (i.e., interface between the $j$-th and $(j+1)$-th fluid layer) are

\[\textrm{verticalfluxes}_{j+1/2} = \int \frac{fโ‚€ยฒ}{g'_{j+1/2} H} (U_j - U_{j+1}) \, v_{j+1} ฯˆ_{j} \frac{๐–ฝx ๐–ฝy}{L_x L_y} , \ j = 1, ..., n-1.\]