# On closures for reduced order models – A spectrum of first-principle to
machine-learned avenues

Shady E. Ahmed  Suraj Pawar  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.  Traian Iliescu  Department
of Mathematics, Virginia Tech, Blacksburg, VA 24061, USA.  Bernd R. Noack
School of Mechanical Engineering and Automation, Harbin Institute of
Technology, Shenzhen 518058, China.  Hermann-Föttinger-Institut für
Strömungsmechanik, Technische Universität Berlin, D-10623 Berlin, Germany

(August 27, 2024)

###### Abstract

For over a century, reduced order models (ROMs) have been a fundamental
discipline of theoretical fluid mechanics. Early examples include Galerkin
models inspired by the Orr-Sommerfeld stability equation and numerous vortex
models, of which the von Kármán vortex street is one of the most prominent.
Subsequent ROMs typically relied on first principles, like mathematical
Galerkin models, weakly nonlinear stability theory, and two- and three-
dimensional vortex models. Aubry et al. [N. Aubry, P. Holmes, J. Lumley, and
E. Stone, Journal of Fluid Mechanics, 192, 115–173 (1988)] pioneered data-
driven proper orthogonal decomposition (POD) modeling. In early POD modeling,
available data was used to build an optimal basis, which was then utilized in
a classical Galerkin procedure to construct the ROM. But data has made a
profound impact on ROMs beyond the Galerkin expansion. In this paper, we take
a modest step and illustrate the impact of data-driven modeling on one
significant ROM area. Specifically, we focus on ROM closures, which are
correction terms that are added to classical ROMs in order to model the effect
of the discarded ROM modes in under-resolved simulations. Through simple
examples, we illustrate the main modeling principles used to construct
classical ROMs, motivate and introduce modern ROM closures, and show how data-
driven modeling, artificial intelligence, and machine learning have changed
the standard ROM methodology over the last two decades. Finally, we outline
our vision on how state-of-the-art data-driven modeling can continue to
reshape the field of reduced order modeling.

Reduced order modeling, projection methods, closure models, statistical
inference, neural networks

##  I Introduction

One of the very first exciting experiences that kids go through is playing
with water; they might throw a stone in a lake, float a rubber duck in a
bathtub, or even stir a straw while enjoying a tasty cup of juice! They like
doing this over and over again because of the magnificent patterns that keep
forming every time. These patterns or _coherent structures_ are ubiquitous in
the world, in general, and in fluid flows, in particular. They attracted
Leonardo da Vinci more than five centuries ago, resulting in some of his
outstanding artwork Marusic and Broomhall (2021). Fluid dynamicists are
especially lucky to enjoy the beauty of these formations on a daily basis. But
other than their aesthetic value, these patterns come with a practical
benefit. In particular, these coherent structures are the cornerstone in the
development of reduced order models (ROMs) for fluid flows. ROMs are built by
using available data to identify and rank these structures, choosing the most
effective few of them, and tracking their dynamical behavior in order to
approximate the evolution of the underlying flow. The computational cost of
the relatively low-dimensional ROMs is dramatically lower than the
computational cost of a direct numerical simulation, which aims at capturing
all the flow scales. Since their introduction to the field of fluid dynamics
more than fifty years ago Lumley (1967), ROMs have witnessed tremendous
changes. Arguably, data-driven modeling has been the main driving force behind
these changes. Over the last two decades, state-of-the-art methods from
machine learning (ML) have reshaped the field of reduced order modeling.

The main objective of this study is to provide an overview of data-driven
reduced order modeling strategies relevant to fluid dynamics applications. The
topic spans a wide spectrum, and there are many review articles on pertinent
discussions, methodologies, and applications in fluids Kou and Zhang (2021);
Beck and Kurz (2021); Brunton, Noack, and Koumoutsakos (2020); Taira  _et al._
(2020); Mendonça, Afonso, and Lau (2019); Yu, Yan, and Guo (2019); Brenner,
Eldredge, and Freund (2019); Yondo  _et al._ (2019); Taira  _et al._ (2017);
Rowley and Dawson (2017); Brunton and Noack (2015); Lassila  _et al._ (2014);
Mezić (2013); Chinesta, Ladeveze, and Cueto (2011); Massarotti  _et al._
(2010); Lucia, Beran, and Silva (2004); Berkooz, Holmes, and Lumley (1993) as
well as closely related fields San, Rasheed, and Kvamsdal (2021); Heinlein
_et al._ (2021); Blechschmidt and Ernst (2021); Kashinath  _et al._ (2021);
Bauer, Stevens, and Hazeleger (2021); Frank, Drikakis, and Charissis (2020);
Zhu and Piotrowski (2020); Tahmasebi  _et al._ (2020); Benner, Gugercin, and
Willcox (2015); Asher  _et al._ (2015); Machairas, Tsangrassoulis, and Axarli
(2014); Mignolet  _et al._ (2013); Bazaz and Janardhanan (2012); Razavi,
Tolson, and Burn (2012); Theodoropoulos (2011); Wagner  _et al._ (2010);
Kleijnen (2009); Pinnau (2008); Ong, Nair, and Keane (2003); Freund (2003);
Bai (2002); Chatterjee (2000); Bonvin and Mellichamp (1982); Elrazaz and Sinha
(1981). Therefore, it is not our intention to include a detailed discussion,
but rather to survey one important ROM research area, closure modeling, and
provide our subjective perspectives on how data-driven modeling has made an
impact in this area. In particular, given the recent interest in ML
applications in fluid dynamics, our survey is intended to encourage cross-
disciplinary efforts between practitioners, physicists, mathematicians, and
data scientists. We hope that our paper will shed light on the new ideas of
integrating both physics in ML models and ML-enabled capabilities in
principled models, a rapidly emerging field that came to be known as physics-
guided ML (PGML). These developments are born to conform with the scientific
foundations that are moving rapidly to the industry to enable the next
generation of digital twin technologies Rasheed, San, and Kvamsdal (2020).

Figure 1: An overview of data-driven reduced order modeling in fluid dynamics.

This paper first aims at identifying imminent practical and mathematical needs
in designing closure approaches for ROMs of nonlinear parameterized fluid
dynamics systems, i.e., complex natural or engineered systems comprising
coupled partial differential equations (PDEs) with variable parameters, and
initial and boundary conditions Karatzas, Ballarin, and Rozza (2020). Such
models usually serve as the inner-workhorse for outer-workflow loops such as
optimal design Aguado  _et al._ (2017), control Aguado  _et al._ (2015);
Zainib  _et al._ (2019), estimation, and discovery Quarteroni and Rozza
(2014); Benner, Gugercin, and Willcox (2015); Wang  _et al._ (2020); Swischuk
_et al._ (2019); Lee and Carlberg (2020). In particular, there has recently
been an increasing interest in ROMs from the fluid dynamics community, where
emerging data-driven methods prevail. This is primarily due to the fact that
data and centralized powerful open-source machine learning and optimization
libraries have become widespread, as indicated in Figure 1. Although we mainly
focus on incompressible flows, we emphasize that there have been inspiring
works done in the compressible case Iollo, Lanteri, and Désidéri (2000);
Rowley, Colonius, and Murray (2004); Bourguet, Braza, and Dervieux (2007);
Gloerfelt (2008); Bourguet  _et al._ (2009); Barone  _et al._ (2009);
Kalashnikova and Barone (2010); Placzek, Tran, and Ohayon (2011); Balajewicz,
Tezaur, and Dowell (2016); Ferrero, Iollo, and Larocca (2018); Yano (2019);
Yu, Wei-Xi, and Chun-Xiao (2019); Rezaian and Wei (2020); Zucatti, Wolf, and
Bergmann (2021); Krath  _et al._ (2021).

To begin with, reduced order modeling can be viewed as the art of converting
existing prior information and collected data into a dramatically more
efficient, yet relatively accurate, surrogate model to be used on demand. For
example, a conceivable strategy of flow control is to put most of the
demanding calculations offline and to keep only low rank updates for fluid
flow evolutions online Hovland, Gravdahl, and Willcox (2008); Ballarin, Rozza,
and Strazzullo (2020). Emerging digital twin infrastructures are one of the
main beneficiaries and driving forces behind efficient surrogate model
development efforts Kapteyn  _et al._ (2020); Kapteyn, Pretorius, and Willcox
(2021). Although the ROM concept is not new in fluid dynamics, there are still
many new fronts and opportunities, mainly due to the recent advances in ML
algorithms and easy-to-use open-source packages that can be utilized in many
control and optimization processes. We also note that, in many fluid dynamics
applications, the typical data sparsity (due to the number of necessarily
resolved degrees of freedom being orders of magnitude larger than the
available sensors) and corruption (e.g., due to signal noise, interference,
and sensor malfunctioning) motivate physics informed data-driven modeling.

Among fluid dynamicists, projection-based linear methods have become popular.
Both proper orthogonal decomposition (POD) Holmes, Lumley, and Berkooz (1998)
and dynamic mode decomposition (DMD)Kutz  _et al._ (2016) enabled approaches
have been exploited. In our work, we mostly focus on POD relevant literature
and refer the reader to Kutz  _et al._ (2016) for the DMD principles. The
mathematical foundations behind the POD-based linear subspace approaches go
back to the principal component analysis (PCA), pioneered by Karl Pierson
Pearson (1901) in 1901 and later demonstrated graphically by Harold Hotelling
Hotelling (1933) in 1933. This powerful statistical approach (also known as
Kosambi-Karhunen-Loève expansionKosambi (2016); Loeve (1955) or empirical
orthogonal functionsLorenz (1956); Monahan  _et al._ (2009)) was first
introduced in the fluid dynamics community by John Lumley Lumley (1967, 1972,
1981), and came to be known as POD. In practice, the _method of snapshots_ ,
established by Lawrence Sirovich Sirovich (1987), was a key enabler to
efficiently determine the POD modes for large scale problems, often
encountered in fluid dynamics. Of particular interest when characterizing the
dynamics of coherent structures in wall bounded flows, the beauty of the POD
modeling approach was demonstrated in seminal works by Nadine Aubry Aubry  _et
al._ (1988); Aubry (1991). An admittedly incomplete chronological evolution of
projection-based ROMs is given in Table 1.

A key advantage of POD is the guaranteed minimal representation error for the
employed snapshots with respect to all other Galerkin expansions with the same
number of modes. Another advantage is the orthogonality of the constructed
modes that naturally leverages the use of a Galerkin type projection onto the
governing PDEs to obtain a system of ordinary differential equations, defining
a dynamical system for the amplitudes of POD modes. These models are often
_intrusive_ in the sense that both the governing equations and the prerecorded
snapshot data are required to build a ROM. Postulating an alternative
approach, DMD based approaches and many other _nonintrusive_ models bypass
this equation dependency in order to construct the ROM solely based on the
prerecorded snapshots Willcox and Peraire (2002); Barrault  _et al._ (2004a);
Bai (2002); Mezić (2005); Grepl  _et al._ (2007); Gugercin, Antoulas, and
Beattie (2008); Chaturantabut and Sorensen (2010); Schmid (2010); Williams,
Kevrekidis, and Rowley (2015); Carlberg, Bou-Mosleh, and Farhat (2011); Tu
_et al._ (2014); Lassila  _et al._ (2014); Benner, Gugercin, and Willcox
(2015); Brunton and Noack (2015).

Table 1: An incomplete chronological list of key contributions to projection-based ROMs in fluid dynamics. Year |  Study |  Key Contribution
---|---|---
1915 |  Galerkin (1915) |  Galerkin method for solving (initial) boundary value problems
1962 |  Saltzman (1962) |  Low-dimensional modeling (with 7 modes, see also RefLakshmivarahan, Lewis, and Hu (2019) for a revisit)
1963 |  Lorenz (1963) |  Low-dimensional modeling (with 3 modes)
1967 |  Lumley (1967) |  Proper orthogonal decomposition (POD)
1987 |  Sirovich (1987) |  Method of snapshots
1988 |  Aubry  _et al._ (1988) |  First POD model: Dynamics of coherent structures and global eddy viscosity modeling
1994 |  Rempfer and Fasel (1994) |  Linear modal eddy viscosity closure
1995 |  Everson and Sirovich (1995) |  Gappy POD
2000 |  Ravindran (2000) |  Galerkin ROM for optimal flow control problems
2001 |  Kunisch and Volkwein (2001a) |  First numerical analysis of Galerkin ROM for parabolic problems
2002 |  Willcox and Peraire (2002) |  Balanced truncation with POD
2003 |  Couplet, Sagaut, and Basdevant (2003) |  Guidelines for modeling unresolved modes in POD–Galerkin models
2004 |  Sirisup and Karniadakis (2004) |  Spectral viscosity closure for POD models
2004 |  Barrault  _et al._ (2004a) |  Empirical interpolation method (EIM)
2005 |  Mezić (2005) |  Spectral decomposition of the Koopman operator
2007 |  Rozza, Huynh, and Patera (2007) |  Reduced basis approximation
2007 |  Cao  _et al._ (2007) |  Galerkin ROM for four‐dimensional variational data assimilation
2008 |  Amsallem and Farhat (2008) |  Interpolation method based on the Grassmann manifold approach
2008 |  Astrid  _et al._ (2008) |  Missing point estimation
2009 |  Rowley  _et al._ (2009) |  Spectral analysis of nonlinear flows
2009 |  Sapsis and Lermusiaux (2009) |  Dynamically orthogonal field equations
2010 |  Schmid (2010) |  A purely nonintrusive perspective: Dynamic mode decomposition (DMD)
2010 |  Chaturantabut and Sorensen (2010) |  Discrete empirical interpolation method (DEIM)
2013 |  Carlberg  _et al._ (2013) |  The Gauss–Newton with approximated tensors (GNAT) method
2013 |  Cordier  _et al._ (2013) |  Proof of global boundedness of nonlinear eddy viscosity closures
2014 |  Östh  _et al._ (2014) |  $\sqrt{K}$-scaled eddy viscosity concept
2015 |  Ballarin  _et al._ (2015) |  Stabilization of POD Galerkin approximations
2015 |  Schlegel and Noack (2015) |  On bounded solutions of Galerkin models
2016 |  Peherstorfer and Willcox (2016) |  Data-driven operator inference nonintrusive ROMs
2016 |  Brunton, Proctor, and Kutz (2016) |  Sparse identification of nonlinear dynamics (SINDy)
2016 |  Sieber, Paschereit, and Oberleithner (2016) |  Spectral POD
2018 |  Towne, Schmidt, and Colonius (2018) |  On the relationship between spectral POD, DMD, and resolvent analysis
2018 |  Reiss  _et al._ (2018) |  Shifted/transported snapshot POD
2018 |  Loiseau, Noack, and Brunton (2018) |  Feature-based manifold modeling
2019 |  Mendez, Balabane, and Buchlin (2019) |  Multi-scale proper orthogonal decomposition
2021 |  Li  _et al._ (2021), Fernex, Noack, and Semaan (2021) |  Cluster-based network models

If we adopt the training and testing terminologies from ML, both POD and DMD
based models require the _training_ snapshot data to build the ROM (i.e.,
data-driven modeling). The fundamental question in practice is how well these
ROMs will perform in testing conditions (i.e., the conditions that are not
included in the training data). The trade-off performance between training and
testing constitutes one of the crucial questions about the credibility of the
proposed ROMs, and motivates more efforts, ideas, and collaborations to push
the frontiers of existing ROM frameworks. Physics informed ML has also made an
impact in reduced order modeling. In this hybrid approach, while the training
data provides a set of global basis functions, the underlying governing
equations (i.e., physics) constrain the evolution within the linear subspace
defined by these POD basis functions. Deep discussions on hybrid approaches
that combine deterministic and statistical modeling can be found elsewhere
Krasnopolsky and Fox-Rabinovitz (2006a, b); Reichstein  _et al._ (2019); de
Bezenac, Pajot, and Gallinari (2019); Karpatne  _et al._ (2017); Childs and
Washburn (2019); Willard  _et al._ (2020); San, Rasheed, and Kvamsdal (2021).

Projection-based ROMs have been explored for decades, and these explorations
have been paying off in many applications. They have great promise for flow
control of industrial processes, for enlarging ensemble size for flow problems
with uncertain data, and even for providing accurate forecasts of fluid
behaviour. The apparent success of low-rank ensemble nonlinear filtering
methods utilized in weather forecasting centers also suggests that there is a
prospect of using a system whose dimension is substantially lower than the
dimension of the state space. Yet, the ROMs’ potential has been only realized
for a small collection of canonical flows. One of the main roadblocks for ROMs
of realistic flows is that they are not accurate models for the dominant
modes. In practice, a closure or correction term is generally added Kalb and
Deane (2007); Amsallem and Farhat (2012); Balajewicz and Dowell (2012); Wang
_et al._ (2012); Cordier  _et al._ (2013); Östh  _et al._ (2014); Protas,
Noack, and Östh (2015); Benosman  _et al._ (2017); Stabile  _et al._ (2019);
Reyes and Codina (2020). In many cases, there are complementary physical,
statistical, and computational challenges that arise in the development of
ROMs and ROM closures, topics that we will systematically survey in this work
toward establishing foundations to close the gap between what ROMs can do and
where they are needed.

##  II Reduced Order Representation

A _reduced order representation_ (ROR) can be viewed as a generalization of
the latent space or manifold, i.e., a simplifying kinematic approximation. For
example, the POD procedure introduced in Section II.2 constitutes a best-fit
linear manifold to establish a ROR. The ROR facilitates a data compression for
an ensemble of snapshot data. Physically interpretable RORs are possible when
dominant coherent structure are present. This is clearly ubiquitous in the
fluid flows that we encounter in our daily life, as introduced in Section I,
as well as large-scale and industrial settings. For spatio-temporal dynamical
systems, the rank-r ROR of the state $u(\bm{x},t)$ can be simply written as

$u(\bm{x},t)=\sum_{i=1}^{r}a_{i}(t)\psi_{i}(\bm{x}),$ (1)

where $\bm{x}$ refers to the spatial coordinates, $t$ is the time, $\psi_{i}$
denotes the $i$-th mode in the ROR, and $a_{i}$ is the corresponding amplitude
or coefficient. Although it is generally assumed that the basis functions
$\psi_{i}$ are time-independent and the dynamical evolution is encapsulated in
the coefficients $a_{i}$, there have been studies that admit a time-evolving
basis functions as well Sapsis and Lermusiaux (2009); Cheng, Hou, and Zhang
(2013); Ramezanian, Nouri, and Babaee (2021); Patil and Babaee (2020).

###  II.1 Eigenfunction expansion

Any set of $n$ linearly independent vectors can serve as a basis for an
$n$-dimensional vector space. Any vector in this space can be expressed as a
linear combination of these linearly independent basis vectors. In an infinite
dimensional vector space of functions, there exists an infinite set of
linearly independent basis functions $\\{\psi_{i}(\bm{x})\\}_{i=1,2,\dots}$
such that a given function $u(\bm{x})$ in this space can be written as a
linear combination of these functions. It is straightforward to show that any
periodic, piecewise continuous function can be written as an infinite sum of
sines and cosines (e.g., Fourier series Haberman (2012)). The eigenfunction
expansion can be viewed as a generalization of the Fourier series expansion
for arbitrary boundary conditions, where the Sturm-Liouville theory provides
an infinite sequence of eigenvalue-eigenfunction pairs. ROR also seeks an
expansion of an arbitrary function in terms of a given set of basis functions.
However, in contrast to the methods mentioned in this section, ROR aims at
finding a low-dimensional basis instead of an infinite dimensional one.

To illustrate these concepts, let us consider a linear advection problem in a
spatial domain $[0,L]$:

$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0,$ (2)

where $c$ is the wave speed. To select the appropriate basis functions, we
consider the boundary conditions. For example, if we have homogeneous
Dirichlet boundary conditions, i.e.,

$u(x=0,t)=0,\quad u(x=L,t)=0,\quad\mbox{for}\quad t\in[0,T],$ (3)

we might choose a set of _orthonormal_ basis functions $\psi_{i}(x)$ defined
as

$\psi_{i}(x)=\sqrt{\frac{2}{L}}\sin\Big{(}\frac{i\pi}{L}x\Big{)},$ (4)

in order to approximate $u$ as follows:

$u(x,t)=\sum_{i=1}^{r}a_{i}(t)\psi_{i}(x).$ (5)

In a more sophisticated scenario with the homogeneous Neumann condition on the
left boundary ($x=0)$ and the homogeneous Drichlet condition on the right
boundary ($x=L=1$), i.e.,

$\frac{\partial u}{\partial x}\Big{|}_{x=0}=0,\quad
u\Big{|}_{x=1}=0,\quad\mbox{for}\quad t\in[0,T],$ (6)

we can define a set of _non-orthogonal_ functions that satisfy Eq. 6 as
follows:

$\phi_{i}(x)=\cos(i\pi x)-(-1)^{i}$ (7)

and apply the Gram-Schmidt _orthonormalization_ process to obtain the
following set of basis functions:

$\displaystyle\psi_{i}(x)$
$\displaystyle=\sqrt{\frac{4i-2}{2i+1}}\Bigg{[}\frac{(-1)^{i+1}}{2i-1}+\cos(i\pi
x)$ (8) $\displaystyle+\frac{2}{2i-1}\sum_{j=1}^{i-1}(-1)^{i+j+1}\cos(j\pi
x)\Bigg{]}.$

Here, we note that these basis functions are orthonormal, i.e.,

$\int_{0}^{1}\psi_{i}(x)\psi_{k}(x)dx=\delta_{ik},$ (9)

and they are derived from the Fourier harmonics that satisfy the boundary
conditions.

Next, we focus on an illustrative example with periodic boundary conditions,
with the domain length $L=2\pi$, maximum time $T=2\pi$, and wave speed $c=1$.
For a given initial condition $u(x,t=0)=\cos(x)$, Eq. 2 admits an analytical
solution in the form of a right travelling wave $u(x,t)=\cos(x-t)$. Let us
approximate $u$ using a modal expansion with only two Fourier harmonics
defined by

$u(x,t)=a_{1}(t)\cos(x)+a_{2}(t)\sin(x).$ (10)

Substituting Eq. 10 into Eq. 2, we get

$\displaystyle\frac{\partial}{\partial
t}\Big{(}a_{1}(t)\cos(x)+a_{2}(t)\sin(x)\Big{)}$
$\displaystyle+\frac{\partial}{\partial
x}\Big{(}a_{1}(t)\cos(x)+a_{2}(t)\sin(x)\Big{)}=0.$ (11)

Once we multiply Eq. II.1 with $\cos(x)$ and integrate over the domain, we
obtain an equation for $a_{1}(t)$,

$\displaystyle\frac{da_{1}}{dt}=-a_{2},$ (12)

and similarly, multiplying Eq. II.1 with $\sin(x)$, the evolution equation for
$a_{2}$ becomes

$\displaystyle\frac{da_{2}}{dt}=a_{1}.$ (13)

Eq. 12 and Eq. 13 constitute the well-known Galerkin system. Using the initial
condition given at $t=0$,

$a_{1}(0)=1,\quad a_{2}(0)=0,$ (14)

we can obtain an analytical solution of the Galerkin system given by Eq. 12
and Eq. 13 as

$a_{1}(t)=\cos(t),\quad a_{2}(t)=\sin(t).$ (15)

Therefore, the two-mode Galerkin model approximation given by Eq. 10 yields a
solution

$u(x,t)=\cos(t)\cos(x)+\sin(t)\sin(x),$ (16)

which can be further written as

$u(x,t)=\cos(x-t).$ (17)

As we illustrated in this example, the two-mode approximation retrieves the
exact solution. One of the key aspects in such a modal scheme is therefore
related to the characteristics of the selected basis functions, which
ultimately provided the best possible expansion in this example. A central
question is how we would know a priori the appropriate $\cos(x)$ and $\sin(x)$
basis functions to approximate $u$.

We also note that the multimodal method utilizes similar arguments to
represent the solution as a superposition of an infinite set of generalized
Fourier bases and time-dependent coefficients. The multimodal method has been
extensively exploited to study the sloshing problem Faltinsen and Timokha
(2009); Lukovsky and Timokha (2017); Gavrilyuk, Lukovsky, and Timokha (2000),
where a set of natural harmonic functions are defined to satisfy the boundary
conditions. This definition is challenging since each individual tank shape
requires a dedicated applied mathematical and physical study. Moreover,
Faltinsen and Timokha (2001) reported that the simple limitation of the
infinite sum by a finite number $r$ can yield either inaccurate or expensive
computations. Thus, the selection of the dominant modes is a non-trivial task.
A truncation based on employing special asymptotic relationships, postulated
following mathematical or physical arguments, has been shown to produce good
results for the sloshing problem Faltinsen  _et al._ (2000); Faltinsen and
Timokha (2002, 2010); Lukovsky, Ovchynnykov, and Timokha (2012); Ansari,
Firouz-Abadi, and Ghasemi (2011); Gómez-Goñi  _et al._ (2013). Nevertheless,
these relations are valid under reasonable assumptions for specific tank
geometries, which poses a fundamental challenge in the multimodal method’s
application in arbitrary settings Faltinsen and Timokha (2001). It is,
therefore, tempting to explore emerging data-driven tools to mitigate such
problems. For example, as discussed in Section II.2, one could consider the
POD procedure, which provides a systematic framework that yields a set of
basis functions (accompanied by a sorting mechanism) from a set of snapshots.

###  II.2 Proper orthogonal decomposition: linear best-fit basis functions

In addition to boundary conditions, we might have archival data (i.e.,
snapshot fields) to help us construct the basis functions. PCA Pearson (1901)
can be used to construct the basis functions that optimally represent the
data. In 1933, a geometric representation of PCA has been proposed by
Hotelling Hotelling (1933), and this concept has later become popular as
empirical orthogonal functions (EOF) Obukhov (1947) in environmental science,
and POD Aubry (1991) in the fluid dynamics community. The _method of
snapshots_ Sirovich (1987) has been instrumental in the development of POD
based approaches Berkooz, Holmes, and Lumley (1993). To compute the POD basis
functions, let us assume that we have access to $m$ snapshots
$\bm{u}(\bm{x},t_{i})$ for $i=1,2,\dots,m$. A Reynolds decomposition-like
expansion can be written as

$\displaystyle{\color[rgb]{0,0,0}\bm{u}(\bm{x},t_{i})=\bar{\bm{u}}(\bm{x})+\bm{\upsilon}(\bm{x},t_{i})},$
(18)

where $\bar{\bm{u}}(\bm{x})$ is a reference (e.g., ensemble mean) field, which
can be obtained as

$\displaystyle{\color[rgb]{0,0,0}\bar{\bm{u}}(\bm{x})=\frac{1}{m}\sum_{i=1}^{m}\bm{u}(\bm{x},t_{i})},$
(19)

and the set of anomaly snapshots $\bm{\upsilon}(\bm{x},t_{i})$ for
$i=1,2,\dots,m$ can be defined as

$\displaystyle{\color[rgb]{0,0,0}\bm{\upsilon}(\bm{x},t_{i})=\bm{u}(\bm{x},t_{i})-\bar{\bm{u}}(\bm{x})}.$
(20)

For clarity, we introduce the POD procedure for a scalar field (POD is
normally applied to the velocity vector field). Let us denote
$\upsilon(\bm{x},t_{j})$ a component of the anomaly velocity vector field
(e.g., the $x$-component). A temporal correlation matrix
$\textbf{A}=[\alpha_{ij}]$ can be constructed from these anomaly snapshots:

$\displaystyle\alpha_{ij}=\int_{\Omega}\upsilon(\bm{x},t_{i})\upsilon(\bm{x},t_{j})d\bm{x},\
\ \bm{x}\in\Omega,$ (21)

where $\Omega$ is the spatial domain, and $i$ and $j$ refer to the snapshot
indices. We define the $L^{2}$ inner product of two functions $f$ and $g$ as

$\displaystyle\Big{(}f(\cdot),g(\cdot)\Big{)}=\int_{\Omega}f(\bm{x})g(\bm{x})d\bm{x},$
(22)

which yields
$\alpha_{ij}=\Big{(}\upsilon(\bm{x},t_{i}),\upsilon(\bm{x},t_{j})\Big{)}$ from
Eq. 21. The data correlation matrix $\textbf{A}=[\alpha_{ij}]$ is a non-
negative, symmetric $m\times m$ matrix, also known as the Gramian matrix of
$\upsilon(\bm{x},t_{1})$, $\upsilon(\bm{x},t_{2})$, $\dots$,
$\upsilon(\bm{x},t_{m})$. If we define the diagonal eigenvalue matrix
$\bm{\Lambda}=\text{diag}[\lambda_{1},....,\lambda_{m}]$ and a right
eigenvector matrix $\bm{\Gamma}=[\bm{\gamma}_{1},....,\bm{\gamma}_{m}]$ whose
columns are the corresponding eigenvectors of A, we can solve the following
eigenvalue problem to obtain the optimal POD basis functions Ravindran (2000):

$\displaystyle\textbf{A}\bm{\Gamma}=\bm{\Gamma}\bm{\Lambda}.$ (23)

In general, most of the subroutines for solving Eq. 23 give $\bm{\Gamma}$ with
all of the eigenvectors normalized to unity. The orthonormal POD basis
functions for the anomaly field, $\upsilon$, can be thus calculated as
follows:

$\displaystyle\psi_{i}(\bm{x})=\frac{1}{\sqrt{\lambda_{i}}}\sum_{k=1}^{m}\gamma_{i}^{k}\upsilon(\bm{x},t_{k}),$
(24)

where $\lambda_{i}$ is the $i^{th}$ eigenvalue, $\gamma_{i}^{k}$ is the
$k^{th}$ component of the $i^{th}$ eigenvector, and $\psi_{i}(\bm{x})$ is the
$i^{th}$ POD mode.

The eigenvalues are often stored in descending order for practical purposes,
i.e., $\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{m}\geq 0$, and the
eigenvectors are normalized in such a way that the basis functions satisfy the
following orthonormality condition:

$\Big{(}\psi_{i},\psi_{j}\Big{)}=\left\\{\begin{aligned} 1,\ \ i=j;\\\ 0,\ \
i\neq j.\end{aligned}\right.$ (25)

Now, we can linearly represent the anomaly field variable $\upsilon(\bm{x},t)$
using the POD modes as follows:

$\displaystyle\upsilon(\bm{x},t)=\sum_{i=1}^{r}a_{i}(t)\psi_{i}(\bm{x}),$ (26)

where $a_{i}$ are the time-dependent (pseudo) modal coefficients, and $r$ is
the total number of retained modes after the truncation, with $r\ll m$. These
$r$ modes with the largest energy content correspond to the largest
eigenvalues ($\lambda_{1},\lambda_{2},...,\lambda_{r}$). In general, adding
more POD modes reduces the POD-ROM error. We note, however, that this is not
always true. For example, adding POD modes that are polluted by numerical
noise can actually decrease the POD-ROM error (see, e.g., the numerical
investigation in Ref.Giere  _et al._ (2015)). Often the value of $r$ is
determined by using the relative information content (RIC) index Gunzburger
(2002), which is defined as

$\displaystyle\mbox{RIC}=\frac{\sum_{i=1}^{r}\lambda_{i}}{\sum_{i=1}^{m}\lambda_{i}},$
(27)

where $\mbox{RIC}=1$ refers to a complete representation of the data
snapshots. For example, if one records a set of snapshots from the field given
by Eq. 17 (i.e., the solution of the advection problem in Eq. 2), let us say
$m=100$ or more equally distributed snapshots between $t=0$ and $t=T$, the POD
analysis could offer a perfect representation with $\mbox{RIC}=1$ using only
two retained modes ($r=2)$ since the underlying dynamics can be constructed by
a linear superposition of two harmonics. Of course, that is not always the
case, and RIC becomes smaller than unity even if we retain a substantial
number of modes, especially for turbulent flows.

This need for a large number of modes is one of the chief motivating factors
for developing closure models to compensate the effects of the truncated modes
in ROMs. However, it is believed that there is no separation of scales in
turbulence, and therefore, most turbulent flow problems cannot be
characterized by a high RIC index. Specifically, if there is no significant
pattern in the evolution dynamics, there is a slow decay rate for the
eigenvalues $\lambda_{k}$, and retaining only a few modes cannot capture the
essential dynamics of turbulence. Thus, it is not surprising that many ROM
practitioners have often demonstrated their proposed methodologies for
problems that show somehow an underlying pattern (e.g., a shedding pattern in
simulating the von Karman street).

This picture can be linked to the Kolmogorov barrier Ahmed and San (2020),
where the linear reducibility (i.e., representing the underlying fluctuation
field as a linear superposition/span of a finite/limited number of basis
functions) is hindered. _Modal expansions have an elliptic nature by
construction, and using such tools for convection-dominated flows with higher
degrees of hyperbolicity might often add another level of complexity when
designing projection ROMs._ Then, a central question might arise about this
data-driven procedure: _Do we really get any benefit using the POD basis
functions generated from prerecorded snapshots?_ Although there might be a
trade-off between storage, accuracy, and efficiency, the answer probably
depends on the problem at hand. It might be a big _yes_ if there is an
underlying pattern (e.g., limit cycles or quasi-periodic oscillations), and
might be a _no_ if the flow is highly turbulent in a statistically non-
equilibrium and chaotic state. In the latter case, one might consider a
standard local discretization (e.g., finite difference/element/volume) method
or a pseudo-spectral method (supported by the harmonics that satisfy the
boundary conditions) without attempting to perform the POD procedure to
compute a set of data-driven global basis functions.

For instance, the fast Fourier transform (FFT) provides an extremely efficient
computational framework for models with such global basis functions without
requiring any additional storage for precomputed or measured snapshots to
generate a set of data-driven basis functions. The trade-off between accuracy
and computational efficiency should always be considered carefully in
generating data-driven models like Galerkin ROMs. The complexity of a typical
right-hand side (RHS) computation of a pseudo-spectral solver becomes slightly
bigger than $\mathcal{O}(n)$, where $n$ refers to the number of grid points.
In contrast, the complexity of a typical Galerkin ROM is $\mathcal{O}(r^{3})$
(there are also additional costs associated with, e.g., collecting and
processing snapshots or solving an eigenvalue problem to generate a set of
basis functions). Therefore, the Galerkin ROM becomes a computationally
feasible approach _if and only if_ a few number of retained modes are
utilized. As a rule of thumb, $r$ should be significantly less than the number
of grid points in each direction for a canonical 3D problem (e.g., $r\ll 256$
for a $256^{3}$ problem). Otherwise, it would be hard to justify that the
model is indeed _reduced order_ , since instead we could simply use the FFT
algorithm to integrate the dynamical system equations in the harmonic space.
Same arguments hold true for using a more flexible and convenient localized
model, especially for problems with more complicated geometries (e.g., with
the finite element, finite difference, or finite volume method, where the RHS
can be obtained in $\mathcal{O}(n)$ computations).

###  II.3 Leveraging ROR

One of the major reasons for the inaccuracy of current ROMs in the numerical
simulation of complex flows is the quality of ROR, which is the ability or
inability of the ROM framework to represent the underlying complex dynamics.
Specifically, in order to determine whether there is a valid ROR of the given
system, we need to answer the following questions: (i) Is the ROM basis able
to accurately approximate the dynamics? (ii) Is the Galerkin projection
yielding an accurate ROM?

Furthermore, the issue of the selection of a convenient domain often comes
into play. If the domain of influence is too small, there might be no good
dynamical prediction. On the other hand, when it is too large, there could be
too many uncorrelated events that have to be lumped in global modes. These
uncorrelated events might work against the modeling accuracy since they often
increase the deformation of modal expansion or degradation of the model
representation. In other words, the domain should be as small as possible and
as large as necessary. The correlation length might be a good initial
reference scale to define the domain of interest. Dynamic mode adaptation,
parameter-space-time domain partitioning, as well as smart clustering ideas
have been explored, although we believe this topic is still in its infancy.

In practice, the construction of _a good_ low-order space is a cornerstone in
projection-based ROM. That said, the representability of POD basis functions
becomes questionable for non-stationary, strongly-evolving, and convection-
dominated flows. Being a linear-based approach, POD might not be sufficient to
describe nonlinear processes. More importantly, using a Galerkin projection
based on elliptical ansatz for a hyperbolic problem could generate numerical
oscillations. Moreover, the POD is optimal _globally_ in the sense that it
minimizes the _averaged_ $L^{2}$ error across all the snapshot data. This
raises the issue of modal deformation by the rapidly varying flow field state
in such a way that the resulting modes are not representative of any of the
system’s states. Furthermore, since the POD modes are ranked based on their
energy content, excursions in state spaces that contain a small amount of
energy can be overlooked by POD even if these excursions might have
significant impact on the dynamical evolution (see Cazemier, Verstappen, and
Veldman (1998) for example). Similar scenarios arise for parameterized systems
spanning a large parameter space when the system’s behavior highly depends on
the parameter value. Therefore, we devote the rest of this section to efforts
aimed at enhancing the basis representability by either improving the offline
construction stages or efficiently updating the ROM during online deployment.

One of the simplest approaches to improve the quality of the POD basis
functions is to enrich the snapshot data matrix with _extra_ information. For
example, in addition to the exact flow field data, the scaled difference
between consecutive snapshots (i.e., the difference quotients) can be utilized
such that the time derivative information is better represented in the
resulting modes Kunisch and Volkwein (2001b); Koc  _et al._ (2021a). Moreover,
instead of collecting snapshot data at arbitrary time intervals and/or
parameter values, more effective sampling techniques should be pursued. In
Ref. Cui, Marzouk, and Willcox (2015), a ROM is integrated into a Markov chain
Monte Carlo (MCMC) framework, where the posterior distribution estimated by
the MCMC algorithm is utilized to adaptively select the parameter values at
which snapshots are evaluated.

In an effort towards the accurate identification of coherent structures from
experimental or numerical data, a spectral POD approach has been developed
Sieber, Paschereit, and Oberleithner (2016), and its relationship to DMD and
resolvent analysis has been established (Towne, Schmidt, and Colonius, 2018).
These studies use a technique that performs several PODs on individual
frequencies obtained from Fourier-transformed windows of snapshot data. Thus,
the modes we get for each frequency correspond to a coherent frequency domain
structure. If an analysis of the eigenvalue spectrum for each frequency
reveals coherent structures, it can indicate that there is a physical process
which is occurring at that frequency. Hence, this becomes a useful data
analysis tool on top of providing orthogonal modes for ROMs. In addition,
transported snapshots POD approaches (Reiss  _et al._ , 2018; Mendible  _et
al._ , 2020) have been introduced for convection-dominated transport systems.
In particular, these studies use a shifting operator on the snapshots
(requiring interpolation on unstructured grids and some knowledge of the
transport speed) to allow POD or DMD to (more) efficiently approximate
advective systems. In their recent works, Mendible  _et al._ (2021) employed
an unsupervised traveling wave identification with shifting and truncation
(UnTWIST) algorithm Mendible  _et al._ (2020) to discover moving coordinate
frames into which the data are shifted, thus overcoming limitations imposed by
the underlying translational invariance and allowing for the application of
traditional dimensionality reduction techniques. Etter and Carlberg (2020)
proposed a novel online adaptive basis refinement mechanism for efficiently
enriching the trial basis in a manner that ensures convergence of the ROM to
the FOM.

Localization methods have been successfully pursued to mitigate the modal
deformation of the POD basis by partitioning the state space Amsallem, Zahr,
and Farhat (2012); Washabaugh  _et al._ (2012); Peherstorfer  _et al._ (2014);
Taddei, Perotto, and Quarteroni (2015), time domain IJzerman (2000); Dihlmann,
Drohmann, and Haasdonk (2011); Borggaard, Hay, and Pelletier (2007); San and
Borggaard (2015); Babaee and Sapsis (2016); Chaturantabut (2017); Ahmed and
San (2018); Ahmed  _et al._ (2019); Ahmed and San (2020), physical domain
Maday and Rønquist (2002); Løvgren, Maday, and Rønquist (2006); Iapichino,
Quarteroni, and Rozza (2012); Eftang and Patera (2013), or parameter space
Eftang and Stamm (2012); Moosavi, Stefanescu, and Sandu (2015); Hess  _et al._
(2019) using multiple local, piecewise affine subspace approximations instead
of a single global approximation. These partitioning or time varying
approaches work by parsing the available snapshot data into a few overlapping
or non-overlapping groups (e.g., based on solution value, time, parameter,
geometry, or component) and applying standard modal decomposition techniques
(e.g., POD) for each region separately. This eventually yields a library of
compact ROMs, each suitable for a specific region and/or dynamics, and
interpolation methods can be utilized when the region of interest does not
exist in the available library.

In this context, clustering techniques can be also utilized to effectively
perform such partitioning. Indeed, cluster-based reduced order models (CROMs)
have been proposed to tackle some of the potential pitfalls of classical GROMs
(e.g., the mismatch between the modal expansion approach and the underlying
dynamics, see Noack (2016)). CROMs start by sorting the snapshot data into a
small number of clusters (e.g., using k-means approach) with centroids being
the representative states in each cluster. Conceptually, this is similar to
coarse-graining the state-space (or generally the feature-space) into
centroidal Voronoi tessellation (CVT) generators Burkardt, Gunzburger, and Lee
(2006). The transition dynamics between these centroids can be modeled as a
probabilistic Markov model Kaiser  _et al._ (2014, 2017); Li and Tan (2020) or
a deterministic–stochastic network model Fernex  _et al._ (2019); Li  _et al._
(2021); Fernex, Noack, and Semaan (2021).

As highlighted in Section II.2, POD provides an efficient way to compress data
and explain the variance of the data better than any other linear combination
Jolliffe (2002). Indeed, from a linear algebra perspective, it can often be
formulated as a singular value decomposition, providing an optimal low-rank
matrix approximation. This can leverage highly performant and scalable
algorithms to handle extremely large datasets, benefiting from the rich legacy
of linear algebra investigations. From a statistical point of view, this
orthogonal projection provides _linearly uncorrelated_ features. However, it
cannot reveal nonlinear correlations in the data. In contrast, _manifold
learning_ (or representation learning) techniques aim at accounting for such
nonlinear correlations to further reduce the dimensionality of the problem.
The generalizations of PCA to nonlinear settings often define a curve in the
latent space which minimizes the mean squared error of all variables. Yet, the
smoothness of the curve can be varied by the method. For example, an
autoassociative or autoencoding neural network model Kramer (1991); Hsieh
(2001, 2009) and a kernel PCA Schölkopf, Smola, and Müller (1998) are two
successful approaches of such a nonlinear PCA (NLPCA) framework. We refer the
reader to recent works Otto and Rowley (2019); Fukami, Nakamura, and Fukagata
(2020); Agostini (2020) for excellent discussions on autoencoder technology in
fluid dynamics (see Figure 2). We also note that other nonlinear
dimensionality reduction techniques, such as principal curves Hastie and
Stuetzle (1989), locally linear embedding Roweis and Saul (2000), isomap
Tenenbaum, De Silva, and Langford (2000) and self-organized map Kohonen (1982)
approaches, can also be regarded as a discrete version of NLPCA.

Figure 2: A schematic diagram for an autoencoder (AE) for latent space
construction, where the input $\mathcal{X}$ is the full field, the output
$\mathcal{Y}$ designates its reconstruction, and the bottleneck $\mathcal{Z}$
represents the latent space (compressed) variables.

It is worth noting that during the modal truncation step (reducing the
dimensionality of the system), dependencies among retained and discarded modes
generally yield inaccurate results if the closure problem is not addressed. In
this regard, persistent homology (PH) Edelsbrunner, Letscher, and Zomorodian
(2000); Edelsbrunner, Harer  _et al._ (2008); Edelsbrunner and Morozov (2014);
Zomorodian and Carlsson (2008, 2005) provides a delicate balance between data
simplification and intrinsic structure extraction. PH is a tool in topological
data analysis that aims at studying and extracting the features that persist
across multiple scales by casting the multi-scale organization into a
mathematical formalism. In particular, PH measures the lifetime of intrinsic
topological features using a filtration process to distinguish between the
long-lived features and the short-lived ones (which are considered topological
noise) Otter  _et al._ (2017).

However, the application of PH has been largely dedicated to qualitative data
classification and analysis, and its utility for quantitative modeling and
prediction (including ROM) is scarce Wang and Wei (2016). PH has been utilized
to characterize the time series from dynamical systems based on topological
features that appear in the solution manifold or attractor Pereira and de
Mello (2015); Garland, Bradley, and Meiss (2016); Maletić, Zhao, and Rajković
(2016). Rieck and Leitte (2015) used PH as an evaluation tool to compare the
performance of different dimensionality reduction algorithms (e.g., PCA, and
isomap mentioned earlier). One of the major challenges of employing PH is its
prohibitive computational cost (for the worst-case scenario). To increase the
PH efficiency, Moitra, Malott, and Wilsey (2018) utilized a clustering
technique to represent similar groups of data points with their cluster
centroid and applied PH onto these clusters.

##  III First Principle Galerkin Method

Finite dimensional low-order models routinely arise when we apply Galerkin
type projection techniques to infinite dimensional PDE models Lakshmivarahan
and Wang (2008a, b); Wang and Lakshmivarahan (2009). We formalize the model
reduction problem for fluid flow systems, considering a generic prognostic
equation as follows:

$\dfrac{\partial\mathbf{u}}{\partial
t}=\mathcal{\bm{F}}(\mathbf{u};\mathbf{x},t),$ (28)

where $\mathbf{u}$ denotes the discrete approximation of a three-dimensional
(3D) dependent variable (e.g., density, velocity, temperature, moisture);
$\mathbf{x}$ denotes the independent spatial variables (e.g., latitude,
longitude, and height); and $\mathcal{\bm{F}}$ defines the model’s dynamical
core (e.g., semi-discretized PDEs representing mass, momentum, and energy
conservation), all written in vector-form. More specifically, we explore the
autonomous dynamics case for a specific quantity of interest $\mathbf{u}$,
with $\mathcal{\bm{F}}$ being decomposed into linear $\mathcal{\bm{L}}$ and
nonlinear $\mathcal{\bm{N}}$ operators as

$\dfrac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=\mathcal{\bm{L}}\mathbf{u}+\mathcal{\bm{N}}(\mathbf{u}),$
(29)

where $\mathbf{u}\in\mathbb{R}^{n}$, where $n$ is the number of degrees of
freedom in the spatial discretization. Considering the Navier-Stokes equations
as a typical mathematical framework for fluid flow modeling, we highlight that
the linear and nonlinear operators often represent the diffusive and
convective effects, respectively.

In order to build the projection-based ROM, the solution $\mathbf{u}$ is
approximated in a low-dimensional affine subspace of dimension $r$ via the
Galerkin ansatz as follows:

$\mathbf{u}(t)\approx\bar{\mathbf{u}}+\bm{\Psi}\bm{a}(t),$ (30)

where $\bar{\mathbf{u}}\in\mathbb{R}^{n}$ is a reference solution representing
the affine offset, $\bm{\Psi}\in\mathbb{R}^{n\times r}$ denotes the trial
basis, and $\bm{a}(t)\in\mathbb{R}^{r}$ is the vector of reduced (generalized)
coordinates, also called modal weights or coefficients. The reference solution
$\bar{\mathbf{u}}$ as well as the basis $\bm{\Psi}$ are constructed during an
offline stage from a collection of FOM evaluations (called snapshots). Without
loss of generality, we suppose that the time-averaged field defines the
reference solution, $\bar{\mathbf{u}}$, and the basis $\bm{\Psi}$ is
constructed using the POD technique. Then, the low-rank approximation given by
Eq. 30 is substituted into Eq. 29 and an inner product with a test basis is
performed to yield a system of ODEs for the unknown modal coefficients,
$\bm{a}(t)$. In Galerkin projection-based ROM (GROM), the test basis is chosen
to be the same as the trial basis. We make use of the orthonormality property
(i.e., $\bm{\Psi}^{\top}\bm{\Psi}=\bm{I}_{r}$, where $\bm{I}_{r}$ is the
$r\times r$ identity matrix and the superscript $\top$ denotes the matrix
transpose, assuming a Euclidean state space) as follows:

$\bm{\Psi}^{\top}\dfrac{\mathrm{d}}{\mathrm{d}t}\left(\bar{\mathbf{u}}+\bm{\Psi}\bm{a}\right)=\bm{\Psi}^{\top}\mathcal{L}\left(\bar{\mathbf{u}}+\bm{\Psi}\bm{a}\right)+\bm{\Psi}^{\top}\mathcal{N}\left(\bar{\mathbf{u}}+\bm{\Psi}\bm{a}\right).$
(31)

Since both $\bar{\mathbf{u}}$ and $\bm{\Psi}$ are considered time-independent,
Eq. 31 reduces to

$\bm{\Psi}^{\top}\bm{\Psi}\dfrac{\mathrm{d}\bm{a}}{\mathrm{d}t}=\bm{\Psi}^{\top}\mathcal{L}\bar{\mathbf{u}}+\bm{\Psi}^{\top}\mathcal{L}\bm{\Psi}\bm{a}+\bm{\Psi}^{\top}\mathcal{N}\left(\bar{\mathbf{u}}+\bm{\Psi}\bm{a}\right).$
(32)

Note that $\bm{\Psi}^{\top}\mathcal{L}\bar{\mathbf{u}}$ and
$\bm{\Psi}^{\top}\mathcal{L}\bm{\Psi}$ can be precomputed during the offline
construction stage, reducing the online computational cost of evaluating the
first two terms on the right-hand side to $O(r)$ independent of the FOM
dimension, $n$. However, generally speaking, computing the third term
representing the system’s nonlinearity depends on $n$, limiting the
computational benefit of ROM. In order to mitigate this limitation,
hyperreduction approaches have been developed to relieve this dependency on
$n$, by _approximating_ , rather than _evaluating_ , the nonlinear term in a
reduced order subspace Ştefănescu, Sandu, and Navon (2014); Dimitriu,
Ştefănescu, and Navon (2017); Yano (2019). Examples of hyperreduction include
the empirical interpolation method (EIM) Barrault  _et al._ (2004a), its
discrete version (DEIM) Chaturantabut and Sorensen (2009, 2010), the gappy POD
Everson and Sirovich (1995); Bui-Thanh, Damodaran, and Willcox (2004);
Carlberg  _et al._ (2013), and the missing point estimation (MPE) Astrid  _et
al._ (2008); Zimmermann and Willcox (2016), where the approximation is
performed using sampling techniques. On the other hand, tensorial ROM can
benefit from the quadratic (or generally polynomial) nonlinearity, which is
ubiquitous in fluid flow systems, to rewrite Eq. 32 as follows:

$\dfrac{\mathrm{d}\bm{a}}{\mathrm{d}t}=\mathbf{\mathfrak{B}}+\mathbf{\mathfrak{L}}\bm{a}+\bm{a}^{\top}\mathbf{\mathfrak{N}}\bm{a},$
(33)

where the vector $\mathbf{\mathfrak{B}}$, the matrix $\mathbf{\mathfrak{L}}$,
and the tensor $\mathbf{\mathfrak{N}}$ are precomputed during the offline
stage, reducing the computational cost of solving the GROM defined in Eq. 33
to $O(r^{3})$ in case of quadratic nonlinearity (which is the case for the
Navier-Stokes equations). When the true underlying dynamics of the system are
non-polynomial, _lifting_ transformations can be exploited to yield a finite-
dimensional coordinate representation in which the system dynamics have
quadratic structure Gu (2011); Benner and Breiten (2015); Benner, Goyal, and
Gugercin (2018); Kramer and Willcox (2019). Although such transformation is
not universally guaranteed, a large class of smooth nonlinear systems that
appear in engineering applications (e.g., elementary functions like
exponential and trigonometric functions or polynomials) can be equivalently
lifted to quadratic form.

###  III.1 POD Galerkin projection: Burgers equation

To illustrate the POD Galerkin approach for flow systems with quadratic
nonlinearity, let us consider the Burgers equation

$\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=\nu\frac{\partial^{2}u}{\partial x^{2}},$ (34)

which is often used as a simplified prototype by fluid dynamicists. Using the
POD procedure outlined in Section II.2, we can define the $u(x,t)$ field as a
linear superposition of the mean field and the POD basis functions,

$\displaystyle u(x,t)=\bar{u}(x)+\sum_{i=1}^{r}a_{i}(t)\psi_{i}(x),$ (35)

and substitute this approximation of our field variable into Eq. 34. Once we
perform an orthonormal Galerkin projection, the resulting dynamical system for
$a_{k}(t)$ can be written as

$\frac{da_{k}}{dt}=\mathfrak{B}_{k}+\sum_{i=1}^{r}\mathfrak{L}^{i}_{k}a_{i}+\sum_{i=1}^{r}\sum_{j=1}^{r}\mathfrak{N}^{ij}_{k}a_{i}a_{j},$
(36)

where

$\displaystyle\mathfrak{B}_{k}=\big{(}\nu\frac{\partial^{2}\bar{u}}{\partial
x^{2}}-\bar{u}\frac{\partial\bar{u}}{\partial x},\psi_{k}\big{)},$
$\displaystyle\mathfrak{L}^{i}_{k}=\big{(}\nu\frac{\partial^{2}\psi_{i}}{\partial
x^{2}}-\bar{u}\frac{\partial\psi_{i}}{\partial
x}-\psi_{i}\frac{\partial\bar{u}}{\partial x},\psi_{k}\big{)},$
$\displaystyle\mathfrak{N}^{ij}_{k}=\big{(}-\psi_{i}\frac{\partial\psi_{j}}{\partial
x},\psi_{k}\big{)}.$ (37)

This tensorial system consists of $r$ coupled ODEs and it is often written as
Eq. 33, where $\bm{a}$ is the vector of unknown coefficients $a_{k}(t)$,
$k=1,2,\dots,r$, $\mathbf{\mathfrak{B}}$ is a scaling vector coming from the
reference mean field with entries $\mathfrak{B}_{k}$, $\mathbf{\mathfrak{L}}$
is an $r\times r$ matrix with entries $\mathfrak{L}^{i}_{k}$ for the
contribution stemming from the linear viscous term, and
$\mathbf{\mathfrak{N}}$ is an $r\times r\times r$ tensor with entries
$\mathfrak{N}^{ij}_{k}$ arising from the nonlinear advection term, $1\leq
i,j,k\leq r$. In this tensorial form, the corresponding model coefficients
$\mathfrak{B}_{k}$, $\mathbf{\mathfrak{L}}$, and $\mathbf{\mathfrak{N}}$ are
_precomputed_ from the available snapshots. Alternatively, there are a number
of _online_ approaches where we can compute the nonlinear part using
hyperreduction or principled sampling strategies to approximate the full
nonlinear state from a small number of measurement or collocation points
Everson and Sirovich (1995); Barrault  _et al._ (2004a); Willcox (2006);
Yildirim, Chryssostomidis, and Karniadakis (2009); Chaturantabut and Sorensen
(2010); Drmac and Gugercin (2016). Ştefănescu, Sandu, and Navon (2014)
performed a comparative study between the direct (online) and tensorial
(precomputed) methods. Moreover, Karasözen, Yıldız, and Uzunca (2021) recently
discussed structure preserving ROMs and compared the direct and tensorial POD
approaches.

###  III.2 Projection based ROMs

Computational models for the Navier-Stokes equations could make a tremendous
impact in critical applications, such as the biomedical and engineering
applications that we describe next. Despite their enormous potential, FOMs
have not fully transitioned to engineering practice. The main roadblock is the
extraordinary computational cost incurred by computational models in many
applications. For example, although preliminary studies of aortic dissections
showed that uncertainties in the geometry and inflow conditions have a
fundamental role, performing an uncertainty quantification study requires a
huge number of computational model runs. Similarly, performing a shape
optimization study to determine the optimal vascular configuration for the
total cavopulmonary connection surgery requires again many computational model
runs. Also, in renewable energy applications, performing data assimilation to
incorporate the available observations in the control of wind-power production
requires numerous model runs.

Since running current computational models hundreds and thousands of times can
take days and weeks on high performance computing (HPC) platforms, a brute
force computational approach for these biomedical and engineering applications
is simply not possible. Therefore, what is needed is a modeling strategy that
allows model runs that take minutes to hours on a laptop.

For structure-dominated systems, ROMs can decrease the FOM computational cost
by orders of magnitude. ROMs are (extremely) low-dimensional models that are
trained (constructed) from available data. As explained in Section III, in an
offline phase, the FOM is run for a few parameters values to construct a low-
dimensional (e.g., $10$-dimensional) ROM basis
$\\{\psi_{1},\ldots,\psi_{10}\\}$, which is used to build the ROM:

$\displaystyle\frac{d\bm{a}}{dt}=\bm{f}(\bm{a}),$ (38)

where $\bm{a}$ is the vector of coefficients in the ROM approximation
$\sum_{i=1}^{10}a_{i}(t)\psi_{i}(\bm{x})$ of the variable of interest and
$\bm{f}$ comprises the ROM operators (e.g., vectors, matrices, and tensors)
that can be preassembled from the ROM basis in the offline phase. In the
online phase, the low-dimensional ROM given by Eq. 38 is then used for
parameters values that are different from those used in the training stage.
Since ROM is low-dimensional (10-dimensional), its computational cost is
orders of magnitude lower than the FOM cost. Thus, for the biomedical and
engineering applications described above, ROMs appear as a natural alternative
to the prohibitively expensive FOMs.

Unfortunately, current ROMs cannot be used in clinical and engineering
practice, since they would require too many modes (degrees of freedom). For
example, to capture all the geometric scales in aortic dissection, one might
need hundreds or even thousands of ROM modes (e.g., see Table 2). Similarly,
to cope with the high Reynolds number in the wind farm optimization, a large
number of ROM modes are necessary. Thus, although ROMs decrease the FOM
computational cost by orders of magnitude, they are still too expensive:
Current ROMs cannot be run in minutes or hours on a laptop and thus cannot be
used easily in clinical and engineering practice.

Table 2: A non-exhaustive list illustrating energy characteristics and range of the number of retained modes $r$ for $m$ given snapshots. Study |  Problem |  $r$ |  $m$ |  Comment
---|---|---|---|---
Östh  _et al._ (2014) |  Ahmed body |  10 – 100 |  2000 |  The first 500 modes resolve 60% of the kinetic energy.
San and Borggaard (2015) |  Marsigli flows |  6 – 30 |  400 |  The first 30 modes resolve 90% of the kinetic energy.
Rahman, Ahmed, and San (2019) |  Quasigeostrophic flows |  10 – 80 |  400 |  The first 50 modes resolve 80% of the kinetic energy.
Ballarin  _et al._ (2016) |  Hemodynamics |  50 |  400 |  O(10–100) modes are required to obtain a reliable approximation.
VerHulst and Meneveau (2014) |  Wind farm |  N/A |  7200 |  430 POD modes are required to capture 80% of the total energy.
Shah and Bou-Zeid (2014) |  Atmospheric boundary layer |  N/A |  2500 |  500 POD modes are required to capture 80% of the total energy.
Zhang and Stevens (2020) |  Atmospheric boundary layer |  N/A |  5000 |  2000 POD modes are required to capture 80% of the total energy.
Figure 3: Evolution of the ROM approaches.

With the evolution of ROM approaches outlined in Figure 3 in mind, the ROM
community is at a _crossroads_. On the one hand, current ROMs can be used for
academic test problems for which a handful of ROM modes can model simple
dynamics with substantial success. On the other hand, realistic, complex flows
require high-dimensional ROMs that cannot be used in clinical and engineering
practice. What is needed is low-dimensional, _efficient_ ROMs that are
_accurate_ so that they can be utilized in such vital applications.

One of the main reasons for the notorious inaccuracy of current ROMs in
complex clinical and engineering settings is the drastic ROM truncation:
Instead of using many (e.g., 100) ROM modes
$\\{\psi_{1},\ldots,\psi_{100}\\}$, current ROMs use only a handful of ROM
modes $\\{\psi_{1},\ldots,\psi_{10}\\}$ to ensure a low computational cost.
This drastic truncation yields acceptable results in simple, academic test
problems, but produces inaccurate results in practical clinical and
engineering settings Östh  _et al._ (2014). Thus, for accurate results, the
ROM closure problem needs to be solved: One needs to model the effect of the
discarded ROM modes $\\{\psi_{11},\ldots,\psi_{100}\\}$ on the ROM dynamics,
i.e., on the time evolution of resolved ROM modes
$\\{\psi_{1},\ldots,\psi_{10}\\}$:

$\displaystyle\frac{d\bm{a}}{dt}=\bm{f}(\bm{a})+\text{Closure}(\bm{a}),$ (39)

where $\text{Closure}(\bm{a})$ is a low-dimensional term that models the
effect of the discarded ROM modes $\\{\psi_{11},\ldots,\psi_{100}\\}$ on
$\\{\psi_{1},\ldots,\psi_{10}\\}$. The closure term is also known as
_unresolved tendency_ , or _model error_ in different disciplines.

The closure problem is prevalent in numerical simulation of complex systems.
For example, classical numerical discretization of turbulent flows (e.g.,
finite element or finite volume methods), inevitably takes place in the under-
resolved regime (e.g., on coarse meshes) and requires closure modeling (i.e.,
modeling the sub-grid scale effects). In computational fluid dynamics (CFD),
e.g., large eddy simulation (LES), there are hundreds (if not thousands) of
closure models Sagaut (2006). This is in stark contrast with reduced order
modeling, where only relatively few ROM closure models have been investigated.
The reason for the discrepancy between ROM closure and LES closure is that the
latter has been mostly built around physical insight stemming from
Kolmogorov’s statistical theory of turbulence (e.g., the concept of eddy
viscosity), which is generally posed in the Fourier setting Pope (2000);
Sagaut (2006). Much of this physical insight is not generally available in a
ROM setting. Thus, current ROM closure models have been deprived of this
powerful methodology that represents the core of most LES closure models. To
construct low-dimensional and efficient ROMs that are _accurate_ , a set of
principled, mathematical and/or data-driven ROM closure modeling strategies
need to be utilized. In Section IV, we survey the main types of closure models
developed in the reduced order modeling community.

##  IV Closure Modeling

Although the solution of Eq. 33 becomes independent of the FOM dimension $n$,
the cubic scaling with respect to $r$ hurts the turnaround of such ROMs. This
is especially true for fluid flows of practical interest (e.g., turbulent and
convection-dominated flows), where the FOM solution manifold is characterized
by a large and slowly decaying Kolmogorov $n$-width Kolmogoroff (1936); Pinkus
(1985). Thus, a large number of modes are required to maintain the solution
accuracy, resulting in excessive computational overhead, which may even exceed
the FOM computational cost. Therefore, in these complex settings, ROM will
always incur a degree of under-resolution by sacrificing some degree of
accuracy for the sake of computational efficiency. This under-resolution has
direct and indirect consequences. The direct outcome is the projection error
affecting the Galerkin ansatz (Eq. 30), where some of the underlying flow
features are lost. The indirect ramifications are related to the nonlinearity
of the system, implying that the discarded modes indeed interact with the
retained ones. By performing severe modal truncation (remember, computational
efficiency is a priority!), we suppress these interactions and Eq. 33 no
longer captures the projected trajectory, decreasing the solution accuracy
Ahmed  _et al._ (2020a).

To illustrate the above discussion, consider a state variable
$\mathbf{u}(t)\in\mathbb{R}^{n}$, which can be _exactly_ written as a
superposition of $n$ basis functions as
$\mathbf{u}(t)=\bm{\Psi}\bm{a}(t)+\bm{\Phi}\bm{b}(t)$, where
$\bm{\Psi}\in\mathbb{R}^{n\times r}$ and
$\bm{\Phi}\in\mathbb{R}^{n\times(n-r)}$ represent the modes to be retained and
truncated, respectively, and $\bm{a}(t)\in\mathbb{R}^{r}$ and
$\bm{b}(t)\in\mathbb{R}^{(n-r)}$ are the corresponding time-dependent
coefficients. A Galerkin projection of the governing equations onto
$\bm{\Psi}$ and $\bm{\Phi}$ yields the following:

$\frac{d}{dt}\begin{bmatrix}\bm{a}\\\
\bm{b}\end{bmatrix}=\begin{bmatrix}\bm{f}_{\bm{a}}(\bm{a},\bm{b})\\\
\bm{f}_{\bm{b}}(\bm{a},\bm{b})\end{bmatrix}.$ (40)

We note that Eq. 40 is an exact representation of the system’s dynamics. In
reduced order modeling, we are only interested in the resolved part of the
dynamics, which can be written as

$\frac{d\bm{a}}{dt}=\bm{f}_{\bm{a}}(\bm{a},\bm{b}).$ (41)

Nevertheless, Eq. 41 is not practical because its solution requires the
knowledge of the unresolved variable, $\bm{b}$. In a classic truncated ROM, it
is often assumed that
$\bm{f}_{\bm{a}}(\bm{a},\bm{b})=\bm{f}_{\bm{a}}(\bm{a},\bm{0})=\bm{f}(\bm{a})$.
However, for nonlinear cases, this relation does not hold (i.e.,
$\bm{f}_{\bm{a}}(\bm{a},\bm{b})\neq\bm{f}_{\bm{a}}(\bm{a},\bm{0})$).

Figure 4: Closure modeling analogy between LES and ROM, where higher values of
$k$ and $m$ refer to smaller scales (adapted from Akhtar  _et al._ (2012); San
and Borggaard (2014); Azadani and Staples (2015)). Shaded areas in energy
spectra represent the discarded scales that must be modeled.

Following the Kolmogorov hypotheses Kolmogorov (1941, 1991) from turbulence
modeling and assuming an analogy between POD and Fourier modes (see Figure 4),
it is commonly agreed in the ROM community that the first POD modes resolve
the large energy-containing flow scales, while the last modes correspond to
the low-energy dissipative scales. Indeed, this analogy has been demonstrated
theoretically and numerically for different flow scenarios (e.g., flow over a
cylinder Noack  _et al._ (2003) and a turbulent flow past a backward-facing
step Couplet, Sagaut, and Basdevant (2003)). Thus, truncating the low-energy
scales is believed to result in a pile-up of energy levels, leading to
solution instability. We also highlight that this argument has been recently
the focus of scientific revisits. For instance, Grimberg, Farhat, and Youkilis
(2020) state, using mathematical arguments and analogies from finite element
analysis, that the solution instabilitiy observed in most studies dealing with
GROM is a byproduct of the Galerkin projection step. Moreover, they show that
a ROM based on Petrov-Galerkin projection, where the test basis differs from
the trial basis, yields more accurate and stable solution than standard
Galerkin projection. However, the test basis needs to be updated at each
iteration and time step, increasing the computational complexity of the
resulting ROM. As a highly promising approach, an adjoint Petrov-Galerkin
method for nonlinear model reduction has been recently put forth by Parish,
Wentland, and Duraisamy (2020). Rather than constructing a low-dimensional
subspace for the entire state space in a monolithic fashion, Hoang, Choi, and
Carlberg (2021) recently proposed a dynamic methodology to construct separate
subspaces for the different subdomains. Although the Petrov-Galerkin
projection Xiao  _et al._ (2013) could mitigate some of the challenges the
Galerkin ROMs have to face in the under-resolved simulation of turbulent
flows, we limit ourselves to Galerkin projection-based ROMs in the current
review, where closure models have been mainly developed to improve the
solution accuracy and stability properties. The major aim of closure models is
to make up for the effects of discarded modes onto the dynamics of resolved
modes. Specifically, the objective is to modify Eq. 33 to correctly resolve
the time dynamics of $\bm{\Psi}$:

$\dfrac{\mathrm{d}\bm{a}}{\mathrm{d}t}=\mathbf{\mathfrak{B}}+\mathbf{\mathfrak{L}}\bm{a}+\bm{a}^{\top}\mathbf{\mathfrak{N}}\bm{a}+\mathbf{\mathfrak{C}},$
(42)

where $\mathbf{\mathfrak{C}}$ represents the closure model that needs to be
determined.

The closure problem has historical roots in CFD, in particular in turbulence
modeling, including the Reynolds averaged Navier-Stokes (RANS) and LES. In
contrast, there are relatively fewer closure models that have been
investigated in a ROM context. In general, the ROM closure modeling approaches
can be classified into (1) functional (phenomenological), which use physical
insights to postulate a model form for the closure term (e.g., a dissipative
term) and (2) structural (mathematical), which often rely on filtering
techniques to reveal the closure term without using any physical assumptions
or additional phenomenological arguments. We emphasize that placing a given
ROM closure model in one of these categories is not always straightforward, as
these categories sometimes overlap. We also refer readers to the unified
exposition of several mean field modeling ideas Johnston, Simpson, and Baker
(2012) as well as other closure techniques for probability density function
(PDF) Brennan and Venturi (2018) or moment closures for kinetic theories
Levermore (1996).

Table 3: A chronological list of key contributions to ROM closure modeling. Year |  Study |  Key Contribution
---|---|---
1988 |  Aubry  _et al._ (1988) |  First closure model: global eddy viscosity modeling
1994 |  Rempfer and Fasel (1994) |  Linear modal eddy viscosity closure
1995 |  Selten (1995) |  Time averaging closure modeling
1997 |  Selten (1997) |  A statistical closure of a barotropic model
1998 |  Cazemier, Verstappen, and Veldman (1998) |  Penalty term closure model based on energy conservation principles
2003 |  Couplet, Sagaut, and Basdevant (2003) |  Guidelines for modeling unresolved modes in POD–Galerkin models
2004 |  Sirisup and Karniadakis (2004) |  Spectral viscosity closure for POD models
2008 |  Noack  _et al._ (2008) |  Finite time thermodynamics and ensemble averaging closure models
2009 |  Bergmann, Bruneau, and Iollo (2009a, b) |  Residual-based variational multiscale POD
2011 |  Borggaard, Iliescu, and Wang (2011) |  First numerical analysis of closure models: artificial viscosity model
2011 |  Akhtar  _et al._ (2012) |  Nonlinear eddy viscosity model based on the Frobenius norm of the Jacobian
2011 |  Wang  _et al._ (2011) |  Two-level discretization model
2012 |  Wang  _et al._ (2012) |  Eddy viscosity variational multiscale and dynamic Smagorinsky closures
2013 |  Balajewicz, Dowell, and Noack (2013) |  Subspace calibration using the Navier-Stokes equations
2013 |  Cordier  _et al._ (2013) |  Proof of global boundedness of nonlinear eddy viscosity closures
2014 |  Östh  _et al._ (2014) |  $\sqrt{K}$-scaled eddy viscosity concept
2014 |  Iliescu and Wang (2014) |  Projection-based eddy viscosity variational multiscale POD
2014 |  San and Iliescu (2014) |  Smagorinsky and Chollet-Lesieur spectral vanishing eddy viscosity models
2015 |  Stinis (2015); Chorin and Lu (2015); Li  _et al._ (2015) |  Mori–Zwanzig (MZ) formalism
2017 |  Gouasmi, Parish, and Duraisamy (2017) |  MZ ROM closures
2017 |  Rebollo  _et al._ (2017) |  Reduced basis methods for the Smagorinsky closure model
2017 |  Xie  _et al._ (2017) |  Approximate deconvolution reduced order modeling
2017 |  Benosman  _et al._ (2017) |  Lyapunov control theory to design learning-based closure models
2018 |  San and Maulik (2018a) |  Extreme learning machine closure model
2018 |  San and Maulik (2018b) |  Neural network closures for ROM
2018 |  Pan and Duraisamy (2018) |  Sparse polynomial regression and neural network for closure model
2019 |  Rahman, Ahmed, and San (2019) |  Dynamic closure model based on a test (secondary) truncation approach
2019 |  Stabile  _et al._ (2019) |  Reduced order variational multiscale approach for turbulent flows
2020 |  Imtiaz and Akhtar (2020) |  Nonlinear closure model based on the Jacobian of the Galerkin model
2020 |  Reyes and Codina (2020) |  Variational multiscale ROMs
2020 |  Xie, Webster, and Iliescu (2020) |  Residual neural network closures
2020 |  Wang, Ripamonti, and Hesthaven (2020) |  Recurrent neural network closures
2021 |  Mou  _et al._ (2021a) |  Data-driven variational multiscale ROMs
2021 |  Gupta and Lermusiaux (2020) |  Neural closure models

###  IV.1 Functional closure models

The functional, or phenomenological, closure modeling investigations have been
largely focused on the concept of eddy viscosity, which is added to the
physical viscosity of the system to drain the excessive energy. This modeling
concept is inspired by Kolmogorov’s ideas Kolmogorov (1941, 1991) about the
energy spectrum and energy cascade. In this section, we outline several
functional (phenomenological) ROM closure modeling strategies centered around
the concept of eddy viscosity.

#### Mixing-length ROM closure.

The first mixing-length ROM closure model was proposed by Aubry  _et al._
(1988), who studied the wall region of the turbulent boundary layer and used a
simple generalization of the Heisenberg spectral model in homogeneous
turbulence to provide the eddy viscosity closure term. Specifically, the
authors assumed that the deviatoric component of the Reynolds stress
$\hat{\bm{\tau}}$ of the unresolved field (represented by truncated modes),
acting on the resolved field (i.e., retained modes), is proportional to the
strain rate $\mathbf{S}$ of the resolved field:

$\hat{\bm{\tau}}=-2\alpha\nu_{e}\mathbf{S},$ (43)

where $\nu_{e}$ is the eddy viscosity parameter and $\alpha$ is a
dimensionless adjustable parameter. Moreover, they expressed the eddy
viscosity term as a function of the eigenvalues and eigenfunctions of the
first neglected modes based on the assumption that the energy decreases
rapidly with increasing mode index. The adequacy of this model was
quantitatively validated using numerical simulations in Ref. Podvin (2001),
and further investigations have been performed by Podvin and Lumley (e.g., for
minimal flow unit Podvin and Lumley (1998) and channel flow Podvin (2009)) and
also in Refs. Wang  _et al._ (2012); Mou  _et al._ (2021b). The two main
drawbacks of the global eddy viscosity modeling approach are:

  * i.

This formulation is equivalent to using Navier-Stokes equations at a lower
Reynolds number.

  * ii.

In this formulation, a linear closure term models nonlinear turbulence
dynamics.

#### Smagorinsky ROM closure.

As an improvement of the mixing-length model in Ref. Aubry  _et al._ (1988),
the celebrated Smagorinsky model Smagorinsky (1963) developed for LES has been
utilized for the ROM closure problem. The eddy viscosity coefficient adapts in
time in the Smagorinsky ROM closure model, but not in the mixing-length ROM
closure model. Thus, the former is expected to be more accurate than the
latter. To our knowledge, the first Smagorinsky ROM closure model was proposed
in Noack, Papas, and Monkewitz (2002) (see also Borggaard  _et al._ (2008)).
Ullmann and Lang (2010) used the same Smagorinsky closure model of the
original FOM simulations for ROMs based on LES snapshots of the turbulent flow
around a circular cylinder. However, the eddy viscosity term does not appear
explicitly in their ROM equation. Instead, the reconstructed velocity field is
utilized to update the (spatially-varying) Smagorinsky eddy viscosity term in
the FOM space, which in turn updates the corresponding model coefficients at
each time step of the time integration of the ROM. Borggaard, Iliescu, and
Wang (2011) proposed the inclusion of an artificial viscosity term in the ROM
equation that resembles the one used in the Smagorinsky model even if the
original FOM (for data generation) does not involve such a term. Rebollo  _et
al._ (2017) investigated the Smagorisnky ROM closure model in a reduced basis
method (RBM) setting. A rigorous numerical analysis of the Smagorinsky ROM
closure model was performed in Borggaard, Iliescu, and Wang (2011), where
error estimates for the time discretization were proven. To our knowledge,
this represents the first numerical analysis of ROM closures. Error estimates
for the time and space discretizations of the Smagorinsky ROM closure model
were later proven in Rebollo  _et al._ (2017) in an RBM context.

#### Dynamic SGS ROM closure.

The dynamic SGS model Germano  _et al._ (1991) is the state-of-the-art closure
model in LES. The main improvement in the dynamic SGS model over the standard
Smagorinsky model is that it uses an eddy viscosity coefficient that is
updated in time by using a secondary filtering operation. The dynamic SGS
closure model was extended for the first time to a ROM setting by Borggaard
_et al._ (2008) and was later investigated by Wang  _et al._ (2012) in the
numerical simulation of a 3D flow past a cylinder, where it yielded
significantly more accurate results than both the mixing-length and the
Smagorinsky ROM closure models. A more efficient numerical discretization of
the dynamic SGS ROM closure model was proposed by Rahman, Ahmed, and San
(2019).

#### Mode-dependent eddy viscosity closure.

Rather than adopting a single global eddy viscosity value $\nu_{e}$ for all
the modes (as in the mixing-length ROM closure model), a mode-dependent eddy
viscosity was proposed by Rempfer (1991) and Rempfer and Fasel (1993) to use a
different amount of dissipation for each scale. In Refs. Rempfer (1991);
Rempfer and Fasel (1993), the effective viscosity is calculated by requiring
the energy variation of different modes in ROM to match the energy variation
of the coherent structures in FOM. A modification to the mixing-length model
can be incorporated by introducing a mode-dependent kernel. The importance of
such a mode dependent kernel was first stressed by Rempfer Rempfer (1991);
Rempfer and Fasel (1993, 1994). Sirisup and Karniadakis (2004) applied a
vanishing viscosity kernel, which adds a small amount of mode-dependent
dissipation that satisfies the entropy condition, yet retains spectral
accuracy. The intrinsic stabilization scheme proposed in Ref. Kalb and Deane
(2007) utilizes information from available snapshots and POD modes to define a
mode-dependent stabilization. In Ref. San and Iliescu (2014), linear,
quadratic, and square-root kernels were investigated for the 1D Burgers
problem.

#### Variational multiscale eddy viscosity ROM closure.

The variational multiscale (VMS) methods developed by Hughes and his
collaborators Hughes  _et al._ (1998); Hughes, Mazzei, and Jansen (2000);
Hughes, Oberai, and Mazzei (2001) have made a significant impact in classical
CFD. The VMS methods center around the principle of locality of energy
transfer, which states that energy is transferred mainly between neighboring
scales or modes. Since ROMs use hierarchical bases in which the large and
small structures are clearly displayed, the VMS framework was naturally
extended to the ROM setting. Next, we present some of the eddy viscosity ROM
closure models developed in a VMS framework. Borggaard  _et al._ (2008)
proposed the first VMS ROM closure model, which was later investigated in Wang
_et al._ (2012) in the numerical simulation of a 3D flow past a circular
cylinder. The VMS ROM in Refs. Borggaard  _et al._ (2008); Wang  _et al._
(2012) used a _three-scale_ decomposition of the flow field into resolved
large, resolved small, and unresolved scales, and employed the Smagorinsky
model to dissipate energy only from the resolved small scales. A _two-scale_
decomposition of the flow field into resolved and unresolved scales was used
by Bergmann, Bruneau, and Iollo (2009b) to develop a VMS ROM closure model
with a residual based eddy viscosity component. Iliescu and Wang (2013, 2014)
put forth a three-scale VMS ROM closure model, in which the ROM projection was
used to construct an eddy viscosity term that acts only on the small resolved
scales.

A similar three-scale VMS ROM was proposed by Roop (2013) for a generalized
Oseen problem. Eroglu, Kaya, and Rebholz (2017) developed a different three-
scale VMS ROM that uses the ROM projection to add an eddy viscosity term
acting only on the small resolved scales in a modular fashion. This VMS ROM
was successfully tested in the numerical simulation of a turbulent channel
flow at $Re_{\tau}=395$ Eroglu, Kaya, and Rebholz (2017) and was extended by
Eroglu, Kaya, and Rebholz (2019) to the Darcy-Brinkman equations with double
diffusive convection. Stabile  _et al._ (2019) proposed a two-scale residual-
based VMS ROM closure model in which the VMS strategy is used at both the FOM
and the ROM levels to ensure model consistency. Reyes and Codina (2020) (see
also Ref. Reyes (2020)) developed a two-scale VMS ROM which is equipped with
time-dependent orthogonal subgrid scales that leverage the orthonormal nature
of the POD basis. Two-scale VMS-ROMs based on orthogonal subgrid scales were
used by Reyes  _et al._ (2018) for thermally coupled low Mach flows and by
Tello, Codina, and Baiges (2020) for a fluid structure interaction problem.
The first numerical analysis of VMS ROMs was performed in Iliescu and Wang
(2013, 2014), where stability and convergence were rigorously proven.
Numerical analysis of VMS ROMs was also performed by Roop (2013) and Eroglu,
Kaya, and Rebholz (2017). We also refer to the recent studies by Rubino and
his coworkersRubino (2020); Azaïez, Rebollo, and Rubino (2019) for multi-stage
ROM stabilization approaches in advection-dominated problems.

#### Finite-time thermodynamics ROM closure.

In the majority of the aforementioned studies, the closure term eventually
appears as a linear term in the GROM (i.e.,
$\mathbf{\mathfrak{C}}=\mathbf{\tilde{\mathfrak{B}}}+\mathbf{\tilde{\mathfrak{L}}}\bm{a}$).
(The Smagorinsky and dynamic SGS ROM closure models are notable exceptions.)
Noack, Morzynski, and Tadmor (2011) highlighted that energy transfer is
actually caused by nonlinear mechanisms. Thus, they introduced a nonlinear
eddy viscosity term $\nu_{e}(\bm{a})$ that is state-dependent. A finite-time
thermodynamics (FTT) Noack  _et al._ (2008) approach was utilized to quantify
the nonlinear eddy viscosity by matching the modal energy transfer effect as
follows:

$\nu_{e}(\bm{a})=\nu_{0}\sqrt{\dfrac{K(t)}{\bar{K}}},$ (44)

where $K(t)=\sum_{i=1}^{r}\dfrac{1}{2}a_{i}(t)^{2}$ represents the total
turbulence kinetic energy resolved by the Galerkin expansion and $\bar{K}$
denotes its time-averaged value. This led to damping levels more consistent
with energy fluctuations than those defined by a linear eddy viscosity model.
The FTT-based nonlinear eddy viscosity with an energy-based scaling model was
successfully applied to a 3D turbulent jet Schlegel  _et al._ (2009) and a 2D
mixing layer Cordier  _et al._ (2013). It was further extended to a mode-
dependent nonlinear eddy viscosity for a high Reynolds number flow over a
square-back Ahmed body Östh  _et al._ (2014).

#### Efficient numerical discretization of ROM closures.

Although the eddy viscosity closure models discussed in this section can
significantly improve the ROM accuracy, their brute-force numerical
discretization can be extremely inefficient. For example, the Smagorinsky ROM
closure model depends on the Frobenius norm of the deformation tensor, which
is a non-polynomial nonlinearity that cannot be preassembled in the offline
stage. Thus, alternative, efficient numerical discretizations have been
proposed over the last decade, which we outline next. Wang  _et al._ (2011)
proposed a two-level method to avoid the brute-force discretization of the
closure term onto the FOM fine mesh. Specifically, the POD bases constructed
from the original fine grid resolution snapshot data were interpolated onto a
coarse grid, and then they were used to efficiently compute the ROM closure
term. To avoid the assembly of the FOM strain rate tensor at each time step,
Akhtar  _et al._ (2012) used the Jacobian of the GROM right-hand side as an
eddy viscosity coefficient. A pre-computed eddy viscosity approach was adopted
by San and Iliescu (2015) by simplifying the nonlinear interaction in the
original Smagorinsky model. San and Iliescu (2014) also explored various
closure approaches including constant, polynomial, and spectral vanishing
viscosity models. Rebollo  _et al._ (2017) were the first to use an efficient
hyper-reduction method Yano (2019) (i.e., EIM Barrault  _et al._ (2004b)) to
discretize the Smagorinsky ROM closure in an RBM setting.

###  IV.2 Structural closure models

Structural closure models are generally derived through mathematical rather
than phenomenological arguments. This often includes a filtering procedure,
where the filtered field is assumed to have larger spatial structures than
those in the original one. Therefore, the filtered flow variables require
fewer modes in the ROM approximation. In other words, for the same number of
modes, ROM is capable of approximating the filtered flow field more accurately
than the unfiltered field. This approach is similar to LES, where the filtered
flow variables can be approximated on the given coarse mesh more accurately
than the original unfiltered flow variables. In this section, we survey ROM
closure models developed by using different types of ROM filtering.

#### Spatial filtering: projection.

Given the hierarchical nature of the ROM basis, not surprisingly, the most
popular type of ROM filtering has been the ROM projection, i.e., the
projection of various (nonlinear) terms living in the $r$-dimensional ROM
space spanned by the first $r$ ROM basis functions onto a smaller,
$s$-dimensional ROM space spanned by the first $s$ ROM basis functions, where
$s<r$. A classical example of ROM closure models constructed by using the ROM
projection is the VMS-ROMs(Stabile  _et al._ , 2019; Reyes and Codina, 2020;
Reyes, 2020; Mou  _et al._ , 2021c), which are discussed in Section IV.1. The
ROM projection, however, has been used to develop other types of ROM closures.
For example, the ROM projection has been utilized to construct parametrized
manifolds ROM closures (Chekroun, Liu, and Wang, 2014, 2015), which are based
on dynamical systems approaches, e.g., approximate inertial manifolds. The ROM
projection has also been used to build ROM closures based on stochastic
dynamical systems ideas (Majda, Harlim, and Gershgorin, 2010; Majda and
Harlim, 2012; Harlim, Mahdi, and Majda, 2014).

#### Spatial filtering: differential filter and approximate deconvolution.

Using the analogy between LES and ROM, we mention that a lot of ideas and
techniques in image and signal processing are also applicable in ROM, and vice
versa! In LES, the approximate deconvolution (AD) represents one of the most
popular techniques in this class. It is based on the deconvolution approaches
developed in the image processing and inverse problems communities to recover
the original signal from a blurred filtered signal.

In stark contrast to the abundance of functional closure studies (beginning in
the 1980s), there are only a few structural closure models in ROM literature.
The AD-ROM was proposed by Xie  _et al._ (2017) for the three-dimensional flow
past a circular cylinder. To construct the AD-ROM, a ROM differential filter
is applied to the Navier-Stokes equations, followed by a Galerkin projection
of the _filtered_ equations. It is usually assumed that the filtering and
differentiation operators commute, and the repercussions of this assumption
are investigated in Ref. Koc  _et al._ (2019). Nonetheless, it is well known
that, in general, nonlinearity and ROM spatial filtering do not commute.
Therefore, the resulting equations include a _filtered_ nonlinear term of the
_unfiltered_ variables (i.e., $\widehat{\mathcal{N}(\mathbf{u})}$, where the
hat operator denotes the filtering process), rather than a nonlinear operator
of the _filtered_ variables (i.e., $\mathcal{N}(\hat{\mathbf{u}}))$. A
regularized deconvolution is adopted to provide the ROM approximation of the
unfiltered flow variables in order to compute the nonlinear term. Thus, the
filtering process increases the accuracy of the ROM in the sense that the
filtered field contains larger spatial structures, and thus can be
sufficiently captured by the ROM approximation. In addition, the AD technique
solved the ROM closure problem by providing an estimate of the unfiltered flow
variables.

######  Remark.

We note that ROM spatial filtering has also been used to develop regularized
ROMs (Reg-ROMs), i.e., ROMs in which spatial filtering is used to smoothen
(regularize) various terms in the Navier-Stokes equations and increase the
numerical stability of the ROM. We emphasize that, while related,
regularization and closure are different: The latter adds a closure term,
whereas the former usually does not. ROM spatial filtering has been used to
develop various types of Reg-ROMs: Wells  _et al._ (2017) proposed, for the
first time, an evolve-then-filter approach in which the GROM (Eq. 33) is
integrated (evolved) for one time step, after which a ROM spatial filter is
applied to filter the intermediate solution obtained in the evolve step. This
filtering reduces the numerical oscillations of the flow variables (i.e., adds
numerical stabilization to the ROM). Gunzburger  _et al._ (2019) proposed an
evolve-filter-relax approach that considers the additional step of relaxation,
which averages the unfiltered and filtered flow variables to control the
amount of numerical dissipation introduced by the filter. Recently, Girfoglio
et al. Girfoglio, Quaini, and Rozza (2021a, b) have investigated the evolve-
filter approach in a finite volume setting. The ROM differential filter has
also been used to develop the Leray Reg-ROM in Refs. Sabetghadam and Jafarpour
(2012); Xie  _et al._ (2018a); Gunzburger, Iliescu, and Schneier (2020).

#### PDF filtering: Mori–Zwanzig formalism and memory effects.

A different type of filtering, based on filtering with respect to the PDF of
the initial conditions, has been instrumental in adding memory effects to ROM
closures, with rationale based on the Mori–Zwanzig (MZ) formalism Mori (1965);
Zwanzig (1980); Chorin, Hald, and Kupferman (2000). More recently, the MZ
formalism has been intensely used to define closures for both LES and ROM
settings. Next, we outline some of these developmemnts. Stinis (2015)
introduced a generalized MZ framework for the construction of ROMs for systems
without scale separation. For example, assume that Eq. 40 defines the
following linear system Zwanzig (2001); Gouasmi, Parish, and Duraisamy (2017);
Pan and Duraisamy (2018); Wang, Ripamonti, and Hesthaven (2020):

$\frac{d}{dt}\begin{bmatrix}\bm{a}\\\
\bm{b}\end{bmatrix}=\begin{bmatrix}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{bmatrix}\begin{bmatrix}\bm{a}\\\ \bm{b}\end{bmatrix}.$ (45)

The evolution of the unresolved state $\bm{b}$ can be evaluated as follows
(assuming that $\bm{a}$ is known) Zwanzig (2001); Gouasmi, Parish, and
Duraisamy (2017); Wang, Ripamonti, and Hesthaven (2020):

$\bm{b}(t)=\int_{0}^{t}e^{A_{22}(t-s)}A_{21}\bm{a}(s)ds+e^{A_{22}t}\bm{b}(0).$
(46)

Therefore, Eq. 41 for the dynamics of $\bm{a}$ can be written as

$\displaystyle\frac{d\bm{a}}{dt}$ $\displaystyle=A_{11}\bm{a}+A_{12}\bm{b}$
$\displaystyle=A_{11}\bm{a}$
$\displaystyle+A_{12}\int_{0}^{t}e^{A_{22}(t-s)}A_{21}\bm{a}(s)ds$
$\displaystyle+A_{12}e^{A_{22}t}\bm{b}(0).$ (47)

Eq. 47 expresses the dynamics of the resolved scales using a Markovian term
(i.e., $A_{11}\bm{a}$, which depends only on the current value of $\bm{a}$), a
memory integral term depending on the history of the resolved scales $\bm{a}$,
and a term describing the contribution of the initial conditions. This
derivation can be extended to nonlinear settings where, for instance, the
nonlinear ordinary differential equation can be written as a linear partial
differential equation using the Liouville operator. The exact evolution
equations for the reduced state can be written as

$\displaystyle\frac{d\bm{a}}{dt}$ $\displaystyle=\bm{f}_{\bm{a}}(\bm{a},0)$
$\displaystyle+\int_{0}^{t}\mathcal{K}(\bm{a}(s),t-s)ds+\mathcal{O}(\bm{a}(0),\bm{b}(0)),$
(48) $\displaystyle\approx\bm{f}(\bm{a})+\text{Closure}(\bm{a}).$ (49)

In Eq. 48, $\mathcal{K}$ is called the memory kernel and $\mathcal{O}$
designates the contribution from the initial conditions. The memory integral
term implies that the accurate resolution of $\bm{a}$ comprises a non-
Markovian contribution. However, the direct computation of Eq. 48 is generally
prohibitive, and estimation of the memory-effect is often sought.

Li  _et al._ (2015) included a great discussion on incorporation of memory
effects in coarse-grained modeling via the MZ formalism. A discrete approach
to stochastic parametrization, dimension reduction, and their connections to
the MZ formalism of statistical physics has been proposed by Chorin and Lu
(2015). In an LES setting, Parish and Duraisamy (2017a) framed the MZ closure
modeling approach by exploiting similarities between two levels of coarse-
graining via the Germano identity of fluid mechanics and by assuming that
memory effects have a finite temporal support. The concept has been also
generalized to provide a mathematically consistent framework for the
construction of ROMs of dynamical systems Gouasmi, Parish, and Duraisamy
(2017). Moreover, Parish and Duraisamy (2017b) established an analogy between
MZ and VMS approaches.

#### Ensemble averaging.

Noack, Morzynski, and Tadmor (2011) first used ensemble averaging to construct
a finite-time thermodynamics (FTT) Noack  _et al._ (2008) framework.
Gunzburger, Jiang, and Schneier (2017) built ensemble-based POD ROMs, where
the nonlinear advection term in the Navier-Stokes equations is replaced by a
linear term in the equations for the resolved scales. This _linearization_ is
performed by using an ensemble of solution trajectories by propagating an
ensemble of ROMs with varying parameters and/or initial conditions and
updating the ensemble average at each time step. Later on, this ensemble-based
approach was equipped with Leray regularization to develop regularized ROMs
for high values of Reynolds number (Gunzburger, Iliescu, and Schneier, 2020).

#### Time averaging.

Selten (1995, 1997) used time averaging to develop ROM closures. In
particular, by estimating the rate at which the ROM trajectory drifts away
from the projection of the FOM solution on the ROM subspace, Selten (1995)
added a linear damping to expand the doubling-time of the error resulting from
the modal truncation. Berselli  _et al._ (2020) developed mathematical support
for eddy viscosity modeling of time-averaged ROM closures. While being
interested in a statistical equilibrium problem exploring possible forward and
backward average transfer of energy among ROM basis functions, they found that
the time-averaged energy exchange from low index POD modes to high index POD
modes is positive for long enough time intervals. This study provides, for the
first time, mathematical support for the ROM eddy viscosity methodology, where
the energy transfer to the truncated modes is modeled by employing extra
viscous dissipation.

#### Calibrating the POD space with a Navier-Stokes based side constraint.

The last modeling approach that we discuss in this section is that proposed by
Balajewicz, Dowell, and Noack (2013). Although this approach does not add a
ROM closure model, it does leverage mathematical arguments to model the effect
of the truncated modes. In this approach, the POD subspace is subjected to a
Navier-Stokes based side constraint. Specifically, the power balance for the
fluctuation energy is required to be satisfied by the attractor data after the
Galerkin projection on the adjusted POD space. This procedure can be
conceptualized as rotating the POD subspace into a more dissipative regime, in
which the extra dissipation is now performed by more dissipative POD modes.

###  IV.3 Stochastic closure models

Although we are mainly focusing on deterministic closure modeling in this
review, we emphasize that the need for stochastic modeling was already
formulated in Aubry  _et al._ (1988) to avoid statistically nonstationary
behaviour for some homoclinic orbits. The dynamics of the unresolved scales,
and hence their interactions with the resolved scales are unknown. Thus, we
can only form an approximate idea of how the truncated modes behave and affect
the ROM solution. Even with the best closure model, we can never be certain
about its accuracy in practical settings. Thus, it is natural to model the
dynamics of the unresolved modes using a random or _stochastic_ process, from
which we can infer the unresolved modes’ contribution to the evolution of
large scales in a statistical sense. We refer to Leith (1996); Chorin and Hald
(2009); Majda, Harlim, and Gershgorin (2010); Majda and Harlim (2012); Harlim,
Mahdi, and Majda (2014); Majda (2015); Resseguier, Mémin, and Chapron (2015);
Chorin and Lu (2015); Lu, Lin, and Chorin (2017); Lu (2020); Sieber,
Paschereit, and Oberleithner (2021); Chekroun, Liu, and Wang (2014) for
detailed discussions on the probabilistic modeling of such random/chaotic
systems as well as the development of statistically accurate ROMs and
stochastic closure models Wilks (2005); Sapsis and Lermusiaux (2009); Sapsis
(2013); Sapsis and Majda (2013a, b); Arnold, Moroz, and Palmer (2013); Pulido
_et al._ (2018); Chattopadhyay, Subel, and Hassanzadeh (2020); Gagne  _et al._
(2020). We also note that nonparametric stochastic modeling approaches have
been proposed for representative stochastic Itô drift diffusion forecast
models Berry, Giannakis, and Harlim (2015). Next, we briefly outline a few of
these strategies.

Stochastic closure approaches seek to account for the effects of the
unresolved scales on the long-term statistics of the resolved scales. In
particular, the closure term is modeled by a stochastic process, usually
represented by Markovian and/or non-Markovian dynamics with a random forcing
(e.g., random noise). For a truncated ROM of the Kuramoto-Sivashinsky system,
Lu, Lin, and Chorin (2017) defined a discrete-time closure term $z^{n}$ at
time $t_{n}$ as follows:

$z^{n}=\Phi^{n}+\xi^{n},$ (50)

where $\xi$ is a sequence of independent identically distributed random
variables, which are sampled from Gaussian distributions and characterize the
stochastic component of the closure, while $\Phi$ is a function of current and
past values of the resolved scales $\bm{a}$ and the forcing $\xi$. The authors
used the nonlinear autoregression moving average with exogenous input (NARMAX)
approach to parameterize $\Phi$. A similar approach was adopted in Ref.Chorin
and Lu (2015) for the Lorenz 96 model and Ref.Lu (2020) for the stochastic
Burgers equation. The multiscale Lorenz 96 model Lorenz (1996), which has been
considered as a non-trivial test problem for stochastic paramerization in
geophysical fluid dynamics studies, can be written as

$\displaystyle\frac{dX_{i}}{dt}$
$\displaystyle=-X_{i-1}(X_{i-2}-X_{i+1})-X_{i}-\frac{hc}{b}\sum_{j=1}^{J}Y_{j,i}+F,$
(51) $\displaystyle\frac{dY_{j,i}}{dt}$
$\displaystyle=-cbY_{j+1,i}(Y_{j+2,i}-Y_{j-1,i})-cY_{j,i}+\frac{hc}{b}X_{i},$
(52)

where Eq. 51 represents the evolution of slow, high-amplitude variables
$X_{i}\leavevmode\nobreak\ (i=1,\dots,I)$, and Eq. 52 describes the evolution
of coupled fast, low-amplitude variables $Y_{j,i}\leavevmode\nobreak\
(j=1,\dots,J)$. In order to investigate different closure approaches, $X$ can
be considered as the resolved scales, while $Y$ can be considered as the
unresolved ones. Therefore, Eq. 52 is assumed to be unknown and is only used
for generating _true_ data. Furthermore, the form of the term
$-\frac{hc}{b}\sum_{j=1}^{J}Y_{j,i}$, representing the contribution of $Y$ to
the dynamics of $X$, is also assumed to be unknown. A closure term is
parameterized as a function of the resolved scales. Although the Lorenz 96
equations are determinstic, Wilks (2005) showed the existence of multiple
closure values that are consistent with any given large-scale variable, i.e.,
that different values of the closure term yield _statistically_ similar
results. Therefore, Wilks (2005) defined the closure term using both
deterministic and stochastic components. Specifically, the author used a
fourth-order polynomial fitting for the deterministic part that represents the
average trend, and a first-order auto-regression model for the stochastic part
that defines the deviation of different realizations from the fitted curve.
Arnold, Moroz, and Palmer (2013) explored several parametrization schemes for
the stochastic component, including additive and multiplicative noise.

Although the distinction between resolved and unresolved variables in the
multiscale Lorenz 96 system is not driven by a Galerkin truncation as is the
case for most projection-based ROMs (which is the focus of the current
review), the same arguments apply in both scenarios. For example, Mémin (2014)
assumed that the flow field is decomposed into a deterministic resolved
component and a generalized random field that models the unresolved flow
component and all the uncertainties in the flow. In other words, the
projection of the flow field onto the truncated space is treated as a
_realization_ or sample of the stochastic component of the flow and is modeled
using Brownian motion. Resseguier, Mémin, and Chapron (2015) provided
numerical investigations of this methodology using POD-Galerkin projection for
flow past cylinder. Nonetheless, we believe that this is an open research area
that needs fresh ideas to translate statistical closure strategies Zhou (2021)
from turbulence modeling to the ROM arena.

##  V Data-driven closure modeling

With the abundant supply of big data, open-source cutting edge and easy-to-use
machine learning libraries, cheap computational infrastructure, and high
quality, readily available training resources, data-driven closure modeling
has become very popular. Since projection-based ROMs are usually constructed
from snapshot data (either collected experimentally or computationally), it is
natural to further exploit this set of data to estimate the closure term
efficiently.

In this section, we survey data-driven closure modeling approaches in which
the closure problem is cast into a regression task. The majority of the recent
data-driven closure studies can be viewed as a regression task, where the
closure model is defined partially or completely as a function of available
information (e.g., resolved scales). We first discuss the early investigations
based on classical least-squares approaches. Then, we introduce ML tools that
perform this regression task using neural networks and Gaussian processes
regression. Finally, we explore legacy data assimilation, parameter
estimation, and system identification tools that can efficiently be used to
solve the closure problem. We note that some subsections are entirely devoted
to closure modeling (Sections V.1, V.2), some only partially address closure
modeling (Sections V.3, V.5), and some do not address closure modeling at all
(Sections V.4, V.6, V.7). Although the approaches in the latter subsections
have not yet been used for closure modeling, we believe that they will soon
make an impact in this dynamic research field.

###  V.1 Trajectory regression vs. model regression

There are two main schools of thought in data-driven ROM closure modeling: (i)
trajectory regression, and (ii) model regression.

The trajectory regression approach aims at finding the ROM closure model
$\mathbf{\mathfrak{C}}$ that yields the best ROM trajectory. In this approach,
the following constrained regression problem is solved:

$\displaystyle\begin{aligned} &\underset{\mathbf{\mathfrak{C}}\text{
parameters}}{\text{minimize}}&&\|\bm{a}^{\small\rm ROM}-\bm{a}^{\small\rm
FOM}\|^{2},\\\ &\text{subject to}&&\bm{a}^{\small\rm
ROM}\text{\leavevmode\nobreak\ solves closed ROM (Eq.\leavevmode\nobreak\
\ref{eq:gromC})}\end{aligned}$ (53)

where $\bm{a}^{\small\rm FOM}$ is the vector of ROM coefficients computed with
the FOM data and $\bm{a}^{\small\rm ROM}$ is the vector of ROM coefficients
yielded by the closed ROM.

The model regression approach aims at finding the ROM closure model
$\mathbf{\mathfrak{C}}$ that yields the best ROM closure model. In this
approach, the following unconstrained regression problem is solved:

$\displaystyle\begin{aligned}
&\underset{\textbf{c}}{\text{minimize}}&&\|\mathbf{\mathfrak{C}}^{\small\rm
ROM}-\mathbf{\mathfrak{C}}^{\small\rm FOM}\|^{2},\end{aligned}$ (54)

where $\mathbf{\mathfrak{C}}^{\small\rm FOM}$ is the ROM closure model
computed with the FOM data, $\mathbf{\mathfrak{C}}^{\small\rm ROM}$ is the
postulated ROM closure model form, and c is the vector of parameters used to
define the ROM closure model form. We note that the model regression approach
is similar in spirit to the a priori testing used in LES Sagaut (2006).

We emphasize that the two approaches are fundamentally different. The
trajectory regression is a black-box approach in which the precise formula for
the ROM closure term, $\mathbf{\mathfrak{C}}$, is not actually known. Instead,
the trajectory regression first postulates a model form for
$\mathbf{\mathfrak{C}}$ (e.g., by using one of the functional models in
Section IV.1), and then solves the constrained optimization problem Eq. 53 to
find the closure model parameters that yield the most accurate ROM trajectory
(i.e., the trajectory $\bm{a}^{\small\rm ROM}$ that is closest to the
projection of the FOM data on the ROM basis, $\bm{a}^{\small\rm FOM}$). In
contrast, the model regression first employs one of the filters described in
Section IV.2 to determine a precise formula for the ROM closure term,
$\mathbf{\mathfrak{C}}^{\small\rm FOM}$. Then, it postulates a model form for
$\mathbf{\mathfrak{C}}^{\small\rm ROM}$. Finally, it solves the unconstrained
optimization problem Eq. 54 to find the closure model parameters that yield
the most accurate ROM closure model (i.e., the closure model
$\mathbf{\mathfrak{C}}^{\small\rm ROM}$ that is closest to the “true” closure
model, $\mathbf{\mathfrak{C}}^{\small\rm FOM}$, computed from FOM data).

There are pros and cons for both approaches. The trajectory regression is
conceptually simpler than the model regression since it does not need to
determine the actual form of the ROM closure term. The trajectory regression
is also more flexible than the model regression since it can model not only
the ROM closure term, but also other sources of ROM uncertainty, such as the
numerical discretization error and the missing data. Finally, according to
Noack’s conjecture Noack, Papas, and Monkewitz (2005), the model regression is
more accurate in the predictive regime (i.e., outside the training interval),
whereas the trajectory regression is more accurate in the reconstructive
regime (i.e., inside the training interval). Noack motivated his conjecture by
noting that using data to match models appears more robust to perturbations
than using data to match trajectories. To our knowledge, Noack’s conjecture
has not been investigated numerically.

###  V.2 Least-squares regression

In this section, we survey the data-driven ROM closure models that use a
least-squares formulation in the optimization problems given by Eq. 53 and Eq.
54.

#### Trajectory regression.

Given its conceptual simplicity, the trajectory regression has been used from
the earliest days of reduced order modeling to develop closures. The general
idea used to develop these closure models is simple: (i) Postulate a ROM
closure model form, either functional (such as the eddy viscosity models
surveyed in Section IV.1) or structural (such as the models surveyed in
Section IV.2). (ii) Use a least-squares problem in Eq. 53 to determine the
various parameters in the postulated model form.

_Functional models_ : Probably the first functional trajectory regression
closure is the mixing-length model proposed in the pioneering work by Aubry
_et al._ (1988), in which trajectory regression is used to determine the
mixing-length constant (see Wang  _et al._ (2012) for related work). A least-
squares trajectory regression was also used by Wang  _et al._ (2012) to
determine the eddy viscosity constants in the Smagorinsky and VMS closure
models. Further improvements to the least-squares trajectory regression of
eddy viscosity ROM closure models were proposed by Östh  _et al._ (2014) and
Protas, Noack, and Östh (2015). The eddy viscosity trajectory regression has
also been used by Stabile  _et al._ (2019); Reyes and Codina (2020), and
Bergmann, Bruneau, and Iollo (2009b).

_Structural models_ : Probably the first structural trajectory regression
closures are those constructed with _calibration_ methods, which have been
introduced to directly modify (or calibrate) the GROM polynomial coefficients
(i.e., $\mathbf{{\mathfrak{B}}}$, $\mathbf{{\mathfrak{L}}}$, and
$\mathbf{{\mathfrak{N}}}$), rather than introducing an additional closure
term. Buffoni  _et al._ (2006) calibrated the constant and linear terms (i.e.,
$\mathbf{{\mathfrak{B}}}$ and $\mathbf{{\mathfrak{L}}}$), while leaving the
nonlinear term, $\mathbf{{\mathfrak{N}}}$, as derived from the Galerkin
projection step. The modified coefficients are then found using a pseudo-
spectral method such that the model prediction is as close as possible to the
actual reference solution Galletti  _et al._ (2007). An extension to calibrate
_all_ the polynomial coefficients (linear and quadratic) was employed by
Couplet, Basdevant, and Sagaut (2005), where the cost function is defined to
penalize the deviation of the calibrated ROM behavior with respect to the
projection of true snapshot data. Moreover, Perret, Collin, and Delville
(2006) considered a cubic polynomial to represent the dynamics of the POD
modal amplitudes for supersonic jet-mixing layer data, and adopted a least-
squares regression to define the coefficients of the polynomial. Baiges,
Codina, and Idelsohn (2015) used a calibration method in a VMS framework to
build a structural trajectory regression closure. Specifically, the authors
postulated a linear model for the unresolved sub-scale term as a function of
the resolved field (i.e.,
$\mathbf{\mathfrak{C}}=\mathbf{\tilde{\mathfrak{B}}}+\mathbf{\tilde{\mathfrak{L}}}\bm{a}$),
and then solved the constrained least-square problem Eq. 53 for the components
of $\mathbf{\tilde{\mathfrak{B}}}$ and $\mathbf{\tilde{\mathfrak{L}}}$.

An assessment of various calibration techniques using the two-dimensional ﬂow
around a cylinder can be found in Ref. Cordier, Abou El Majd, and Favier
(2010). A vital merit of the calibration methods is that the computational
costs of these methods are reasonable since they employ the temporal part of
the POD information for the regression task, while the Galerkin projection
method exploits the much more voluminous spatial POD information to construct
the ROM polynomial.

#### Model regression.

These closure models are constructed as follows: (i) Use the filters surveyed
in Section IV.2 to determine a mathematical formula for the “true” closure
model, $\mathbf{\mathfrak{C}}^{\small\rm FOM}$, computed from FOM data. (ii)
Postulate a ROM closure model form for the ROM closure model,
$\mathbf{\mathfrak{C}}^{\small\rm ROM}$. (iii) Use the least-squares problem
in Eq. 54 to determine the various parameters in the postulated model form. To
our knowledge, the vast majority of ROM closure model forms that have been
proposed in this direction are of structural type. A notable exception is the
model proposed by Hijazi  _et al._ (2020), which uses finite volume RANS data
in conjunction with a model regression approach to determine the eddy
viscosity component of the ROM.

A model regression closure that uses the ROM projection as a spatial filter is
the data-driven VMS-ROM model proposed by Mou  _et al._ (2021c). A two-scale
version was investigated in Xie  _et al._ (2018b); Mou  _et al._ (2020), and a
three-scale version was proposed in Mou  _et al._ (2021c). Linear Mou  _et
al._ (2020), quadratic Xie  _et al._ (2018b); Mou  _et al._ (2021c), and even
cubic Mou (2018) terms were used for the model form:
$\mathbf{\mathfrak{C}}=\mathbf{\tilde{\mathfrak{B}}}+\mathbf{\tilde{\mathfrak{L}}}\bm{a}+\bm{a}^{\top}\mathbf{\tilde{\mathfrak{N}}}\bm{a}+\bm{a}^{\top}(\bm{a}^{\top}\mathbf{\tilde{\mathfrak{N}}}\bm{a})$.
A data-driven VMS-ROM to increase the pressure accuracy was proposed in
Ivagnes (2021). The verifiability of the data-driven VMS-ROM was proven by Koc
et al Koc  _et al._ (2021b). Other model regression closures that use the ROM
projection as a spatial filter were employed to build parameterized manifold
closures by Liu and his collaborators Chekroun, Liu, and Wang (2014) , and by
Lu and his coworkers Lu (2020); Lin and Lu (2021); Lu, Lin, and Chorin (2017)
to construct stochastic ROM closures. A model regression closure that uses the
ROM differential filter was proposed by Koc  _et al._ (2019). The resulting
data-driven LES-ROM uses the ROM differential filter to determine a
mathematical formula for the “true” closure model,
$\mathbf{\mathfrak{C}}^{\small\rm FOM}$, computed from FOM data. To our
knowledge, the data-driven LES-ROM in Koc  _et al._ (2019) is the only model
regression approach that utilizes a spatial filter (i.e., the differential
filter) instead of the commonly used ROM projection.

A model regression closure that uses the PDF filtering in an MZ setting was
proposed by Duraisamy, Parish, and their collaborators Gouasmi, Parish, and
Duraisamy (2017); Parish, Wentland, and Duraisamy (2018, 2020). In the MZ
framework, PDF filtering is used to express the “true” closure model,
$\mathbf{\mathfrak{C}}^{\small\rm FOM}$, as a memory term, which is then
approximated by using the FOM data and solving the least-squares problem Eq.
54. Moreover, a model regression closure that uses time filtering was proposed
by Selten (1997) and utilized in the numerical simulation of a barotropic
model.

Mohebujjaman, Rebholz, and Iliescu (2019) used physical constraints in the
model regression closure to improve the stability and accuracy of the data-
driven VMS-ROM model Mou  _et al._ (2021c). Specifically, they equipped the
least-squares problem Eq. 54 with physical constraints to enforce the
regressed matrix and tensor to have similar characteristics as the GROM
operators (e.g., $\mathbf{\tilde{\mathfrak{L}}}$ being negative semi-definite
and $\mathbf{\tilde{\mathfrak{N}}}$ being energy conserving). A similar
quadratic formula was adopted in Ref. Pawar  _et al._ (2020a) to recover the
hidden physical processes (e.g., source terms) for system with incomplete
governing equations.

###  V.3 Neural network regression

The introduction of neural network regression into ROM was highly motivated by
the desire to construct purely data-driven nonintrusive ROM (NIROM) frameworks
Yu, Yan, and Guo (2019); Fresca, Dede, and Manzoni (2021), which solely rely
on data to learn the dynamics of the relevant solution manifold without the
need to access the governing equations. Non-intrusive approaches are
attractive due to their portability since they do not necessarily require the
exact form of the equations or the methods used to generate the data. In
addition, non-intrusive models offer a unique advantage in multidisciplinary
collaborative environments, where only data can be shared without revealing
the proprietary or sensitive information. Non-intrusive approaches are also
useful when the detailed governing equations of the problem are unknown. This
modeling approach can benefit from the enormous amount of data collected from
experiments, sensor measurements, and large-scale simulations to build NIROMs
based on the assumption that data is a manifestation of _all_ the underlying
dynamics and processes.

Machine learning tools, in particular artificial neural networks (ANNs)
equipped with the universal approximation theorem Csáji (2001), have been
widely exploited in this regard. A typical feed-forward neural network is
depicted in Figure 5, where a mapping $\bm{\mathcal{M}}$ from the input
$\bm{\mathcal{X}}$ to the output
$\bm{\mathcal{Y}}=\bm{\mathcal{M}}(\bm{\mathcal{X}})$ is inferred through a
learning algorithm. For transient flows, a single-layer feed-forward neural
network was proposed by San, Maulik, and Ahmed (2019) to provide accurate
predictions of the ROM coefficients with varying control parameter values,
using sequential and residual approaches. In the sequential approach, a
mapping from the current values of $\bm{a}$ to their future values is
approximated. Moreover, the input layer is augmented with the acting Reynolds
number and the time. That is,
$\bm{\mathcal{X}}=\\{\text{Re},t_{n},\bm{a}(t_{n})\\}$, while
$\bm{\mathcal{Y}}=\\{\bm{a}(t_{n+1})\\}$. On the other hand, the residual
implementation relies on learning the deviation of the future state from the
current values (i.e., $\bm{\mathcal{Y}}=\\{\bm{a}(t_{n+1})-\bm{a}(t_{n})\\}$).
Pawar  _et al._ (2019) employed deep neural networks (DNN) to bypass the
Galerkin projection step and build a fully NIROM for the two-dimensional
Boussinesq equations with a differentially heated cavity flow setup at various
Rayleigh numbers. In particular, the evolution of the POD modal amplitudes
$\bm{a}(t_{n+1})$ was predicted from their past values using residual and
backward difference scheme formulas. The application of variants of ANNs as
regression models for the dynamics of low-order states (e.g., POD amplitudes)
has gained substantial popularity Hesthaven and Ubbiali (2018); Wang,
Hesthaven, and Ray (2019); Wu  _et al._ (2020a), owing to the availability of
open-source and user-friendly ML libraries. This is a hot topic and dozens of
new papers appear every week in different journals and conferences all over
the world, dealing with different aspects of NIROM based on ANNs (e.g.,
different architectures, test bed problem, and error bounds).

Figure 5: A schematic diagram of a typical feed-forward neural network with an
input layer, hidden layers, and an output layer. We note that these general-
purpose dense deep network architectures have been evolving to more specific
neural network designs Goodfellow, Bengio, and Courville (2016). We discuss
some of them briefly in Section VI.

Recurrent neural networks (RNNs) are very effective for sequence predictions
in numerous applications, e.g., speech recognition and translation. RNNs
contain loops that allow them to retain information from one time step to
another so as to enforce the temporal dependencies. A deep residual recurrent
neural network was utilized by Kani and Elsheikh (2017) for the model
reduction of nonlinear dynamical systems. However, one of the limitations of
RNNs is vanishing (or exploding) gradient to capture the long-term
dependencies, stemming from the repetitive multiplication of gradient with
potentially ill-conditioned weight matrices during the back-propagation
learning algorithm. The long short-term memory (LSTM) neural networks mitigate
this issue by employing a gating mechanism that allows information to be
forgotten. Vlachas  _et al._ (2018) trained an LSTM to predict the derivative
of $\bm{a}$ with respect to time from a short history of $\bm{a}$ values,
where a first-order forward difference scheme was adopted to represent the
temporal derivative. They also combined the LSTM with a mean stochastic model
to cope with attractor regions that are not captured in the training set.
Mohan and Gaitonde (2018) explored the bidirectional variant of LSTM,
employing two LSTM networks: one in the forward and the other in the reverse
direction, for NIROM of forced isotropic turbulence and magneto-hydrodynamic
turbulence using the Johns Hopkins turbulence database Graham  _et al._
(2016). Rahman  _et al._ (2019) utilized LSTMs for the NIROM implementation of
the two-dimensional single-layer quasi-geostrophic ocean circulation model. A
sliding window approach was adopted to predict the evolution of the POD
amplitudes. In another interesting article, Wang, Ripamonti, and Hesthaven
(2020) utilized a conditioned LSTM for the memory term in the GROM equations,
representing the closure model, for parametric systems.

Instead of using ML to entirely replace the GROM with NIROM, data-driven ML
can be utilized along with the physics-based GROM to construct the closure
model. This hybrid approach was adopted in Ref. San and Maulik (2018a), which
used a dissipative term employing an eddy viscosity coefficient and utilized a
single layer extreme learning machine (ELM) to estimate a modal $\nu_{e}$ as a
function of the mode index, GROM right-hand side (RHS), and modal amplitudes.
Furthermore, a clipping procedure was carried out by discarding negative
values of $\nu_{e}$. ML tools can be exploited to provide the numerical value
of the closure term, without constraints on the functional form of the closure
model. San and Maulik (2018c, d) utilized an ELM Huang, Zhu, and Siew (2006)
to learn the value of the closure term as a function of the GROM RHS, i.e.,
$\mathbf{\mathfrak{C}}=f(\mathbf{\mathfrak{B}}+\mathbf{\mathfrak{L}}\bm{a}+\bm{a}^{\top}\mathbf{\mathfrak{N}}\bm{a})$.
For training purposes, the true closure term is computed from the projection
of the pure evolution of the PDE onto the reduced space. In other words, Eq.
29 is first evaluated in the FOM space, then projected onto the basis
functions, while the GROM RHS (Eq. 33) is computed directly in the reduced
space.

Wan  _et al._ (2018) utilized LSTM architectures to learn the mismatch
resulting from the imperfect GROM RHS as a function of the sequence of past
values of the resolved state. The ML model is exploited to assist the
imperfect model whenever data is available, while for locations with sparse
data, the GROM still provides an acceptable baseline for the prediction of the
system state. The estimation of the closure term as a function of the time
history of the resolved scales has roots in the Mori-Zwanzig formalism Chorin,
Hald, and Kupferman (2002); Chorin and Hald (2009) and the memory embedding of
LSTM implementation for closure modeling is supported by the Takens embedding
theorem Takens (1981). Gupta and Lermusiaux (2020) employed RNNs to learn the
non-Markovian term appearing in Eq. 48 using neural delay differential
equations with discrete-delays as follows:

$\frac{d\bm{a}}{dt}\approx\bm{f}(\bm{a})+g_{RNN}(\bm{a}(t),\bm{a}(t-\tau_{1}),\dots,\bm{a}(t-\tau_{K});\theta),$
(55)

where $K$ is the number of discrete-delays and $\theta$ represents the neural
network weights. For distributed-delays, the ROM dynamics is written as

$\frac{d\bm{a}}{dt}\approx\bm{f}(\bm{a})+g\bigg{(}\bm{a}(t),\int_{t-\tau_{2}}^{t-\tau_{1}}h(\bm{a}(\tau),\tau)d\tau,t\bigg{)},$
(56)

where the delay is distributed over past time periods $t-\tau_{2}$ and
$t-\tau_{1}$. Gupta and Lermusiaux (2020) approximated the $g$ and $h$
functions defining the delay term using two different coupled neural networks
as follows:

$\displaystyle\frac{d\bm{a}}{dt}$
$\displaystyle\approx\bm{f}(\bm{a})+g_{NN}(\bm{a}(t),y(t),t;\theta_{g})$ (57)
$\displaystyle\frac{dy}{dt}$ $\displaystyle\approx
h_{NN}(\bm{a}(t-\tau_{1}),t-\tau_{1};\theta_{h})-h_{NN}(\bm{a}(t-\tau_{2}),t-\tau_{2};\theta_{h}),$
(58)

where the memory effect can be embedded without the use of any specific
recurrent neural network. The authors showed that the non-Marovian closure
outperforms its Marovian counterpart (with no time delays).

###  V.4 Kernel regression

Gaussian process regression (GPR) Rasmussen (2003) has the advantage of the
simultaneous prediction of system’s dynamics and the associated uncertainty.
Although GPR provides a powerful tool for probabilistic inference that enables
modelers to strike a balance between model complexity and data fitting Raissi
and Karniadakis (2018), its use is often limited to relatively small training
data sets due to the well-known cubic scaling characteristics of the Gaussian
processes. However, fast algorithms using approximate matrix-vector products
can be utilized for large data sets Raykar, Duraiswami, and Zhao (2010);
Chalupka, Williams, and Murray (2013). GPR has gained prominence providing
surrogate models for complex and multidimensional systems Forrester, Sóbester,
and Keane (2007); Perdikaris  _et al._ (2015); Feldstein  _et al._ (2020). In
ROM applications, these confidence measures can be particularly informative
when the ROM dimension is lower than the intrinsic dimension of the system.
Wan and Sapsis (2017) formulated a stochastic model based on GPR dynamics and
utilized a Monte-Carlo framework for the forecast of the system’s state and
corresponding uncertainty. Maulik  _et al._ (2020) utilized a Gaussian process
emulator for the dynamical evolution of the latent space state variables,
obtained from POD and autoencoder compression, for the shallow water
equations. We also note that Raissi and Karniadakis (2018) developed a
Gaussian process framework to learn PDEs from relatively small quantities of
data. Xiao  _et al._ (2015a) applied a second order Taylor series scheme and a
Smolyak sparse grid collocation method to calculate the POD modal coefficients
at each time step from their values at earlier time steps. A radial basis
function (RBF) multi‐dimensional interpolation was used in Ref. Xiao  _et al._
(2015b) for similar purposes. RBF interpolations have been also utilized for
parameterized problems as mappings from the parameter space to the ROM space
Rajaram  _et al._ (2020).

###  V.5 Data assimilation and error correction

An excellent discussion of the similarities between ML and data assimilation
(DA) tools has been recently provided by Alan Geer (Geer, 2021). The
synergistic integration of ML and DA is essential for developing improved
approaches Tang and Hsieh (2001); Liaqat, Fukuhara, and Takeda (2003); Brajard
_et al._ (2021, 2020); Bonavita and Laloyaux (2020); Farchi  _et al._ (2020);
Bocquet  _et al._ (2020); Farchi  _et al._ (2021); Pawar and San (2021); Pawar
_et al._ (2021a). Looking at the ROM framework from a data assimilation point
of view opens up innovative avenues to tackle the closure problem. Ideas from
optimal control theory, data assimilation, and parameter estimation were
proven to be valuable in this regard. Data assimilation is a generic framework
combining the available observations with the underlying dynamical principles
governing the system to estimate the physical quantities of interest. This is
usually accomplished by starting from a background solution and computing an
optimal estimation of the true state of the system that minimizes the
discrepancy between model predictions and collected observations. This
minimization problem is solved while taking into account the respective
statistical confidence of different observations, background solution, and
model uncertainty Ahmed, Pawar, and San (2020). In order to emphasize the
model’s error (due to truncation), D’Adamo  _et al._ (2007) added a Gaussian
variable to the ROM equation. This is similar to the weak variational data
assimilation framework implemented in Ref. Artana  _et al._ (2012) for ROMs
using real experimental conditions with noisy particle image velocimetry data.