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, with the bottom-most layer is the $n$-th layer.
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 both the relative vorticity in each layer $\nabla^2 \psi_j$ and the vortex stretching terms, i.e.,
\[\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 the $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 where $\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 S and S⁻¹ respectively. Additionally, the background PV gradients $\partial_x Q$ and $\partial_y Q$ are also included in Params as Qx and Qy respectively.
One can get the $\widehat{\psi}_j$ from $\widehat{q}_j$ via streamfunctionfrompv!. The inverse, i.e., to obtain $\widehat{q}_j$ from $\widehat{\psi}_j$, is done via pvfromstreamfunction!.
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 operator $L$ includes only (hyper)-viscosity and is computed via hyperviscosity(params, grid).
The nonlinear term is computed via 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!.
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.
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.
GeophysicalFlows.MultiLayerQG.updatevars! — Function
updatevars!(vars, params, grid, sol)
updatevars!(prob)Update all problem variables using sol.
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 .\]
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. (For a single fluid layer, i.e., 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.\]
Examples
examples/multilayerqg_2layer.jl: Simulate the growth and equilibration of baroclinic eddy turbulence in the Phillips 2-layer model.