# Frame invariant neural network closures for Kraichnan turbulence

Suraj Pawar  School of Mechanical & Aerospace Engineering, Oklahoma State
University, Stillwater, OK 74078, USA.  Omer San<EMAIL_ADDRESS>School of
Mechanical & Aerospace Engineering, Oklahoma State University, Stillwater, OK
74078, USA.  Adil Rasheed  Department of Engineering Cybernetics, Norwegian
University of Science and Technology, N-7465, Trondheim, Norway.  Prakash
Vedula  School of Aerospace & Mechanical Engineering, The University of
Oklahoma, Norman, OK 73019, USA.

###### Abstract

Numerical simulations of geophysical and atmospheric flows have to rely on
parameterizations of subgrid scale processes due to their limited spatial
resolution. Despite substantial progress in developing parameterization (or
closure) models for subgrid scale (SGS) processes using physical insights and
mathematical approximations, they remain imperfect and can lead to inaccurate
predictions. In recent years, machine learning has been successful in
extracting complex patterns from high-resolution spatio-temporal data, leading
to improved parameterization models, and ultimately better coarse grid
prediction. However, the inability to satisfy known physics and poor
generalization hinders the application of these models for real-world
problems. In this work, we propose a frame invariant closure approach to
improve the accuracy and generalizability of deep learning-based subgrid scale
closure models by embedding physical symmetries directly into the structure of
the neural network. Specifically, we utilized specialized layers within the
convolutional neural network in such a way that desired constraints are
theoretically guaranteed without the need for any regularization terms. We
demonstrate our framework for a two-dimensional decaying turbulence test case
mostly characterized by the forward enstrophy cascade. We show that our frame
invariant SGS model (i) accurately predicts the subgrid scale source term,
(ii) respects the physical symmetries such as translation, Galilean, and
rotation invariance, and (iii) is numerically stable when implemented in
coarse-grid simulation with generalization to different initial conditions and
Reynolds number. This work builds a bridge between extensive physics-based
theories and data-driven modeling paradigms, and thus represents a promising
step towards the development of physically consistent data-driven turbulence
closure models.

Neural network, Subgrid scale closure model, Equivariant network, Galilean
invariance, Kraichnan turbulence

##  I Introduction

Computational modeling of turbulent flows remains a key issue in many
engineering and geophysical applications. Turbulence involves a wide range of
spatio-temporal scales that makes the direct numerical simulation (DNS)
computationally infeasible for many complex systems. Coarse-graining
approaches like large eddy simulation (LES) alleviate the computational burden
by resolving only large and intermediate scales of the flow. The non-linearity
of the Navier-Stokes equations introduces a subgrid scale (SGS) closure
problem in LES which can be addressed via modeling of unresolved scales on the
resolved flow quantities. The choice of the SGS model directly affects the
accuracy of LES-based solution, and, therefore, the SGS modeling has been an
active area of research for the past few decades Pope and Pope (2000); Sagaut
(2006); Berselli  _et al._ (2006); Meneveau and Katz (2000). The development
of SGS models has largely been driven by physical insights, mathematical
considerations, and often problem-specific intuition. More recently, the
availability of data from observations and high-resolution simulation along
with advances in hardware and algorithms has fuelled interest in the
development of data-driven turbulence models Duraisamy  _et al._ (2019); Beck
and Kurz (2021); Brunton  _et al._ (2019); Brenner  _et al._ (2019).

The initial efforts towards data-driven SGS modeling include training a neural
network to predict computationally expensive SGS model for channel flow
Sarghini  _et al._ (2003) with the aim to speed-up LES. Similar frameworks
includes applying neural network to determine the eddy-viscosity of the
dynamic Smagorinsky model Pal (2019); Pawar  _et al._ (2020), SGS model
classification and blending Maulik  _et al._ (2019a), and reinforcement
learning to predict SGS dissipation coefficient Novati  _et al._ (2021). Deep
learning (DL) has been applied to discovering new SGS models from the DNS data
without any assumption of prior structural or functional form of the model
Gamahara and Hattori (2017); Maulik  _et al._ (2019b); Xie  _et al._ (2019);
Wang  _et al._ (2021a); Rasp  _et al._ (2018); Gagne  _et al._ (2020); Pawar
_et al._ (2021). The data-driven approach that employs convolutional neural
network for learning the SGS model has also been used for different problems
like two-dimensional decaying turbulence Pawar  _et al._ (2020); Pawar and San
(2021); Guan  _et al._ (2021), three-dimensional decaying homogeneous
isotropic turbulence Beck  _et al._ (2019), momentum forcing in ocean models
Bolton and Zanna (2019), and subgrid-scale scalar flux modeling Frezat  _et
al._ (2021). Moreover, neural networks have also been utilized to learn the
optimal map between filtered and unfiltered variables in the approximate
deconvolution framework for SGS modeling Maulik  _et al._ (2018); Yuan  _et
al._ (2020). Apart from SGS closure modeling, machine learning (ML) and in
particular DL is being increasingly applied for different problems in fluid
mechanics, like superresolution of turbulent flows Fukami  _et al._ (2019);
Kim  _et al._ (2021), Reynolds-Average Navier-Stokes (RANS) closure modeling
Wu  _et al._ (2018); Parish and Duraisamy (2016); Srivastava and Duraisamy
(2021), and reduced-order modeling Xiao  _et al._ (2019); Murata  _et al._
(2020); Kaptanoglu  _et al._ (2021); Ahmed  _et al._ (2021).

Despite their early success, ML models are faced with an array of challenges,
such as poor generalization, lack of interpretability, and in some cases,
violation of the known governing laws of the physical systems. For example,
the SGS model derived through supervised learning may be numerically unstable,
and diverge from the original trajectory, and this issue is exposed in many
studies Maulik  _et al._ (2019b); Wu  _et al._ (2019); Beck  _et al._ (2019);
Nadiga and Livescu (2007). These issues can be addressed by leveraging our
prior knowledge about the physical systems into an ML model. Readers are
referred to recent review articles on physics-informed machine learning
Karniadakis  _et al._ (2021); Kashinath  _et al._ (2021) that detail different
methods of incorporating physics into ML models and discuss the capabilities
and limitations of these methods. In the context of SGS modeling, there are
many ways to embed physical constraints into the ML model. One such method for
constructing a robust and generalizable SGS model is through the selection of
suitable non-dimensionalized input and output quantities of the ML model to
ensure that the known symmetries are respected Prakash  _et al._ (2021).
Another class of methods pertains to the customized neural network
architectures that encode the prior physical or mathematical knowledge as hard
constraints. Some of the examples of this methods applied in fluid dynamics
are tensor basis neural network Ling  _et al._ (2016), transformation
invariant neural network Frezat  _et al._ (2021), physics-embedded neural
network Mohan  _et al._ (2020), spatial transformer Chattopadhyay  _et al._
(2021); Kneer  _et al._ (2021), and equivariant networks Wang  _et al._
(2021b); Siddani  _et al._ (2021).

In this work, we address the challenges associated with data-driven turbulence
modeling by introducing a frame invariant convolutional neural network (FI-
CNN) for SGS closure model discovery. Specifically, we select model inputs
that are Galilean invariant, and replace the convolution operation with group
convolutions Cohen and Welling (2016a); Weiler  _et al._ (2018a) to embed
rotation invariance. Therefore, the FI-CNN preserves various symmetries,
including translation, Galilean, and rotation both during training and
inference. This makes the FI-CNN framework physically consistent and robust to
extrapolation, and consequently, it produces accurate and stable results in
their a posteriori deployment. We demonstrate our framework for two-
dimensional turbulence which is often used as a prototypical test case for
large-scale geophysical flows Boffetta and Ecke (2012); Bouchet and Venaille
(2012). Although we focus on SGS closure model development in this study, this
framework has a promising application for many scientific problems where
physical symmetries are very common. For example, there are several invariant
finite-difference schemes based on equivariant moving frames that preserve Lie
symmetries that have been developed for the solution of partial differential
equations (PDEs) via consideration of modified forms of the underlying PDEs
Ozbenli and Vedula (2020, 2017a, 2017b). These symmetries can be exploited
along with data-driven discretization Bar-Sinai  _et al._ (2019) to design
numerical schemes that are more accurate than their non-invariant counterpart.

This paper is organized as follows. In Sec II, the symmetries of Navier-Stokes
equations and the SGS closure modeling problem for two-dimensional turbulence
is introduced. The detailed procedure on how to embed frame symmetries,
including translation, Galilean, and rotation invariance within the CNN is
provided in Sec III. In Sec IV, the details on data generation and training
are discussed. The performance of the FI-CNN in the a priori and a posteriori
settings along with a detailed discussion of the results are presented in Sec
V. Finally, the concluding remarks and summary of the work are given in Sec
VI.

##  II SGS Closure Modeling

###  II.1 Symmetries of Navier-Stokes Equations

The Navier-Stokes equations governing incompressible fluid flows can be
written in primitive variable (velocity–pressure) form as

$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0,$ (1)
$\displaystyle\frac{\partial\bm{u}}{\partial t}+\bm{u}\cdot\nabla\bm{u}$
$\displaystyle=-\frac{1}{\rho}\nabla p+\nu\nabla^{2}\bm{u},$ (2)

where $\bm{u}$ is the velocity, $p$ is the pressure, $\rho$ is the density,
and $\nu$ is the kinematic viscosity of the fluid. The governing equations for
LES (also called as the filtered Navier-Stokes equations) are obtained by
applying a spatial filter operation and it can be written as follows

$\displaystyle\nabla\cdot\bar{\bm{u}}$ $\displaystyle=0,$ (3)
$\displaystyle\frac{\partial\bar{\bm{u}}}{\partial
t}+\bar{\bm{u}}\cdot\nabla\bar{\bm{u}}$
$\displaystyle=-\frac{1}{\rho}\nabla\bar{p}+\nu\nabla^{2}\bar{\bm{u}}+\nabla\cdot\underbrace{(\bar{\bm{u}}\bar{\bm{u}}-\mkern
1.5mu\overline{\mkern-1.5mu\bm{u}\bm{u}\mkern-1.5mu}\mkern
1.5mu)}_{\bm{\tau}(\bm{u},\bm{u})},$ (4)

where the overbar is used to denote the filtered variables and
$\bm{\tau}(\bm{u},\bm{u})$ is subgrid-scale stress tensor. The problem of
determining subgrid-scale stress tensor $\bm{\tau}$ using the filtered
variables is called the subgrid scale closure problem in LES.

There are many possible SGS closure models and any mathematical, physical
constraints will lead to a specific type of SGS model. Requirements such as
frame invariance, realizability, finite kinetic energy can act as guiding
principles for a satisfactory SGS closure model, and readers are referred to
Berselli  _et al._ (2006) for more details. The frame invariance constraint on
the SGS model is derived by enforcing the symmetry of the original Navier-
Stokes equations Sagaut (2006); Pope and Pope (2000); Oberlack (1997);
Berselli  _et al._ (2006) upon the filtered Navier-Stokes equations with the
SGS closure model. Let $G$ denote a group of transformation acting on space-
time functions $\bm{u}(\bm{x},t)$. We say that the group $G$ is a symmetry
group of the Navier-Stokes equations if, for all $\bm{u}$ which are solutions
of the Navier-Stokes equations, and all $g\in G$, the function $g\bm{u}$ is
also a solution Frisch and Kolmogorov (1995). The frame invariance constraint
involves preservation of the symmetry property of the original Navier-Stokes
equations to translation, Galilean, and rotation transformations, and they can
be written as follows

  * •

Space-translation:
$g_{\bm{\delta}}^{\text{space}}\bm{u}(\bm{x},t)=\bm{u}(\bm{x}-\bm{\delta},t),~{}\forall\bm{\delta}\in\mathbb{R}^{3}$,
where $g_{\bm{\delta}}^{\text{space}}$ is the space-translation operator with
the arbitrary displacement $\bm{\delta}$.

  * •

Galilean transformation:
$g_{\bm{\alpha}}^{\text{Gal}}\bm{u}(\bm{x},t)=\bm{u}(\bm{x}-\bm{\alpha}t,t)+\bm{\alpha},~{}\forall{\bm{\alpha}}\in\mathbb{R}^{3}$,
where $g_{\bm{\alpha}}^{\text{Gal}}$ is the Galilean operator and
$\bm{\alpha}$ is a fixed but arbitrary constant vector.

  * •

Space-rotations:
$g_{\bm{A}}^{\text{rot}}\bm{u}(\bm{x},t)=\bm{A}\bm{u}(\bm{A}^{-1}\bm{x},t)$,
where $g_{\bm{A}}^{\text{rot}}$ is the rotation operation and $\bm{A}\in
SO(3)$.

Imposing the symmetry preservation constraint give some structure to the the
SGS model, and this insights have been extensively used in turbulence models
Speziale (1985); Sagaut (2006); Berselli  _et al._ (2006). We exploit these
symmetries as physical constraints while building a frame invariant data-
driven SGS model.

###  II.2 Two-dimensional Turbulence

In this work, we are interested in the SGS modeling for two-dimensional
turbulence that is usually applied for modeling geophysical flows in the
atmosphere and ocean Bouchet and Venaille (2012); Boffetta and Ecke (2012)
where rotation and stratification dominate, and the most efficient way to
model it is using the vorticity transport equation. Taking the curl of Eq. 2
yields the Navier-Stokes equations in vorticity-velocity formulation, and, for
incompressible fluid flows, it can be written as follows

$\displaystyle\frac{\partial\omega}{\partial t}+(\bm{u}\cdot\nabla)\omega$
$\displaystyle=\nu\nabla^{2}\omega,$ (5)

where $\omega$ is the vorticity, and for two-dimensional flows, we have
$\omega={\partial v}/{\partial x}-{\partial u}/{\partial y}$. A scalar
function called the streamfunction is defined in such a way that the
continuity equation is satisfied if the velocity expressed in terms of the
streamfunction is substituted in the continuity equation. This leads to the
definition of velocity in terms of the streamfunction as follows

$u=\frac{\partial\psi}{\partial y},\quad v=-\frac{\partial\psi}{\partial x},$
(6)

where $\psi$ is the streamfunction. The Poisson equation relating the
vorticity and streamfunction is obtained by substituting the above velocity
components in the definition of vorticity. Thus, we have

$\nabla^{2}\psi=-\omega.$ (7)

It is convenient to write Eq. 5 in the vorticity-streamfunction formulation as
follows

$\displaystyle{\frac{\partial\omega}{\partial
t}+J(\omega,\psi)=\frac{1}{\text{Re}}\nabla^{2}\omega},$ (8)
$\displaystyle{J(\omega,\psi)=\frac{\partial\omega}{\partial
x}\frac{\partial\psi}{\partial y}-\frac{\partial\omega}{\partial
y}\frac{\partial\psi}{\partial x},}$ (9)

where $J(\cdot,\cdot)$ is the Jacobian (or the nonlinear term), and Re is the
Reynolds number of the flow. The above equation is also called the vorticity
transport equation.

The filtered Navier-Stokes equations for two-dimensional turbulence is
obtained by applying a spatial filtering operation to Eq 8 as follows

$\frac{\partial\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}{\partial t}+\mkern
1.5mu\overline{\mkern-1.5muJ({\omega},{\psi})\mkern-1.5mu}\mkern
1.5mu=\frac{1}{\text{Re}}\nabla^{2}\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern 1.5mu.$ (10)

The above equation can be rewritten as

$\frac{\partial\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu}{\partial t}+{J(\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern 1.5mu,\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu)}=\frac{1}{\text{Re}}\nabla^{2}\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern 1.5mu+\Pi,$ (11)

where the overbar quantities represent filtered variables and are evolved on a
grid that is significantly coarse compared to the DNS resolution. The effect
of the unresolved scales due to truncation of high wavenumber flow scales is
encapsulated in a subgrid scale (SGS) source term $\Pi$ and must be modeled
solely based on the resolved variables $(\mkern
1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern 1.5mu,\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu)$. Mathematically,
the true SGS source term $\Pi$ can be expressed as

${\Pi}=J(\mkern 1.5mu\overline{\mkern-1.5mu\omega\mkern-1.5mu}\mkern
1.5mu,\mkern 1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu)-\mkern
1.5mu\overline{\mkern-1.5muJ(\omega,\psi)\mkern-1.5mu}\mkern 1.5mu.$ (12)

The functional and structural models are the most commonly used approaches for
modeling the SGS closure term in LES of turbulent flows Sagaut (2006). The
functional models are based on the concept of eddy viscosity where the effect
of unresolved scales are approximated by artificial dissipation Smagorinsky
(1963); Leith (1971). The functional models can be further improved by dynamic
adaptations of the coefficients that control the dissipation of the model and
are determined adaptively by the use of a low-pass spatial test filter Germano
_et al._ (1991); Lilly (1992); Frederiksen and Kepert (2006). Although the
dynamic formulation allows for spatial and temporal variation of coefficients
in the eddy viscosity model, the ensemble averaging procedure does not allow
for true back-scattering in order to limit the growth of numerical
instabilities during the a posteriori testing Kirkil  _et al._ (2012); Iliescu
and Fischer (2004). The structural models on the other hand aim at obtaining
an accurate approximation of the SGS term and are based on the approximate
deconvolution procedure Stolz  _et al._ (2001); San  _et al._ (2013) and
scale-similarity arguments Bardina  _et al._ (1980). Scale-similarity models
address the SGS closure term by extrapolation from the smallest resolved
scales to unresolved scales and have found to be the most accurate in a priori
testing Sagaut (2006); Sarghini  _et al._ (1999). However, numerical
instabilities have been reported with scale-similarity models, and this has
led to development of many mixed models with additional eddy viscosity term
for stability reasons Layton and Lewandowski (2003); Liu  _et al._ (1994);
Maulik and San (2017). More recently, data-driven methods are emerging as a
new paradigm to build turbulence closure models by extracting information from
the data, and are seen as the potential applications to address the
limitations of existing SGS models Duraisamy  _et al._ (2019); Beck and Kurz
(2021); Duraisamy (2021).

##  III Frame invariant SGS closure model

In this work, we consider the frame invariance property that must be satisfied
by any SGS model and demonstrate how to include them within a neural network
as hard constraints. The SGS source term $\Pi$ is approximated using a neural
network as shown below

$\widetilde{\Pi}\approx\mathcal{M}(\bar{\omega},\bar{\psi}),$ (13)

where $\mathcal{M}$ is a neural network-based model, and $\widetilde{\Pi}$ is
the approximation of true SGS source term $\Pi$. We remark here that the
vorticity is defined using the spatial derivative of the velocity field, and,
therefore it is invariant to Galilean transformations. Additionally, the
streamfunction is computed using the vorticity, and therefore both the inputs
to our model are Galilean invariant. We now discuss how to embed the
translation and rotational invariance/symmetry properties into the neural
network-based model.

###  III.1 Translation invariance

In this work, we employ the convolutional neural network (CNN) for learning
the SGS closure model based on filtered vorticity and streamfunction as the
model inputs. The CNN is an attractive choice for high-dimensional data and it
does not suffer from the curse of dimensionality due to its weight-sharing
feature. The CNN is composed of many convolutional layers and each of the
layers is parameterized by filters, also called kernels, that has to be
learned through training. Let
$\underaccent{\bar}{f},\underaccent{\bar}{k}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{N_{c}}$
be vector-valued two-dimensional features and kernel, i.e.,
$\underaccent{\bar}{f}=(f_{1},\cdots,f_{N_{c}})$ and
$\underaccent{\bar}{k}=(k_{1},\cdots,k_{N_{c}})$, then the convolutional
operation can be defined as

$(\underaccent{\bar}{k}\star\underaccent{\bar}{f})(\mathbf{x})=\sum_{c=1}^{N_{c}}\int_{\mathbb{R}^{2}}k_{c}(\mathbf{x}-\mathbf{x^{\prime}})f_{c}(\mathbf{x^{\prime}})d\mathbf{x^{\prime}},$
(14)

where $\mathbf{x^{\prime}}$ is a dummy variable spanning over $\mathbb{R}^{2}$
space. The convolutional layer maps a feature vector
$\underaccent{\bar}{f}^{(l-1)}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{N_{l-1}}$
with $N_{l-1}$ channels to feature vector
$\underaccent{\bar}{f}^{(l)}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{N_{l}}$
using a set of $N_{l}$ kernels
$\mathbf{k}^{(l)}:=(\underaccent{\bar}{k}_{1}^{(l)},\cdots,\underaccent{\bar}{k}_{N_{l}}^{(l)})$
and this operation can be defined as

$\underaccent{\bar}{f}^{(l)}=\zeta(\mathbf{k}^{(l)}\star\underaccent{\bar}{f}^{(l-1)}):=\zeta(\underaccent{\bar}{k}_{1}^{(l)}\star\underaccent{\bar}{f}^{(l-1)},\cdots,\underaccent{\bar}{k}_{N_{l}}^{(l)}\star\underaccent{\bar}{f}^{(l-1)}),$
(15)

where $\zeta$ is an activation function. The parameters of the kernel are
shared for the whole image as the kernel is convolved relative to the position
about $\mathbf{x}$ and this aspect of the relative motion makes the CNN
translation invariant. Although we present the convolution operation with
continuous kernels, convolutional layers are equipped with discretized-
filtering operations in their practical implementation. From here on, we refer
to the model build using convolutional layers and nonlinear activation
function as $\mathcal{M}_{\text{CNN}}$. The inputs to our model are the
vorticity and streamfunction and the output is the SGS source term. Therefore,
the learning map for $\mathcal{M}_{\text{CNN}}$ can be expressed as

$\mathcal{M}_{\text{CNN}}:\\{\bar{\omega},\bar{\psi}\\}\in\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}\mapsto\\{\widetilde{\Pi}\\}\in\mathbb{R}^{2}\rightarrow\mathbb{R}^{1},$
(16)

where $\widetilde{\Pi}$ is the predicted SGS source term.

###  III.2 Rotation invariance

The rotational invariance of the SGS model requires that it maps as a tensor
under the coordinate rotation Sagaut (2006). As discussed in Section III.1,
the CNN is often invariant to only translation and not for other groups of
transformations. However, there are recent developments on this front to
exploit polar mapping of input images to convert rotation to translation Kim
_et al._ (2020). In this work, we apply the group equivariant convolutions
within the _E(2)-CNN_ framework Weiler and Cesa (2019) for embedding
rotational symmetry. The first roto-translation equivariant CNN was called the
group convolutional neural network (GCNN) and it considered the rotations by
multiples of $\pi/2$ Cohen and Welling (2016a). The GCNN was further augmented
by defining filters in terms of the steerable basis that are equivariant to
rotations by multiples of $2\pi/N$, with $N>4$ Weiler  _et al._ (2018b). The
_E(2)-CNN_ library is based on the framework of steerable CNNs Cohen and
Welling (2016b); Weiler  _et al._ (2018a) and it has different options for the
group that takes the form of the semi-direct group $H=\mathbb{R}^{2}\rtimes G$
where the group $G\leq O(2)$ (here $O(2)$ is the group of continuous rotations
and reflections). For example, the group $H=\mathbb{R}^{2}\rtimes SO(2)=SE(2)$
is the semi-direct product of the group of planar translations
$\mathbb{R}^{2}$ and continuous rotations $SO(2)$. In this work, we utilize
the cyclic group $G=C_{N}$ containing the discrete rotations of $2\pi/N$
(i.e., $H=\mathbb{R}^{2}\rtimes C_{N}$). For a large value of $N$, the
difference between continuous rotations and discrete rotations is
indistinguishable due to space discretization.

A full understanding of the steerable CNNs requires some knowledge of the
group representation theory, but the implementation of the steerable CNNs is
similar to ordinary CNNs. Readers are suggested to read Weiler and Cesa Weiler
and Cesa (2019) and references therein for a more comprehensive discussion on
the general framework of steerable CNNs. Here, we briefly explain the
$G$-equivariant convolutions. A $G$-convolution between a vector-valued two-
dimensional image
$\underaccent{\bar}{f}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{N_{c}}$ and a
filter $\underaccent{\bar}{k}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{N_{c}}$
where $\underaccent{\bar}{f}=(f_{1},\cdots,f_{N_{c}})$ and
$\underaccent{\bar}{k}=(k_{1},\cdots,k_{N_{c}})$ can be expressed as follows

$(\underaccent{\bar}{k}\tilde{\star}\underaccent{\bar}{f})(g)=\sum_{c=1}^{N_{c}}\int_{\mathbb{R}^{2}}k_{c}(g^{-1}\mathbf{x^{\prime}})f_{c}(\mathbf{x^{\prime}})d\mathbf{x^{\prime}},$
(17)

where $g=(\mathbf{x},\theta)\in H=\mathbb{R}^{2}\rtimes C_{N}$,
$\mathbf{x^{\prime}}\in\mathbb{R}^{2}$, and $\tilde{\star}$ denotes the group
correlation operation under joint translation and rotation. This operation
corresponds to lifting of the data on two-dimensional space to the data that
lives on a three-dimensional position orientation space $H$. The first layer
maps a two-dimensional image
$\underaccent{\bar}{f}^{(l-1)}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{N_{l-1}}$
with $N_{l-1}$ channels at $(l-1)$th layer to $H$ vector image
$\underaccent{\bar}{F}^{(l)}:H\rightarrow\mathbb{R}^{N_{l}}$ using a set of
$N_{l}$ kernels
$\mathbf{k}^{(l)}:=(\underaccent{\bar}{k}_{1}^{(l)},\cdots,\underaccent{\bar}{k}_{N_{l}}^{(l)})$
as follows

$\underaccent{\bar}{F}^{(l)}=\zeta(\mathbf{k}^{(l)}\tilde{\star}\underaccent{\bar}{f}^{(l-1)}):=\zeta(\underaccent{\bar}{k}_{1}^{(l)}\tilde{\star}\underaccent{\bar}{f}^{(l-1)},\cdots,\underaccent{\bar}{k}_{N_{l}}^{(l)}\tilde{\star}\underaccent{\bar}{f}^{(l-1)}).$
(18)

Since the $\underaccent{\bar}{F}$ is a function on $H$, the filters from the
second layer onward should also be functions on $H$. The subsequent group
convolutions are defined as Weiler  _et al._ (2018b); Bekkers  _et al._ (2018)

$(\underaccent{\bar}{K}\tilde{\star}\underaccent{\bar}{F})(g)=\sum_{c=1}^{N_{c}}\int_{H}K_{c}(g^{-1}h)F_{c}(h)dh.$
(19)

A group convolution layer is defined by a set of $H$ kernels
$\mathbf{K}:=(\underaccent{\bar}{K}_{1}^{(l)},\cdots\underaccent{\bar}{K}_{N_{l}}^{(l)})$
that maps $\underaccent{\bar}{F}^{(l-1)}$ with $N_{l-1}$ channels to
$\underaccent{\bar}{F}^{(l)}$ with $N_{l}$ channels as shown below

$\underaccent{\bar}{F}^{(l)}=\zeta(\mathbf{K}^{(l)}\tilde{\star}\underaccent{\bar}{F}^{(l-1)}):=\zeta(\underaccent{\bar}{K}_{1}^{(l)}\tilde{\star}\underaccent{\bar}{F}^{(l-1)},\cdots,\\\
\underaccent{\bar}{K}_{N_{l}}^{(l)}\tilde{\star}\underaccent{\bar}{F}^{(l-1)}).$
(20)

Finally, the feature field at the last layer can be synthesized from $H$ space
to $\mathbb{R}^{2}$ space. The user interface of the _E(2)-CNN_ library Weiler
and Cesa (2019) hides most of the intricacies of group theory, solutions of
the steerable kernels space constraints, and requires users to specify only
the transformation laws of the feature spaces. We use the regular
representation for all hidden layers and the action of regular representation
is given by permutation matrices (Appendix B in Weiler and Cesa (2019)). From
here on, the model built using the equivariant CNN is called as
$\mathcal{M}_{\text{FI-CNN}}$. The learning map for $\mathcal{M}_{\text{FI-
CNN}}$ is same as the $\mathcal{M}_{\text{CNN}}$ given in Eq. 16.

##  IV Data Generation and Training

The parameters of the neural network based SGS models are learned through
supervised training that requires a set of labeled inputs and outputs, usually
obtained from direct numerical simulation (DNS). The dataset should encompass
a range of dynamics that is expected to be reproduced by the SGS model. The
data for training is generated from DNS of two-dimensional Kraichnan
turbulence in a doubly periodic square domain with $L_{x}\times
L_{y}=[0,2\pi]\times[0,2\pi]$, and the domain is discretized using $2048^{2}$
degrees of freedom. Our DNS solver is based on a second-order accurate energy-
conserving Arakawa scheme Arakawa (1997) for the nonlinear Jacobian and
second-order accurate finite-difference scheme for the Laplacian of the
vorticity. The elliptic equation for the relationship between the
streamfunction and vorticity is solved using a second-order accurate FFT-based
Poisson solver, and the time integration is performed with a third-order
accurate Runge-Kutta method. The vorticity distribution at the start of the
simulation is initialized based on the energy spectrum given by Orlandi (2000)

$E(k)=Ak^{4}\text{exp}\bigg{(}-\bigg{(}\frac{k}{k_{p}}\bigg{)}^{2}\bigg{)},$
(21)

where $A=4k_{p}^{-5}/3\pi$ and $k=|\mathbf{k}|=\sqrt{k_{x}^{2}+k_{y}^{2}}$.
For our numerical experiments, we use $k_{p}=10$. The initial vorticity
distribution in Fourier space is obtained through the introduction of random
phase as follows

$\tilde{\omega}(\mathbf{k})=\sqrt{\frac{k}{\pi}E(k)}~{}\text{e}^{\mathbf{i}\xi(\mathbf{k})},$
(22)

where the phase function is given by
$\xi(\mathbf{k})=\phi(\mathbf{k})+\eta(\mathbf{k})$. Here, $\phi(\mathbf{k})$
and $\eta(\mathbf{k})$ are independent random values chosen in $[0,2\pi]$ at
each grid point in the first quadrant of the $k_{x}-k_{y}$ plane (i.e.,
$k_{x},k_{y}\geq 0$). The phase function for other quadrants is obtained
through conjugate relations as follows

$\displaystyle\xi(\mathbf{k})$
$\displaystyle=-\phi(\mathbf{k})+\eta(\mathbf{k})~{}\text{for}~{}k_{x}<0~{}\text{and}~{}k_{y}\geq
0,$ (23) $\displaystyle\xi(\mathbf{k})$
$\displaystyle=-\phi(\mathbf{k})-\eta(\mathbf{k})~{}\text{for}~{}k_{x}<0~{}\text{and}~{}k_{y}<0,$
(24) $\displaystyle\xi(\mathbf{k})$
$\displaystyle=\phi(\mathbf{k})-\eta(\mathbf{k})~{}\text{for}~{}k_{x}\geq
0~{}\text{and}~{}k_{y}<0,$ (25)

Further details on the problem setup and the numerical schemes can be found in
our previous work Maulik and San (2017). Different realizations of the initial
vorticity field can be obtained by using different phase functions with a
different seed for random value generation.

The DNS is performed from time $t=0$ to $t=4$ with the time step $\Delta
t=5\times 10^{-4}$. In the Kraichnan turbulence problem, the initial vorticity
field is dominated by a population of vortices and small-scale structure
starts appearing as the flow evolves. The initial spin-up time from $t=0$ to
$0.5$ is neglected and we start collecting the data for training from time
$t=0.5$. From time $t\approx 0.5$, the flow has started following
Kraichnan–Batchelor–Leith (KBL) theory Kraichnan (1967); Batchelor (1969);
Leith (1971) of energy cascade where energy is transferred from the smaller
scales to the larger scales. From time $t\approx 0.5$ onward, large coherent
vortices start emerging through vortex merging mechanism and viscous
dissipation of small-scale structures. The vorticity field and angle-averaged
energy spectrum are displayed in Fig. 1 and we can see that the energy
spectrum has started exhibiting $k^{-3}$ scaling from approximately $t=0.5$.

Figure 1: Visualization of the vorticity field and energy spectrum at
different time instances for Re = 16000 with grid resolution $2048\times
2048$.

The data for training a neural network-based SGS model is stored at every
$20\Delta t$, i.e., we have 350 snapshots of the vorticity and streamfunction
between $t=0.5$ to $4.0$. We emphasize here that the neural network-based SGS
model is trained only for a single Reynolds number Re = 16000 and we assess
the performance of the model for Reynolds number up to Re = 128000. The
filtered DNS data for training is obtained by first applying a Gaussian filter
transfer function to the DNS data and then coarse-graining the filtered
solution to the LES grid Zanna and Bolton (2020); Guan  _et al._ (2021). The
Gaussian filter provides a smooth transition between resolved and subgrid
scales and is also positive definite in physical and wave space Piomelli  _et
al._ (1991); Rogallo and Moin (1984). Additionally, our numerical solver is in
physical space, and therefore we select the Gaussian filter instead of a
spectral cut-off filter. The coarse-grid level for LES is $256^{2}$ which
corresponds to 64 times fewer spatial degrees of freedom compared to DNS.

We do not pre-process the filtered DNS data before training as the DNS data is
generated from a non-dimensionalized vorticity transport equation. The total
data is divided into 80% of the data for training and 20% for the validation
set. While the input and output of both $\mathcal{M}_{\text{CNN}}$ and
$\mathcal{M}_{\text{FI-CNN}}$ are the same, the user needs to specify the type
of representation for intermediate feature field while constructing an FI-CNN
Weiler and Cesa (2019), similar to the number of channels for CNN. We use the
kernel size of $5\times 5$, six hidden layers and ReLU activation function for
the $\mathcal{M}_{\text{CNN}}$ and $\mathcal{M}_{\text{FI-CNN}}$. The number
of channels for the CNN and FI-CNN models is set to 30 and 16, respectively.
With these hyperparameters, the number of trainable parameters is roughly the
same around $O(1.1\times 10^{5})$ for both models. Both the models are trained
for 100 iterations using an Adam optimizer. Fig. 2 shows the history of
training loss versus iterations for both neural network-based SGS models and
we can observe that the loss for $\mathcal{M}_{\text{FI-CNN}}$ is almost one
order magnitude less than the loss for $\mathcal{M}_{\text{CNN}}$. This can be
attributed to rotational invariances incorporated in the
$\mathcal{M}_{\text{FI-CNN}}$ against $\mathcal{M}_{\text{CNN}}$, which is
only invariant to translation and Galilean transformation. For both neural
network-based SGS models, we use the parameters (i.e., weights) corresponding
to minimum validation loss obtained while training the neural network.

Figure 2: History of the training loss versus iterations for
$\mathcal{M}_{\text{CNN}}$ and $\mathcal{M}_{\text{FI-CNN}}$.

##  V Numerical Results

We first outline the numerical results of our framework in the a priori
settings where the neural network-based models are utilized in predicting the
SGS source term. We analyze the capability of $\mathcal{M}_{\text{CNN}}$ and
$\mathcal{M}_{\text{FI-CNN}}$ in incorporating the frame-invariance property
over the testing data. Then, we present the results of a posteriori LES
coupled with neural network-based SGS models and evaluate their performance
using numerous statistical metrics.

###  V.1 A Priori investigation

Here, we assess the performance of neural network-based models in predicting
the SGS source term compared to the true SGS source term for the out-of-
training data. The out-of-training data is obtained for a different initial
condition and corresponds to 70 snapshots stored randomly between time $t=0.5$
to $t=4.0$. We remark here that the initial energy spectrum for the testing
data is also given by Eq. 21 and the difference is due to a different phase
function. Fig. 3 displays the probability distribution function and cumulative
distribution function for the test data. There is a very good agreement
between the true SGS source term and the predicted SGS source term from both
models. However, we notice that the $\mathcal{M}_{\text{FI-CNN}}$ is more
accurate near the tails of the distribution (Fig 3, left) compared to
$\mathcal{M}_{\text{CNN}}$. This difference is also observable in the
cumulative distribution function of true and predicted SGS source terms and is
highlighted in the zoom-in portion (Fig 3, right). Based on these results, we
may conclude that both neural network-based SGS model has learned the
relationship between filtered quantities and the SGS source term. Both models
are able to produce viable physical results for the completely unseen data
with similar physics.

Figure 3: Probability distribution function (left) and cumulative distribution
function (right) of the SGS source term over the entire testing dataset. The
testing dataset corresponds to 70 snapshots selected randomly between time
$t=0.5$ to $t=4.0$ for the initial condition different from the one used in
training.

Next, we evaluate neural network-based models in respecting rotational
symmetry on the test data. Specifically, we perturbed the test data based on
the rotation transformation, and generate multiple test datasets. Then, we
compute the root mean squared error (RMSE) for each dataset, and calculate the
expected value and variance for all test datasets. The rotation matrix
$\mathbf{A}$ is sampled uniformly between $[0^{\circ},360^{\circ}]$ in the
multiple of $90^{\circ}$. If the rotation symmetry constraint is satisfied
strictly, then the RMSE will be the same for each test dataset leading to zero
variance for the RMSE metric. The results in Table 1 demonstrate the
remarkable ability of $\mathcal{M}_{\text{FI-CNN}}$ to respect the rotation
symmetry in contrast to $\mathcal{M}_{\text{CNN}}$ which violates this
symmetry. Furthermore, the RMSE for $\mathcal{M}_{\text{FI-CNN}}$ is one order
of magnitude lower than $\mathcal{M}_{\text{CNN}}$ and is consistent with the
training loss (Fig. 2). The performance of both models is substantially
accurate in terms of the Pearson’s cross-correlation coefficient, with
$\mathcal{M}_{\text{FI-CNN}}$ slightly better than $\mathcal{M}_{\text{CNN}}$.
We note here that it is relatively straightforward to embed Galilean
invariance constraint within neural network-based SGS model through
intelligent selection of model inputs, and translation invariance through
simple CNN. However, incorporating rotational symmetry in a neural network-
based SGS model is more complex and requires special consideration. Although a
relatively simple method like data augmentation can be utilized to impose the
rotation symmetry as a soft constraint, it does not satisfy rotation
invariance strictly Frezat  _et al._ (2021). The strict enforcement of
rotation symmetry is challenging and requires the use of tailored neural
network architecture, such as equivariant CNN.

Table 1: Evaluation of the rotational symmetry constraints provided by $\mathcal{M}_{\text{CNN}}$ and $\mathcal{M}_{\text{FI-CNN}}$. The expected value and variance of the root mean squared error on the SGS source term predicted by both models is computed from many realizations (20 ensembles) on the testing data. The testing dataset corresponds to 70 snapshots selected randomly for the initial condition different from the one used for training. The rotational angle $\mathbf{A}$ is sampled uniformly between $[0^{\circ},360^{\circ}]$ in the multiple of $90^{\circ}$ and is used in the rotational operator $g_{\mathbf{A}}^{\text{rot}}$. The Pearson’s cross-correlation coefficient between the predicted SGS source term and the filtered DNS solution is computed as $\mathcal{P}(X,Y)=\text{cov}(X,Y)/\sigma_{X}\sigma_{Y}$. Metric | $\mathcal{M}_{\text{CNN}}$ | $\mathcal{M}_{\text{FI-CNN}}$
---|---|---
$\text{E}[\mathcal{L}]$ |  10.6941 |  7.3462
$\sigma[\mathcal{L}]$ |  4.2442 $\times 10^{-2}$ |  5.6587 $\times 10^{-8}$
$\mathcal{P}(X,Y)$ |  0.9600 |  0.9776

###  V.2 A Posteriori deployment

We now evaluate the performance of neural network-based SGS models in the LES
of Kraichnan turbulence. The spatial resolution for LES is reduced by a factor
of eight in each direction and this gives us $256^{2}$ degrees of freedom. The
time step for LES simulation is ten times larger compared to the DNS, i.e.,
$\Delta t_{\text{LES}}=5\times 10^{-3}$. The performance of neural network-
based SGS models is compared with the widely used dynamic Smagorinsky model
(DSM) Germano  _et al._ (1991); Lilly (1992). The a posteriori deployment is a
rigorous task for any data-driven SGS model due to the presence of numerical
instabilities, and the challenges and remedies have been highlighted in many
studies Maulik  _et al._ (2019b, a, 2020); Guan  _et al._ (2021); Beck  _et
al._ (2019); Stoffer  _et al._ (2021); Zhou  _et al._ (2019). For example,
Maulik  _et al._ Maulik  _et al._ (2019b) and Zhou  _et al._ Zhou  _et al._
(2019) achieved the stable LES results by truncating SGS source term
corresponding to negative eddy viscosity. Stoffer  _et al._ Stoffer  _et al._
(2021) attained stable a posteriori results by resorting to artificially
introducing additional dissipation (via eddy-viscosity models). Guan  _et al._
Guan  _et al._ (2021) provided sufficient amount of data during training to
obtain a stable a posteriori results. While the exact reason for this behavior
is unknown, several issues such as error accumulation, aliasing errors,
numerical instability, extrapolation beyond the training data, chaotic nature
of turbulence, presence of multiple attractors might be responsible for
unstable a posteriori simulation Beck  _et al._ (2019); Stoffer  _et al._
(2021); Beck and Kurz (2021); Nadiga and Livescu (2007).

From our a posteriori simulation, it is revealed that
$\mathcal{M}_{\text{CNN}}$ is unstable, while $\mathcal{M}_{\text{FI-CNN}}$ is
able to produce a stable and physical solution without any kind of clipping or
by adding artificial dissipation. We note here that, perhaps
$\mathcal{M}_{\text{CNN}}$ can also achieve stable a posteriori simulation,
provided there is sufficient data available for training or some kind of post-
processing is carried out for the predicted SGS source term. However, our main
motivation in this work is to construct a physically consistent data-driven
SGS model that can be trained in a data-sparse regime and is also stable in
the a posteriori simulation. We assess the performance of our a posteriori
simulation using several statistical metrics and compare it with the
statistics from filtered DNS solution. The turbulent kinetic energy at time
$t_{k}$ is computed as follows

$TKE=\lambda(u_{f}^{2}+v_{f}^{2}(t_{k})),$ (27)

where $u_{f}$ and $v_{f}$ are the fluctuating components of velocity given by

$\displaystyle u_{f}=\bar{u}-\lambda(\bar{u}),$ (28) $\displaystyle
v_{f}=\bar{v}-\lambda(\bar{v}),$ (29)

where $\lambda(a)$ represents the spatial average of the field variable $a$.
The velocity $\bar{u}$, and $\bar{v}$ are computed by spectral differentiation
of streamfunction as shown below

$\bar{u}=\frac{\partial\bar{\psi}}{\partial
y};\quad\bar{v}=-\frac{\partial\bar{\psi}}{\partial x}.$ (30)

The vorticity variance at each time step is computed as

$\sigma^{2}=\lambda((\bar{\omega}-\lambda(\bar{\omega}))^{2}).$ (31)

We compare the kinetic-energy spectra and the vorticity structure function at
intermediate time $t=2.0$ and at final time $t=4.0$ with the $k^{-3}$ scaling
which is observed in two-dimensional turbulence based on the classical KBL
theory. The vorticity structure function is calculated using the formula given
by Grossmann and Mertens (1992) for two-dimensional turbulence as follows

$S_{\omega}=<|\bar{\omega}(\mathbf{x}+\mathbf{r})-\bar{\omega}(\mathbf{x})|^{2}>,$
(32)

where $<>$ indicates ensemble averaging, $\mathbf{x}$ is the position on the
grid, and $\mathbf{r}$ is certain distance from this location. The PDF of the
vorticity increment is utilized to assess the capability of SGS models in
predicting the coherent vortices in the flow. The vorticity increments at
different separations $\mathbf{r}$ is defined as

$\delta\omega(\mathbf{r})=\omega(\mathbf{x}+\mathbf{r})-\omega(\mathbf{r}).$
(33)

We reiterate here that neural network-based SGS models are trained using the
data for Reynolds number Re = 16000 and a single initial condition. Once the
models are trained, the LES coupled with SGS models is performed for Reynolds
number up to Re = 128000 and for five different initial conditions. Fig. 4
shows the evolution of turbulent kinetic energy and vorticity variance for LES
runs with five different initial conditions and for several Reynolds numbers.
For all the LES runs, we initialize the vorticity field at $t=0.5$ after the
initial spin-up period using the filtered DNS solution. We can observe that
the model $\mathcal{M}_{\text{CNN}}$ is stable only for short time and quickly
becomes unstable after $t\approx 2.0$ even for Reynolds number Re = 16000
which was included in the training. In contrast to $\mathcal{M}_{\text{CNN}}$,
model $\mathcal{M}_{\text{FI-CNN}}$ is stable for all test cases conducted
here without any post-processing of the predicted SGS source term. The
ensemble averaging procedure in DSM leads to highly dissipative results and is
noticeable in the overprediction of the energy decay rate. The results of the
LES runs with $\mathcal{M}_{\text{FI-CNN}}$ have the best agreement with
filtered DNS solution for both turbulent kinetic energy and the vorticity
variance.

Figure 4: The time evolution of the turbulent kinetic energy $TKE(t_{k})$
normalized by the initial turbulent kinetic energy $TKE(t_{0})$ (top row) and
vorticity variance (bottom row) for different Reynolds numbers at $256^{2}$
grid resolution. The solid line shows the mean from LES runs for five
different initial conditions and the shaded area corresponds to one standard
deviation. The LES simulation starts at $t=0.5$ after the initial spin-up time
(i.e., once the turbulence has set in). The CNN and FI-CNN models are trained
using the data generated from a single initial condition at Reynolds number Re
= 16000.

Fig. 5 displays the kinetic-energy spectra at intermediate time $t=2.0$ and at
final time $t=4.0$ obtained from LES runs with five different initial
conditions for multiple Reynolds number. Although the LES runs coupled with
$\mathcal{M}_{\text{CNN}}$ is stable at $t=2.0$, the solution is unphysical as
seen by the energy pile up near grid cutoff wavenumbers. This behavior is also
illustrated in Fig. 6 through a large value of vorticity structure function at
$t=2.0$ across all Reynolds numbers. The LES runs with
$\mathcal{M}_{\text{CNN}}$ has diverged around $t\approx 2.5$ (see Fig. 4),
and, therefore the kinetic-energy spectra and vorticity structure function are
missing at $t=4.0$ in Fig. 5 and Fig. 6, respectively. There is a very good
agreement between the kinetic-energy spectra for LES runs with
$\mathcal{M}_{\text{FI-CNN}}$ and filtered DNS solution, especially in the
inertial subrange and $k^{-3}$ theoretical scaling is captured accurately.
From Fig. 6, we can see that the model $\mathcal{M}_{\text{FI-CNN}}$ is
successful in producing the $r^{3/2}$ scaling Kramer  _et al._ (2011) for the
vorticity structure function at small scales and it gradually flattens near
the large scales. The excessive dissipation of DSM is also illustrated in Fig.
5 and Fig. 6 via mismatch between kinetic-energy spectra and vorticity
structure function between DSM and filtered DNS solution. The successful
performance of LES runs with $\mathcal{M}_{\text{FI-CNN}}$ demonstrates that
incorporating frame symmetries as hard constraints has been effective in
stabilizing the coarse-grid simulation and in ensuring generalized learning
across different initial conditions and Reynolds numbers.

Figure 5: A posteriori kinetic-energy spectra for different Reynolds numbers
at $t=2.0$ (top row) and $t=4.0$ (bottom row). These results are obtained from
LES runs with five different initial conditions and only mean kinetic energy
spectrum is shown. Note here that the CNN model has diverged and the kinetic-
energy spectra for the CNN model is missing at the final time $t=4.0$ (bottom
row).

Figure 6: A posteriori second-order vorticity structure for different Reynolds
numbers at $t=2.0$ (top row) and $t=4.0$ (bottom row). These results are
obtained from LES runs with five different initial conditions and the solid
line shows the mean vorticity structure and the shaded area corresponds to one
standard deviation. Note here that the CNN model has diverged and the
vorticity structure for the CNN model is not present at the final time $t=4.0$
(bottom row).

Figs. 7-10 provides the visualization of vorticity field and probability
density function (PDF) of vorticity increments for Reynolds number Re = 16000
to Re = 128000 computed using the filtered DNS solution, LES with DSM model,
and LES with $\mathcal{M}_{\text{FI-CNN}}$ at final time $t=4.0$. We remark
here that these results correspond to only one initial condition that is
different from the one used for training the neural network-based SGS models.
Even though the LES with DSM model is successful in capturing large-scale
structures in the flow, it fails to capture the small-scale structure due to
excessive dissipation. The LES with $\mathcal{M}_{\text{FI-CNN}}$ is able to
capture both large- and small-scale structures in the flow, and this can be
ascertained to the stabilizing property of $\mathcal{M}_{\text{FI-CNN}}$ in
the a posteriori deployment without any post-processing of the predicted SGS
source term. Qualitatively, the vorticity field obtained from LES with model
$\mathcal{M}_{\text{FI-CNN}}$ is very similar to the filtered DNS solution.
The similarity in the shape of the PDF of vorticity increments as shown in
Figs. 7-10 suggests the scale-invariant statistics of turbulence at all
Reynolds numbers investigated in this study. The shape of the PDF of vorticity
increments predicted by the LES with $\mathcal{M}_{\text{FI-CNN}}$ matches
with the shape of the filtered DNS solution, and the heavy exponential tails
in the PDF are related to the presence of coherent vortices in the flow. These
heavy exponential tails are missing in the PDF of the solution obtained from
LES with DSM, and it follows the Gaussian distribution.

Figure 7: Snapshots of the vorticity distribution for Re = 16000 taken at
final time $t=4.0$ (top row) from different models and compared qualitatively
against the FDNS solution. The bottom row displays the probability density
function $P(\delta\omega)$ of the vorticity increments $\delta\omega$ for
separations
$r=2\pi/256,2\pi/128,2\pi/64,2\pi/32,2\pi/16,2\pi/8,\text{and}2\pi/4$ computed
from different models at Re = 16000. A Gaussian distribution is given in gray
dashed line for comparison. We also note that the plain vanilla CNN becomes
numerically unstable and unbounded before $t=4.0$.

Figure 8: Snapshots of the vorticity distribution for Re = 32000 taken at
final time $t=4.0$ (top row) from different models and compared qualitatively
against the FDNS solution. The bottom row displays the probability density
function $P(\delta\omega)$ of the vorticity increments $\delta\omega$ for
separations
$r=2\pi/256,2\pi/128,2\pi/64,2\pi/32,2\pi/16,2\pi/8,\text{and}2\pi/4$ computed
from different models at Re = 32000. A Gaussian distribution is given in gray
dashed line for comparison. We also note that the plain vanilla CNN becomes
numerically unstable and unbounded before $t=4.0$.

Figure 9: Snapshots of the vorticity distribution for Re = 64000 taken at
final time $t=4.0$ (top row) from different models and compared qualitatively
against the FDNS solution. The bottom row displays the probability density
function $P(\delta\omega)$ of the vorticity increments $\delta\omega$ for
separations
$r=2\pi/256,2\pi/128,2\pi/64,2\pi/32,2\pi/16,2\pi/8,\text{and}2\pi/4$ computed
from different models at Re = 64000. A Gaussian distribution is given in gray
dashed line for comparison. We also note that the plain vanilla CNN becomes
numerically unstable and unbounded before $t=4.0$.

Figure 10: Snapshots of the vorticity distribution for Re = 128000 taken at
final time $t=4.0$ (top row) from different models and compared qualitatively
against the FDNS solution. The bottom row displays the probability density
function $P(\delta\omega)$ of the vorticity increments $\delta\omega$ for
separations
$r=2\pi/256,2\pi/128,2\pi/64,2\pi/32,2\pi/16,2\pi/8,\text{and}2\pi/4$ computed
from different models at Re = 128000. A Gaussian distribution is given in gray
dashed line for comparison. We also note that the plain vanilla CNN becomes
numerically unstable and unbounded before $t=4.0$.

Next, we examine the robustness of the neural network-based SGS model by
training an ensemble of neural networks using randomization-based approaches
where different random initialization of weights are utilized for generating
ensembles. Specifically, we train five neural networks for both models
$\mathcal{M}_{\text{CNN}}$ and $\mathcal{M}_{\text{FI-CNN}}$ using the same
dataset as discussed in Section IV. This method is also applied to quantify
the model-form uncertainty in deep learning Lakshminarayanan  _et al._ (2017).
Fig. 11 shows the time evolution of turbulent kinetic energy and vorticity
variance at different Reynolds numbers. The time evolution of the TKE in Fig.
11 implies that the weights of the neural networks for $\mathcal{M}_{\text{FI-
CNN}}$ are learned in such a way that the final models are overall dissipative
in nature (as indicated by the solid line for mean from different LES runs).
The vorticity variance predicted by the model $\mathcal{M}_{\text{FI-CNN}}$ is
more accurate compared to DSM and is very close to the filtered DNS solution.

Figure 11: The time evolution of the turbulent kinetic energy $TKE(t_{k})$
normalized by the initial turbulent kinetic energy $TKE(t_{0})$ (top row) and
vorticity variance (bottom row) for different Reynolds numbers at $256^{2}$
grid resolution. The LES simulation starts at $t=0.5$ after the initial spin-
up time (i.e., once the turbulence has set in). For the CNN and FI-CNN models,
an ensemble of neural networks is trained using the data generated from a
single initial condition at Reynolds number Re = 16000 with different weights
initialization. The solid line shows the mean from LES runs for a single
initial condition with different trained networks and the shaded area
corresponds to one standard deviation.

Fig. 12 depicts the kinetic-energy spectra at intermediate time $t=2.0$ and at
final time $t=4.0$ obtained from LES runs for a single initial condition with
different network-based SGS models for multiple Reynolds numbers. We observe
the energy pile up near grid cutoff wavenumbers for the LES runs coupled with
$\mathcal{M}_{\text{CNN}}$ and this suggests that the solution is unphysical.
This behavior is also demonstrated in Fig. 13 through a large value of
vorticity structure function at $t=2.0$ across all Reynolds numbers. The LES
runs with $\mathcal{M}_{\text{CNN}}$ has diverged around $t\approx 2.5$ (as
seen by large $TKE$ in Fig. 11), and, therefore the kinetic-energy spectra and
vorticity structure function are missing at $t=4.0$ in Fig. 12 and Fig. 13,
respectively. The kinetic-energy spectra for LES runs with
$\mathcal{M}_{\text{FI-CNN}}$ is highly accurate and shows an excellent
agreement with the filtered DNS solution, especially in the inertial subrange.
The small uncertainty band also suggests that an ensemble of neural networks
have produced very similar statistics for $\mathcal{M}_{\text{FI-CNN}}$. Fig.
13 shows that the model $\mathcal{M}_{\text{FI-CNN}}$ is successful in
capturing the $r^{3/2}$ scaling for the vorticity structure function at small
scales and flattening near large scales. With this numerical experiment, we
can establish that the $\mathcal{M}_{\text{FI-CNN}}$ is robust, trustworthy,
and stable in the LES, and it also ensures generalizable learning across
different initial conditions and Reynolds numbers.

Figure 12: A posteriori kinetic-energy spectra for different Reynolds numbers
at $t=2.0$ (top row) and $t=4.0$ (bottom row). The solid line shows the mean
from LES runs for a single initial condition with different trained networks
and the shaded area corresponds to one standard deviation. Note here that the
CNN model has already diverged and the kinetic-energy spectra for the CNN
model is missing at the final time $t=4.0$ (bottom row).

Figure 13: A posteriori second-order vorticity structure for different
Reynolds numbers at $t=2.0$ (top row) and $t=4.0$ (bottom row). The solid line
shows the mean from LES runs for a single initial condition with different
trained networks and the shaded area corresponds to one standard deviation.
Note here that the CNN model has diverged and the vorticity structure for the
CNN model is not present at the final time $t=4.0$ (bottom row).

##  VI Concluding remarks

Closure modeling in fluid dynamics simulations refers to parameterizing the
interactions between high-fidelity and coarse-fidelity descriptions. In this
study, we explore data-driven closure modeling strategies to improve both the
accuracy and generalizability of such residual models. The motivation behind
data-driven closure modeling stems from the fact that most of the existing SGS
models are derived based on physical and mathematical considerations, and
might not account for the important transfer of kinetic energy from small
scales to large scales (i.e., back-scatter) Hewitt  _et al._ (2020). We
introduce a frame invariant neural network architecture aiming at embedding
physical symmetries directly into the structure of the convolutional neural
networks. Thus, our model theoretically guarantees the frame symmetries,
including translation, Galilean, and rotation invariance both during training
and inference. The embedding of physical symmetries as hard constraints not
only improves the accuracy of the model but notably improves the
generalization of the model, and eventually makes the model stable in their a
posteriori deployment without any clipping.

We test the proposed framework for subgrid-scale modeling of Kraichnan
turbulence in a priori and a posteriori settings. The performance of the
proposed framework is evaluated using several metrics like kinetic energy
spectra, vorticity structure, and vorticity increments. Based on our analysis,
we concluded that symmetry preservation has the potential to improve the
accuracy, generalizability, and stability of the SGS model, besides embedding
important geometric properties of the underlying PDEs into deep learning
models. This work also illustrates a broader lesson on how to combine machine
learning with physics for scientific computing. It may be argued that two-
dimensional turbulence is far from reality. However, it is generally
considered as a canonical testbed for geophysical turbulence in the atmosphere
and oceans. Our future development will be focused on scaling up the proposed
frame invariant closure modeling framework to solve more realistic three-
dimensional turbulent flows. Another interesting avenue is to apply this
framework for learning parameterization models for geophysical flows, paving
the way for improved weather and climate prediction.

## Data availability

The data that supports the findings of this study are available within the
article. Implementation details and Python scripts can be accessed from the
Github repository Pawar (2021).

###### Acknowledgements.

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research under
Award Number DE-SC0019290. O.S. gratefully acknowledges their support. Direct
numerical simulations and neural network training for this project was carried
out using resources of the Oklahoma State University High-Performance
Computing Center. Disclaimer: This report was prepared as an account of work
sponsored by an agency of the United States Government. Neither the United
States Government nor any agency thereof, nor any of their employees, makes
any warranty, express or implied, or assumes any legal liability or
responsibility for the accuracy, completeness, or usefulness of any
information, apparatus, product, or process disclosed, or represents that its
use would not infringe privately owned rights. Reference herein to any
specific commercial product, process, or service by trade name, trademark,
manufacturer, or otherwise does not necessarily constitute or imply its
endorsement, recommendation, or favoring by the United States Government or
any agency thereof. The views and opinions of authors expressed herein do not
necessarily state or reflect those of the United States Government or any
agency thereof.

## References

  * Pope and Pope (2000) S. B. Pope and S. B. Pope,  _Turbulent flows_ (Cambridge university press, 2000).
  * Sagaut (2006) P. Sagaut,  _Large eddy simulation for incompressible flows: an introduction_ (Springer Science & Business Media, 2006).
  * Berselli  _et al._ (2006) L. C. Berselli, T. Iliescu, and W. J. Layton,  _Mathematics of Large Eddy Simulation of Turbulent Flows_ (Scientific Computation, Springer-Verlag, Berlin, 2006).
  * Meneveau and Katz (2000) C. Meneveau and J. Katz, Annual Review of Fluid Mechanics 32, 1 (2000).
  * Duraisamy  _et al._ (2019) K. Duraisamy, G. Iaccarino, and H. Xiao, Annual Review of Fluid Mechanics 51, 357 (2019).
  * Beck and Kurz (2021) A. Beck and M. Kurz, GAMM-Mitteilungen 44, e202100002 (2021).
  * Brunton  _et al._ (2019) S. L. Brunton, B. R. Noack, and P. Koumoutsakos, Annual Review of Fluid Mechanics 52 (2019).
  * Brenner  _et al._ (2019) M. Brenner, J. Eldredge, and J. Freund, Physical Review Fluids 4, 100501 (2019).
  * Sarghini  _et al._ (2003) F. Sarghini, G. De Felice, and S. Santini, Computers & Fluids 32, 97 (2003).
  * Pal (2019) A. Pal, arXiv preprint arXiv:1905.12765 (2019).
  * Pawar  _et al._ (2020) S. Pawar, O. San, A. Rasheed, and P. Vedula, Theoretical and Computational Fluid Dynamics , 387 (2020).
  * Maulik  _et al._ (2019a) R. Maulik, O. San, J. D. Jacob, and C. Crick, Journal of Fluid Mechanics 870, 784 (2019a).
  * Novati  _et al._ (2021) G. Novati, H. L. de Laroussilhe, and P. Koumoutsakos, Nature Machine Intelligence 3, 87 (2021).
  * Gamahara and Hattori (2017) M. Gamahara and Y. Hattori, Physical Review Fluids 2, 054604 (2017).
  * Maulik  _et al._ (2019b) R. Maulik, O. San, A. Rasheed, and P. Vedula, Journal of Fluid Mechanics 858, 122 (2019b).
  * Xie  _et al._ (2019) C. Xie, J. Wang, K. Li, and C. Ma, Physical Review E 99, 053113 (2019).
  * Wang  _et al._ (2021a) Y. Wang, Z. Yuan, C. Xie, and J. Wang, AIP Advances 11, 055216 (2021a).
  * Rasp  _et al._ (2018) S. Rasp, M. S. Pritchard, and P. Gentine, Proceedings of the National Academy of Sciences 115, 9684 (2018).
  * Gagne  _et al._ (2020) D. J. Gagne, H. M. Christensen, A. C. Subramanian, and A. H. Monahan, Journal of Advances in Modeling Earth Systems 12, e2019MS001896 (2020).
  * Pawar  _et al._ (2021) S. Pawar, O. San, A. Rasheed, and I. M. Navon, GEM - International Journal on Geomathematics (2021).
  * Pawar and San (2021) S. Pawar and O. San, Physical Review Fluids 6, 050501 (2021).
  * Guan  _et al._ (2021) Y. Guan, A. Chattopadhyay, A. Subel, and P. Hassanzadeh, arXiv preprint arXiv:2102.11400 (2021).
  * Beck  _et al._ (2019) A. Beck, D. Flad, and C.-D. Munz, Journal of Computational Physics 398, 108910 (2019).
  * Bolton and Zanna (2019) T. Bolton and L. Zanna, Journal of Advances in Modeling Earth Systems 11, 376 (2019).
  * Frezat  _et al._ (2021) H. Frezat, G. Balarac, J. Le Sommer, R. Fablet, and R. Lguensat, Physical Review Fluids 6, 024607 (2021).
  * Maulik  _et al._ (2018) R. Maulik, O. San, A. Rasheed, and P. Vedula, Physics of Fluids 30, 125109 (2018).
  * Yuan  _et al._ (2020) Z. Yuan, C. Xie, and J. Wang, Physics of Fluids 32, 115106 (2020).
  * Fukami  _et al._ (2019) K. Fukami, K. Fukagata, and K. Taira, Journal of Fluid Mechanics 870, 106 (2019).
  * Kim  _et al._ (2021) H. Kim, J. Kim, S. Won, and C. Lee, Journal of Fluid Mechanics 910 (2021).
  * Wu  _et al._ (2018) J.-L. Wu, H. Xiao, and E. Paterson, Physical Review Fluids 3, 074602 (2018).
  * Parish and Duraisamy (2016) E. J. Parish and K. Duraisamy, Journal of Computational Physics 305, 758 (2016).
  * Srivastava and Duraisamy (2021) V. Srivastava and K. Duraisamy, arXiv preprint arXiv:2103.16042 (2021).
  * Xiao  _et al._ (2019) D. Xiao, F. Fang, C. Heaney, I. Navon, and C. Pain, Computer Methods in Applied Mechanics and Engineering 354, 307 (2019).
  * Murata  _et al._ (2020) T. Murata, K. Fukami, and K. Fukagata, Journal of Fluid Mechanics 882 (2020).
  * Kaptanoglu  _et al._ (2021) A. A. Kaptanoglu, K. D. Morgan, C. J. Hansen, and S. L. Brunton, Physical Review E 104, 015206 (2021).
  * Ahmed  _et al._ (2021) S. E. Ahmed, S. Pawar, O. San, A. Rasheed, T. Iliescu, and B. R. Noack, Physics of Fluids 33, 091301 (2021).
  * Wu  _et al._ (2019) J. Wu, H. Xiao, R. Sun, and Q. Wang, Journal of Fluid Mechanics 869, 553 (2019).
  * Nadiga and Livescu (2007) B. Nadiga and D. Livescu, Physical Review E 75, 046303 (2007).
  * Karniadakis  _et al._ (2021) G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Nature Reviews Physics 3, 422 (2021).
  * Kashinath  _et al._ (2021) K. Kashinath, M. Mustafa, A. Albert, J. Wu, C. Jiang, S. Esmaeilzadeh, K. Azizzadenesheli, R. Wang, A. Chattopadhyay, A. Singh, _et al._ , Philosophical Transactions of the Royal Society A 379, 20200093 (2021).
  * Prakash  _et al._ (2021) A. Prakash, K. E. Jansen, and J. A. Evans, arXiv preprint arXiv:2106.13410 (2021).
  * Ling  _et al._ (2016) J. Ling, A. Kurzawski, and J. Templeton, Journal of Fluid Mechanics 807, 155 (2016).
  * Mohan  _et al._ (2020) A. T. Mohan, N. Lubbers, D. Livescu, and M. Chertkov, arXiv preprint arXiv:2002.00021 (2020).
  * Chattopadhyay  _et al._ (2021) A. Chattopadhyay, M. Mustafa, P. Hassanzadeh, E. Bach, and K. Kashinath, Geoscientific Model Development Discussions , 1 (2021).
  * Kneer  _et al._ (2021) S. Kneer, T. Sayadi, D. Sipp, P. Schmid, and G. Rigas, arXiv preprint arXiv:2111.02893 (2021).
  * Wang  _et al._ (2021b) R. Wang, R. Walters, and R. Yu, in  _International Conference on Learning Representations_ (2021).
  * Siddani  _et al._ (2021) B. Siddani, S. Balachandar, and R. Fang, Physics of Fluids 33, 103323 (2021).
  * Cohen and Welling (2016a) T. Cohen and M. Welling, in  _International conference on machine learning_ (PMLR, 2016) pp. 2990–2999.
  * Weiler  _et al._ (2018a) M. Weiler, M. Geiger, M. Welling, W. Boomsma, and T. Cohen, arXiv preprint arXiv:1807.02547 (2018a).
  * Boffetta and Ecke (2012) G. Boffetta and R. E. Ecke, Annual Review of Fluid Mechanics 44, 427 (2012).
  * Bouchet and Venaille (2012) F. Bouchet and A. Venaille, Physics Reports 515, 227 (2012).
  * Ozbenli and Vedula (2020) E. Ozbenli and P. Vedula, Physical Review E 101, 023303 (2020).
  * Ozbenli and Vedula (2017a) E. Ozbenli and P. Vedula, Journal of Computational Physics 349, 376 (2017a).
  * Ozbenli and Vedula (2017b) E. Ozbenli and P. Vedula, Physical Review E 96, 063304 (2017b).
  * Bar-Sinai  _et al._ (2019) Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner, Proceedings of the National Academy of Sciences 116, 15344 (2019).
  * Oberlack (1997) M. Oberlack, Annual Research Briefs , 3 (1997).
  * Frisch and Kolmogorov (1995) U. Frisch and A. N. Kolmogorov,  _Turbulence: the legacy of AN Kolmogorov_ (Cambridge university press, 1995).
  * Speziale (1985) C. G. Speziale, Journal of Fluid Mechanics 156, 55 (1985).
  * Smagorinsky (1963) J. Smagorinsky, Monthly Weather Review 91, 99 (1963).
  * Leith (1971) C. Leith, Journal of the Atmospheric Sciences 28, 145 (1971).
  * Germano  _et al._ (1991) M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, Physics of Fluids A: Fluid Dynamics 3, 1760 (1991).
  * Lilly (1992) D. K. Lilly, Physics of Fluids A: Fluid Dynamics 4, 633 (1992).
  * Frederiksen and Kepert (2006) J. S. Frederiksen and S. M. Kepert, Journal of the Atmospheric Sciences 63, 3006 (2006).
  * Kirkil  _et al._ (2012) G. Kirkil, J. Mirocha, E. Bou-Zeid, F. K. Chow, and B. Kosović, Monthly Weather Review 140, 266 (2012).
  * Iliescu and Fischer (2004) T. Iliescu and P. Fischer, Computers & Fluids 33, 783 (2004).
  * Stolz  _et al._ (2001) S. Stolz, N. A. Adams, and L. Kleiser, Physics of fluids 13, 997 (2001).
  * San  _et al._ (2013) O. San, A. E. Staples, and T. Iliescu, Ocean Modelling 63, 1 (2013).
  * Bardina  _et al._ (1980) J. Bardina, J. Ferziger, and W. Reynolds, in  _13th fluid and plasmadynamics conference_ (1980) p. 1357.
  * Sarghini  _et al._ (1999) F. Sarghini, U. Piomelli, and E. Balaras, Physics of Fluids 11, 1596 (1999).
  * Layton and Lewandowski (2003) W. Layton and R. Lewandowski, Applied Mathematics Letters 16, 1205 (2003).
  * Liu  _et al._ (1994) S. Liu, C. Meneveau, and J. Katz, Journal of Fluid Mechanics 275, 83 (1994).
  * Maulik and San (2017) R. Maulik and O. San, Computers & Fluids 158, 11 (2017).
  * Duraisamy (2021) K. Duraisamy, Physical Review Fluids 6, 050504 (2021).
  * Kim  _et al._ (2020) J. Kim, W. Jung, H. Kim, and J. Lee, arXiv preprint arXiv:2007.10588 (2020).
  * Weiler and Cesa (2019) M. Weiler and G. Cesa, arXiv preprint arXiv:1911.08251 (2019).
  * Weiler  _et al._ (2018b) M. Weiler, F. A. Hamprecht, and M. Storath, in  _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ (2018) pp. 849–858.
  * Cohen and Welling (2016b) T. S. Cohen and M. Welling, arXiv preprint arXiv:1612.08498 (2016b).
  * Bekkers  _et al._ (2018) E. J. Bekkers, M. W. Lafarge, M. Veta, K. A. Eppenhof, J. P. Pluim, and R. Duits, in  _International conference on medical image computing and computer-assisted intervention_ (Springer, 2018) pp. 440–448.
  * Arakawa (1997) A. Arakawa, Journal of Computational Physics 135, 103 (1997).
  * Orlandi (2000) P. Orlandi,  _Fluid flow phenomena: a numerical toolkit_ , Vol. 55 (Springer Science & Business Media, 2000).
  * Kraichnan (1967) R. H. Kraichnan, The Physics of Fluids 10, 1417 (1967).
  * Batchelor (1969) G. K. Batchelor, The Physics of Fluids 12, II (1969).
  * Zanna and Bolton (2020) L. Zanna and T. Bolton, Geophysical Research Letters 47, e2020GL088376 (2020).
  * Piomelli  _et al._ (1991) U. Piomelli, W. H. Cabot, P. Moin, and S. Lee, Physics of Fluids A: Fluid Dynamics 3, 1766 (1991).
  * Rogallo and Moin (1984) R. S. Rogallo and P. Moin, Annual review of fluid mechanics 16, 99 (1984).
  * Maulik  _et al._ (2020) R. Maulik, O. San, and J. D. Jacob, Physica D: Nonlinear Phenomena 406, 132409 (2020).
  * Stoffer  _et al._ (2021) R. Stoffer, C. M. van Leeuwen, D. Podareanu, V. Codreanu, M. A. Veerman, M. Janssens, O. K. Hartogensis, and C. C. van Heerwaarden, Geoscientific Model Development 14, 3769 (2021).
  * Zhou  _et al._ (2019) Z. Zhou, G. He, S. Wang, and G. Jin, Computers & Fluids 195, 104319 (2019).
  * Grossmann and Mertens (1992) S. Grossmann and P. Mertens, Zeitschrift für Physik B Condensed Matter 88, 105 (1992).
  * Kramer  _et al._ (2011) W. Kramer, G. Keetels, H. Clercx, and G. van Heijst, Physical Review E 84, 026310 (2011).
  * Lakshminarayanan  _et al._ (2017) B. Lakshminarayanan, A. Pritzel, and C. Blundell, in  _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17 (Curran Associates Inc., Red Hook, NY, USA, 2017) p. 6405–6416.
  * Hewitt  _et al._ (2020) H. T. Hewitt, M. Roberts, P. Mathiot, A. Biastoch, E. Blockley, E. P. Chassignet, B. Fox-Kemper, P. Hyder, D. P. Marshall, E. Popova, _et al._ , Current Climate Change Reports 6, 137 (2020).
  * Pawar (2021) S. Pawar, “FI-CNN,” https://github.com/surajp92/FI-CNN (2021).