A posteriori learning of quasi-geostrophic turbulence
parametrization: an experiment on integration steps

Hugo Frezat
Univ. Grenoble Alpes, CNRS UMR LEGI, Grenoble, France
Univ. Grenoble Alpes, CNRS UMR IGE, Grenoble, France
IMT Atlantique, CNRS UMR Lab-STICC, Brest, France
hugo.frezat@univ-grenoble-alpes.fr
Julien Le Sommer
Univ. Grenoble Alpes, CNRS UMR IGE, Grenoble, France
Ronan Fablet
IMT Atlantique, CNRS UMR Lab-STICC, Brest, France
Guillaume Balarac
Univ. Grenoble Alpes, CNRS UMR LEGI, Grenoble, FranceInstitut Universitaire de France (IUF), Paris, France
Redouane Lguensat
Learning, Data and Robotics Lab, ESIEA, Paris, France
LOCEAN-IPSL, Sorbonne Université, Institut Pierre Simon Laplace, Paris, France

Abstract
Modeling the subgrid-scale dynamics of reduced models is a long standing open
problem that finds application in ocean, atmosphere and climate predictions where
direct numerical simulation (DNS) is impossible. While neural networks (NNs)have already been applied to a range of three-dimensional flows with success, twodimensional flows are more challenging because of the backscatter of energy fromsmall to large scales. We show that learning a model jointly with the dynamicalsolver and a meaningful a posteriori-based loss function lead to stable and realistic
simulations when applied to quasi-geostrophic turbulence.

1

Introduction

Understanding and predicting the evolution of various natural systems would not be possible without
simulations of turbulent flows. However, solving all the spatial and temporal scales of the associatedpartial differential equation (PDE), i.e., the Navier-Stokes equations remains computationally prohibitive. One popular solution (e.g. [14, 8, 4]) is to resolve only the largest scales and use subgrid
closures (or physical parametrizations) to represent the smaller ones.
Recently, neural networks (NNs) have been proposed as a promising alternative to algebraic
parametrizations in three-dimensional incompressible turbulence [6, 15, 5]. Two-dimensional problems, however, are more challenging due to the inverse cascade of energy that leads to negativeviscosities, and it has been demonstrated that numerical stability of the trained model in decaying
turbulence requires either the removal of negative eddy viscosities [13] or a large training dataset [9].
Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

y(t)f (y(t))y(t + ∆t)···y(t + N ∆t)f (y(t+N ∆t))y(t + (N + 1)∆t)

T (y(t))g(ȳ(t))Laposteriori···

Laposteriorig(ȳ(t + N ∆t))

Laposteriori

ȳ(t)

+

ȳ(t + ∆t)···ȳ(t + N ∆t)+ȳ(t + (N + 1)∆t)

MNN (ȳ(t))

MNN (ȳ(t + N ∆t))

LaprioriR(ȳ(t))Figure 1: Sketch of one learning step for the presented strategies. The a priori loss is computed at
instantaneous time t (dashed, red), while the a posteriori loss is evaluated on states forward in time
(dashed, blue).
2D turbulence is also a simple model that well approximates more complex rotating stratified flows(found in atmosphere and ocean dynamics).
We show that on a forced configuration of 2D turbulence, stability can be obtained when the NN
is trained end-to-end on a temporal horizon, i.e. jointly with the forward solver (see Fig. 1).
At equivalent computational complexity, the standard learning strategy does not lead to a stableparametrization and is thus not usable in practice. Additionally, the presented strategy allows us
to constrain the calibration of the model on a posteriori metrics, i.e. the evolution of quantities ofinterest over time, which are the absolute performance measures in turbulence modeling. Finally, weshow that the number of iterations performed during the training process have a significant impact on
the long-time statistics of the model.

2

A priori and a posteriori learning

In the machine learning community, it is well known that models can take advantage of additional
informations during training. For instance, the basic, purely data-based approach has been extended
to train models end-to-end using an entire system. In practice, some examples of those systems
include autonomous vehicles [2], objects detection [19] and super-resolution [3]. This idea has also
recently found many applications in the physical sciences community (e.g. [10, 11]).
In the context of partial differential equation parametrization, we are given with a high-resolution
system f (y) and a low resolution system g(ȳ) describing the evolution of high and low resolution
variables y(t) and ȳ(t) respectively by,∂y= f (y),
y∈Ω



 ∂t
∂ ȳ
(1)
= g(ȳ) + R(ȳ),
ȳ ∈ Ω̄


∂t




T (y) = ȳwhere Ω̄ is a smaller domain than Ω and T is a known projection operator that goes from high to low
resolution variables.
In this framework, we are interested in modeling the residual term, or parametrization R(ȳ) with aneural network MNN (ȳ), based on high resolution data coming from a direct numerical simulationof f (y), collected in a dataset
D := {ȳ(t)} → {R(ȳ(t))}.
(2)
If we connect the machine learning naming to fluid dynamics concepts, we could say that the databased approach only gives the ability to optimize a model on a priori metrics, i.e., instantaneouspredictions of the targeted term. With the end-to-end approach, however, we have access to simulation
quantities over time, which can be used to train the model on a posteriori metrics (See [17] for2

more details). Using the data D, a standard learning strategy based on the direct prediction of theresidual term can be written as arg minθ L(R(ȳ), MNN (ȳ|θ)), which has already been applied to
fluid dynamics parametrization, and corresponds to an a priori minimization problem. Now, if wehave a differentiable solver for g(ȳ), the standard minimization problem performed during learning
can be rewritten as a series of a posteriori losses on low resolution variables that are advanced intime, such thatarg min L(y(t), ȳ(t0 ) + Φθ (ȳ(t0 ), t0 , t)), t ∈ [0, T ](3)θ

where L is the loss function, [0, T ] is the temporal horizon on which the forward solver is integrated
and Φ the flow operator that advances the reduced model in time, expressed as
Z tΦθ (ȳ, t0 , t) =g(ȳ(τ )) + MNN (ȳ(τ )|θ) dτ.
(4)t0

In practice, the forward solver performs N discrete timesteps, which defines the temporal horizon
T = N ∆t. A sketch of the a posteriori learning strategy is shown in Fig. 1.

3

Application to quasi-geostrophic turbulence parametrization

The described strategies are applied to a classical two-dimensional single-layer barotropic quasigeostrophic (QG) potential vorticity equation [12], defined by∂t ω + J(ψ, ω) = ν∇2 ω − µω + F(5)

u = (−∂y ψ, ∂x ψ)

(6)

with
2

ω=∇ ψ

(7)

where ω is the vorticity, ψ the stream function, u the velocity vector and J(ψ, ω) = ∂x ψ∂y ω −∂y ψ∂x ω is the Jacobian operator. The model is parametrized by a viscosity ν, a linear drag µ and asource term F .
Projecting to the reduced system. In turbulence modeling, a common approach is to suppose thatthe reduced system has been filtered so that the smallest scales are removed, and instead predictedby a model. In practice, the projection operator T is a spatial filter Λδ (x, y) at spatial scale δ > 0
defined such that
T (ω) := ω ∗ Λδ = ω̄(x, y, t)(8)and the reduced system becomes∂t ω̄ + J(ψ̄, ω̄) = ν∇2 ω̄ − µω̄ + F̄ + J(ψ̄, ω̄) − J(ψ, ω)
|
{z
}

(9)

R(ψ,ω)

where R(ψ, ω) is the subgrid-scale (SGS) term or parametrization that needs to be modeled.
Configuration. In our experiment, we discretize the full solutions on 2048 × 2048 grid points usinga pseudo-spectral method and use δ = 16 with a spectral cut-off filter Λδ (k) = 0, ∀k > π/δ∆ where∆ is the filter width so that our reduced solution lies on a 128 × 128 grid. The source term F is atime-dependent wind-forcing acting at large scales k = 4,
F = CF (t) [cos(4y + π sin(1.4t)) − cos(4x + π sin(1.5t))]

(10)

2

with steady enstrophy rate CF such that hF i/2 = 3. Following the procedure given by [8], we
spin-up multiple simulations with 1024 × 1024 grid initialized with random modes at k = 4 forapproximately 1300 days. The simulations are then integrated for 35 days on the high-resolution
grid to obtain a turbulent state that will be used as initial conditions for training and testing. Finally,
training data is extracted from 10 independent ≈ 65 days long (48000 iterations) trajectories. The
parameters of the simulations are summarized in dimensional and non-dimensional units in Table 3.
Our model 1 is implemented with PyTorch [16] in order to benefit from automatic differentiation,
which is required when computing the gradient of the flow operator (7).
1https://github.com/hrkz/torchqg

3

Table 1: Parameters of the forced configuration in dimensional and non-dimensional units, with
Ld (x) = 1.2 × 106 s and Td (x) = 504 × 104 /πm.
Parameter
Domain length (L)Spatial resolution (∆x)
Time step (∆t)
Linear drag (µ)
Viscosity (ν)
Reynolds number (Re)Dimensional
Numerical
104 km2π
10 km
3.00 × 10−3
120 s
1.00 × 10−4−8−1
1.25 × 10 m
2.00 × 10−2
2 −1
22.0 m s
1.02 × 10−5
≈ 220000

Figure 2: Vorticity fields at the end of the testing simulation with the filtered DNS (leftmost) and thediscussed models. Note that the model trained with a priori learning was numerically unstable andits states diverged to infinity.

Learning components. The neural network MNN used in this experiment for both a priori and aposteriori strategies is realized with a fully convolutional architecture. It consist of 10 convolutional
layers with 64 features each, interspresed with ReLU activation functions using kernel sizes of 5 × 5.
The parameters are optimized for a fixed number of epochs with Adam and a learning rate of 10−4 .
In particular, the a posteriori version is trained for a different number of iterations N = {1, 5, 30},
with a loss function that minimizes the mean squared error (MSE) of vorticity ω̄ over time, such thatN

Laposteriori :=1 X(ω(i∆t) ∗ Λδ − ω̄(i∆t))2
N i=1

(11)

Results. We run a new simulation at δ = 16, i.e. on a 128 × 128 grid, comparing the NNs to thereference direct numerical simulation DNS and the dynamic Smagorinsky model [18, 7], a purely
dissipative algebraic model which is a common baseline in the fluid dynamics community. The run
is performed for 3000 and 48000 iterations of the low resolution and direct numerical simulation
systems, respectively, since the timestep for low resolution systems can be δ times larger than the
DNS ∆t.
First, we can see in Fig. 2 the vorticity fields ω at the end of the new simulation. As expected, thedynamic Smagorinsky model is only dissipative and thus the smallest scales of the simulation are notvisible anymore. On the other side, the model learnt from the a priori strategy accumulated too muchenergy on the smallest scales which led to numerical instabilities, as expected from similar studies[13, 9]. The 3 models trained with the a posteriori strategy remain stable at the end of the simulation.
Moreover, these models are preserving the small-scale features of the flow. Now, we look at quadratic
invariants of 2D turbulence, which develops as a double cascade in statistically stationary conditions[1]. In particular, the enstrophy wavenumber spectrum Z(k) and flux ΠZ (k) are relevant measuresof the success of turbulence parametrizations, as proposed by [8].
In Fig. 3 (left), the enstrophy spectra confirm a large deviation at the smallest scales (highest
wavenumbers k) for the Smagorinsky model, and a visible deviation at the largest scales (smallest
wavenumbers k) for the a priori model and a posteriori model trained with N = 1. Overall statistical
performance of the model is verified by the enstrophy flux in Fig. 3 (right), where the NN trainedwith the a posteriori strategy at N = 30 is the closest to the exact flux predicted by the DNS.
4

0.8

100

Π Z ( k)Z ( k)

0.6
10-1
DNS
M DynSmagorinsky
M apriori= 30
MN
aposteriori(states)N=5
M aposteriori(states)N=1
M aposteriori(states)

10-2

100

k

0.4

DNS
M DynSmagorinsky
M apriori= 30
MN
aposteriori(states)N=5
M aposteriori(states)N=1
M aposteriori(states)

0.2

0.0
101

100

k

101

Figure 3: Enstrophy spectrum (left) and time-averaged enstrophy flux (right) of the different models
w.r.t the reference DNS (dashed). For the NNs models, only those trained with the a posteriori
strategy are stable and thus have valid statistics.

4

Conclusion

We show that a NN trained jointly with the forward solver on a posteriori informations lead to stable
simulations and outperforms the classical data-driven strategy and algebraic models based on different
statistical metrics for the prediction of quasi-geostrophic turbulence parametrization. However, it isimportant to consider that this requires a differentiable solver for the reduced order system, which is
not the case for most global circulation models developed nowadays.
Acknowledgments
The authors thank Laure Zanna, Olivier Pannekoucke and Corentin Lapeyre for the helpful discussions.
This research was supported by the CNRS through the 80 PRIME project and the ANR throughthe Melody, OceaniX, and HRMES projects. Additional support was provided by Schmidt Futures,a philanthropic initiative founded by Eric and Wendy Schmidt, as part of its Virtual Earth System
Research Institute (VESRI). Computations were performed using GPU resources from GENCIIDRIS.

References[1] Guido Boffetta. Energy and enstrophy fluxes in the double cascade of two-dimensional turbulence. Journalof Fluid Mechanics, 589:253–260, 2007.
[2] Mariusz Bojarski, Davide Del Testa, Daniel Dworakowski, Bernhard Firner, Beat Flepp, Prasoon Goyal,
Lawrence D Jackel, Mathew Monfort, Urs Muller, Jiakai Zhang, et al. End to end learning for self-drivingcars. arXiv preprint arXiv:1604.07316, 2016.
[3] Yu Chen, Ying Tai, Xiaoming Liu, Chunhua Shen, and Jian Yang. Fsrnet: End-to-end learning facesuper-resolution with facial priors. In Proceedings of the IEEE Conference on Computer Vision and Pattern
Recognition, pages 2492–2501, 2018.
[4] Anurag Dipankar, Bjorn Stevens, Rieke Heinze, Christopher Moseley, Günther Zängl, Marco Giorgetta,
and Slavko Brdar. Large eddy simulation using the general circulation model icon. Journal of Advances inModeling Earth Systems, 7(3):963–986, 2015.
[5] Hugo Frezat, Guillaume Balarac, Julien Le Sommer, Ronan Fablet, and Redouane Lguensat. Physicalinvariance in neural networks for subgrid-scale scalar flux modeling. Physical Review Fluids, 6(2):024607,
2021.
[6] Masataka Gamahara and Yuji Hattori. Searching for turbulence models by artificial neural network.
Physical Review Fluids, 2(5):054604, 2017.
[7] Massimo Germano, Ugo Piomelli, Parviz Moin, and William H Cabot. A dynamic subgrid-scale eddy
viscosity model. Physics of Fluids A: Fluid Dynamics, 3(7):1760–1765, 1991.

5

[8] Jonathan Pietarila Graham and Todd Ringler. A framework for the evaluation of turbulence closures usedin mesoscale ocean large-eddy simulations. Ocean Modelling, 65:25–39, 2013.
[9] Yifei Guan, Ashesh Chattopadhyay, Adam Subel, and Pedram Hassanzadeh. Stable a posteriori les of 2dturbulence using convolutional neural networks: Backscattering analysis and generalization to higher re
via transfer learning. arXiv preprint arXiv:2102.11400, 2021.
[10] Philipp Holl, Nils Thuerey, and Vladlen Koltun. Learning to control pdes with differentiable physics. In
International Conference on Learning Representations, 2019.
[11] Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer.
Machine learning–accelerated computational fluid dynamics. Proceedings of the National Academy of
Sciences, 118(21), 2021.
[12] Andrew Majda and Xiaoming Wang. Nonlinear dynamics and statistical theories for basic geophysicalflows. Cambridge University Press, 2006.
[13] Romit Maulik, Omer San, Adil Rasheed, and Prakash Vedula. Subgrid modelling for two-dimensional
turbulence using neural networks. Journal of Fluid Mechanics, 858:122–144, 2019.
[14] Charles Meneveau and Joseph Katz. Scale-invariance and turbulence models for large-eddy simulation.
Annual Review of Fluid Mechanics, 32(1):1–32, 2000.
[15] Jonghwan Park and Haecheon Choi. Toward neural-network-based large eddy simulation: application toturbulent channel flow. Journal of Fluid Mechanics, 914, 2021.
[16] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen,
Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deeplearning library. In Advances in Neural Information Processing Systems, volume 32, pages 8026–8037,
2019.
[17] Stephen B Pope. Turbulent flows, 2001.
[18] Joseph Smagorinsky. General circulation experiments with the primitive equations: I. the basic experiment.
Monthly Weather Review, 91(3):99–164, 1963.
[19] Yin Zhou and Oncel Tuzel. Voxelnet: End-to-end learning for point cloud based 3d object detection. In
Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4490–4499, 2018.

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes] Yes, thoroughly.
(b) Did you describe the limitations of your work? [Yes] Yes, see Conclusion.
(c) Did you discuss any potential negative societal impacts of your work? [N/A](d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes] Yes, no potential negative societal impacts have been identified related tothis paper.
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A](b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] The code isavailable at https://github.com/hrkz/torchqg.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they werechosen)? [Yes] Yes, we did describe the training details in the "Learning components"
part of Section 3.
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [No] No
6(d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [No] No, but runtime is not an importantfactor for this paper.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] Yes, we used PyTorch.
(b) Did you mention the license of the assets? [Yes] We cited the paper, which then refers
to the licence.
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A](d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]7