GitHub

MultiLayerQG

Basic Equations

This module solves the layered Boussinesq 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}\]

with

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

where

\[b_{j} = - g \frac{\delta \rho_j}{\rho_0}\]

is the Boussinesq buoyancy in each layer, with $\rho = \rho_0 + \delta \rho$ the total density, $\rho_0$ a constant reference density, and $|\delta \rho| \ll \rho_0$ the perturbation density.

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}\]

where

\[\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 ,\]

with

\[\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}\]

the background PV gradient components in each layer and with $\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)$ is the two-dimensional Jacobian. On the right hand side, $\mu$ is linear bottom drag, and $\nu$ is hyperviscosity of order $n_\nu$. Plain old viscosity corresponds to $n_\nu = 1$.

Implementation

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!.

source
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 .\]

source

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 .\]

source

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.

source

Diagnostics

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)
energies(prob)

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 .\]

source

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)
fluxes(prob)

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.\]

source

Examples