# 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!`

.

`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. (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.