MultiLayerQG

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

with

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

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

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_q! — Function

set_q!(sol, params, vars, grid, q)
set_q!(prob, q)

Set the solution prob.sol as the transform of q and update variables.

source

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

GeophysicalFlows.MultiLayerQG.updatevars! — Function

updatevars!(vars, params, grid, sol)
updatevars!(prob)

Update all problem variables using sol.

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

examples/multilayerqg_2layer.jl: Simulate the growth and equilibration of baroclinic eddy turbulence in the Phillips 2-layer model.