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

We present a simple parameterization for coarse-resolution ocean models. To replace computationally expensive high-resolution ocean models, we develop a computationally cheap parameterization for coarse-resolution models based solely on the modification of the viscosity term in advection equations. It is meant to reproduce the mean quantities like pressure, velocity, or vorticity computed from a high-resolution reference solution or using observations. We test this new parameterization on a double-gyre quasi-geostrophic model in the eddy-permitting regime. Our results show that the proposed scheme improves significantly the energy statistics and the intrinsic variability on the coarse mesh. This method shall serve as a deterministic basis model for coarse-resolution stochastic parameterizations in future works.


Introduction
Ocean general circulation models used at climatic scales are limited for evident computational reasons to too coarse horizontal resolutions to solve correctly ocean mesoscale and sub-mesoscale eddies, even with large computational infrastructure. The horizontal resolution of the most recent climatic ocean models is of the order of the Rossby radius of deformation. These models are hence in the so-called eddy-permitting regime and they can solve partially the mesoscale (i.e. 10-100 km) eddy field. These models however suffer from strong limitations. In particular, they are unable to reproduce accurately large-scale structures such as the eastward turbulent jet in an idealized double-gyre configuration.
Recent parameterizations have shown significant improvements in coarse-resolution models compared to high-resolution reference solutions [3]. However, it remains an important topic of research, as the actual generation of parametrizations is not completely able to resolve the effects of the unresolved scales on the large-scale flow structures.
A wide range of subgrid parametrizations relies on eddy viscosity such as Laplacian and biharmonic schemes [17,11,5,4]. It has been shown in [10] that including only these (hyper)viscosity in coarseresolution models often causes too much dissipation and results in an artificial energy sink at large scales. In general, even eddy-permitting models are not energetic enough and as a result, the long-time average of any coarse model's variable of interest departs completely from the long-time average of high-resolution models subsampled at the same scale. This becomes the main motivation of the present work. In particular, we would like to answer the following question: how can we reduce the excessive resolved kinetic energy loss due to the viscosity while simultaneously ensuring numerical stability?
We propose a simple affine parameterization of (hyper)viscosity. The (bi)laplacian operator Δ is replaced by Δ ( − ), where is a field of same dimension as that does not depend upon time. We interpret this method as a mathematical regularization technique to guide the solutions towards prior information. We frame as the solution of an optimal control problem to reproduce statistics computed from a reference solution or observations. We present a method to solve this optimal control problem.
We test the proposed method with an idealized double-gyre configuration. For that purpose, we release with this article a fast, concise, and CPU-GPU portable Pytorch implementation of a multi-layer quasi-geostrophic model on a rectangular domain. We implement and test our optimization procedure within this setting.
This article is organized as follows: we present in section 2 the double gyre quasi-geostrophic model we use and detail its implementation, we present in section 3 our modified viscosity parameterization and we show and discuss numerical results in section 4.

Governing equations
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): The forced and damped quasi-geostrophic (QG) equations can be then written as where Δ = 2 + 2 denotes the horizontal Laplacian, Δ 2 the bi-laplacian operator, ( , ) = − stands for the Jacobi operator, 0 + ( − 0 ) is the Coriolis parameter under betaplane approximation with the meridional axis center 0 , 2 and 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 ì = ( , ) on the surface and by a linear Ekman stress at the bottom. In that case, the forcing term can be specified by where 0 is the magnitude of surface wind, is the background thickness of layer , and ek is the bottom Ekman layer thickness. The vertical stratification level of such a model is described by the term − 2 0 p in Equation (2) where +0.5 is the reduced gravity defined across the interface between layers and + 1. A multilayered 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 ( 0 p/ ì 2 ) ≈ − 0 p approximated by finite differences, and in which ì denotes the buoyancy (or Brunt-Vaisala) frequency.

Pytorch implementation
To facilitate numerical developments and benefit from built-in automatic differentiation, we develop a Pytorch [13] implementation of the above-described multilayer QG model . For this purpose, we follow rigorously the strategy of [8]: Available at https://github.com/louity/qgm_pytorch 1. we use a regular numerical grid with finite differences 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. 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. We update the boundary values of the potential vorticity 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.

Eddy-resolving and eddy-permitting regimes
We consider two spatial settings for our simulations : 1. The eddy-resolving regime, our high-resolution reference with a 5 km resolution.
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 eddyresolving 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.

Motivation
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.

Modified viscosity
Here we propose a simple affine modification parameterization of hyperviscosity. We add a bias to the term Δp in eq. (1), which becomes Δ (p − p ) where p is a dimensional field that does not depend upon time. The PV advection equation with hyperviscosity becomes  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 p HR can be downsampled to the targeted coarse grid resolution in p HR ↓, and we want the average of the modified low-resolution p LR model to be as close as possible to the high-resolution reference p HR ↓.
We face here an optimal control problem, as the low-resolution average is a function of the control parameter p . We state it with the following least-square formulation 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 heuristicp of the optimal solution p opt .
Computationally, the implementation of this modified hyperviscosity is simple and computationally cheap. We precompute Δp and subtract it from Δp at each time-integration step. It increases the integration time of the advection equation (1) by less than 1% on CPUs and GPUs.

Modified viscosity regularization
The continuously stratified QG equations can be rewritten in a variational formulation [9] with a Hamiltonian J defined as 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 added to the Hamiltonian J (p) to produce a smooth solution. The Gradient norm penalization of Laplacian 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 2 and 4 quantify the strength of these regularization constraints.
Here, we simply propose to replace it with the following penalization We now penalize (p − p ) instead of p, meaning that we guide the solution to a possibly non-smooth reference p that will produce the correct large scale behavior.

Iterative procedure
Here we present a method to find a solution to the optimization problem (4). A natural guess for p opt is p HR ↓. We solve the equations and compute the average pressure p LR . Results are shown in Fig. 4. It is a good first-guess, but the difference p HR ↓ −p LR is still large. We propose the following iterative procedure to find a better guess for p opt . In the following we assume that we are in low resolution, i.e. p = p LR and p = p LR unless explicitly written.
• We set p 0 and we compute the average pressure p 0 solving standard equations (1,2) without modified viscosity.
• Evolve the ensemble for years and compute the corresponding average pressure p 1 with ensemble average. • For = 1 . . . : -Evolve the ensemble for years and compute new average pressure p +1 .

• return p and p
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.

Statistics
We use ensemble averages to compute the statistics. To create ensembles of size , 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 snapshots out of these 5000 snapshots. The ensemble averages are simply average over these ensemble members that we evolve in parallel. Such ensemble averages are denoted with • in the following, i.e. the average pressure is denoted by p, average velocity by , etc.

Iterative procedure
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 = 10 years to evolve the ensemble after each iterate. We compute the reference pressure average p HR with the same model in the eddy-resolving regime. Fig. 3 shows the relative square error p − p ↓ 2 / p ↓ 2 at iterations of the procedure with = 1 and = 0.7. The procedure converges with = 0.7 and oscillates with = 1. Fig. 4 shows the output average pressure p of the iterative procedure, the reference p HR and the difference between the two, as well as for zonal velocity 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.

Conclusion
We presented a simple modified-viscosity scheme for coarse resolution ocean modeling that we derived and tested on a double-gyre multi-layer quasi-geostrophic model. We interpret it as a modified regularization technique that will guide the solution to a reference rather than producing a too smooth solution in the eddy-permitting regime. The technique requires solving an optimization problem, and we presented a procedure to find a good guess for the solutions. We showed that it converges to a reasonable solution that fairly reproduces the input reference.
If this method mimics the average of the high-resolution, it only reproduces partially the variability and higher-order statistics of the high-resolution. We see in Fig. 5 our model's snapshots resemble the averages. In future works, we consider using this method as a deterministic basis for stochastic parameterizations such as Location-Uncertainty [12].