Modified (hyper-)viscosity for coarse-resolution ocean models

We use the same multi-layer quasi-geostrophic model in a non-periodic rectangular domain as in [7]. Here, we only give a brief review of this system. The quasi-geostrophic pressure and potential vorticity (PV) are stacked in three isopycnal layers. We adopt vector forms to denote the layered pressure and potential vorticity (PV):

$$\mathbf{p}=\begin{bmatrix}p_{1}\\ p_{2}\\ p_{3}\end{bmatrix},\quad\mathbf{q}=\begin{bmatrix}q_{1}\\ q_{2}\\ q_{3}\end{bmatrix}\ .$$

The forced and damped quasi-geostrophic (QG) equations can be then written as

	$\displaystyle\partial_{t}\mathbf{q}=\frac{1}{f_{0}}J(\mathbf{q},\mathbf{p})+f_{0}B\mathbf{e}+\frac{1}{f_{0}}\big{(}a_{2}\Delta-a_{4}\Delta^{2}\big{)}(\Delta\mathbf{p}),$		(1)
	$\displaystyle\big{(}\Delta-f_{0}^{2}A\big{)}\mathbf{p}=f_{0}\mathbf{q}-f_{0}\beta(y-y_{0}),$		(2)

where $\Delta=\partial_{xx}^{2}+\partial_{yy}^{2}$ denotes the horizontal Laplacian, $\Delta^{2}$ the bi-laplacian operator, $J(a,b)=\partial_{x}a\partial_{y}b-\partial_{x}b\partial_{y}a$ stands for the Jacobi operator, $f_{0}+\beta(y-y_{0})$ is the Coriolis parameter under beta-plane approximation with the meridional axis center $y_{0}$, $a_{2}$ and $a_{4}$ are the Laplacian and biharmonic viscosity coefficients. Parameters of the configuration are listed in the Tables 1 and 2 in Appendix. Besides, the second term on the right-hand side of equation ((1)) represents the external forcing applied on different layers. In this work, we only consider an idealized case in which the ocean basin is driven by a stationary and symmetric wind stress $\vec{\tau}=(\tau^{x},\tau^{y})$ on the surface and by a linear Ekman stress at the bottom. In that case, the forcing term can be specified by

$$B=\begin{bmatrix}\frac{1}{H_{1}}&\frac{-1}{H_{1}}&0&0\\ 0&\frac{1}{H_{2}}&\frac{-1}{H_{2}}&0\\ 0&0&\frac{1}{H_{3}}&\frac{-1}{H_{3}}\\ \end{bmatrix},\ \quad\mathbf{e}=\begin{bmatrix}\partial_{x}\tau^{y}-\partial_{y}\tau^{x}\\ 0\\ 0\\ \frac{\delta_{\rm ek}}{2|f_{0}|}\Delta p_{3}\end{bmatrix},\ \quad\vec{\tau}=\tau_{0}\begin{bmatrix}-\cos(2\pi y/L_{y})\\ 0\end{bmatrix},$$

where $\tau_{0}$ is the magnitude of surface wind, $H_{k}$ is the background thickness of layer $k$, and $\delta_{\rm ek}$ is the bottom Ekman layer thickness. The vertical stratification level of such a model is described by the term $-f_{0}^{2}A\mathbf{p}$ in Equation (2) with

$$A=\begin{bmatrix}\frac{1}{H_{1}g_{1.5}^{\prime}}&\frac{-1}{H_{1}g_{1.5}^{\prime}}&0\\ \frac{-1}{H_{2}g_{1.5}^{\prime}}&\frac{1}{H_{2}}\left(\frac{1}{g_{1.5}^{\prime}}+\frac{1}{g_{2.5}^{\prime}}\right)&\frac{-1}{H_{2}g_{2.5}^{\prime}}\\ 0&\frac{-1}{H_{3}g_{2.5}^{\prime}}&\frac{1}{H_{3}g_{2.5}^{\prime}}\end{bmatrix},$$

where $g_{k+0.5}$ is the reduced gravity defined across the interface between layers $k$ and $k+1$. A multi-layered generalization of this model can be found in [6]. Note also that such a multi-layered model can be considered as a vertical discretized approximation of the continuously stratified QG system [18] with $\partial_{z}(f_{0}\partial_{z}\mathbf{p}/\vec{N}^{2})\approx-f_{0}A\mathbf{p}$ approximated by finite differences, and in which $\vec{N}$ denotes the buoyancy (or Brunt-Vaisala) frequency.

To facilitate numerical developments and benefit from built-in automatic differentiation, we develop a Pytorch [13] implementation of the above-described multilayer QG model ¹¹1Available at https://github.com/louity/qgm_pytorch. For this purpose, we follow rigorously the strategy of [8]:

1. 1.

   we use a regular numerical grid with finite differences

2. 2.

   We solve the PV advection equation (1) on the whole domain except the boundaries. We use a standard 5-point finite difference scheme for the (bi-)laplacian and the energy-enstrophy conservative Arakawa-Lamb scheme for the Jacobian [1].

3. 3.

   We apply a vertical change of coordinate to equation (2) which becomes a set of three inhomogeneous Helmholtz equations. We solve these equations with the spectral Discrete Sine Transform (DST) method, and we add corresponding homogeneous Helmholtz equation solutions to ensure mass conservation.

4. 4.

   We update the boundary values of the potential vorticity $\mathbf{q}$ using equation (2).

Detailed equations and numerical routine design choices can be found in [8]. We use a Heun-Runge-Kutta 2 time-stepping instead of the Leap-Frog time scheme used by [8].

For sake of numerical efficiency, we follow the recommendation of [15]: we compile computationally demanding routines and simplify finite difference calculations by reducing as much as possible the number of multiplications. We end up with a very concise code (less than 300 lines) that only depends upon Numpy and Pytorch libraries. This implementation will be open-sourced at the time of the publication.

We consider two spatial settings for our simulations :

1. 1.

   The eddy-resolving regime, our high-resolution reference with a 5 km resolution.

2. 2.

   The eddy-permitting regime, our low-resolution setting with a 40 km resolution.

Parameters for these two different regimes are written in table 2 in Appendix.

[16] studied the resulting flows’ differences between these two regimes. The high-resolution eddy-resolving model shows a well-pronounced eastward jet fuelled by mesoscale eddies circulating while the low-resolution eddy-permitting model does not induce a proper eastward jet as shown on Fig. 1. Temporal statistics significantly differ between high- and low-resolution simulations.

In both resolutions, we use biharmonic viscosity as in [17, 11, 5, 4] essentially because it is less dissipating at large scales than a Laplacian. Compared to the usual Laplacian viscosity, it preserves large-scale structures. However, hyperviscosity remains much too dissipative in the “eddy-permitting” regime [10]. This too strong dissipation kills the eastward jet that is present in the high-resolution and that we expect to see in such a double-gyre quasi-geostrophic model. Fig. 2 shows a sequence of snapshots of the low-resolution models where we input a downsampled snapshot of the high-resolution (see Appendix for details on downsampling). After as few as three years, the eastward jet has almost disappeared, showing that the model is too dissipating. Lowering the hyper-viscosity coefficient by a factor of 10 does not solve this problem, and creates spurious gradients in the potential vorticity as shown in Fig. 2. These numerical artifacts are due to a bad representation of the direct enstrophy cascade, causing a piling up of the small-scale vorticity gradients at the cut-off frequency together with aliasing effects.

Here we propose a simple affine modification parameterization of hyperviscosity. We add a bias to the term $\Delta\mathbf{p}$ in eq. (1), which becomes $\Delta\left(\mathbf{p}-\mathbf{p}^{\prime}\right)$ where $\mathbf{p}^{\prime}$ is a dimensional field that does not depend upon time. The PV advection equation with hyperviscosity becomes

$\displaystyle\partial_{t}\mathbf{q}=\frac{1}{f_{0}}J(\mathbf{q},\mathbf{p})+f_{0}B\mathbf{e}+\frac{1}{f_{0}}\big{(}a_{2}\Delta-a_{4}\Delta^{2}\big{)}\left(\Delta(\mathbf{p}-\mathbf{p}^{\prime})\right).$

(3)

The elliptic equation (2) remains unchanged.

The goal of this additional term is to reproduce a relevant time-average pressure field relying on observations or high-resolution solutions. For example the high-resolution average $\overline{\mathbf{p}_{\rm HR}}$ can be downsampled to the targeted coarse grid resolution in $\overline{\mathbf{p}_{\rm HR}}\downarrow$, and we want the average of the modified low-resolution $\overline{\mathbf{p}_{\rm LR}}$ model to be as close as possible to the high-resolution reference $\overline{\mathbf{p}_{\rm HR}}\downarrow$.

We face here an optimal control problem, as the low-resolution average is a function of the control parameter $\mathbf{p}^{\prime}$. We state it with the following least-square formulation

	$\displaystyle\mathbf{p}^{\prime}_{\rm opt}$	$\displaystyle=\underset{\mathbf{p}^{\prime}}{\text{argmin }}\mathcal{F}\left(\mathbf{p}^{\prime}\right)$		(4)
	$\displaystyle\mathcal{F}\left(\mathbf{p}^{\prime}\right)$	$\displaystyle=\big{\|}\ \overline{\mathbf{p}_{\rm LR}}\left(\mathbf{p}^{\prime}\right)-\overline{\mathbf{p}_{\rm HR}}\downarrow\big{\|}^{2}$		(5)

This optimization problem is a priori non-convex and we shall not expect to find a global optimum. In the following, we propose a numerical procedure to find a heuristic $\hat{\mathbf{p}^{\prime}}$ of the optimal solution $\mathbf{p}^{\prime}_{\rm opt}$.

Computationally, the implementation of this modified hyperviscosity is simple and computationally cheap. We precompute $\Delta\mathbf{p}^{\prime}$ and subtract it from $\Delta\mathbf{p}$ at each time-integration step. It increases the integration time of the advection equation (1) by less than 1% on CPUs and GPUs.

The continuously stratified QG equations can be rewritten in a variational formulation [9] with a Hamiltonian $\mathcal{J}$ defined as

$$\mathcal{J}(\mathbf{p})=\frac{1}{2}\int_{\Omega}\frac{1}{f_{0}}|\nabla\mathbf{p}|^{2}+\frac{f_{0}}{N^{2}}(\partial_{z}\mathbf{p})^{2}.$$

Our model is a discretized version of the continuous stratification. Since we add an external wind forcing term and we use an energy conservative Arakawa advection scheme, we need to add some viscosity or hyperviscosity to dissipate energy. In a variational formulation, these (hyper-)viscous terms become the following penalization

$$\frac{1}{2}\int_{\Omega}a_{2}|\Delta\mathbf{p}|^{2}+a_{4}\left|\nabla\left(\Delta\mathbf{p}\right)\right|^{2},$$

added to the Hamiltonian $\mathcal{J}(\mathbf{p})$ to produce a smooth solution. The Gradient norm penalization of Laplacian $\mathbf{p}$ guides the minimization toward solutions of smooth Laplacian. Hyperviscosity corresponds to the Laplacian norm penalization and enforces a solution of minimum Laplacian norm. The parameters $a_{2}$ and $a_{4}$ quantify the strength of these regularization constraints.

Here, we simply propose to replace it with the following penalization

$$\frac{1}{2}\int_{\Omega}a_{2}|\Delta(\mathbf{p}-\mathbf{p}^{\prime})|^{2}+a_{4}\left|\nabla\left(\Delta(\mathbf{p}-\mathbf{p}^{\prime})\right)\right|^{2}\ .$$

We now penalize $(\mathbf{p}-\mathbf{p}^{\prime})$ instead of $\mathbf{p}$, meaning that we guide the solution to a possibly non-smooth reference $\mathbf{p}^{\prime}$ that will produce the correct large scale behavior.

Here we present a method to find a solution to the optimization problem (4). A natural guess for $\mathbf{p}^{\prime}_{\rm opt}$ is $\overline{\mathbf{p}_{\rm HR}}\downarrow$. We solve the equations and compute the average pressure $\overline{\mathbf{p}_{\rm LR}}$. Results are shown in Fig. 4. It is a good first-guess, but the difference $\overline{\mathbf{p}_{\rm HR}}\downarrow-\overline{\mathbf{p}_{\rm LR}}$ is still large.

We propose the following iterative procedure to find a better guess for $\mathbf{p}^{\prime}_{\rm opt}$. In the following we assume that we are in low resolution, i.e. $\mathbf{p}=\mathbf{p}_{\rm LR}$ and $\overline{\mathbf{p}}=\overline{\mathbf{p}_{\rm LR}}$ unless explicitly written.

• •

  We set $\mathbf{p}^{\prime}_{0}$ and we compute the average pressure $\overline{\mathbf{p}}_{0}$ solving standard equations (1,2) without modified viscosity.

• •

  Choose $k\in\left]0,1\right]$.

• •

  Start with $\mathbf{p}^{\prime}_{1}=\overline{\mathbf{p}_{\rm HR}}\downarrow$.

• •

  Evolve the ensemble for $n$ years and compute the corresponding average pressure $\overline{\mathbf{p}}_{1}$ with ensemble average.

• •

  For $n=1\dots$:

  • –

    Set $\mathbf{p}^{\prime}_{n+1}=\mathbf{p}^{\prime}_{n}+k\left(\overline{\mathbf{p}_{\rm HR}}\downarrow-\overline{\mathbf{p}}_{n}\right)$.

  • –

    Evolve the ensemble for $n$ years and compute new average pressure $\overline{\mathbf{p}}_{n+1}$.

• •

  return $\mathbf{p}^{\prime}_{n}$ and $\overline{\mathbf{p}}_{n}$

There is no theoretical guarantee that this procedure converges, but we observe in the next section that it converges with the double-gyre QG model that we use.

We use ensemble averages to compute the statistics. To create ensembles of size $N$, we start from a zero solution and spin up the models for 100 years with a timestep of 1200 seconds to reach statistically steady states as in [14]. Then we run the models for 500 years and save 10 snapshots a year to get 5000 snapshots, and we randomly select $N$ snapshots out of these 5000 snapshots. The ensemble averages are simply average over these $N$ ensemble members that we evolve in parallel. Such ensemble averages are denoted with $\overline{\bullet}$ in the following, i.e. the average pressure is denoted by $\overline{\mathbf{p}}$, average velocity by $\overline{u}$, etc.

We test the iterative procedure described in section 4.2 with the double-gyre model presented in section 2 in the eddy-permitting regime. We use $n=10$ years to evolve the ensemble after each iterate. We compute the reference pressure average $\overline{\mathbf{p}}_{\rm HR}$ with the same model in the eddy-resolving regime.

Fig. 3 shows the relative square error $\|\overline{\mathbf{p}}_{n}-\overline{\mathbf{p}_{HR}}\downarrow\|^{2}/\|\overline{\mathbf{p}_{HR}}\downarrow\|^{2}$ at iterations of the procedure with $k=1$ and $k=0.7$. The procedure converges with $k=0.7$ and oscillates with $k=1$.

Fig. 4 shows the output average pressure $\overline{\mathbf{p}}_{n}$ of the iterative procedure, the reference $\overline{\mathbf{p}}_{\rm HR}$ and the difference between the two, as well as for zonal velocity $\mathbf{u}$ . Our model can reproduce the eastward jet produced by the high-resolution reference model. Kinetic energy spectra shown on Fig. 5 shows also the improvement of our model compared to low-resolution. Finally, Fig. 6 shows high-resolution and low-resolution snapshots as well as a snapshot of the proposed model at low-resolution. Our model effectively produces the eastward jet and a re-circulation zone around it where eddies are created. Artifacts can be also observed on the zonal velocity and potential vorticity on the right of Fig. 6. They can likely be Rossby waves created by the harmonic regularization terms, which remain an artificial constraint, but this needs to be studied further.

Downsampling the high-resolution solution on a low-resolution grid consists of interpolating the high-resolution (769$\times$ 961) streamfunction on the low-resolution (97$\times$121) grid. Then we can compute the potential vorticity using equation (2). Because of the no-flow constraint, the downsampled streamfunction should be constant on the boundaries and should satisfy a mass conservation constraint [8]. We also want to preserve the frequency information and prevent aliasing.

We use the following procedure:

1. 1.

   we apply a Gaussian filter and downsample the streamfunction on the domain except on the boundaries.

2. 2.

   we adding homogeneous solutions of the Elliptic equation (2) to the streamfunction in order to satisfy the mass conservation as in [8].