1.5-layers Quasi-Geostrophic model

Equations solved

We model the ocean by the superposition of one dynamical thin upper layer on top of a denser and thicker static layer, with flat bathymetry (see figure).

../_images/model_one_layer.png

The adiabatic Quasi-Geostrphic (QG) dynamics in the 1.5-layer model are governed by the conservation of Potential Vorticity (PV), \(q\), whose prognostic equation without atmospheric forcing is:

\[\frac{\partial q}{\partial t} + \mathbf{u_g} \cdot \nabla q = 0\]

where \(\mathbf{u_g}\) is the geostrophic velocity vector diagnosed from streamfunction \(\psi\), proportional to SSH \(\eta\):

\[\mathbf{u_g}=-\mathbf{k}\times \mathbf{\psi}=-\frac{g}{f} \mathbf{k}\times \mathbf{\nabla} \eta\]

\(\mathbf{k}\) denotes the vertical direction, \(g\) the gravity constant, and \(f\) the Coriolis parameter. The PV is linked to the streamfunction by the elliptical equation:

\[q=\nabla^{2}\psi - {\frac {f^2}{{c_1}^{2}}}\psi\]

where \(c_1\) is the first baroclinic phase velocity.

Formulation with Sea Level Anomalies

In reality, an altimeter only provides accurate observations of the time-fluctuating part of SSH, called Sea Level Anomaly (SLA). The time-averaged of SSH, called Mean Dynamical Topography (MDT) is computed with the combination of data and other satellite observations.

To formulate the QG dynamics with the SLA variable, we decompose the geostrophic velocity and the PV using the Reynolds decomposition:

\[\begin{split}\begin{aligned} \mathbf{u_g}&=\overline{\mathbf{u_g}}+\mathbf{u_g}' \\ q &= \overline{q}+q' \end{aligned}\end{split}\]

where \(\overline{\mathbf{u_g}}\) and \(\overline{q}\) stand for the mean components (SLA-independent) and \(\mathbf{u_g'}\) and \(q'\) the time-fluctuating components (diagnosed from SLA).

The prognostic equation for the potential vorticity fluctuation is then:

\[\frac{\partial q'}{\partial t} + \mathbf{\overline{u_g}} \cdot \nabla q' + \mathbf{u_g'} \cdot \nabla \overline{q} + \mathbf{u_g'} \cdot \nabla q' = 0\]

Numerical implementations

The numerical codes of the model are located in the mapping/models/model_qg1l directory. Two different implementations of the model are provided:

  • qgm.py with the numpy python library;

  • jqgm.py with jax python library, enabling automatic differentiation for 4Dvar data assimilation technique.

The time stepping from time \(t_i\) to \(t_{i+1}\) is implemented as follows.

  • The PV \(q_i\) is diagnosed from \(\eta_i\) (see Qgm.h2pv function).

  • \(q_{i+1}\) is computed by integrating the prognostic equation (see Qgm.qrhs function).

  • \(\psi_{i+1}\) (and therefore \(\eta_{i+1}\)) is obtained by inverting the elliptical equation using a gradient conjugate technique (see Qgm.pv2h function).

Both forward Euler and fourth-order Runge-Kutta time schemes are implemented (see Qgm.step_euler and Qgm.step_rk4 functions respectively).

The advection can be performed either with first, second, or third upwind spatial schemes (see Qgm._rq1, Qgm._rq2 and Qgm._rq3 functions respectively).

Both numpy and jax implementations allow Reynolds decomposition, with the prescription of a mean field such as the one provided by the MDT CMEMS products (see Example ???). In addition, jax impelmentation allows to consider multiscale QG formulation, closed to the SLA formulation, but with a time-variable ‘mean’ flow, that we can associate to the large spatial scales (see jqgm.Qgm.step_multiscales function). This will be useful for providing background data in data assimilation experiments (see Example ???).

Configuration parameters

Here are the parameters specific to the model to be prescribed in the configuration file. If some are not prescribed, then the default values shown hereafter will be used.

# Name of the model. Either:
# 'QG1L': numpy implementation of the model
# 'JAX-QG1L': JAX implementation of the model
# 'QG1LM': JAX implementation of the multiscale model
name_model = 'QG1L'

# Time scheme of the model (e.g. Euler,rk4)
qg_time_scheme = 'Euler'

# Order of the upwind scheme for PV advection (either 1,2 or 3)
upwind = 3

# Order of the upwind scheme for adjoint PV advection (either 1,2 or 3).
# Only relevent for 4Dvar assimilation scheme.
# If None, default value is set equal to upwind.
upwind_adj = None

# Number of iterations to perform the gradient conjugate algorithm in the QG elliptical equation
qgiter = 20

# Number of iterations to perform the gradient conjugate algorithm in the adjoint QG elliptical equation
# Only relevent for 4Dvar assimilation scheme.
# If None, default value is set equal to qgiter.
qgiter_adj = None

# First baroclinic phase velocity
c0 = 2.7

# Path of the auxilliary file to be used as phase velocity field.
# Spatial interpolation is handled inline
# c0 has to be set to None
filec_aux = None

# Variable names for the phase velocity auxilliary file
name_var_c = {'lon':'','lat':'','var':''}

# Whether to use Reynolds decomposition or not (True or False).
# Observed variable has to be SLA and MDT has to be provided
Reynolds = False

# Path of the auxilliary file to be used as MDT
# Spatial interpolation is handled inline
path_mdt = None

# Variable names for the MDT auxilliary file
# If 'mdu' and 'mdv' are None, the mean currents will be computed thanks to MDT and the geostrophy assumption
name_var_mdt = {'lon':'','lat':'','mdt':'','mdu':None,'mdv':None}

# Diffusion coefficient
# Note: here the diffusion is done on PV (instead of SSH in the Diffusion model, previous section)
Kdiffus = 0

# Whether to use only diffusion or not.
# Useful for quick tests
only_diffusion = False