QuasiGeostrophic Equations

Quasi-Geostrophic Equations Formulation

J. Emmanuel Johnson
CNRS
MEOM

In this section, we look at how we can solve the QG equations using elements from this package.

Generalized QG Equations

Note: This formational is largely based on [Amraoui et al., 2023]. See this paper for more information including a derivation about the energetics.

We can define two fields as via the potential vorticity (PV) and stream function (SF).

Stream Function:ψ=ψ(x⃗,t)x⃗∈Ω⊂RDs,t∈T⊂R+Potential Vorticity:q=q(x⃗,t)x⃗∈Ω⊂RDs,t∈T⊂R+\begin{aligned} \text{Stream Function}: && \psi & =\boldsymbol{\psi}(\vec{\mathbf{x}},t) && && \vec{\mathbf{x}}\in\Omega\sub\mathbb{R}^{D_s} , && t\in\mathcal{T}\sub\mathbb{R}^{+} \\ \text{Potential Vorticity}: && q& =\boldsymbol{q}(\vec{\mathbf{x}},t) && && \vec{\mathbf{x}}\in\Omega\sub\mathbb{R}^{D_s}, && t\in\mathcal{T}\sub\mathbb{R}^{+} \\ \end{aligned}
(1)

Let's assume we have a stack of said fields, i.e. qk,ψkq_k,\psi_k. So we have a state defined as:

Stream Function:ψ⃗=[ψ1,ψ2,…,ψK]⊤Potential Vorticity:q⃗=[q1,q2,…,qK]⊤\begin{aligned} \text{Stream Function}: && \vec{\boldsymbol{\psi}} = [\psi_1, \psi_2, \ldots, \psi_K]^\top \\ \text{Potential Vorticity}: && \vec{\boldsymbol{q}} = [q_1, q_2, \ldots, q_K]^\top \\ \end{aligned}
(2)

We can write the QG PDE to describe the spatiotemporal relationship between the PV and the SF which is defined as:

∂tqk+u⃗k⋅∇qk=Fk+Dk\begin{aligned} \partial_t q_k + \vec{\boldsymbol{u}}_k\cdot\boldsymbol{\nabla}q_k = F_k + D_k \end{aligned}
(3)

where at each layer kk, we have u⃗k\vec{\boldsymbol{u}}_k is the velocity vector, qkq_k is the potential vorticity, FkF_k are the forcing term(s), and DkD_k are dissipation terms. The forcing term can be made up of any external forces we deem necessary for the PDE. For example, we could have wind stress that affects the top layer or bottom friction that affects the bottom layer. The dissipation term represents all of the diffusion terms. For example, we could have some lateral friction or hyper-viscosity terms.

The advection term in equation (3) includes the velocity vector, u⃗k\vec{\boldsymbol{u}}_k, which is defined in terms of the stream function, ψ\psi. This is given by:

u⃗=[u,v]⊤=[−∂yψ,∂xψ]⊤\vec{\boldsymbol{u}} = [u, v]^\top = \left[ -\partial_y \psi, \partial_x \psi\right]^\top
(4)

The PV is the culmination of the potential energy contained within the dynamical forces, the thermal/stretching forces and the planetary forces.

Potential Vorticity Forces=Dynamical+Thermal+Planetary+…\text{Potential Vorticity Forces} = \text{Dynamical} + \text{Thermal} + \text{Planetary} + \ldots
(5)

Concretely, we can define this the PV and its relationship to the SF and other forces as:

qk=Δψk+(Mψ)k+fk\begin{aligned} q_k = \boldsymbol{\Delta}\psi_k + (\mathbf{M}\psi)_k + f_k \end{aligned}
(6)

where at each layer, kk, we have the dynamical (relative) vorticity, Δψk\boldsymbol{\Delta}\psi_k, the thermal/stretching vorticity, (Mψ)k(\mathbf{M}\psi)_k, and the planetary vorticity, fkf_k.

Idealized QG Model

We can use the above formulation to describe an idealized QG model. This idealized model will be fairly abstract with parameters that represent real things but do not actually correspond to anything tangible. However, these types of flows are very well-studied and are very useful for having a controlled setting which removes all uncertainties. We will use the notation used in [Frezat et al. (2022)Srinivasan et al. (2023)] with some explanations taken from [Gupta & Brandstetter, 2022].

The domain is simplified to be periodic and defined on a 2π2\pi domain, i.e. Lx=Ly=2πLx = L_y = 2\pi. We will also assume a periodic domain. This domain is topologically equivalent to a two-dimensional torus.

We will define a final form for the idealized QG model of the vorticity-stream functions on this idealized, periodic domain as

∂tω+det⁡J(ψ,ω)=−μω+ν∇2ω−β∂xψ+F\partial_t \omega + \det\boldsymbol{J}(\psi,\omega) = - \mu\omega + \nu\boldsymbol{\nabla}^2\omega - \beta\partial_x\psi + F
(7)

where ω\omega is the vorticity and ψ\psi is the stream function. The βψx=βv\beta\psi_x =\beta v-term is geophysical system approximation that captures the effect of the differential rotation. This force is experienced within the Earth system in a tangent plane approximation, i.e. β=∂yf\beta = \partial_y f. This β\beta-term is important and it allows the flow to manifest different turbulent regimes. For example, Rossby waves are common within the planetary systems and can appear when β=0\beta=0. The determinant Jacobian term encapsulates the advection term seen in (3). It is defined as:

Determinant Jacobian:det⁡J(ψ,ω)=∂xψ∂yω−∂yψ∂xω\begin{aligned} \text{Determinant Jacobian}: && && \det \boldsymbol{J}(\psi,\omega) &= \partial_x\psi\partial_y \omega - \partial_y\psi\partial_x \omega \end{aligned}
(8)

We see from that it is directly related to the advection term seen in (3) but written in a different way.

Determinant Jacobian

Like most things in physics, there are often many ways to express the same expression. Ultimately, they are all advection expressions. See the example on wikipedia for more details. The detemrinant Jacobian term can be written in many ways. Let's rewrite it as:

Advection=det⁡J(ψ,ω)\text{Advection} = \det \boldsymbol{J}(\psi,\omega)
(9)

We can expand the full expression which gives us

Advection=∂xψ∂yω−∂yψ∂xω\text{Advection} = \partial_x\psi\partial_y \omega - \partial_y\psi\partial_x \omega
(10)

We can plug int the velocity components of the stream function definition (4) into the above equation

Advection=v∂yω−u∂xω\text{Advection} = v\partial_y \omega - u\partial_x \omega
(11)

The partial derivative operator is commutable so we can take out the operator of both terms

Advection=∂y(vω)−∂x(uω)\text{Advection} = \partial_y (v\omega) - \partial_x (u\omega)
(12)

Alternatively, we can write the velocity and partial derivative operators in the vector format

Advection=u⃗⋅∇ω\text{Advection} = \vec{\boldsymbol{u}} \cdot \boldsymbol{\nabla} \omega
(13)

and we see that we arrive at formulation in (7). I personally prefer this way of writing it as it is more general. Furthermore, it exposes the many ways we can express this term like the determinant Jacobian or simple partial derivatives.

Note: It is important to take note of the many ways to express this as it can be useful for numerical schemes. For example, an upwind scheme might benefit from advection term where the velocity components are multiplied with the partial derivatives. Alternatively the determinant Jacobian on an Arakawa C-grid is a well known formulation for dealing with this.

The forcing term is typically chosen to be a constant wind forcing

Fω(x⃗)=kf[cos⁡(kfx)+cos⁡(kfy)]\boldsymbol{F}_\omega(\vec{\mathbf{x}}) = k_f \left[ \cos (k_f x) + \cos (k_f y)\right]
(14)
Relation to Navier-Stokes Equations

This derivation and explanation was largely taken from [Srinivasan et al., 2023] which has one of the best high-level explanation of the derivation without too much mathematical detail.

Taking a velocity field, u⃗=[u,v]⊤\vec{\boldsymbol{u}} = [u,v]^\top, we can write the non-dimensional form of the Navier-Stokes (N-S) equations as

∂tu⃗+u⃗⋅∇u⃗+f(k×u⃗)=−μu⃗+1Re∇2u⃗+Fu⃗(x⃗)\partial_t \vec{\boldsymbol{u}} + \vec{\boldsymbol{u}} \cdot \boldsymbol{\nabla}\vec{\boldsymbol{u}} + f (k \times \vec{\boldsymbol{u}}) = -\mu \vec{\boldsymbol{u}} + \frac{1}{Re}\boldsymbol{\nabla}^2\vec{\boldsymbol{u}} + \boldsymbol{F}_{\vec{\boldsymbol{u}}}(\vec{\mathbf{x}})
(15)

The, ff, is the Coriolis parameter which is the local rotation rate of the Earth and/or other planetary atmospheric forcings. The example in showcases an example with beta-plane forcing. k\boldsymbol{k} is the unit-vector normal to the (x,y)(x,y)-plane. The μ\mu is the linear drag coefficient which represents the bottom friction. The ReRe is the Reynolds number measuring the strength of the non-linear advection term, relative to the viscous term. In otherwords, the relationship give by:

u⃗⋅∇u⃗∝1Re∇2u⃗\vec{\boldsymbol{u}} \cdot \boldsymbol{\nabla}\vec{\boldsymbol{u}}\propto \frac{1}{Re}\boldsymbol{\nabla}^2\vec{\boldsymbol{u}}
(16)

This Reynolds number is indirectly proportional to the viscosity and proportional to the absolute velocity [Gupta & Brandstetter, 2022]

[Srinivasan et al., 2023] choose a Reynolds number of Re=2,500Re=2,500. The forcing is given by

Forcing:Fu⃗(x⃗)=[−sin⁡(kfy),sin⁡(kfx)]⊤\begin{aligned} \text{Forcing}: && \boldsymbol{F}_{\vec{\boldsymbol{u}}}(\vec{\mathbf{x}}) &= [-\sin(k_f y), \sin(k_f x)]^\top \end{aligned}
(17)

which is a sinusoidal time-invariant forcing field that continuously drives the flow. In [Guan et al., 2023], the wavenumber kf=4k_f=4, was chosen.

The velocity field in equation (15) is required to satisfy the mass conservation principal given by the continuity equation

Continuity Equation:∇⋅u⃗=∂yu+∂xv=0\begin{aligned} \text{Continuity Equation}: && \boldsymbol{\nabla}\cdot\vec{\boldsymbol{u}} &= \partial_yu + \partial_xv = 0 \end{aligned}
(18)

One can satisfy this by defining a stream function, ψ(x⃗)\psi(\vec{\mathbf{x}}), which is a scalar field that is defined as

u=−ψy,v=ψx\begin{aligned} u = - \psi_y, && && v =\psi_x \end{aligned}
(19)

The stream function and the continuity equations can be expressed as an evolution equation of a single scalar field, the vorticity. This scalar field is defined as the two-dimensional curl of the velocity field

ω=∇×u⃗=∂xv−∂yu=∇2ψ\omega = \boldsymbol{\nabla}\times \vec{\boldsymbol{u}}=\partial_x v - \partial_y u = \boldsymbol{\nabla}^2\psi
(20)

This equation captures the local rotation of a fluid parcel. The final result is the incompressible 2D Navier-Stokes equations in the scalar vorticity stream function form. In other words, the Quasi-Geostrophic equations.

∂tω+det⁡J(ψ,ω)−βv=−μω+1Re∇2ω+F\partial_t \omega + \det\boldsymbol{J}(\psi,\omega) - \beta v= - \mu\omega + \frac{1}{Re}\boldsymbol{\nabla}^2\omega + F
(21)

Note: There are some small differences between this equation and . The first is the coefficient in front of the diffusion term, ∇2ω\boldsymbol{\nabla}^2\omega. Here, we have the Reynolds number, 1/Re1/Re instead of the viscosity term, ν\nu, as shown in (7). In addition, we have the β\beta term. In this formulation, it is βv\beta v whereas in (7) it is expressed as β∂x\beta \partial_x. However, these are equivalent because the first component of the stream function velocities (19) is defined as v=∂xψv=\partial_x\psi. So we can plug this into the equation above.

NOTE: I am not sure about the sign issue of the β\beta-term in (21). I think it is a mistake and that it should be positive which would match the equation in (7) along with various other formulations [Frezat et al., 2022]

Case Studies

Flow Regimes

In [Frezat et al., 2022], they were looking at how the QG model could be help train surrogate models to fill in missing dynamics. They did the whole training regime online.

Parameter Details

Below are some experimental parameters found in which showcase 3 different flow regimes based on the parameter scheme.

Table 1:Table with idealized configuration

NameSymbolUnitsDecay FlowForced Flowβ\beta-Plane Flow
ResolutionNx×NyN_x\times N_y2,048×2,0482,048\times 2,0482,048×2,0482,048\times 2,0482,048×2,0482,048\times 2,048
DomainLx×LyL_x\times L_ykm10e3×10e310e3 \times 10e310e3×10e310e3 \times 10e310e3×10e310e3 \times 10e3
Time StepΔt\Delta_ts120120120120120120
Linear Drag Coefficientμ\mum−1^{-1}001.25e−81.25e-81.25e−81.25e-8
Viscosityν\num2^2s−1^{-1}67.067.022.022.022.022.0
Beta-Termβ\betam−1^{-1}s−1^{-1}0.00.00.00.01.14e−111.14e-11
Reynolds NumberReRe32e332e322e422e434e434e4

Toy QG

In this section, we will briefly outline an idealized QG model with terms that correspond to real quantities in nature. This was taken from the lab that's available online (PDF | Course | Notebook | Code). It features a great step-by-step introduction to a complete toy problem.

∂tq+u⃗⋅q=1ρH∇H×τ⃗−κ∇Hψ+a4ΔH2−β∂xψ\partial_t q + \vec{\boldsymbol{u}}\cdot q = \frac{1}{\rho H}\boldsymbol{\nabla}_H\times\vec{\boldsymbol{\tau}} - \kappa\boldsymbol{\nabla}_H\psi + a_4\boldsymbol{\Delta}_H^2 - \beta\partial_x\psi
(22)

where q=∇2ψq=\boldsymbol{\nabla}^2\psi is the PV. This PDE describes the transport of the PV with respect to the following

∂tq+Advection=Wind Stress+Bottom Friction+Dissipation+Planetary Forces\partial_t q + \text{Advection} = \text{Wind Stress} + \text{Bottom Friction} + \text{Dissipation} + \text{Planetary Forces}
(23)

which correspond to all of the terms in equation (22).

Horizontal Operators. Apart from the partial derivative operators, we have some horizontal operators which means that they only operate on the horizontal axis, i.e. (x,y)(x,y). ∇H\boldsymbol{\nabla}_H is the horizontal gradient operator defined as ∇H=[∂x,∂y]⊤\boldsymbol{\nabla}_H = [\partial_x, \partial_y]^\top. ΔH\boldsymbol{\Delta}_H is the horizontal Laplacian operator defined as ΔH=[∂x2,∂y2]⊤\boldsymbol{\Delta}_H = [\partial^2_x, \partial^2_y]^\top. ΔH2\boldsymbol{\Delta}_H^2 is the hyper-Laplacian operator which is normally applied sequentially, i.e. ΔH∘ΔH\boldsymbol{\Delta}_H\circ\boldsymbol{\Delta}_H.

Advection. We have the advection term defined as the dot product velocity vector (equation (4)) defined via the stream function. As mentioned above, there are many ways to write this advection term, for example the determinant Jacobian.

Wind Stress. We have defined a wind stress term which is the cross product of the horizontal gradient and the wind stress, τ⃗\vec{\boldsymbol{\tau}}. As shown above, this term is a 2D vector field for the x and y directions.

τ⃗:=τ⃗(x,y)=[τx(x,y),τy(x,y)]⊤\vec{\boldsymbol{\tau}} := \vec{\boldsymbol{\tau}}(x,y) = \left[ \tau_x(x,y), \tau_y(x,y)\right]^\top
(24)

An "ideal" example os to assume the ocean is a "box" that extends from the equator to a latitude of about 60∘^\circ. The winds are typically constant and easterly near the equation and turn westerly at middle latitudes. In ideal setting, this is prescribed as a constant cosine forcing τ⃗=τ⃗max(−cos⁡y,0)\vec{\boldsymbol{\tau}}=\vec\tau_{\text{max}}(-\cos y,0) which is similar to the idealized QG model listed above. This product is then scaled according to the inverse of the density of water, ρ\rho, and the height of the layer, HH.

Bottom Friction. The bottom friction parameter is a scaling constant that dissipates the energy. The constant, κ\kappa, is defined as inversely proportional to the product of the height, HH, and the square-root of the vertical Eddy viscosity constant and the Coriolis parameter. This constant is multiplied with the Laplacian of the stream function.

Dissipation. There is a dissipation term which features the hyper-viscosity term. The horizontal hyper-Laplacian of the stream function is multiplied by some constant νH\nu_H which represents the lateral Eddy viscosity that may occur in the system.

Planetary Forcing The planetary forcing within this toy model is given by the β\beta-plane approximation for the planetary vorticity. This is given by the β\beta term which is defined as

β≈∂yf\beta \approx \partial_y f
(25)

where ff is the planetary forcing function given by equation (28). The only term in this function that depends upon yy is the second term so we approximate this as β\beta.

Parameter Details

There are many constants that are displayed in equation (22) and they each mean something. Below are some tables which outline each of the constants and parameters seen in the equation. The subsequent tables showcase some derived quantities. The last table showcases the nondimensional versions which might be useful for the non-dimensional version of the PDE.

The first table showcases the constants along with the range of values that it can take for the toy model.

Table 2:Table with constants

NameSymbolUnitsValue
Earth's RadiusRRm6400e36400e3
Earth's Angular FrequencyΩ\Omegas−1^{-1}7.25e−57.25e-5
Depth of Active LayerHHm100→4,000100 \rightarrow 4,000
Length of Ocean BasinLLm1000e3→5000e31000e3 \rightarrow 5000e3
Density of Waterρ\rhokg m−3^{-3}1e31e3
Lateral Eddy ViscosityνH\nu_Hm2^2s−1^{-1}00 or 1e1→1e41e1 \rightarrow 1e4
Vertical Eddy ViscosityνZ\nu_Zm2^2s−1^{-1}1e−4→1e11e-4 \rightarrow 1e1
Maximum Wind Stressτmax\tau_{\text{max}}kg m−1^{-1}s−2^{-2}1e−2→11e-2\rightarrow 1
Mean/Mid Latitudeθ0\theta_00→π/30\rightarrow \pi/3

The second table lists some of the derived constants that we can calculate given the above constants.

Table 3:Table of derived quantities.

NameEquationUnitsValue
β\beta-Planeβ=2Ωcos⁡θ0/R\beta = 2\Omega\cos\theta_0/R1.1→2.3e111.1 \rightarrow 2.3e11
Coriolis Parameterf0=2Ωsin⁡θ0f_0 = 2\Omega\sin\theta_0s−1^{-1}0.0→1.3e−40.0 \rightarrow 1.3e-4
Velocity ScaleU0=τmax/(βρHL)U_0 = \tau_{\text{max}}/(\beta\rho H L)1e−5→1e−11e-5 \rightarrow 1e-1
Bottom Frictionκ=1H(νZf0)/2\kappa = \frac{1}{H\sqrt{(\nu_Z f_0)/2}}0.0→1e−50.0 \rightarrow 1e-5

The last table lists the non-dimensional derived quantities which are useful for the non-dimensional version of the toy QG PDE (22).

Table 4:Table of non-dimensional quantities.

NameEquationValue
ϵ/Vorticity Ratio=U0/(βL2)\epsilon/\text{Vorticity Ratio} = U_0/(\beta L^2)
τmax/(ϵβ2ρHL3)\tau_{\text{max}}/(\epsilon\beta^2\rho HL^3)1e−12→1e−141e-12 \rightarrow 1e-14
κ/(βL)\kappa/(\beta L)4e−4→6e14e-4 \rightarrow 6e1
νH/(βL3)\nu_H/(\beta L^3)1e−7→1e−41e-7 \rightarrow 1e-4

QG for Sea Surface Height

In this section, we will briefly outline how the QG equations can be used to represent the

There is a known relationship between sea surface height (SSH) and the stream function. This is given by

η=f0gψ\eta = \frac{f_0}{g}\psi
(26)

where f0f_0 and gg is the Coriolis and gravity constant respectively.

This relationship has been exploited in various works in the literature, especially within the data assimilation literature. For example [Guillou et al. (2021)Guillou et al. (2023)Amraoui et al. (2023)] used a 1 layer QG model to assimilate NADIR alongtrack SSH observations. A related example was in [Guillou et al., 2021] where they used a simple 1 layer QG model and a 1 layer linear SW model to model the balanced motions and the internal tides respectively. This was also used to assimilate NADIR and SWOT alongtrack SSH observations.

To write out the PDE, we will use the same formulation listed in equation (3). However, we will change this to reflect the relationship between SSH, the SF and now the PV.

q=(∇2−f02c2)ψ+fq = (\boldsymbol{\nabla}^2 - \frac{f_0^2}{c^2})\psi + f
(27)

where ff is the planetary forcing given by equation (28), f0f_0 is the Coriolis parameter and cc is the phase speed. This formulation is very similar to the PV listed in (6) except for the coefficient listed in front of the ψ\psi component. This can be summarized as the barotrophic deformation wavenumber

k=f0c=1LRk = \frac{f_0}{c}=\frac{1}{L_R}
(28)

where LRL_R is the Rossby radius of deformation. In an idealized setting, this is equal to LR=gH/f0L_R = \sqrt{gH}/f_0 It is a helpful parameter that can be used for the QG model of choice. For example, in [Guillou et al., 2021], they used a Rossby radius of 3030km for a 10∘×10∘10^\circ\times 10^\circ domain.

Case Studies

The most common case study I have seen for this would be applied to data assimilation. In particular, the QG model has been used to assimilate SSH satellite observations [Amraoui et al., 2023].

In [Amraoui et al., 2023], the authors did a free-run of the QG model for an assimilation task of 21 days. Their target area was over the North Atlantic (NA) domain. A slight change to the above formulation is that they used a forcing term as the constant wind forcing term outlined above. Note: They needed to do ∼\sim10,000 spinup steps to get a good running simulation.

Below is an outline of the configuration.

Parameter Details

Note: They used a LeapFrog scheme for the integrator.

Table 1:Fixed for their Eddy resolving simulations.

NameSymbolUnitsValue
Domain SizeLx×LyL_x\times L_ykm2_052 x 3_099
Resolutiondx×dydx\times dykm18 x 18
Grid SizeNx×NyN_x\times N_y113 x 170
Time StepΔt\Delta ts600
Phase Speedccms−2^{-2}[2.5, 1.0, 1.0]
Assimilation WindowTTdays21
Time Splitterμ\mu0.2

Stacked QG Model

In this section, we will look at the formulation from the stacked QG model. The majority of this section comes from papers by [Thiry et al. (2022)Thiry et al. (2023)] with some inspiration from the Q-GCM numerical model.

In this case, we are similar to the generalized QG model outlined in equation (3). We will change the notation slightly. Let's have a stacked set of PV and SF vectors denoted as q⃗\vec{\boldsymbol{q}} and ψ⃗\vec{\boldsymbol{\psi}} respectively. We consider the stream function and the potential vorticity to be NZN_Z stacked isopycnal layers. This is already outlined in equation (2) Now, we will write the stacked model in vector format as

∂tq⃗+u⃗⋅∇q⃗=BF+D\partial_t \vec{\boldsymbol{q}} + \vec{\boldsymbol{u}}\cdot\boldsymbol{\nabla} \vec{\boldsymbol{q}} = \mathbf{BF} + \mathbf{D}
(29)

where FF is a vector of forcing terms to be applied at each layer, kk, B\mathbf{B} is a matrix of interlayer height scaling factors and D\mathbf{D} is a vector of dissipation terms to be applied at each layer, kk. This is analogous to equation (3) but in vector format instead of a component-wise.

We can also define the PV equation in vectorized format (analogous to equation (6)) as

q⃗=(∇2−f02M)ψ+βy\vec{\boldsymbol{q}} = \left(\boldsymbol{\nabla}^2 - f_0^2\mathbf{M}\right)\psi + \beta y
(30)

Below, we will go over some key components of these two formulations.

Inter-Connected Layers

We define B\mathbf{B} as the matrix which connects the inter-layer interactions. This is a bi-diagonal matrix which only maps the corresponding layer kk and k−1k-1.

B=[1H1−1H10……01H1−1H2………………………−1HN−11HN−10………−1HN1HN]\mathbf{B} = \begin{bmatrix} \frac{1}{H_1} & \frac{-1}{H_1} & 0 & \ldots & \ldots \\ 0 & \frac{1}{H_1} & \frac{-1}{H_2} & \ldots & \ldots \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ \ldots & \ldots & \frac{-1}{H_{N-1} } & \frac{1}{H_{N-1}} & 0 \\ \ldots & \ldots& \ldots & \frac{-1}{H_N} & \frac{1}{H_N} \\ \end{bmatrix}
(31)

We also have M\mathbf{M} as the matrix connecting the layers of the stream function

M=[1H1g1′−1H1g2′………−1H2g1′1H1(1g1′+1g2′)−1H2g2′………………………−1Hn−1gn−2′1Hn−1(1gn−2′+1gn−1′)−1Hn−1gn−2′………−1Hngn−1′1Hngn−1′]\mathbf{M} = \begin{bmatrix} \frac{1}{H_1 g_1'} & \frac{-1}{H_1 g_2'} & \ldots & \ldots & \ldots \\ \frac{-1}{H_2 g_1'} & \frac{1}{H_1}\left(\frac{1}{g_1'} + \frac{1}{g_2'} \right) & \frac{-1}{H_2 g_2'} & \ldots & \ldots \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ \ldots & \ldots & \frac{-1}{H_{n-1} g_{n-2}'} & \frac{1}{H_{n-1}}\left(\frac{1}{g_{n-2}'} + \frac{1}{g_{n-1}'} \right) & \frac{-1}{H_{n-1} g_{n-2}'} \\ \ldots & \ldots& \ldots & \frac{-1}{H_n g_{n-1}'} & \frac{1}{H_n g_{n-1}'} \\ \end{bmatrix}
(32)

This matrix is a tri-diagonal matrix as each layer, kk, within the stacked layers will influence or be influenced by the neighbours, k−1,k+1k-1,k+1

Forcing

We have a forcing vector/matrix F\mathbf{F} which defines the forcing we apply at each layer. For example, the top layer could have wind forcing and/or atmospheric forcing. The bottom layer could have bottom friction due to the topography.

F=[F0,0,…,0,FN]⊤\mathbf{F} = \left[ F_0, 0, \ldots, 0, F_N\right]^\top
(33)

where we have the wind stress and the bottom friction. The wind forcing could be

F0=∂xτy−∂yτxτ⃗=τ0ρ0H1[−cos⁡(2πy/Ly),0]⊤\begin{aligned} F_0 =& \partial_x\tau_y - \partial_y\tau_x \\ \vec{\boldsymbol{\tau}} =& \frac{\tau_0}{\rho_0 H_1} \left[ -\cos(2\pi y/L_y), 0\right]^\top \end{aligned}
(34)

and the bottom friction could be described as

FN=δek2HN∇ψN\begin{aligned} F_N =& \frac{\delta_{ek}}{2H_N} \boldsymbol{\nabla}\psi_N \end{aligned}
(35)

where δek\delta_{ek} is the Ekman coefficient.

Dissipation

We have a very similar dissipation stategy as listed above in .... Here, we define this as

Viscosity:D1=a2Δ2ψHyperViscosity:D2=−a4Δ3ψ\begin{aligned} \text{Viscosity}: && && D_1 &= a_2 \boldsymbol{\Delta}^2\psi \\ \text{HyperViscosity}: && && D_2 &= - a_4 \boldsymbol{\Delta}^3\psi \end{aligned}
(36)

Case Studies

Q-GCM

There is an open-source QG GCM model (Q-GCM) that is available. It has a coupled model for the atmosphere and ocean where the atmospheric component is a stacked QG model and the oceanic component is a stacked QG model. They describe in detail two configurations for the double-gyre (North Atlantic) and the southern ocean. We outline in detail their

Parameter Details

Below is a table with the parameter configuration for their experiments.

Table 6:Variable parameters for their Northern and Southern Ocean experiments.

NameSymbolUnitsDouble GyreSouthern Ocean
Domain SizeLx×LyL_x \times L_ykm3_840 x 4_80023_040 x 2_880
Resolutiondx×dydx \times dykm10 x 1010 x 10
Grid SizeNx×NyN_x \times N_y384 x 4802_304 x 288
Time StepΔt\Delta tmin3010
Mean Layer ThicknessHkH_km[300, 1_100, 2_600][300, 1_100, 2_600]
Reduced Gravitygkg_kms−2^{-2}[0.05, 0.025][0.05, 0.025]
Bottom Ekman Layer Thicknessδek\delta_{ek}m12
Ocean Densityρ0\rho_0kgm−3^{-3}1_0001_000
Baroclinic Rossby RadiiLdL_dkm[51, 32][42, 26]
Laplacian Viscosity Coefficienta2a_2m2^2s−1^{-1}00
BiHarmonic Viscosity Coefficienta4a_4m4^4s1e103e10
Mean Coriolis Parameterf0f_0s−1^{-1}1e-4-1.1947e-4
Coriolis Parameter Gradientβ\betam−1^{-1}s−1^{-1}2e-41.313e-11

Dissipation Studies

In the paper [Thiry et al., 2023], they were investigating the impact of implicit dissipation via numerical methods or explicit dissipation via a hyper-viscosity parameter. For their experiments, they use the stacked QG model that was listed above to do a canonical double gyre experiment.

∂tq⃗+u⃗⋅∇Hq⃗=a2ΔH2ψ−a4ΔH3ψ+τ0ρ0H1[∂xτy−∂yτx,0⋯ ,0]−δek2HNZ[0,⋯ ,0,ΔψN]\partial_t \vec{\boldsymbol{q}} + \vec{\boldsymbol{u}}\cdot\boldsymbol{\nabla}_H \vec{\boldsymbol{q}} = a_2\boldsymbol{\Delta}_H^2\psi -a_4\boldsymbol{\Delta}_H^3\psi + \frac{\tau_0}{\rho_0H_1}\left[\partial_x\tau_y - \partial_y\tau_x, 0\cdots,0\right] - \frac{\delta_{ek}}{2H_{N_Z}} \left[0,\cdots,0,\Delta\psi_N\right]
(37)

Using this model, they modified the configuration of the spatial resolution and the hyper-viscosity which coincide with Eddy non-resolving, Eddy permitting, and Eddy resolving cases. They showcased how the numerical scheme they used was better than explicitly prescribing the dissipation.

Parameter Details

Below is a table with the parameter configuration for their experiments.

Table 7:Variable parameters for their Eddy resolving simulations.

NameSymbolUnitsNon-Eddy ResolvingEddy PermittingEddy Resolving
Grid SizeNx×NyN_x \times N_y129 x 129161 x 161193 x 193257 x 257385 x 385513 x 5131,025 x 1,025
Resolutiondx×dydx \times dykm40 x 4032 x 3226 x 2620 x 4013.3 x 13.310 x 105 x 5
Time StepΔt\Delta ts8_0006_0005_4004_0002_7002_0001_000
Munk Scaleδ\delta11111.251.52
Hyperviscositya4a_4m4^4s−1^{-1}1.8e125.9e112.4e115.6e102.6e101.3e101.7e9

Below is a table with the fixed parameters that stayed constant throughout the simulation.

Table 8:Fixed for their Eddy resolving simulations.

NameSymbolUnitsValue
Domain SizeLx×LyL_x\times L_ykm5120 x 5120
Mean Layer ThicknessHkH_km[400, 1_100, 2_600]
Reduced Gravitygkg_kms−2^{-2}[0.025, 0.0125]
Bottom Ekman Layer Thicknessδek\delta_{ek}m1
Wind Stress Magnitudeτ0\tau_0Nm−1^{-1} | Pa0.08
Ocean Densityρ0\rho_0kgm−3^{-3}1_000
Mean Coriolis Parameterf0f_0s−1^{-1}9.375e-5
Coriolis Parameter Gradientβ\betams−1{-1}1.754e-11
Baroclinic Rossby RadiiLdL_dkm[41, 25]

Forcing

They used a stationary symmetric wind stress forcing.

τx=−τ0ρ0cos⁡[2πyLy],τy=0\begin{aligned} \tau_x = - \frac{\tau_0}{\rho_0} \cos \left[ \frac{2\pi y}{L_y} \right], && && \tau_y = 0 \end{aligned}
(38)

Boundaries

  • For the velocity, u⃗\vec{\boldsymbol{u}}, they use free-slip boundaries, i.e., the tangential velocity for the North-South extent is zero.
  • For the relative vorticity, qq, they use zero boundaries.

Note: they used a 60 year spin-up period to get started.

Data Assimilation Benchmark

In [Laloyaux et al., 2020], they were exploring the effectiveness of a data assimilation method (4DVar) when applied to observation data.

They used a simple 2-Layer QG model with the stream function ψk\psi_k and the potential vorticity, qkq_k, as shown in equation (29). In the equation, they don't have any explicit forcing but they do mention some optional constant wind forcing.

They have a two layer system so their M-equation will be

M=[f02H1g−f02H1g−f02Hngf02Hng]\mathbf{M} = \begin{bmatrix} \frac{f_0^2}{H_1 g} & \frac{-f_0^2}{H_1 g} \\ \frac{-f_0^2}{H_n g} & \frac{f_0^2}{H_n g} \\ \end{bmatrix}
(39)

which is similar to equation (32) except they put the constant Coriolis parmater inside. In addition, they don't have any reduced gravities, just a constant gravity for each term.

This was designed for speed, stability, and convenience instead of accuracy and conservation.

Parameter Details

The authors have written the formulation completely out as.

q1=∇2ψ1−F1(ψ1−ψ2)+βyq2=∇2ψ2−F2(ψ2−ψ1)+βy+Rs\begin{aligned} q_1 &= \nabla^2\psi_1 - F_1(\psi_1 - \psi_2) + \beta y \\ q_2 &= \nabla^2\psi_2 - F_2(\psi_2 - \psi_1) + \beta y + R_s \end{aligned}
(40)

where β\beta is the (non-dimensionalized) northward derivative, ff is the Coriolis parameter, yy is the vertical coordinate and RsR_s represents the orography or heating. The parameters, F1,F2F_1,F_2, are the parameters that couple the laters together.

Below are some experimental parameters for the experimental setup

Table 8:Parameters

NameSymbolUnitsValue
Basin LengthLLm1e6
VelocityUUms−1^{-1}10
Coriolis Parameterf0f_0s−1^{-1}1e-4
Northward Derivativeβ0\beta_0m−1^{-1}s−1^{-1}1.5e-11
Layer DepthsD1,D2D_1, D_2m[6_000, 4_000]
Mean Potential Temperatureθˉ\bar{\theta}......
Layer Difference in Potential TemperatureΔθ\Delta\theta......
Ratio PTΔθθˉ\frac{\Delta\theta}{\bar{\theta}}...0.1
Mean Wind - Upper...ms−1^{-1}40
Mean Wind - Lower...ms−1^{-1}10

Note: the Coriolis paramater is at the southern boundary.

Table 8:Non-Dimensionalization

NameSymbolUnitsTransformation
Timettst~UˉL\tilde{t} \frac{\bar{U}}{L}
X-Coordinatexxmx~L\frac{\tilde{x}}{L}
Y-Coordinateyymy~L\frac{\tilde{y}}{L}
u-Velocityuums−1^{-1}u~U\frac{\tilde{u}}{U}
v-Velocityvvms−1^{-1}v~U\frac{\tilde{v}}{U}
Northward Derivativeβ\beta...β0L2U\beta_0\frac{L^2}{U}
Coupling Term IF1F_1...f02L2D1gΔθθˉ\frac{f_0^2L^2}{D_1 g \frac{\Delta\theta}{\bar{\theta}}}
Coupling Term IIF2F_2...f02L2D2gΔθθˉ\frac{f_0^2L^2}{D_2 g \frac{\Delta\theta}{\bar{\theta}}}
Rossby Numberϵ\epsilon...Uˉf0L=0.1\frac{\bar{U}}{f_0 L} = 0.1

Experimental Details

  • Time Stepping - first order upstream
  • PV Advection - Semi-Lagrangian advection
  • Interpolation of upstream PV - bi-cubic
  • Advection Outside Domain - edge values
  • North-South PV Values - user supplied
  • Advection Wind - Inverting PV
  • Domain - Cyclic in Zonal Direction

Initial Condition

They place a Gaussian hill centered at point (10,15) with a dimensional height of 2_000m and an e-folding width of 1_000km

Errors

Table 8:Parameters

Covariance MatrixStandard Deviation (σ\sigma)Horizontal Correlation (chc_h)Vertical Correlation (cvc_v)Scales
B\mathbf{B}0.80.6e60.2Short
Qs\mathbf{Q}_s0.0055550.6e60.2Short
Ql\mathbf{Q}_l0.0055551.6e60.8Long
Qki\mathbf{Q}_{k_i}0.81.6e60.8Long
R\mathbf{R}0.20.00.8Grid Pint
References
  1. Amraoui, S., and Didier Auroux, Blum, J., & and, E. C. (2023). Back-and-forth nudging for the quasi-geostrophic ocean dynamics with altimetry: Theoretical convergence study and numerical experiments with the future SWOT observations. Discrete and Continuous Dynamical Systems - S, 16(2), 197–219. 10.3934/dcdss.2022058
  2. Frezat, H., Sommer, J. L., Fablet, R., Balarac, G., & Lguensat, R. (2022). A posteriori learning for quasi-geostrophic turbulence parametrization. 10.48550/ARXIV.2204.03911
  3. Srinivasan, K., Chekroun, M. D., & McWilliams, J. C. (2023). Turbulence closure with small, local neural networks: Forced two-dimensional and β-plane flows. arXiv. 10.48550/ARXIV.2304.05029
  4. Gupta, J. K., & Brandstetter, J. (2022). Towards Multi-spatiotemporal-scale Generalized PDE Modeling. arXiv. 10.48550/ARXIV.2209.15616
  5. Guan, Y., Subel, A., Chattopadhyay, A., & Hassanzadeh, P. (2023). Learning physics-constrained subgrid-scale closures in the small-data regime for stable and accurate LES. Physica D: Nonlinear Phenomena, 443, 133568. 10.1016/j.physd.2022.133568
Quasi-Geostrophic Model
Overview
Quasi-Geostrophic Model
SSH Free-Run QG