# Multiscale model reduction for incompressible flows

Jared L. Callaham1∗, Jean-Christophe Loiseau2, and Steven L. Brunton1
1 Department of Mechanical Engineering, University of Washington, Seattle, WA
98195, USA
2 Arts et Métiers Institute of Technology, CNAM, DynFluid, HESAM Université,
F-75013 Paris, France

###### Abstract

Many unsteady flows exhibiting complex dynamics are nevertheless characterized
by emergent large-scale coherence in space and time. Reduced-order models
based on Galerkin projection of the governing equations onto an orthogonal
modal basis approximate the flow as a low-dimensional dynamical system with
linear and quadratic terms. However, these Galerkin models often fail to
reproduce the true dynamics, in part because they ignore important nonlinear
interactions with unresolved flow scales. Here, we use a separation of time
scales between the resolved and subscale variables to derive a reduced-order
model with cubic closure terms for the truncated modes, generalizing the
classic Stuart-Landau equation. The leading-order cubic terms are determined
by averaging out fast variables through a perturbation series approximation of
the action of a stochastic Koopman operator. We show analytically that this
multiscale closure model can capture both the effects of mean-flow deformation
and the energy cascade before demonstrating improved stability and accuracy in
models of chaotic lid-driven cavity flow and vortex pairing in a mixing layer.
This approach to closure modeling establishes a general theory for the origin
and role of cubic nonlinearities in low-dimensional models of incompressible
flows.

Keywords: Reduced-order modeling, incompressible fluid dynamics, stochastic
Koopman operator, closure modeling

##  1 Introduction

Despite the complex and often chaotic dynamics exhibited by many unsteady
fluid flows, they are often dominated by energetic coherent structures
evolving on relatively long length and time scales (Holmes et al., 1996). This
realization made it possible to study fluid flows as high-dimensional
dynamical systems evolving on a low-dimensional manifold, providing answers to
longstanding questions on topics such as the route to turbulence (Landau,
1944; Hopf, 1948; Ruelle & Takens, 1971; Swinney & Gollub, 1981) and the role
of nonlinear interactions (Stuart, 1958; Landau & Lifshitz, 1959). The
persistent nature of these coherent structures eventually also raised the
possibility of using low-dimensional surrogate models for optimization and
control objectives (Noack et al., 2011; Brunton & Noack, 2015; Rowley &
Dawson, 2017).

The fundamental challenges in constructing such reduced-order models may be
broken into two categories: approximating the structures themselves, and
approximating their evolution. We refer to these as the kinematic and dynamic
approximations, respectively. A modal separation of variables assumption is
often employed for the former task (Taira et al., 2017, 2020), but although
this may be suitable for many closed or diffusion-dominated flows, it is not a
natural representation of traveling waves or advection-dominated flows (Rowley
& Marsden, 2000; Reiss et al., 2018; Rim et al., 2018; Grimberg et al., 2020;
Mendible et al., 2020). This issue is inextricably linked to the problem of
modeling the coherent structure dynamics; as is well known in many domains, a
proper choice of coordinates can greatly simplify the modeling task (Champion
et al., 2019).

The multiscale nature of fluid flows further complicates both the kinematic
and dynamic aspects of low-dimensional modeling. For instance, the effective
dimensionality of a chaotic or turbulent flow may be orders of magnitude
greater than that of a laminar flow with one or two instability modes.
Meanwhile, the “triadic” structure of the nonlinear interactions in the
wavenumber or frequency domain ensures that the dynamics of all scales across
the flow are linked. Thus, even if the large, energetic coherent structures
can be approximated with a low-dimensional basis, models of their dynamics
that do not account for the role played by the unresolved degrees of freedom
are often unstable or physically inconsistent (Noack et al., 2011; Callaham et
al., 2021a). Similar considerations impact the development of numerical
methods (Bazilevs et al., 2007), self-consistent mean flow modeling (Meliga,
2017), and resolvent analysis (Padovan et al., 2020; Rigas et al., 2021;
Barthel et al., 2021, 2022).

The need for subscale modeling was apparent even in early work combining
empirical modal approximations, such as the proper orthogonal decomposition
(POD), with physics-based model reduction, such as Galerkin projection. For
example, low-dimensional models of vortex shedding in the globally unstable
cylinder wake could accurately predict the dynamics over short times, but were
subject to structural instability over a longer time horizon (Deane et al.,
1991; Ma & Karniadakis, 2002). Eventually Noack et al. (2003) showed that this
was a result of the failure of the standard post-transient POD basis to
resolve the Stuart-Landau mechanism of mean flow deformation associated with
nonlinear interactions between the fluctuations (Landau, 1944; Stuart, 1958),
an insight that enabled low-dimensional modeling of natural and actuated flows
with increasingly more complex dynamics (Luchtenburg et al., 2009; Deng et
al., 2020; Sieber et al., 2020; Deng et al., 2021; Callaham et al., 2021b).

Contemporary work on POD-Galerkin modeling of turbulent shear flows also
recognized the need for closure models that could approximate the effect of
unresolved scales (Aubry et al., 1988; Rempfer & Fasel, 1994; Ukeiley et al.,
2001). These efforts targeted two distinct physical mechanisms: mean flow
deformation and subscale dissipation. For flows that are either parallel
(Aubry et al., 1988) or weakly nonparallel under the assumption of Taylor’s
frozen turbulence hypothesis (Ukeiley et al., 2001), these authors developed a
Boussinesq eddy viscosity relationship between the resolved scales and a
slowly varying parallel mean flow, leading to stabilizing cubic terms
consistent with the Stuart-Landau description. To capture dissipation due to
the unresolved scales, these models also adopted a linear mixing length
approximation.

Although these studies remain landmark explorations of low-dimensional
coherent structure modeling, there are several opportunities to improve the
proposed closure strategies. First, it is difficult to generalize the Reynolds
stress models applied to parallel shear flows by Aubry et al. (1988) and
Ukeiley et al. (2001) to fully inhomogeneous flows. While the “shift mode”
approach to mean flow modeling via an augmented POD basis introduced by Noack
et al. (2003) is agnostic to the geometry of the flow, it requires computation
of the unstable steady state of the Navier-Stokes equations, which is not
generally experimentally accessible. Second, in the Richardson-Kolmogorov
energy cascade description, energy is transferred from large to small scales
through _nonlinear_ interactions, where it is finally dissipated. This is not
consistent with a linear mixing length model, a fact exploited by later work
investigating nonlinear models of subscale dissipation (Wang et al., 2012;
Cordier et al., 2013; Östh et al., 2014), particularly the finite time
thermodynamics approach, which is centered on modeling unresolved nonlinear
energy transfers (Noack et al., 2008).

Recent years have seen a surge in interest in data-driven and machine
learning-based methods (Brenner et al., 2019; Duraisamy et al., 2019; Brunton
et al., 2020), including a number of proposed closure and stabilization
schemes for reduced-order models, either through regression to additional
linear-quadratic terms (Xie et al., 2018; Mohebujjaman et al., 2017, 2018) or
by adding a deep learning model to approximate the residual (San & Maulik,
2018b, a; Menier et al., 2022). Alternative work has explored interpretable
system identification methods that forego the projection-based model
altogether (Brunton et al., 2016; Loiseau & Brunton, 2018; Peherstorfer &
Willcox, 2016; Qian et al., 2020; Callaham et al., 2021b), but may incorporate
physical constraints derived from the Galerkin system (Loiseau et al., 2018b;
Deng et al., 2020; Kaptanoglu et al., 2021a). If an accurate, nonintrusive
model is more important than interpretability or satisfying physical
constraints, then the traditional projection-based framework can be eliminated
altogether with black-box neural network forecasting methods (Hesthaven &
Ubbiali, 2018; Wan et al., 2018).

These empirical closure models are constructed on the assumption that the
influence of the unresolved variables can be approximated based on information
from the resolved variables alone. The success of the Reynolds stress-mean
flow models (Aubry et al., 1988; Ukeiley et al., 2001; Mantic̆-Lugo et al.,
2014), invariant or center manifold reductions (Guckenheimer & Holmes, 1983;
Coullet & Spiegel, 1983; Noack et al., 2003; Carini et al., 2015), and weakly
nonlinear analysis (Stuart, 1958; Sipp & Lebedev, 2007; Meliga et al., 2009;
Meliga & Chomaz, 2011) suggests that there may be circumstances where this
relationship may be derived analytically via traditional analysis. The
unifying thread between these methods is the assumption of a scale separation
between the resolved and unresolved variables that can be exploited to develop
an asymptotically correct closure model.

Figure 1:  Multiscale closure model applied to a mixing layer. The
visualization of the network of average quadratic energy transfer between the
leading harmonic modes shows the cascade of energy to higher-order modes. The
multiscale model reduction (MMR) method approximates the effects of unresolved
higher-order modes via stochastic averaging, which leads to a generalized
Stuart-Landau-type equation with cubic nonlinear interactions. The network
visualization of quadratic energy transfers (left) is computed from the modal
coefficients and Galerkin model for the leading harmonics, while that of of
the multiscale approximation (right) is a notional illustration of the origin
of the cubic terms in Eq. (49b). See Sec. 6.4 for details on the construction
of this figure and the low-dimensional model of the mixing layer.

Beyond the field of fluid dynamics, the method of adiabatic elimination
(Haken, 1983; Risken, 1996) has long been used to discard the fast variables
in systems with emergent large-scale coherence when there is a separation in
time scales, while heterogeneous multiscale methods have played an important
role in simulating physical systems with widely separated scales (E &
Engquist, 2003; E et al., 2007; E, 2011). A similar stochastic averaging
approach has been successful in climate modeling (Majda et al., 2001), where
the primitive equations have the same quadratic nonlinearity as the usual
Navier-Stokes equations without rotation, buoyancy, topography, etc. This
method, also called homogenization in the multiscale modeling literature, has
its roots in the theory of singular perturbations of Markov processes (Kurtz,
1973; Papanicolaou, 1976) and is rigorously supported for stochastic systems
with asymptotic scale separation; see Majda et al. (2001); Givon et al.
(2004); E (2011), or Pavliotis & Stuart (2012) for in-depth presentations.
With some assumptions on ergodicity, a similar approach can also be taken with
deterministic systems, even in a regime where the scale separation is not in
the asymptotic limit (Majda et al., 2006).

This work explores the application of multiscale stochastic averaging methods
developed by Majda et al. (2001), Givon et al. (2004); E (2011); Pradas et al.
(2012), Pavliotis & Stuart (2012), and others to the closure problem in
reduced-order models of incompressible flow or other systems reducible to
linear-quadratic dynamics (Rowley et al., 2004; Qian et al., 2020; Kaptanoglu
et al., 2021b), including introducing an approximation to the form of the fast
dynamics that allows for computation of the averaged dynamics in closed form.
The procedure is summarized in Fig. 1 and shown schematically in Fig. 2. We
refer to this secondary dimensionality reduction of the Galerkin system as
multiscale model reduction (MMR).

Since fluid flows generally do not have a true scale separation away from the
threshold of instability, we demonstrate via numerical simulations that this
method is a robust and systematic approach to stabilizing low-dimensional
models. This extends the work of Majda et al. (2006) exploring the application
of this class of methods to systems beyond the parameter regimes where their
validity can be rigorously proven. Mltiscale model reduction is a unified
framework for understanding the origin and importance of cubic terms in
reduced-order models of the linear-quadratic Navier-Stokes equations, capable
of capturing both mean flow deformation and subscale dissipation. Throughout
this work we also highlight connections to other modeling methods, including
Koopman theory and weakly nonlinear analysis.

Figure 2:  Schematic of nonlinear interactions in the multiscale closure
scheme. The dynamics of one variable ($a_{1}$) in a system with two slow
variables involve quadratic interactions between fast and slow variables that
would be neglected in a standard truncation. Instead, the proposed method
averages over the fast scales, ultimately generating effective cubic
nonlinearities in the closed equations.

##  2 Reduced-order modeling for incompressible flows

For more than three decades, proper orthogonal decomposition (POD) and
Galerkin projection have been the foundational tools in low-order modeling for
nonlinear incompressible fluid flow. POD analysis identifies dominant coherent
structures in the flow and provides an energy-optimal linear modal basis.
Galerkin projection approximates both the state of the flow and its time
derivative in this finite-dimensional subspace, resulting in a reduced system
of ordinary differential equations in terms of the POD mode amplitudes. These
topics have been covered extensively elsewhere (see e.g. Holmes et al. (1996);
Noack et al. (2011); Benner et al. (2015); Rowley & Dawson (2017); Taira et
al. (2017)); here we provide a brief overview.

We assume the flow is governed by the unsteady incompressible Navier-Stokes
equations

$\displaystyle\frac{\partial\bm{u}}{\partial
t}+\nabla\cdot\left(\bm{u}\otimes\bm{u}\right)$ $\displaystyle=-\nabla
p+\frac{1}{\real}\nabla^{2}\bm{u}$ (1) $\displaystyle\nabla\cdot\bm{u}$
$\displaystyle=0,$

where $\bm{u}(\bm{x},t)$ is the velocity field, $p$ is the pressure field, and
$\real$ is the Reynolds number, $\real=UL\nu$, based on the kinematic
viscosity $\nu$ and suitable length and velocity scales $L$ and $U$. The
nonlinear term is expressed in divergence form with the tensor product
$\otimes$. In this work, we only consider two dimensional examples for which
$\bm{u}=\begin{bmatrix}u&v\end{bmatrix}^{T}$, although all of the following
applies equally well to three-dimensional flows.

More broadly, POD can be used for any collection of data for which a suitable
norm can be defined. The specific method described here for reducing the
Navier-Stokes equations to a linear-quadratic system of ODEs via Galerkin
projection relies on the quadratic nonlinearity in Eq. (1), although with
certain restrictions this procedure can also be applied to systems such as
compressible flow (Rowley et al., 2004), magnetohydrodynamics (Kaptanoglu et
al., 2021b), or systems with non-polynomial nonlinearities that have been
“lifted” to a coordinate system in which the dynamics are linear-quadratic
(Qian et al., 2020).

###  2.1 Proper orthogonal decomposition

In order to reduce the system of partial differential equations (1) to a
finite-dimensional system of ordinary differential equations, we assume that
the space and time dependence can be separated, so that an arbitrary velocity
field $\bm{u}(\bm{x},t)$ can be approximated with the $r$-dimensional modal
representation

$\bm{u}(\bm{x},t)\approx\bm{\psi}_{0}(\bm{x})+\sum_{i=1}^{r}\bm{\psi}_{i}(\bm{x})a_{i}(t),$
(2)

where $\bm{\psi}_{0}(\bm{x})$ is a fixed base flow around which the expansion
is performed and each mode individually satisfies the divergence-free
constraint $\nabla\cdot\bm{\psi}_{i}=0$ for $i=0,2,\dots,r$. Typically in this
work we will take $\bm{\psi}_{0}$ to be the mean flow $\bar{\bm{u}}$ as
estimated with a time average, although another common choice is to let
$\bm{\psi}_{0}$ be the steady-state solution to the Navier-Stokes equations.

The proper orthogonal decomposition (POD) identifies the set of modes
$\\{\bm{\psi}_{i}\\}_{i=1}^{r}$ that provide the optimal rank-$r$
approximation for an average velocity field in the energy norm induced by the
velocity inner product

$\langle\bm{u},\bm{u}^{\prime}\rangle=\int_{\Omega}\bm{u}^{\prime}(\bm{x})\cdot\bm{u}^{*}(\bm{x})\differential\,\Omega$
(3)

defined on the spatial domain $\Omega$, where the star indicates complex
conjugation. Since our data may be represented on a nonuniform mesh, we
approximate this integral with a weighted sum

$\langle\bm{u},\bm{v}\rangle\approx\mathsfbi{u}^{T}\mathsfbi{W}\mathsfbi{v},$
(4)

where $\mathsfbi{u}$ and $\mathsfbi{v}$ are the discrete approximations to
$\bm{u}$ and $\bm{v}$, and $\mathsfbi{W}$ is a diagonal weight matrix
containing the cell volumes. The transpose in Eq. (4) is understood to include
complex conjugation in the case that the fields or modes are complex valued.

The POD modes themselves are eigenfunctions of the two-point spatial
correlation tensor, although in practice they are typically approximated
numerically with the method of snapshots (Sirovich, 1987) or the singular
value decomposition (Brunton & Kutz, 2019). In the method of snapshots, the
POD modes are estimated from a weighted sum of velocity field snapshots, where
the weights are determined from eigenvectors of the two-point temporal
correlation matrix. The construction from a linear combination of divergence-
free snapshots also ensures that the velocity modes are divergence-free.

The POD expansion (2) is often only applied to the velocity fields. This can
pose a problem for model reduction methods approximating the pressure gradient
term in the Navier-Stokes equations. While this term vanishes for closed flows
and is often negligible for flows with a localized global instability, it is
important in open shear flows, such as the mixing layer examined in Sec. 6
(Noack et al., 2005). There are various methods for approximating the pressure
term (see e.g. Caiazzo et al. (2014)), but in this work we use a velocity-
pressure expansion with the same inner product (3). Defining
$\bm{q}(\bm{x},t)=\begin{bmatrix}\bm{u}&p\end{bmatrix}^{T}$, we replace (2)
and (3) with

$\displaystyle\bm{q}(\bm{x},t)$
$\displaystyle\approx\begin{bmatrix}\bm{\psi}^{\bm{u}}_{0}(\bm{x})\\\
\psi^{p}_{0}(\bm{x})\end{bmatrix}+\sum_{i=1}^{r}\begin{bmatrix}\bm{\psi}^{\bm{u}}_{i}(\bm{x})\\\
\psi^{p}_{i}(\bm{x})\end{bmatrix}a_{i}(t)$ (5a)
$\displaystyle=\bm{\psi}_{0}(\bm{x})+\sum_{i=1}^{r}\bm{\psi}_{i}(\bm{x})a_{i}(t)$
(5b)

and energy inner product
$\langle\bm{q},\bm{q}^{\prime}\rangle_{\bm{u}}\equiv\langle\bm{u},\bm{u}^{\prime}\rangle$.
Numerically, we simply carry the pressure fields through the method of
snapshots by setting the diagonal entries of $\mathsfbi{W}$ corresponding to
pressure to zero. The pressure components of the POD modes are therefore
estimated from the same linear combination of snapshots as the velocity
components.

Besides providing an optimal linear modal approximation, the POD
representation has the advantage that the velocity modes are orthonormal and
the temporal coefficients are uncorrelated in time:

$\displaystyle\langle\bm{\psi}_{i},\bm{\psi}_{j}\rangle_{\bm{u}}=\delta_{ij}$
(6a) $\displaystyle\overline{a_{i}a_{j}}=\lambda_{i}\delta_{ij},$ (6b)

where the constant of proportionality $\lambda_{i}$ for $\overline{a_{i}^{2}}$
is the POD eigenvalue, with the physical meaning of average fluctuation
kinetic energy in the $i^{\mathrm{th}}$ mode. The subspace dimension $r$ is
typically selected based on residual energy content of the truncated modes.
Once the modes $\\{\bm{\psi}_{i}\\}_{i=1}^{r}$ are calculated, the temporal
coefficients $\bm{a}(t)$ can be computed by projecting the time series onto
the span of the modes, i.e.
$a_{i}(t)=\langle\bm{\psi}_{i}(\bm{x}),\bm{u}(\bm{x},t)\rangle$.

Although POD analysis is a powerful tool, it suffers from several limitations.
First, it assumes a space-time separation of variables, which is typically
appropriate for closed flows, but is often a poor approximation of
inhomogeneous open flows dominated by traveling waves. Second, the condition
that the temporal coefficients be linearly uncorrelated does not rule out
nonlinear correlations. Finally, the optimality guarantee of the POD
approximation is purely based on the instantaneous flow field and does not
account for any time dynamics. The relationship between the first two issues
was examined at length by Callaham et al. (2022), where it was shown that
leveraging nonlinear correlations can stabilize reduced-order models of
advection-dominated flows while reducing the dimensionality of the model. The
third issue motivates alternative modal representations, such as dynamic mode
decomposition (Rowley et al., 2009; Schmid, 2010; Kutz et al., 2016; Schmid,
2022), spectral POD (Towne et al., 2018; Schmidt & Colonius, 2020), balanced
POD Willcox & Peraire (2002); Rowley (2005), resolvent modes (McKeon & Sharma,
2010; Sharma & McKeon, 2013; Zare et al., 2017; Padovan et al., 2020), exact
coherent structures (Graham & Floryan, 2020; Sharma et al., 2016; Rosenberg &
McKeon, 2019), and trajectory-based optimization (Otto et al., 2022). These
drawbacks notwithstanding, we use POD in this work due to its numerical
robustness and orthonormality properties.

###  2.2 Galerkin projection

Galerkin projection is a widely used method in computational and applied
mathematics for discretizing continuous problems, in particular partial
differential equations (PDEs). In the context of reduced-order modeling, this
approach typically refers to projecting the infinite-dimensional PDE onto a
finite-dimensional subspace spanned by global, orthogonal basis functions.
While such a basis can by constructed analytically for simple geometries, the
energy optimality and orthonormality of POD makes it a natural candidate for
statistically stationary (i.e. post-transient) flows with more complex
geometries. Moreover, it is a data-driven approach that tailors the basis
based on observed flow kinematics.

The derivation of the Galerkin system begins by substituting the velocity-
pressure POD approximation (5) into the Navier-Stokes equations (1). In order
for the dynamics to be an optimal continuous-time approximation, the residual
error should be orthogonal to the POD subspace. Using the inner product (3)
and orthonormality properties (6), this optimality condition leads to the
linear-quadratic system of ODEs (Holmes et al., 1996; Noack et al., 2011)

$\dot{a}_{i}=F_{i}+L_{ij}a_{j}+Q_{ijk}a_{j}a_{k},\qquad i,j,k=1,2,\dots,r$ (7)

with constant, linear, and quadratic terms given by

$\displaystyle F_{i}$
$\displaystyle=\langle\bm{\psi}_{i},-\nabla\cdot(\bm{\psi}_{0}^{\bm{u}}\otimes\bm{\psi}_{0}^{\bm{u}})-\nabla\psi_{0}^{p}+\mbox{{Re}}^{-1}\nabla^{2}\bm{\psi}_{0}^{\bm{u}}\rangle_{\bm{u}}$
(8a) $\displaystyle L_{ij}$
$\displaystyle=\langle\bm{\psi}_{i},-\nabla\cdot(\bm{\psi}_{0}^{\bm{u}}\otimes\bm{\psi}^{\bm{u}}_{j}+\bm{\psi}^{\bm{u}}_{j}\otimes\bm{\psi}_{0}^{\bm{u}})-\nabla\psi^{p}_{j}+\mbox{{Re}}^{-1}\nabla^{2}\bm{\psi}^{\bm{u}}_{j}\rangle_{\bm{u}}$
(8b) $\displaystyle Q_{ijk}$
$\displaystyle=\langle\bm{\psi}_{i},-\nabla\cdot(\bm{\psi}_{j}^{\bm{u}}\otimes\bm{\psi}_{k}^{\bm{u}})\rangle_{\bm{u}}.$
(8c)

Throughout most of this work, we take the base mode to be the time-averaged
mean flow. However, if the flow is instead expanded about the steady-state
solution, the constant forcing term $\mathsfbi{F}$ vanishes, ensuring that the
origin is a fixed point.

###  2.3 Symmetries and energy conservation in the Galerkin system

The combination of proper orthogonal decomposition and Galerkin projection is
a popular choice for model reduction of incompressible flows, in part because
it provides a semi-empirical “gray-box” framework that combines some of the
advantages of data-driven analysis with physics-based modeling. In particular,
the modal basis is tailored from data for a particular flow configuration,
while the dynamics are obtained by directly projecting the governing equations
onto this basis. One consequence is that the dynamical system Eq. (7) inherits
certain properties from the Navier-Stokes equations. For example, although we
will mostly take the base flow $\bm{\psi}_{0}(\bm{x})$ to be the time-averaged
mean flow, if the the flow is instead expanded about the steady-state solution
the constant forcing term $\mathsfbi{F}$ vanishes, ensuring that the origin is
a fixed point. If, in addition, an energy inner product such as Eq. (3) is
used for projection, the reduced-order model also preserves Lyapunov stability
of the origin, although not necessarily the stability of limit cycles or other
fixed points (Rowley et al., 2004).

A key feature of the quadratic nonlinearity of the Navier-Stokes equations is
that it is energy-preserving, in the sense that it has no net contribution to
the evolution equation for kinetic energy in the spectral domain (Kraichnan &
Chen, 1989), with some restrictions on the boundary conditions. This is the
foundation of the energy cascade picture of turbulence, in which the role of
the nonlinearity is to transfer energy from the large scales to the small,
dissipative scales. For the Galerkin system, the energy preservation condition
implies that (Schlegel & Noack, 2015)

$Q_{ijk}+Q_{ikj}+Q_{jik}+Q_{jki}+Q_{kij}+Q_{kji}=0.$ (9)

This is often approximately true numerically for POD-Galerkin systems, but
enforcing it explicitly can improve the stability of the model (Cordier et
al., 2013). Based on the structure of the linear-quadratic system (7), the
quadratic tensor can also be symmetrized in the last two indices:

$Q_{ijk}=Q_{ikj}.$ (10)

Finally, Schlegel & Noack (2015) showed that the long-term boundedness of the
Galerkin system can be guaranteed if a criterion based on generalized Lyapunov
functions is satisfied. Although this is a dynamical systems analysis and not
an _a priori_ theoretical guarantee, it can also be used as a constraint in
sparse model identification to construct empirical Galerkin-type models with
guaranteed stability (Kaptanoglu et al., 2021a).

##  3 Generators of deterministic and stochastic processes

In the study of mechanics there are often several equivalent representations
of the same system (e.g. Lagrangian, Hamiltonian, etc.), each of which may be
useful depending on the application. The operator-theoretic perspective will
be especially useful in the development of the multiscale closure model as a
framework for abstract formal manipulations of the dynamical systems model. In
particular, the closure model will be derived via a perturbation series
approximation to the action of a stochastic Koopman operator. Critically, the
infinite-dimensional operator itself need not be computed or directly
approximated; instead, the closed model appears as a solvability condition
analogous to the derivation of the amplitude equation in weakly nonlinear
analysis (Sipp & Lebedev, 2007).

This section gives a brief overview of topics necessary for the development of
the closure models in Section 4; for a more in-depth presentation of these
topics see, e.g. Risken (1996); E (2011); Pavliotis & Stuart (2012); Klus et
al. (2018, 2020); Brunton et al. (2022). Koopman theory is perhaps the most
widely-discussed operator-theoretic method in recent work on the analysis of
fluid flows and other large-scale nonlinear dynamics (Mezić, 2005; Rowley et
al., 2009; Mezić, 2013; Klus et al., 2018; Brunton et al., 2022), making it a
convenient place to begin the discussion.

Suppose the state $\bm{x}(t)$ is governed by an autonomous ordinary
differential equation (ODE)

$\dot{\bm{x}}=\bm{f}(\bm{x}),\qquad\bm{x}\in\mathcal{X}.$ (11)

We will assume, unless otherwise specified, that states are real-valued, e.g.
$\mathcal{X}\subseteq\mathbb{R}^{N}$, and that all other functions are
$L^{2}$-integrable with inner product

$\langle
f(\bm{x}),g(\bm{x})\rangle=\int_{\mathcal{X}}f(\bm{x})g(\bm{x})\differential
x,$ (12)

where the usual complex-conjugation is omitted since $g$ is real-valued.

The Koopman operator $\mathcal{K}^{t}$ acts on the $L^{2}$-integrable space of
scalar observables $g(\bm{x}):\mathcal{X}\rightarrow\mathbb{R}$, advancing
them forward time $t$. More precisely, let $\bm{x}$ be the solution to the
initial value problem of (11) with $\bm{x}(0)=\bm{x}_{0}$. Then the Koopman
operator is defined as

$\mathcal{K}^{t}g(\bm{x}_{0})=g(\bm{x}(t)).$ (13)

Although the Koopman operator is generally difficult to either represent
explicitly or approximate in a useful finite-dimensional subspace, analysis of
its spectral properties has drawn great interest in recent work. For our
purposes, the infinitesimal generator $\mathcal{L}$, defined by

$\mathcal{L}g=\lim_{t\rightarrow 0}\frac{\mathcal{K}^{t}g-g}{t}$ (14)

and sometimes called the Lie operator, is more theoretically useful. If we
consider $g(\bm{x}(t))$ to be an explicit function of time, then $\mathcal{L}$
can be derived by applying the chain rule to $g$. With a slight abuse of
notation, if $g(\bm{x},t)$ is taken to be a function of both time and initial
state $\bm{x}$ (i.e. in the Lagrangian frame of reference), then this can be
extended to a PDE over all of state space:

$\partialderivative{g}{t}=\mathcal{L}g=\bm{f}(\bm{x})\cdot\nabla g.$ (15)

Thus the generator $\mathcal{L}$ is a linear advection operator governing the
evolution of scalar observables of the system described by (11). Similarly,
the adjoint of (15),

$\partialderivative{\rho}{t}=\mathcal{L}^{\dagger}\rho=-\nabla\cdot(\rho\bm{f}(\bm{x})),$
(16)

is a continuity equation governing the evolution of densities $\rho(\bm{x},t)$
in phase space. As a point of reference, equation (16) reduces to the
Liouville equation from classical mechanics under the incompressibility
condition $\nabla\cdot\bm{f}=0$, and $\mathcal{L}^{\dagger}$ is also known as
the generator of the Perron-Frobenius operator (Froyland, 2005; Froyland &
Padberg, 2009; Froyland et al., 2010; Klus et al., 2016).

This description of dynamics can readily be extended to systems governed by
stochastic differential equations (SDEs) of the form

$\dot{\bm{x}}=\bm{f}(\bm{x})+\bm{\Sigma}\bm{w}(t),$ (17)

where the deterministic component $\bm{f}(\bm{x})$ is known as the drift
function and the diffusion matrix $\bm{\Sigma}$ modifies a vector-valued
Gaussian white noise process $\bm{w}(t)$. In this case the evolution of the
probabilitiy distribution $\rho(\bm{x},t)$ is governed by the stochastic
analog of (16), known as the forward Kolmogorov, or Fokker-Planck, equation:

$\partialderivative{\rho}{t}=\mathcal{L}^{\dagger}\rho=-\nabla\cdot(\rho\bm{f})+\nabla\nabla^{T}:(\rho\bm{D}),$
(18)

where $\bm{D}=\bm{\Sigma}\bm{\Sigma}^{T}/2$ is the diffusion tensor and the
colon denotes tensor contraction. The associated backwards Kolmogorov equation
is the adjoint of (18):

$\partialderivative{g}{t}=\mathcal{L}g=\bm{f}\cdot\nabla
g+\bm{D}:\nabla\nabla^{T}g.$ (19)

To interpret (19), consider the expectation of a scalar observable
$g(\bm{x}):\mathcal{X}\rightarrow\mathbb{R}$ defined by the inner product of
$g$ with the probability distribution $\rho_{0}(\bm{x})$:

$\mathop{\mathbb{E}}[g]=\langle\rho_{0},g\rangle=\int_{\mathcal{X}}\rho_{0}(\bm{x})g(\bm{x})\differential\bm{x}.$
(20)

The probability distribution is advanced in time by
$\rho(\bm{x},t)=e^{\mathcal{L}^{\dagger}t}\rho_{0}(\bm{x})$. Using the
definition of the adjoint,
$\langle\rho,\mathcal{L}g\rangle=\langle\mathcal{L}^{\dagger}\rho,g\rangle$,
the expectation of $g$ at time $t$ can be written equivalently as

$\mathop{\mathbb{E}}[g](t)=\langle\rho,g\rangle=\int_{\mathcal{X}}\rho_{0}(\bm{x})e^{\mathcal{L}t}g(\bm{x})\differential\bm{x}.$
(21)

Therefore, just as the Koopman operator advances an observable in time, the
backwards Kolmogorov equation describes the evolution of the expectation of an
observable.

As a simple example relevant for the closure modeling presented below,
consider the Ornstein-Uhlenbeck process defined by the SDE

$\dot{x}=\nu(\mu-x)+\sigma w(t)$ (22)

for positive constants $\nu$, $\mu$, and $\sigma$. The stationary distribution
$\rho^{\infty}(x)$ is given by the steady-state Fokker-Planck equation

$0=\mathcal{L}^{\dagger}\rho^{\infty}=-\derivative{x}\rho^{\infty}(x)\nu(\mu-x)+\frac{\sigma^{2}}{2}\derivative[2]{x}\rho^{\infty}(x),$
(23)

along with the usual normalization condition on $\rho^{\infty}$. This can be
solved analytically by a Gaussian distribution with mean $\mu$ and variance
$\sigma^{2}/2\nu$. Defining an observable that is the state itself $g(x)=x$,
the evolution is given by the backwards Kolmogorov equation (19), which in
this case reduces to the initial value problem

$\derivative{g}{t}=\nu(\mu-g),\qquad g(0)=x.$ (24)

The expectation at time $t$ is the solution $g(t)$, which in this case is
exponential decay towards the mean $\mu$:

$\mathop{\mathbb{E}}[x](t)=g(t)=\mu-(x-\mu)e^{-\nu t}.$ (25)

Given that the model reduction procedure described in Section 2 aims to reduce
the physics model from an infinite-dimensional partial differential equation
to a finite-dimensional system of ordinary differential equations, it may be
counterintuitive that it would be helpful to return to an infinite-dimensional
function space and represent the dynamics with a partial differential
equation. The primary advantage in doing so is that these generators are
linear operators, which in some cases can be amenable to approaches that are
unavailable for the nonlinear dynamics of the ODE (11) or SDE (17).

Aside from simple one-dimensional or linear examples like the Ornstein-
Uhlenbeck process, it is typically very difficult to solve these PDEs
analytically, motivating the approximation in terms of an asymptotic expansion
in Sec. 4. As will be seen, this expansion can be used to derive a solvability
condition based on the Fredholm alternative. Using the manipulations reviewed
in this chapter, particularly in terms of the generators and their adjoints,
this condition can be interpreted as an average over subscale variables, which
can be computed explicitly under suitable assumptions. The averaged slow
dynamics then take the form of a cubic generalized Stuart-Landau equation that
accounts for the leading-order effect of the unresolved fast variables.

##  4 Multiscale closure modeling

One of the primary difficulties of reduced-order modeling for fluid flows is
that the nonlinear interactions in the Navier-Stokes equations transfer energy
between all scales of the flow. Thus, restricting the dynamics to a low-
dimensional subspace can lead to significant approximation errors. Generally
speaking, the goal of closure modeling is to augment this truncated model with
additional terms that account for the effect of the unresolved scales.

The multiscale approach to closure modeling accomplishes this by first
partitioning the dynamics into fast and slow variables, and then approximating
the solution to the associated backwards Kolmogorov equation with a
perturbation series expansion. The result can be interpreted as a Koopman
generator of the form of Eq. (15) corresponding to the coarse-grained dynamics
for the slow variables. As we will show, when applied to the linear-quadratic
Galerkin system (7), this leads to a generalized Stuart-Landau equation
including cubic terms.

An interesting feature of cubic Stuart-Landau-type models of fluid flow is
that they are often more accurate than linear-quadratic models, even though
the nonlinearity in the underlying governing equations is quadratic (Noack et
al., 2003; Loiseau & Brunton, 2018). One reason for this is the artificial
space-time separation of variables introduced by the modal representation (2),
which can introduce spurious degrees of freedom that represent phase-locked
harmonics of traveling waves, for instance (Callaham et al., 2022). As a
result, POD-Galerkin models are vulnerable to decoherence and can fail to
improve with increasing rank, even when the kinematic approximation of the
flow field becomes nearly perfect. As illustrated by the examples in Sec. 6,
the introduction of cubic terms can mitigate this, either by modifying the
dynamics to resemble phase-locked nonlinear oscillators or by eliminating
spurious degrees of freedom altogether.

Although multiscale methods can be rigorously justified when there is a strict
separation of timescales (Majda et al., 2001; E & Engquist, 2003; E, 2011;
Pavliotis & Stuart, 2012), there is no such spectral gap in most fluid flows
of practical interest. The proposed method for model reduction should
therefore properly be viewed as an approximate closure model motivated by the
asymptotic limit in the spirit of the variational multiscale approach to
turbulence modeling (Bazilevs et al., 2007). Nevertheless, it provides a
systematic method for stabilizing low-dimensional models and also highlights a
generic mechanism by which higher-order nonlinearities can arise from the
quadratic term in the governing equations.

###  4.1 Averaging over unresolved variables

Beginning with the Galerkin dynamics (7) for the state consisting of POD
coefficients $\bm{a}(t)\in\mathbb{R}^{r}$, where we assume that $r$ is large
enough for an accurate kinematic reconstruction of a typical flow field, we
partition the system into slow variables
$\bm{x}(t)\in\mathcal{X}=\mathbb{R}^{r_{0}}$ and fast variables
$\bm{y}(t)\in\mathcal{Y}=\mathbb{R}^{r-r_{0}}$, so that
$\bm{a}=\begin{bmatrix}\bm{x}^{T}&\bm{y}^{T}\end{bmatrix}^{T}$.

For concise notation we will use the Einstein convention that repeated indices
imply summation and we will omit explicit summation unless not doing so would
lead to ambiguity. We will index the slow variables with Roman subscripts $i$,
$j$, $k$, $\ell$ ranging from $1$ to $r_{0}$ and the fast variables with Greek
subscripts $\alpha$, $\beta$, $\gamma$ ranging from $1$ to $r-r_{0}$. Without
loss of generality we will also assume that the quadratic term has been
symmetrized in the last two indices, so that $Q_{ijk}=Q_{ikj}$. Then the
partitioned Galerkin system is

$\displaystyle\dot{x}_{i}$
$\displaystyle=f^{x}_{i}(\bm{x},\bm{y})=F^{1}_{i}+L^{11}_{ij}x_{j}+L^{12}_{i\beta}y_{\beta}+Q^{111}_{ijk}x_{j}x_{k}+2Q^{112}_{ij\beta}x_{j}y_{\beta}+Q^{122}_{i\beta\gamma}y_{\beta}y_{\gamma}$
(26a) $\displaystyle\dot{y}_{\alpha}$
$\displaystyle=f^{y}_{\alpha}(\bm{x},\bm{y})=F^{2}_{\alpha}+L^{21}_{\alpha
j}x_{j}+L^{22}_{\alpha\beta}y_{\beta}+Q^{211}_{\alpha
jk}x_{j}x_{k}+2Q^{212}_{\alpha
j\beta}x_{j}y_{\beta}+Q^{222}_{\alpha\beta\gamma}y_{\beta}y_{\gamma}.$ (26b)

Standard truncation of this system is equivalent to retaining only the terms
$\mathsfbi{F}^{1}$, $\mathsfbi{L}^{11}$, and $\mathsfbi{Q}^{111}$. Here we
will attempt to approximate the terms involving the fast variables in
$\bm{f}^{x}(\bm{x},\bm{y})$ in an average sense in order to derive a closed
system $\dot{\bm{x}}=\hat{\bm{f}}(\bm{x})$, largely following the approach of
Pavliotis & Stuart (2012).

We assume that the time scales of the fast and slow dynamics are separated by
a parameter $\epsilon\ll 1$, so that we may define
$\bm{f}^{x}(\bm{x},\bm{y})\equiv\epsilon\tilde{\bm{f}}^{x}(\bm{x},\bm{y})$.
Furthermore, we note that the role of the fast self-interaction term
$\mathsfbi{Q}^{222}(\bm{y},\bm{y})$ is to transfer energy between the
unresolved scales. Since this mechanism is of secondary importance to the
transfers between slow and fast scales, we apply a version of the “working
assumption of stochastic modeling” (Majda et al., 2001):

$Q^{222}_{\alpha\beta\gamma}y_{\beta}y_{\gamma}\approx\sigma_{\alpha}w_{\alpha}(t),$
(27)

where $\bm{w}(t)$ is a Wiener process and $\bm{\sigma}$ is an as-yet-undefined
constant forcing amplitude. Finally, we coarse-grain the dynamics on the time
scale $\tau=\epsilon t$, leading to the slow/fast system

$\displaystyle\partialderivative{x_{i}}{\tau}$
$\displaystyle=\tilde{f}^{x}_{i}(\bm{x},\bm{y})$ (28a)
$\displaystyle\partialderivative{y_{\alpha}}{\tau}$
$\displaystyle=\frac{1}{\epsilon}f^{y}_{\alpha}(\bm{x},\bm{y})+\frac{1}{\sqrt{\epsilon}}\sigma_{\alpha}w_{\alpha}(\tau),$
(28b)

where we have also used the scaling property of Wiener processes that
$\bm{w}(\epsilon t)=\sqrt{\epsilon}\bm{w}(t)$.

Defining an arbitrary scalar-valued observable $g^{\epsilon}(\bm{x},\bm{y})$,
the backwards Kolmogorov equation (19) associated with Eq. (28) is

$\partialderivative{g^{\epsilon}}{\tau}=\frac{1}{\epsilon}\underbrace{\left[f^{y}_{\alpha}(\bm{x},\bm{y})\partialderivative{y_{\alpha}}+\frac{\sigma_{\alpha}^{2}}{2}\partialderivative[2]{y_{\alpha}}\right]}_{\mathcal{L}_{0}}g^{\epsilon}+\underbrace{\tilde{f}^{x}_{i}(\bm{x},\bm{y})\partialderivative{x_{i}}}_{\mathcal{L}_{1}}g^{\epsilon}.$
(29)

Note that $\mathcal{L}_{0}$ can be viewed as the generator of a stochastic
process in $\bm{y}$ with $\bm{x}$ as a fixed parameter. Consequently, we make
the following assumptions related to the ergodicity of $\bm{y}$:

  1. 1.

The operator $\mathcal{L}_{0}$ has a one-dimensional nullspace spanned by
constants in $\bm{y}$:

$\mathcal{L}_{0}g(\bm{x})=0.$ (30)

  2. 2.

The Fokker-Planck operator $\mathcal{L}_{0}^{\dagger}$ has a one-dimensional
nullspace corresponding to the stationary distribution
$\rho^{\infty}_{\bm{x}}$, where again $\bm{x}$ is treated as a fixed
parameter:

$\mathcal{L}_{0}^{\dagger}\rho^{\infty}_{\bm{x}}(\bm{y})=0,$ (31)

along with the usual normalization condition
$\int_{\mathcal{Y}}\rho^{\infty}_{\bm{x}}\differential\bm{y}=1$.

These assumptions obviate the need to express the fast scale $\bm{y}$ as an
instantaneous function of $\bm{x}$, as in an invariant manifold model
(Guckenheimer & Holmes, 1983; Pavliotis & Stuart, 2012). Instead, the fast
variable is modeled with a simplified distribution that can be averaged over
as follows. We assume that the backwards Kolmogorov equation (29) can be
approximated by means of an asymptotic expansion

$g^{\epsilon}(\bm{x},\bm{y},\tau)=g_{0}+\epsilon
g_{1}+\mathcal{O}(\epsilon^{2}).$ (32)

Substituting into Eq. (29) and equating powers of $\epsilon$ gives the
consistency conditions

$\displaystyle 0$ $\displaystyle=\mathcal{L}_{0}g_{0},$ (33a)
$\displaystyle\partialderivative{g_{0}}{\tau}$
$\displaystyle=\mathcal{L}_{0}g_{1}+\mathcal{L}_{1}g_{0}.$ (33b)

By virtue of the assumption of a one-dimensional nullspace, Eq. (33a) is
satisfied if $g_{0}$ is not a function of $\bm{y}$, or
$g_{0}=g_{0}(\bm{x},\tau)$. As expected, the leading-order solution does not
depend on the fast variable.

An effective evolution equation for $g_{0}$ can be derived by considering Eq.
(33b) as a linear equation for $g_{1}(\bm{x},\bm{y},\tau)$. The Fredholm
alternative specifies that for any equation of the form
$\mathcal{L}_{0}g_{1}=b$ to have a unique solution, all functions $\rho$ in
the nullspace of the adjoint operator $\mathcal{L}_{0}^{\dagger}$ must be
orthogonal to $b$. Since we have assumed that the nullspace of
$\mathcal{L}_{0}^{\dagger}$ is one-dimensional and spanned by the stationary
distribution $\rho_{\bm{x}}^{\infty}(\bm{y})$, this implies that
$\langle\rho_{\bm{x}}^{\infty},b\rangle=0$. In other words, the solvability
condition for

$\mathcal{L}_{0}g_{1}=-\partialderivative{g_{0}}{\tau}+\mathcal{L}_{1}g_{0}$
(34)

is that the right-hand side of Eq. (34) has zero mean with respect to the
stochastic process generated by $\mathcal{L}_{0}$:

$\int_{\mathcal{Y}}\rho_{\bm{x}}^{\infty}(\bm{y})\left[-\partialderivative{\tau}+\mathcal{L}_{1}\right]g_{0}(\bm{x},\tau)\differential\bm{y}=0.$
(35)

Using the normalization condition for $\rho_{x}^{\infty}$, this simplifies to
a closed evolution equation for $g_{0}(\bm{x},\tau)$:

$\partialderivative{g_{0}}{\tau}=\left[\int_{\mathcal{Y}}\rho_{\bm{x}}^{\infty}(\bm{y})\tilde{f}_{i}^{x}(\bm{x},\bm{y})\differential\bm{y}\right]\partialderivative{g_{0}}{x_{i}}.$
(36)

Again, the observable $g_{0}$ itself is arbitrary, but Eq. (36) has the form
of the generator of a Koopman operator corresponding to the coarse-grained
dynamics in $\bm{x}$ alone. Undoing the $\epsilon$ scaling in $\tau$ and
$\tilde{\bm{f}}^{x}$, Eq. (36) can be written as

$\partialderivative{g_{0}}{t}=\hat{\bm{f}}(\bm{x})\cdot\nabla g_{0},$ (37)

corresponding to the averaged dynamics

$\displaystyle\dot{\bm{x}}$ $\displaystyle=\hat{\bm{f}}(\bm{x}),$ (38a)
$\displaystyle\hat{\bm{f}}(\bm{x})$
$\displaystyle=\int_{\mathcal{Y}}\rho_{\bm{x}}^{\infty}(\bm{y})f_{i}^{x}(\bm{x},\bm{y})\differential\bm{y}.$
(38b)

The distribution $\rho_{\bm{x}}^{\infty}(\bm{y})$ specifies the probability
distribution of the fast variables $\bm{y}$ and is stationary on the fast
timescale, but implicitly time-varying since it is parameterized by the state
$\bm{x}$ of the slow variables. In the simplest case $\rho_{\bm{x}}^{\infty}$
might be proportional to a delta function in $\bm{x}$, indicating that the
fast variables are a direct function of the slow variables. Physically this
might correspond to the case where $\bm{y}$ represents a slow amplitude-
dependent deformation of the base flow or phase-locked higher harmonics, for
instance.

More generally the functional form of this distribution might be complicated,
but it is not necessary to specify it in closed form provided the integral in
Eq. 38b can be evaluated. For instance, when the dynamics are linear-quadratic
then Eq. 38b simplifies to first and second moments of the distribution. Then
the key to the MMR closure for the POD-Galerkin is deriving an approximation
of the fast dynamics that is consistent with the original system but allows
for computation of these moments in closed form.

As an aside, although the disappearance of the fictitious small parameter
$\epsilon$ is necessary for consistency with the original Galerkin system, it
does call into question the validity of the perturbation series approximation
(32). In this work we do not attempt to make this more rigorous, but instead
demonstrate by example that it is a useful heuristic capable of resolving
important features of the flow physics. We justify this based on the
observation that the underlying linear-quadratic Galerkin systems often poorly
approximate the true dynamics, so it would not be useful to perfectly match
the original model even if it was possible to do so.

###  4.2 Stochastic averaging for a model system

Before applying the formal solution (38) to the linear-quadratic Galerkin
dynamics, we first illustrate the stochastic averaging procedure on a simple
two-dimensional fast-slow system:

$\displaystyle\derivative{x}{\tau}$ $\displaystyle=\lambda x-xy$ (39a)
$\displaystyle\epsilon\derivative{y}{\tau}$ $\displaystyle=-y+\mu
x^{2}+\sqrt{\epsilon}\sigma w(t).$ (39b)

This SDE system is associated with the following backwards Kolmogorov equation
for an arbitrary scalar observable $g^{\epsilon}(x,y,\tau)$:

$\partialderivative{g^{\epsilon}}{\tau}=\epsilon^{-1}\underbrace{\left[(-y+\mu
x^{2})\partialderivative{y}+\frac{\sigma^{2}}{2}\partialderivative[2]{y}\right]}_{\mathcal{L}_{0}}g^{\epsilon}+\underbrace{\left[(\lambda
x-xy)\partialderivative{x}\right]}_{\mathcal{L}_{1}}g^{\epsilon}.$ (40)

Here $\mathcal{L}_{0}$ is the backwards Kolmogorov operator associated with an
Ornstein-Uhlenbeck process in $y$, which we assume satisfies the ergodicity
requirements (30) and (31). By comparison of $\mathcal{L}^{\dagger}_{0}$ with
the steady-state Fokker-Planck equation for an Ornstein-Uhlenbeck process
(23), the stationary distribution $\rho^{\infty}_{x}$ is a Gaussian with mean
and variance

$\int_{-\infty}^{\infty}y\rho^{\infty}_{x}(y)\differential y=\mu
x^{2},\qquad\qquad\int_{-\infty}^{\infty}y^{2}\rho^{\infty}_{x}(y)\differential
y=\frac{\sigma^{2}}{2}.$ (41)

In this case, the solvability condition (35) reduces to

$\partialderivative{g_{0}}{\tau}=\left[\int_{-\infty}^{\infty}\rho_{x}^{\infty}(y)(\lambda
x-xy)\differential y\right]\partialderivative{g_{0}}{x}=0.$ (42)

Comparing with Eq. (36), this condition corresponds to a Koopman operator for
the averaged dynamics. Using the normalization condition on the probability
distribution $\rho_{x}^{\infty}$ and the mean given by Eq. (41), the fast
variable can be integrated analytically, giving the coarse-grained dynamics

$\displaystyle\derivative{x}{\tau}$
$\displaystyle=\int_{-\infty}^{\infty}\rho_{x}^{\infty}(y)(\lambda
x-xy)\differential y$ (43) $\displaystyle=\lambda x-\mu x^{3}.$ (44)

Finally, we note that in this case the same solution could have been reached
by neglecting the diffusion term in the $y$ dynamics from the beginning, as in
an invariant manifold-type model (Guckenheimer & Holmes, 1983; Noack et al.,
2003; Pavliotis & Stuart, 2012), since due to the form of the slow dynamics
the effect of the noise disappears on averaging. This would give the
approximation $y=\mu x^{2}$, which could then be substituted into the
$x$-dynamics to recover Eq. (43). The difference between that approach and the
one described here lies in Eq. (42); rather than assuming the fast variables
can be expressed as an instantaneous function of the slow variables, they are
assumed to fluctuate and their effect is accounted for by averaging. Later we
will discuss some general conditions under which the two are equivalent, but
in this example if there was a term proportional to $y^{2}$ in the $x$
dynamics, the solvability condition would include the second moment of $y$ and
the diffusion $\sigma$ would appear in the coarse-grained dynamics.

###  4.3 Application to the Galerkin system

As illustrated by the previous example, in order to make practical use of the
averaging procedure for the partitioned Galerkin system (26), we must be able
to perform the integral over the distribution $\rho_{\bm{x}}^{\infty}(\bm{y})$
of the fast variables. This distribution is the solution for fixed $\bm{x}$ to
the steady-state Fokker-Planck equation

$\mathcal{L}^{\dagger}_{0}\rho=-\partialderivative{(\rho
f_{\alpha}^{y})}{y_{\alpha}}+\frac{\sigma_{\alpha}^{2}}{2}\partialderivative[2]{\rho}{y_{\alpha}}=0,$
(45)

corresponding to the linear stochastic process

$\dot{y}_{\alpha}=\left[F^{2}_{\alpha}+L^{21}_{\alpha j}x_{j}+Q^{211}_{\alpha
jk}x_{j}x_{k}\right]+\left[L^{22}_{\alpha\beta}+2Q^{212}_{\alpha
j\beta}x_{j}\right]y_{\beta}+\sigma_{\alpha}w_{\alpha}(t).$ (46)

The stationary distribution is a multivariate Gaussian with mean and
covariance that can be determined analytically from the solution of a Lyapunov
equation (Risken, 1996), but this would need to be done at each value of
$\bm{x}$. Instead, we propose the diagonal drift approximation

$\dot{y}_{\alpha}\approx\nu_{\alpha}(\bm{x})\left(\mu_{\alpha}(\bm{x})-y_{\alpha}\right)+\sigma_{\alpha}w_{\alpha}(t),$
(47)

for which the stationary distribution is a product of univariate Gaussians
solving the 1-dimensional Fokker-Planck equation (23), i.e.
$y_{\alpha}\sim\mathcal{N}(\mu_{\alpha},\sigma_{\alpha}^{2}/2\nu_{\alpha})$.
By comparison with Eq. (46), the conditional mean is naturally defined as

$\mu_{\alpha}(\bm{x})=\sum_{j,k=1}^{r}\nu^{-1}_{\alpha}\left[L^{21}_{\alpha
j}x_{j}+Q^{211}_{\alpha jk}x_{j}x_{k}\right].$ (48)

Here we have omitted the contribution of the constant forcing term
$\mathsfbi{F}^{2}$ so that the approximate fast process $\bm{y}$ preserves the
zero-mean property of the POD coefficients for $\bm{x}=\bm{0}$. In this work
we will make the simplifying assumption that the effective damping
coefficients $\nu_{\alpha}$ are constant, although more generally they could
be functions of $\bm{x}$. Appropriate values for $\nu_{\alpha}$ can be
determined from energy balance, as described below.

With this approximation, the average in Eq. (38) is a straightforward
calculation resulting in the generalized Stuart-Landau model

$\displaystyle\hat{f}_{i}(\bm{x})$
$\displaystyle=F^{1}_{i}+L^{11}_{ij}x_{j}+L^{12}_{i\alpha}\mu_{\alpha}(\bm{x})+Q^{111}_{ijk}x_{j}x_{k}+2Q^{112}_{ij\alpha}x_{j}\mu_{\alpha}(\bm{x})+\nu_{\alpha}^{-1}Q^{122}_{i\alpha\alpha}\sigma_{\alpha}\sigma_{\alpha}/2$
(49a)
$\displaystyle\equiv\hat{F}_{i}+\hat{L}_{ij}x_{j}+\hat{Q}_{ijk}x_{j}x_{k}+\hat{C}_{ijk\ell}x_{j}x_{k}x_{\ell},$
(49b)

with the following closed quantities denoted by a hat:

$\displaystyle\hat{F}_{i}$
$\displaystyle=F^{1}_{i}+\nu_{\alpha}^{-1}Q^{122}_{i\alpha\alpha}\sigma_{\alpha}\sigma_{\alpha}/2$
(50a) $\displaystyle\hat{L}_{ij}$
$\displaystyle=L^{11}_{ij}+\nu^{-1}_{\alpha}L^{12}_{i\alpha}L^{21}_{\alpha j}$
(50b) $\displaystyle\hat{Q}_{ijk}$
$\displaystyle=Q^{111}_{ijk}+\nu^{-1}_{\alpha}\left(L^{12}_{i\alpha}Q^{211}_{\alpha
jk}+2Q^{112}_{ij\alpha}L^{21}_{\alpha k}\right)$ (50c)
$\displaystyle\hat{C}_{ijk\ell}$
$\displaystyle=2\nu^{-1}_{\alpha}Q^{112}_{ij\alpha}Q^{211}_{\alpha k\ell}.$
(50d)

We will refer to the secondary reduction of a standard rank-$r$ POD-Galerkin
model to the cubic model (49b) as the multiscale model reduction (MMR)
approach to closure modeling.

##### Interpretation of the cubic closure model

The terms in the closure model are computed by summing the Galerkin tensors
over the fast variables. As a mnemonic for these averaged quantities, the
superscripts for the partitioned system could be thought of as tensor
contractions with the modification of the damping $\nu_{\alpha}^{-1}$. For
example, the cubic term $Q^{112}_{ij\alpha}\nu^{-1}_{\alpha}Q^{211}_{\alpha
k\ell}$ is suggestive of the slow variables “filtering through” the fast
dynamics via the quadratic interactions $\mathsfbi{Q}^{211}$ and
$\mathsfbi{Q}^{112}$.

Although the constant, linear, and quadratic terms are all modified as a
result of the stochastic averaging procedure, the appearance of this cubic
term is perhaps the most noteworthy. It represents the leading-order
contribution of the fast-slow nonlinear interaction in the slow dynamics due
to the slow-slow interaction in the fast dynamics. We will explore this in
more detail in the following section, but it is a generalization of the weakly
nonlinear Stuart-Landau mechanism (Stuart, 1958; Landau & Lifshitz, 1959) that
is capable of resolving the stabilizing influences of both mean flow
deformation and the energy cascade.

The final ingredient in the multiscale closure model (49b) is the
determination of the damping and diffusion coefficients $\nu_{\alpha}$ and
$\sigma_{\alpha}$. A reasonable approximation for the diffusion could be the
mean of the neglected term
$\sigma_{\alpha}=\sum_{\beta}Q^{222}_{\alpha\beta\beta}\overline{y_{\beta}^{2}},$
although in our numerical examples we find little difference with neglecting
the diffusion altogether. For the damping term we apply an energy balance
condition on the closed model, derived from the assumption that the system is
statistically stationary so that the average variation of kinetic energy
vanishes (Noack et al., 2011):

$\derivative{t}\overline{x_{i}^{2}}=\overline{\hat{f}_{i}(\bm{x})x_{i}}=0.$
(51)

When the closed tensors (50) are substituted for $\hat{\bm{f}}(\bm{x})$, this
simplifies to a linear system of equations for $\bm{\nu}^{-1}$. The energy
balance approximation does not require statistics of the fast variables, but
due to the cubic term it does require well-converged fourth moments of the
slow variables. In particular, the vector $\bm{\nu}_{\mathrm{i}nv}$ which is
the element-wise inverse of $\bm{\nu}$ is the solution to the linear system

$\mathsfbi{A}\bm{\nu}_{\mathrm{i}nv}+\bm{b}=0,$ (52a)
$\displaystyle\begin{split}A_{i\alpha}&=\sum_{j=1}^{r_{0}}L^{12}_{i\alpha}L^{21}_{\alpha
j}\overline{a_{i}a_{j}}+\sum_{j,k=1}^{r_{0}}L^{12}_{i\alpha}Q^{211}_{\alpha
jk}\overline{a_{i}a_{j}a_{k}}\\\ &\hskip
56.9055pt+2\sum_{j,k=1}^{r_{0}}Q^{112}_{ij\alpha}L^{21}_{\alpha
k}\overline{a_{i}a_{j}a_{k}}+2\sum_{j,k,\ell=1}^{r_{0}}Q^{112}_{ij\alpha}Q^{211}_{\alpha
k\ell}\overline{a_{i}a_{j}a_{k}a_{\ell}}\end{split}$ (52b) $\displaystyle
b_{i}$
$\displaystyle=L^{11}_{ii}\overline{a_{i}^{2}}+\sum_{j,k=1}^{r_{0}}Q^{111}_{ijk}\overline{a_{i}a_{j}a_{k}}$
(52c)

In general the damping should be positive since the linear stochastic
approximation (46) and its diagonal simplification is only plausible if the
mean $\bm{\mu}(\bm{x})$ is linearly stable for all $\bm{x}$. This is
physically consistent with the energy cascade picture, in which we expect that
energy will primarily flow “downhill” from the slow variables, representing
large coherent structures and global instabilites, to the fast variables,
representing the smaller, dissipative scales of the flow. On this basis,
negative values of the damping coefficients can also be set to zero to avoid
introducing unphysical instabilities in the cubic term, for example.

We conclude the discussion of the application to Galerkin-type systems with a
note on the computational scaling of the closure model. While the simulation
of the original linear-quadratic model is dominated by the quadratic term,
which requires $\mathcal{O}(r^{3})$ operations to evaluate, the cubic term in
the closure model (49b) requires $\mathcal{O}(r_{0}^{4})$ operations to
evaluate. In some cases this may mean that the cubic model is actually more
expensive to simulate than the original Galerkin system, although whether or
not this is the case will depend on the specific values of $r$, the dimension
of the linear-quadratic system, and $r_{0}$, the dimension of the closed
model. For the examples given in Sec. 6 either $r_{0}^{4}<r^{3}$, as in Secs.
6.1, 6.2 and 6.3, or the two are the same order of magnitude, as in Sec. 6.4.
As we will show, the primary advantage of this approach is not necessarily
that it is faster to simulate than the usual POD-Galerkin system, but that it
is more stable and physically faithful with many fewer modes.

###  4.4 Relationship to invariant manifold modeling

Invariant manifold modeling is a method in which the fast variables are
assumed to depend algebraically on the slow variables. It is often applied
near the threshold of bifurcation, when one may assume that there are only a
small number of weakly unstable directions in the modal subspace and that the
modal coefficients are typically small.

Under these assumptions, the fast variables can be approximated with a power
series expansion with coefficients determined by a dynamical consistency
condition (Guckenheimer & Holmes, 1983; Pavliotis & Stuart, 2012). For
example, if the dynamics are close to a supercritical Hopf bifurcation then
the leading-order approximation is a parabolic manifold that depends
quadratically on the slow variables; replacing the fast variables with this
approximation results in a cubic evolution equation for the slow variables,
similar to Eq. (49b). This approach is closely related to the technique of
center manifold reduction, which is typically applied to a system at the
bifurcation point, yielding normal form-type amplitude equations for the
evolution of the instability modes (Coullet & Spiegel, 1983; Carini et al.,
2015). Here we present the parabolic invariant manifold approximation to give
a comparison with the stochastic averaging method, but we will conclude by
showing that a different manifold model gives an equivalent result to Eq. (50)
for $\bm{\sigma}=\bm{0}$.

We first shift the origin of the Galerkin system (7) to a fixed point
$\bm{a}^{0}$ defined by $\bm{f}(\bm{a}^{0})=\bm{0}$. Defining fluctuations
$\bm{a}^{\prime}(t)=\bm{a}(t)-\bm{a}^{0}$, the Galerkin system is

$\displaystyle\dot{a}_{i}^{\prime}$
$\displaystyle=L_{ij}^{\prime}a_{j}^{\prime}+Q_{ijk}a_{j}^{\prime}a_{k}^{\prime}$
(53a) $\displaystyle L_{ij}^{\prime}$
$\displaystyle=L_{ij}+(Q_{ijk}+Q_{ikj})a^{0}_{k}.$ (53b)

We furthermore assume that the linear term $\mathsfbi{L}^{\prime}$ is
diagonalizable with an invertible set of eigenvectors $\mathsfbi{V}$ and
eigenvalues $\bm{\lambda}=\mathrm{diag}(\bm{\Lambda})$, such that
$\bm{a}^{\prime}=\mathsfbi{V}\bm{b}$. The similarity transform defined by
$\mathsfbi{V}$ diagonalizes $\mathsfbi{L}^{\prime}$ so that Eq. (53a) becomes

$\displaystyle\dot{\bm{b}}$
$\displaystyle=\mathsfbi{V}^{-1}\bm{f}^{\prime}(\mathsfbi{V}\bm{b})\equiv\bm{f}^{\Lambda}(\bm{b})$
(54a)
$\displaystyle=\bm{\Lambda}\bm{b}+\mathsfbi{Q}^{\Lambda}(\bm{b},\bm{b}),$
(54b)

where $Q^{\Lambda}_{ijk}=V^{-1}_{i\ell}Q_{\ell mn}V_{mj}V_{nk}$ and the
superscripts $\Lambda$ indicates that this tensor corresponds to the
diagonalized coordinates $\bm{b}$. As with the multiscale averaging method
described in Sec. 4.1, we then partition the state $\bm{b}$ into slow
variables $\bm{x}$ and fast variables $\bm{y}$, so that
$\bm{b}=\begin{bmatrix}\bm{x}^{T}&\bm{y}^{T}\end{bmatrix}^{T}$ and
$\bm{f}^{\Lambda}=\begin{bmatrix}\left(\bm{f}^{\bm{x}}\right)^{T}&\left(\bm{f}^{\bm{y}}\right)^{T}\end{bmatrix}^{T}$.
Similarly, the eigenvalues $\bm{\lambda}$ and quadratic tensor
$\mathsfbi{Q}^{\Lambda}$ can be partitioned into $\mathsfbi{Q}^{111}$,
$\mathsfbi{Q}^{112}$, etc., as before.

The central assumption of the invariant manifold model is that the fast
variables can be expressed as an instantaneous function
$\bm{y}\approx\bm{h}(\bm{x})$ of the slow variables. Dynamical consistency for
the fast variables then requires that, to leading order,

$\displaystyle\derivative{t}\bm{y}=\bm{f}^{\bm{y}}(\bm{x},\bm{h}(\bm{x}))$
$\displaystyle=\derivative{t}\bm{h}(\bm{x})$ (55a)
$\displaystyle=\mathbfcal{D}\bm{h}(\bm{x})\bm{f}^{\bm{x}}(\bm{x},\bm{h}(\bm{x}))$
(55b)

where $\mathbfcal{D}\bm{h}(\bm{x})$ is the Jacobian of $\bm{h}$ evaluated at
$\bm{x}$ (Guckenheimer & Holmes, 1983).

Assuming all of the state variables in $\bm{b}$ are typically small, the
manifold equation can be expanded as a power series in $\bm{x}$. Preservation
of the fixed point at the origin requires $\bm{h}(\bm{0})=\bm{0}$, and in the
case of a center manifold reduction where the real parts of the slow
eigenvalues vanish, the condition that the $\bm{h}(\bm{x})$ be tangent to the
center manifold at the origin also implies that
$\mathbfcal{D}\bm{h}(\bm{0})=\bm{0}$. Under these conditions the first
nontrivial terms in the expansion are quadratic and the expansion coefficients
can be determined by matching powers of $\bm{x}$ in the consistency condition
(55b).

The connection between this method and stochastic averaging can be made clear
by viewing $\bm{y}(t)=\bm{h}(\bm{x},t)$ as an explicit function of time and
rewriting the dynamical consistency condition as a Koopman generator acting on
the $\alpha^{\mathrm{th}}$ component of the manifold equation:

$\partialderivative{h_{\alpha}}{t}=\bm{f}^{\bm{x}}(\bm{x},\bm{h}(\bm{x}))\cdot\nabla
h_{\alpha}.$ (56)

Since the observable in such an equation is arbitrary, this is equivalent to
the averaged backwards Kolmogorov equation (37) if

$f_{i}^{\bm{x}}(\bm{x},\bm{h}(\bm{x}))=\int_{\mathcal{Y}}\rho_{\bm{x}}^{\infty}(\bm{y})f_{i}^{x}(\bm{x},\bm{y})\differential\bm{y}.$
(57)

In the limit of zero diffusion in the fast process, the stationary
distribution approaches a Dirac delta function
$\rho_{\bm{x}}^{\infty}=\delta^{r-r_{0}}(\bm{y}-\bm{\mu}(\bm{x}))$ and the two
are equivalent with $\bm{h}(\bm{x})=\bm{\mu}(\bm{x})$. With the diagonal drift
approximation (47) and constant damping $\bm{\nu}$, the zero-diffusion limit
of the proposed Galerkin closure in Sec. 4.3 can be seen as an invariant
manifold model with $\bm{h}(\bm{x})$ defined as the conditional mean (48).
Because we derived the stochastic averaging model with the assumption of real-
valued coefficients and the parabolic model uses the complex-valued
eigenvector coordinates, it is less obvious that the parabolic invariant
manifold corresponds to a particular choice of $\bm{\nu}$ and $\bm{\mu}$ in
Eq. (47).

Although the invariant manifold reduction is more theoretically
straightforward than the stochastic averaging closure, it generally requires
stronger assumptions on the system. In particular, the development of the
parabolic manifold model requires that the origin be a fixed point and that
the real part of the slow eigenvalues is zero. This does not raise any
difficulties for simple weakly unstable configurations where the Galerkin
model preserves fixed points of the fluid flow, or if the steady state is used
as the base of the modal expansion. However, when the flow is far from a
bifurcation point, such as for a fully turbulent flow, the mean flow is often
a more relevant base of the expansion. In this case the resulting Galerkin
models may not necessarily have fixed points that coincide with meaningful
states of the underlying system. In addition, the invariant manifold reduction
relies on smallness of the coefficients (justifying the power series
expansion) and the assumption that the higher-order modes can be expressed as
instantaneous algebraic functions of the active coefficients. For chaotic or
turbulent dynamics neither of these is necessarily justified.

###  4.5 Finite-time thermodynamics

Before concluding the discussion of the multiscale closure model, we give a
brief comparison to the finite-time thermodynamics (FTT) framework introduced
by Noack et al. (2008). This approach begins from the linear-quadratic POD-
Galerkin system, Eq. (7), but rather than assuming that the unresolved
variables are the higher-order modal coefficients, the FTT model decomposes
the modal coefficients $\bm{a}(t)\in\mathbb{R}^{r}$ into an “ensemble-
averaged” value $\bm{x}(t)\in\mathbb{R}^{r}$ and a “fluctuating” value
$\bm{y}(t)\in\mathbb{R}^{r}$ in the spirit of the Reynolds decomposition. As a
result, both the resolved and unresolved variables correspond to the _same_
spatial modes, and the modal expansion (2) becomes

$\bm{u}(\bm{x},t)\approx\bm{\psi}_{0}(\bm{x})+\sum_{i=1}^{r}(x_{i}(t)+y_{i}(t))\bm{\psi}_{i}(\bm{x}).$
(58)

The second major difference is the assumption of the behavior of the
unresolved variables $\bm{y}$ on ensemble averaging. The ensemble average is
with respect to the probability distribution $\rho$ of the fluctuations
$\bm{y}$. The MMR method assumes that this distribution is parameterized
explicitly by the slow variables $\bm{x}$, so that
$\rho=\rho(\bm{y};\bm{x}(t))$ (written as $\rho_{\bm{x}}^{\infty}(\bm{y})$ in
Sec. 4.1). On the other hand, the Reynolds-averaging analogy in FTT leads to
the assumption that the distribution $\rho$ may be time-varying but does not
depend on $\bm{x}$, so that $\rho=\rho(\bm{y};t)$.

| MMR | FTT
---|---|---
$\rho$ | $\rho(\bm{y};\bm{x})$ | $\rho(\bm{y};t)$
$\mathop{\mathbb{E}}[y_{i}]$ | $\mu_{i}(\bm{x})$ | 0
$\mathop{\mathbb{E}}[y_{i}^{2}]$ | $\sigma_{i}^{2}$ | $2E_{i}(t)$
Table 1:  Statistics of unresolved variables for the multiscale and finite-
time thermodynamics models. MMR assumes the distribution is parameterized by
the slow variables $\bm{x}$, while FTT assumes it is independent and zero-
mean, although with possibly fluctuating variance. We have found the variance
$\sigma_{i}^{2}$ to be negligible in MMR if it is assumed constant, but more
generally it could be also modeled as a function of $\bm{x}$. Note that the
unresolved variables $\bm{y}$ are also defined differently between the two
methods (see Sec. 4.5).

By definition, the fluctuating amplitudes $\bm{y}$ have zero mean with respect
to this distribution, so that

$\mathop{\mathbb{E}}[y_{i}]\equiv\int_{\mathcal{Y}}y_{i}\rho(\bm{y};t)\differential\bm{y}=0,\qquad
i=1,2,\dots r.$ (59)

Here we use the notation $\mathop{\mathbb{E}}[\cdot]$ to denote the ensemble
average with respect to the unknown distribution $\rho(\bm{y},t)$ to
distinguish it from the usual time average of observed values used in the rest
of this work.

As a result of Eq. (59), the mean and fluctuating amplitudes are linearly
uncorrelated with respect to this average:

$\mathop{\mathbb{E}}[x_{i}y_{j}]=0,\qquad i,j=1,2,\dots r.$ (60)

The normalization of the probability distribution also implies that
$\mathop{\mathbb{E}}[x_{i}]=x_{i}$. Finally, FTT assumes that the fluctuations
themselves are uncorrelated, with the usual variance-kinetic energy
relationship for POD coefficients. Combining these assumptions, the second-
order statistics of the POD coefficients with respect to the ensemble average
defined by $\rho(\bm{y};t)$ are

$\mathop{\mathbb{E}}[a_{i}a_{j}]=x_{i}x_{j}+2E_{i}(t)\delta_{ij},$ (61)

where $E_{i}(t)=\mathop{\mathbb{E}}[y_{i}^{2}]/2$. A comparison of the
assumptions on the distributions is shown in Table 1. Note that since the
variance $E_{i}$ must be non-negative, Eq. (61) implies that the statistics of
$\bm{a}$ and $\bm{x}$ are not identical unless $\bm{E}(t)=0$ for all time.

Expanding the POD coefficients into mean and fluctuating component and taking
the ensemble average gives a set of $2r$ evolution equations for the mean
values $\bm{x}$ and the fluctuation kinetic energy $\bm{E}$:

$\displaystyle\dot{x}_{i}$
$\displaystyle=\sum_{j=1}^{r}L_{ij}x_{j}+\sum_{j,k=1}^{r}Q_{ijk}x_{j}x_{k}+2\sum_{j=1}^{r}Q_{ijj}E_{j}$
(62a) $\displaystyle\dot{E}_{i}$
$\displaystyle=2\sum_{j=1}^{r}\left(L_{ii}+2Q_{iij}x_{j}\right)E_{j}+\sum_{j,k=1}^{r}Q_{ijk}\mathop{\mathbb{E}}[y_{i}y_{j}y_{k}].$
(62b)

The remainder of the FTT method centers on deriving a nonlinear closure model
for the third-order statistics $\mathop{\mathbb{E}}[y_{i}y_{j}y_{k}]$; see
Noack et al. (2008) for details.

The main similarity between the FTT and multiscale methods is that both seek
to approximate the effect of unresolved fluctuations on the resolved
variables, resulting in nonlinear closure models. However, the two are
philosophically different in that MMR assumes the unresolved variables
correspond to different spatial structures and evolve on faster time scales,
while FTT models them as linearly uncorrelated fluctuations associated with
the same spatial modes as the resolved variables.

Noack et al. (2008) showed that the FTT model is quite successful at
predicting post-transient energy levels in several flow configurations,
including the cylinder wake and a turbulent homogeneous shear flow. This
method also has the advantage that it does not assume any separation of scales
in amplitude or time. On the other hand, there are some theoretical
difficulties in working with the FTT framework, including that Eq. (61)
muddies the relationship between the second-order statistics of the modal
coefficients $\bm{a}$ and $\bm{x}$, which are foundational to the POD method.
FTT also increases the state dimension by a factor of two and cannot be used
to derive the generalized Stuart-Landau equations that directly account for
effects such as mean flow deformation and the energy cascade. However,
Schlegel et al. (2009), Noack et al. (2011), and Cordier et al. (2013) showed
that a simplified nonlinear closure model based on Eq. (62) has been shown to
stabilize POD-Galerkin models of a mixing layer similar to that examined in
Sec. 6.4 without increasing the state dimension (see App. C).

##  5 Flow configurations

In this work we consider models of three fluid flows: the canonical flow past
a cylinder at Reynolds number 100, a lid-driven cavity flow, and an
incompressible mixing layer with two domain extents. We construct projection-
based models of these flows based on direct numerical simulation (DNS) using
the open-source Nek5000 spectral element solver (Fischer et al., 2008). The
semi-implicit time-stepping integrates diffusive terms with third order
backwards differentiation and convective terms with a third order
extrapolation.

Once the DNS is complete, we estimate the POD modes and coefficients using the
modred library, which provides a set of parallelized, high-performance
algorithms in Python for linear modal decompositions and model reduction
(Belson et al., 2014). We then compute the components of the Galerkin system
(7) by extracting the weight matrix and necessary gradients from the POD modes
with the built-in Nek5000 post-processor. Both the linear-quadratic Galerkin
model (7) and the cubic closure (49b) are simulated using the LSODA time-
stepping algorithm available from scipy with analytic Jacobian matrices for
each system.

###  5.1 Flow past a circular cylinder

The vortex shedding in the wake behind a cylinder at Reynolds number 100,
based on the free-stream velocity and cylinder diameter, is a canonical flow
configuration for reduced-order modeling (Noack et al., 2003) and it is shown
in Fig. 3. We simulate this flow on a domain of 2600 sixth-order spectral
elements on $x,y\in(-5,15)\times(-5,5)$ refined close to the cylinder wall;
further details and analysis can be found in Loiseau et al. (2018a, b). We
first perform a global stability analysis of the the unstable steady state,
determined by solving the stationary Navier-Stokes equations with selective
frequency damping (Åkervik et al., 2006), using a Krylov-Schur time-stepping
algorithm (Loiseau et al., 2019). The transient simulation is then initialized
with the unstable steady state perturbed by the least-stable global eigenmode,
normalized so that its energy is $10^{-4}$.

Figure 3:  Vortex shedding in the wake of a circular cylinder at
$\mbox{{Re}}=100$. En route to the periodic post-transient flow, Reynolds
stresses deform the background flow from the unstable steady state to the
marginally stable mean flow. This can be approximated in a model by augmenting
the POD basis with the “shift mode” $\bm{\psi}_{\Delta}$. The usual POD modes
show the typical structure of spatial and temporal harmonics describing
periodically advecting flow features.

The flow regime at Reynolds number 100 is beyond the Hopf bifurcation leading
to vortex shedding but below the onset of three-dimensional instability (Noack
& Eckelmann, 1994; Barkley & Henderson, 1996; Ma & Karniadakis, 2002). The
flow is also far enough from the threshold of bifurcation that weakly
nonlinear analyses fail to predict the amplitude and frequency of the vortex
shedding limit cycle, although a variety of techniques have been used to
derive accurate models for the transient and post-transient wake (Noack et
al., 2003, 2011; Mantic̆-Lugo et al., 2014; Loiseau & Brunton, 2018; Loiseau
et al., 2018b).

Following Noack et al. (2003) and Loiseau et al. (2018a), we compute POD modes
from 100 snapshots of the mean-subtracted velocity field, equally spaced over
one period of vortex shedding; the pressure gradient term can shown to be
negligible for Galerkin models. In order to construct a basis that can span
the deformation between the unstable steady state and mean flow, we augment
the POD basis with the “shift mode”, visualized in Fig. 3. This additional
mode is computed from the difference between the base and mean flows and then
orthonormalized with respect to the rest of the POD modes with the Gram-
Schmidt process. We denote the coefficient of this mode by $a_{\Delta}(t)$ to
distinguish it from the usual POD modes.

The full expansion basis then consists of eight POD modes capturing the limit-
cycle vortex shedding, plus the shift mode. We take the unstable steady state
to be the base flow of the modal expansion so that the origin is a fixed point
of the POD-Galerkin system. This basis accurately represents the kinematics of
the post-transient flow, with a normalized energy residual of
$\mathcal{O}(10^{-4})$ (Loiseau et al., 2018a).

###  5.2 Lid-driven cavity flow

Figure 4:  Regularized lid-driven cavity flow at $\mbox{{Re}}=20,000$. The
singular value spectrum is slow to converge (left), indicating that relatively
many modes are necessary for an accurate kinematic approximation. The leading
POD modes themselves appear roughly in pairs with similar spatial structure
and frequency content. As for the cylinder wake shown in Fig. 3, this usually
indicates oscillatory dynamics, although the evolution of the temporal
coefficients is much more irregular in this chaotic flow than for the periodic
wake.

Lid-driven cavity flow is another idealized geometry that serves as a
benchmark problem for numerical methods and model reduction (Arbabi & Mezić,
2017; Arbabi & Sapsis, 2019; Lee et al., 2019; Rubini et al., 2020). The flow,
shown in Fig. 4, takes place on a square domain with
$(x,y)\in(-L,L)\times(-L,L)$ with $L=1$. A velocity profile is imposed on the
upper boundary, given by

$u(x,L)=U(1-x^{2})^{2}.$ (63)

This driving velocity profile roughly approximates the shear-driven flow over
an open cavity at a much lower computational cost, although the dynamics of
the shear- and lid-driven cavities are different. The regularized fourth-order
velocity profile on the lid avoids numerical problems where the lid meets the
no-slip boundaries on the other walls.

The standard Reynolds number for this flow is defined as $\mbox{{Re}}=UL/\nu$
based on the cavity half-height and forcing velocity. The flow follows the
stereotypical Ruelle-Takens route to chaos parameterized by Reynolds number,
undergoing a Hopf bifurcation at $\mbox{{Re}}_{c}^{(2)}\approx 10,250$, a
secondary Hopf bifurcation to quasiperiodic flow at
$\mbox{{Re}}_{c}^{(2)}\approx 15,500$, and a final bifurcation to chaos near
$\mbox{{Re}}_{c}^{(3)}\approx 18,000$ (Lee et al., 2019). We simulate the flow
in the chaotic regime at $\mbox{{Re}}=20,000$.

The numerical configuration is similar to that described by Arbabi & Sapsis
(2019). The domain is discretized with 50 seventh-order spectral elements in
each direction, refined towards the walls. The flow is integrated for 3000
time units, saving snapshots at a sampling rate $\Delta t=0.1$. Again, the
pressure gradient term vanishes for this closed flow (Noack et al., 2005), so
the velocity POD modes are computed from the first 5000 snapshots, with the
remaining 25000 snapshots retained for the statistical comparison in Sec. 6.3.
Figure 4 shows a snapshot of the DNS along with the leading POD modes and the
singular value spectrum.

###  5.3 Incompressible mixing layer

\begin{overpic}[width=433.62pt]{figs/mix_modes.pdf} \end{overpic} Figure 5:
Proper orthogonal decomposition (POD) applied to the mixing layer. The primary
shear layer instability forms Kelvin-Helmholtz waves that roll up into
vortices, which in turn undergo successive vortex pairing (top). Modes
computed on both the short and long domains reveal modes related to the
dominant flow features: the shear layer instability and the downstream vortex
pairing. Although the vortex pairing is secondary to the upstream instability,
on the longer domain it accounts for a larger proportion of the fluctuation
kinetic energy. Higher-order modes not pictured here include harmonics,
nonlinear crosstalk, and modes related to irregularity in the location of the
vortex pairing events.

As a more difficult test of the proposed MMR closure method, we also examine
models of an incompressible mixing layer, pictured in Fig. 11. The mixing
layer is a canonical example of a free shear flow exhibing convective
instability (Monkewitz & Huerre, 1982; Huerre & Monkewitz, 1985, 1990). Inlet
forcing is used to excite Kelvin-Helmholtz instability waves, which roll up
into approximately discrete vortices. These vortices are subject to a
secondary instability (Brancher & Chomaz, 1997) and undergo successive helical
pairing at the subharmonic of the primary Kelvin-Helmholtz frequency, a
process that drives the linear growth of the mixing layer (Winant & Browand,
1974). At higher Reynolds numbers and in three dimensions, the turbulent
mixing layer is dominated by the linear growth of coherent structures similar
to those seen in 2D simulation (Brown & Roshko, 1974; Stanley & Sarkar, 1997),
which play a major role in mixing, transport, and entrainment of the turbulent
flow (Hussain, 1981). From a global perspective, the flow behaves as an
amplifier, where strong non-normality in the linear operator results in
transient algebraic energy growth of small perturbations (Chomaz, 2005).

By analogy to numerical methods, modeling an advection-dominated flow with a
POD-Galerkin approach is generally ill-advised in much the same way that
pseudospectral methods are not ideal for solving hyperbolic problems. Although
the vortices shed in the flow past a cylinder are also traveling waves in a
sense, the Kelvin-Helmholtz instability in the mixing layer is convective in
nature as opposed to the absolute instability in the cylinder wake. For these
reasons, Galerkin projection onto a global modal basis such as POD is not a
natural way to approximate shear flows like the mixing layer. Nevertheless,
both free and wall-bounded shear layers are present in a wide variety of
flows, so in order to employ a low-dimensional dynamical systems model it is
important that it be able to capture these features. The challenge of modeling
advection-driven phenomena in a global reduced-order modeling framework
motivates the use of similar configurations in several studies of
stabilization methods for Galerkin-type systems (Ukeiley et al., 2001; Noack
et al., 2005; Cordier et al., 2013; Balajewicz et al., 2013).

Since the flow becomes increasingly complex downstream, the difficulty of
modeling the flow can be tuned by varying the streamwise extent of the domain.
We construct models on both short ($L_{x}=150$) and long ($L_{x}=250$)
domains, as shown in Fig. 5. Over the extent of the short domain, the flow is
mainly characterized by the initial instability and vortex rollup, with the
earliest signs of vortex pairing at $x\gtrsim 130$. On the full domain, the
vortex pairing accounts for a much larger fraction of the fluctuation kinetic
energy, as shown by the POD analysis below.

The inlet profile is a standard hyperbolic tangent:

$U(y)=\bar{U}+(\Delta U/2)\tanh 2y,$ (64)

where $\bar{U}=(U_{1}+U_{2})/2\equiv 1$ and $\Delta U=U_{1}-U_{2}$ in terms of
the free-stream velocities $U_{1}$ and $U_{2}$ above and below the layer,
respectively, and the length scale is nondimensionalized by the initial
vorticity thickness $\delta_{\omega}={\Delta U}/{U^{\prime}(0)}.$ We define a
Reynolds number $\mbox{{Re}}=\delta_{\omega}(0)\Delta U/\nu$ based on intial
vorticity thickness, velocity difference $\Delta U$ across the layer, and
kinematic viscosity $\nu$ and set $\mbox{{Re}}=500$ in the simulation. The
domain consists of 5200 eleventh-order spectral elements on
$x,y\in(-10,300)\times(-50,50)$, but the full streamwise extent is masked to
$x,y\in(0,250)\times(-20,20)$ for modeling in order to discount any boundary
effects. The domain is shown in Fig. 11.

This configuration roughly follows that described by Balajewicz et al. (2013);
Cordier et al. (2013), with two key differences. First, we simulate
incompressible flow, while these authors simulate isothermal compressible flow
at low Mach numbers ($\mathrm{Ma}\sim 0.1$). Consequently, we replace the
“sponge” boundary conditions for far-field acoustic waves with a region of
increased viscous damping corresponding to $\mbox{{Re}}=50$ for $x>250$ to
avoid numerical instability at the outflow. Second, rather than disturbing the
inlet with random solenoidal perturbations, we use eigenfunction forcing
described in App. A, similar to that employed by Colonius et al. (1997).

The DNS is run until a final time of $t=2820$, corresponding to fifty periods
of the lowest-frequency component of the eigenfunction forcing, saving at
every tenth time step with $\Delta t=0.00705$. The POD modes are then computed
from the first 10% of these snapshots (4000 fields), with the remainder
reserved for statistical comparison with the models. The domain truncation can
be accomplished by simply setting elements of the weight matrix $\mathsfbi{W}$
in Eq. (4) to zero for mesh locations with $x>L_{x}$ so that these locations
do not contribute to the statistics in the correlation matrix. For the mixing
layer we compute modes with both velocity and pressure fields, as described in
Sec. 2.1, since the pressure gradient term is not necessarily negligible in
Galerkin models of free shear flows (Noack et al., 2005).

Figure 5 shows the leading POD eigenvalues and modes for both domains. In both
cases the bulk of the fluctuation kinetic energy is contained in the first
four modes, with the remainder of the spectrum decaying relatively slowly. As
expected, the dominant mode pair for the shorter domain corresponds to the
Kelvin-Helmholtz waves, while the second pair represents the onset of vortex
pairing at roughly half the spatial wavenumber. Higher-order modes are either
harmonics (such as the third mode pair) or less regular downstream structures
related to slight variations in the vortex pairing.

The directed graph of energy transfers in Fig. 1 is constructed from the modes
computed on the short domain. Specifically, the leading mode pair and its
first three harmonics are shown for a straightforward visualization of the
energy cascade; the nodes labeled $a_{1}-a_{8}$ thus actually correspond to
mode pairs $(\bm{\varphi}_{1},\bm{\varphi}_{2})$,
$(\bm{\varphi}_{5},\bm{\varphi}_{6})$,
$(\bm{\varphi}_{17},\bm{\varphi}_{18})$, and
$(\bm{\varphi}_{24},\bm{\varphi}_{25})$. The thickness of the edge from node
$a_{i}$ to $a_{j}$ is proportional to
$\sum_{k}Q_{ijk}\overline{a_{i}a_{j}a_{k}}$ for the quadratic Galerkin tensor
$\mathsfbi{Q}$, which is roughly related to average nonlinear energy transfer
from one mode to the other, although it should not necessarily be interpreted
in any rigorous way beyond the purposes of conceptual visualization of the
energy cascade.

On the longer domain the vortex pairing makes up a larger proportion of the
fluctuation kinetic energy. Consequently, the leading four modes are spatially
localized downstream, representing vortex pairing, while the upstream linear
instability is not evident until the third mode pair, even though the vortex
pairing is secondary to the shear layer instability from a physical
perspective.

The slowly decaying POD spectrum and wavelike spatial structure of the modes
are a consequence of the advection-driven nature of this flow. It is well
known that the space-time separation of variables assumed by POD is not a
natural representation of traveling waves. This is one reason that free shear
flows have long posed difficulties for POD-Galerkin modeling, making this a
challenging test case for the MMR approach.

##  6 Results

In this section we evaluate the proposed multiscale model reduction (MMR)
closure method on several numerical test problems. Before modeling the flows
described in Sec. 5, we first explore the capability of MMR to resolve
stereotypical subscale stabilization mechanisms in simple analytically
tractable model problems in Sec. 6.1, namely for the mean-field model of Noack
et al. (2003) and for Burgers’ equation.

In Sec. 6.2 we analyze vortex shedding behind a circular cylinder, showing
that the method reproduces the classic invariant manifold model due to Noack
et al. (2003), which must account for both mean flow deformation and the
energy cascade in order to accurately model the transient dynamics. As a more
challenging test in Sec. 6.3, we show that the multiscale closure stabilizes a
reduced-order model of the chaotic lid-driven cavity flow, while closely
matching the spectral energy content of the DNS. Finally, in Sec. 6.4 we
develop a model for convective instability and subharmonic vortex pairing in
the spatially evolving mixing layer.

By means of these examples, we show that the proposed approach for multiscale
closure is a systematic method for eliminating degrees of freedom from
reduced-order models that can resolve critical subscale phenomena and produce
stable models of dynamically complex unsteady flows. More broadly, the
multiscale framework represents a general description of the means by which
effective cubic nonlinearity arises in low-dimensional models.

###  6.1 Analytically tractable model problems

There are two principal mechanisms by which nonlinear interactions with
higher-order modes effectively stabilize the flow: Reynolds stress-induced
mean flow deformation and the energy cascade. Before demonstrating the MMR
closure introduced in Sec. 4 on DNS of unsteady fluid flow, we first confirm
that it can resolve both of these effects by examining simple representative
analytic models.

####  6.1.1 Nonlinear stability in a mean-field model

Nonlinear stability analysis was born out of the observation that the
exponential energy growth predicted by linear theory when the flow has at
least one unstable eigenvalue is untenable at long times. Furthermore, the key
assumption of infinitesimal perturbation amplitudes is violated once the
instability grows to finite amplitude. It is conceivable that the energy
growth eventually leads to a breakdown into turbulence, but it was observed
experimentally (e.g. Taylor (1923)) that this is not always the case.

Instead, if the perturbations are assumed to be of small but finite amplitude,
the linear stability problem can be considered the first term in an asymptotic
expansion. The next leading-order terms suggest that the Reynolds stresses
from the nonlinear self-interaction of the perturbation tends to deform the
slowly-varying background flow from the unstable steady-state towards a
neutrally stable mean flow (Stuart, 1958). This type of weakly nonlinear
analysis has since been applied to a wide range of flows and bifurcation
scenarios (Sipp & Lebedev, 2007; Meliga et al., 2009; Meliga & Chomaz, 2011;
Rigas et al., 2017). Although this has proven to be a highly successful
approach, it is only strictly applicable near the threshold of bifurcation and
requires invasive access to the various linear operator “building blocks” of
the governing equations. The self-consistent mean flow modeling framework
(Mantic̆-Lugo et al., 2014; Mantic̆-Lugo & Gallaire, 2016; Meliga, 2017; Rigas
et al., 2021) is more general in terms of parametric dependence, but it does
not lead to a low-dimensional dynamical system.

Resolving mean flow deformation has been shown to be a key ingredient in
stabilizing reduced-order models. Models that only include POD modes derived
from expansion around the mean give accurate approximations on short time
horizons, but they eventually tend towards structural instability (Deane et
al., 1991; Ma & Karniadakis, 2002). However, explicitly incorporating the
stabilizing feedback mechanism via Reynolds stress modeling (Aubry et al.,
1988; Holmes et al., 1996), an augmented POD basis (Noack et al., 2003;
Luchtenburg et al., 2009; Noack et al., 2011), or a normal form ansatz (Deng
et al., 2020) in models of globally unstable flows, has significantly improved
their robustness. Here we show that if the POD basis spans the mean flow
deformation, the proposed multiscale model reduction method can resolve its
stabilizing effect without introducing additional degrees of freedom. This is
consistent with the invariant manifold reduction approach of Noack et al.
(2003), but does not require treating the mean flow deformation differently
than any of the other unresolved fast variables.

A classic illustrative model of the nonlinear stability mechanism is the
three-dimensional “mean-field” system proposed by Noack et al. (2003). The
model consists of a pair of unstable slow variables $a_{1}$ and $a_{2}$ along
with the fast variable $a_{\Delta}$:

$\displaystyle\dot{a}_{1}$ $\displaystyle=\lambda a_{1}-\omega
a_{2}-a_{1}a_{\Delta}$ (65a) $\displaystyle\dot{a}_{2}$ $\displaystyle=\omega
a_{1}-\lambda a_{2}-a_{2}a_{\Delta}$ (65b) $\displaystyle\dot{a}_{\Delta}$
$\displaystyle=-a_{\Delta}+\mu(a_{1}^{2}+a_{2}^{2}).$ (65c)

The slow variables represent amplitudes of a pair of global instability modes,
while $a_{\Delta}$ models a background deformation that tends to reduce the
linear growth rate. This system is a simplified model with similar behavior to
the cylinder wake discussed in Sec. 6.2.

In many cases the growth rate $\lambda$ is relatively small, so that the
standard approach is to use the adiabatic elimination approximation
$\dot{a}_{\Delta}=0$. The fast variable is then given by an algebraic
relationship that pins the system to the parabolic invariant manifold
$a_{\Delta}\approx\mu(a_{1}^{2}+a_{2}^{2})$. Replacing $a_{\Delta}$ in Eq.
(65) results in a cubic Stuart-Landau equation for the complex amplitude
$A(t)=a_{1}+ia_{2}$:

$\dot{A}=(\lambda+i\omega)A-\mu A|A|^{2}.$ (66)

We recapitulate these results in order to show that this mechanism is also
resolved by the multiscale closure (49b). Here the fast dynamics are already
linear in $a_{\Delta}$ with damping constant $\nu=1$. The fast-slow linear
blocks $\mathsfbi{L}^{12}$ and $\mathsfbi{L}^{21}$ are also zero, so the only
modification due to the closure model is the cubic term $\hat{\mathsfbi{C}}$.
The only nonzero entries in the quadratic tensor are
$Q_{11\Delta}=Q_{22\Delta}=-1$ and $Q_{\Delta 11}=Q_{\Delta 22}=\mu$, leading
to cubic terms
$\hat{C}_{1111}=\hat{C}_{1122}=\hat{C}_{2211}=\hat{C}_{2222}=-\mu$. Replacing
the fast-slow nonlinear interactions in Eq. (65) with the cubic closure, the
three-dimensional model reduces to

$\displaystyle\dot{a}_{1}$ $\displaystyle=\lambda a_{1}-\omega a_{2}-\mu
a_{1}(a_{1}^{2}+a_{2}^{2})$ (67a) $\displaystyle\dot{a}_{2}$
$\displaystyle=\omega a_{1}-\lambda a_{2}-\mu a_{2}(a_{1}^{2}+a_{2}^{2}),$
(67b)

identical to the real and imaginary parts of the invariant manifold reduction
(66).

####  6.1.2 Energy transfer in the Burgers model

Along with mean flow interactions, one of the primary challenges in applying
model reduction to the Navier-Stokes equations is the under-resolution of the
fine scales of the flow, which typically results in underestimating the energy
dissipation. In the classic spectral dynamics picture of isotropic turbulence,
although the bulk of the energy is contained in the larger scales of the flow,
the viscous energy dissipation scales with the square of the wavenumber, so
that even very small scales play a key role in the energy balance of the flow.
The problem of truncating the cascade of energy from large to small scales
thus applies equally well to “model reduction” in the sense of discretizing
the governing equations on a mesh that cannot resolve the finest scales of the
flow.

Since POD modes are ordered hierarchically by average energy content, this
problem often manifests in reduced-order modeling when attempting to severely
truncate a POD representation; the higher-order modes are often much more
dissipative, even though they may not be of practical interest (Noack et al.,
2011). Further complicating this is the issue that advection-dominated flows
may require many more modes for a good kinematic approximation than for a good
dynamic approximation, and these additional degrees of freedom can lead to
decoherence and instability in the reduced-order model (Callaham et al.,
2022). The approach of partitioning the system and averaging over the “fast”
variables can address both of these issues by retaining only the most
energetic degrees of freedom while accounting for the average effect of
unresolved dissipation at fine scales.

To illustrate the role of nonlinear advection, Burgers (1948) and Hopf (1948)
introduced simple one-dimensional toy models of turbulence, now generally
studied as Burgers‘ equation. Here we apply the multiscale closure approach to
the Fourier-Galerkin representation of a forced, viscous Burgers‘ equation on
a periodic domain (Noack et al., 2008):

$\partialderivative{u}{t}+u\partialderivative{u}{x}=\eta\partialderivative[2]{u}{x}+\zeta
g(x,t),\qquad x\in(0,2\pi),$ (68)

where $\eta$ and $\zeta$ are the viscous and forcing amplitudes, respectively,
and the forcing $g(x,t)$ is constructed to act on only the leading
wavenumbers. Majda et al. (2006) used a similar model to explore a similar
stochastic model reduction procedure on large deterministic systems without
asymptotic scale separation.

Since the domain is periodic, we forego the empirical POD basis and expand the
solution in a Fourier basis

$u(x,t)=\sum_{k=-\infty}^{\infty}a_{k}(t)e^{ikx}.$ (69)

The standard Galerkin projection gives the coupled system of ODEs

$\dot{a}_{k}=-\eta k^{2}a_{k}-i\ell
a_{\ell}a_{k-\ell}+\zeta\delta_{1|k|}a_{k},$ (70)

where the forcing reduces to the Kronecker delta function $\delta_{1|k|}$ by
construction. The quadratic nonlinearity would typically result in a double
sum over wavenumbers, but the orthogonality property of the Fourier basis
restricts the interaction to the “triads” $(k,\ell,k-\ell)$. As for the full
Navier-Stokes equations, this nonlinearity satisfies the conservative property
discussed in Sec. 2.3. The system is forced at large scales $|k|=1$ and is
dissipative at smaller scales $|k|\gg 1$, and the triadic interactions mediate
energy transfers between these scales. Although the closure method described
in Sec. 4 was derived for real-valued systems, Eq. (70) could be transformed
to be purely real by separating the real and imaginary parts of the complex-
conjugate pairs, as was done for instance by Noack et al. (2008). In this case
the system is simple enough that it does not matter, so we work with the
complex representation for compact notation.

We can apply the results of Sec. 4 directly by rewriting the quadratic term
with the sparse tensor $\mathsfbi{Q}$ defined by

$Q_{k\ell m}=\begin{cases}-i\ell,\qquad&\ell+m=k\\\
0,\qquad&\mathrm{otherwise}.\end{cases}$ (71)

The closure modeling is simplified by the diagonal structure of the linear
term, so that the only modification to Eq. (70) is the cubic term
$\hat{\mathsfbi{C}}$ given by Eq. (50d). Based on the sparse structure of
$\mathsfbi{Q}$, the sum over fast scales disappears so that

$\hat{C}_{k\ell
mn}=\frac{1}{\nu_{k-\ell}}\begin{cases}-km,\qquad&\ell+m+n=k\\\
0,\qquad&\mathrm{otherwise}.\end{cases}$ (72)

Note that these entries are also nonzero only for $k$, $\ell$, $m$, and $n$ in
the slow scales and $k-\ell$ in the fast scales. For example, if we model only
the forced scales $|k|=1$ as slow variables, the only nonvanishing
contributions to $\hat{\mathsfbi{C}}$ are from the wavenumber sets
$(k,\ell,m,n)=(1,-1,1,1)$ and $(-1,1,-1,-1)$. Since the slow quadratic tensor
$\hat{\mathsfbi{Q}}$ is also zero for $|k|=1$, the closed ODE system
simplifies to the one-dimensional Stuart-Landau equation

$\dot{a}_{1}=(\zeta-\eta)a_{1}-\nu_{2}^{-1}a_{1}|a_{1}|^{2}.$ (73)

The cubic term in Eq. (73) captures the leading-order influence of the energy
transfers to smaller scales. This effect is naturally stabilizing; as in the
physical mechanism, if the slow variables are excited, the amplitude term
$|a_{1}|^{2}$ representing the unresolved scales increases and the cubic
damping term removes energy from the resolved scales.

In this case the multiscale closure directly accounts for only the next set of
wavenumbers $|k|=2$ because the purely triadic structure of the quadratic
interactions does not admit contributions from any smaller scales. The damping
constant $\nu_{2}$ must therefore be selected so that the explicitly modeled
transfers from $|k|=1$ to $|k|=2$ are representative of the entire cascade.
This can be accomplished with the energy balance method described in Sec. 4,
for example.

While traditional closure models for Galerkin systems approximate the
dissipation at higher orders with a linear linear eddy-viscosity-type closure
(Aubry et al., 1988; Rempfer & Fasel, 1994), the mechanism of energy transfer
to smaller scales is fundamentally nonlinear in nature. As clearly illustrated
by this model problem, the multiscale closure model preserves the nonlinearity
of the energy transfer in its approximation of the effects of the unresolved
small scales. This is consistent with the findings of Östh et al. (2014), that
a nonlinear term modeling subscale turbulence can outperform the linear eddy
viscosity closure in POD-Galerkin models.

###  6.2 Vortex shedding in the cylinder wake

\begin{overpic}[width=346.89731pt]{figs/cyl_results.pdf} \put(2.0,52.5){a)}
\put(73.0,52.5){b)} \put(2.0,25.0){c)} \end{overpic} Figure 6:  Reduced-order
models of the cylinder wake. The standard 9-mode Galerkin model accurately
estimates the stable limit cycle (a), but the transient dynamics deviate from
the slow manifold, leading to an energy overshoot (b). Both the two-
dimensional invariant manifold reduction (Noack et al., 2003) and MMR closure
models prevent this by eliminating the rapidly-equilibrating variable
associated with mean flow deformation, but MMR more accurately estimates the
limit cycle amplitude. The full flow field can be reconstructed by
approximating the harmonics with polynomial regression (c).

In this section we revisit the canonical two-dimensional flow past a circular
cylinder to show that the MMR closure can naturally resolve the stabilizing
effects, both of mean flow deformation and of energy transfer to higher-order
modes, without assuming proximity to a bifurcation.

In keeping with previous work demonstrating that the cylinder wake can be
modeled accurately with as few as two degrees of freedom, we reduce the 9-mode
POD-Galerkin system to a 2-mode generalized Stuart-Landau model (49b) with the
closure model (50) and damping derived from the energy-balance relation (51).
In this case, the elimination of the higher-order modes does not require
sacrificing kinematic resolution, since the unresolved coefficients can be
approximated as sparse polynomial functions of the resolved variables (Loiseau
et al., 2018a; Callaham et al., 2022).

Figure 6a compares the energy of the first mode pair,
$E(t)=a_{1}^{2}(t)+a_{2}^{2}(t)$, predicted by the POD-Galerkin model and the
2-mode closure model to DNS of the transient flow initialized from the
unstable steady state and perturbed by the least-stable eigenmode. Also
included for comparison is the invariant manifold model introduced by Noack et
al. (2003) and described in Sec. 4.4. Although the Galerkin model is stable
and accurately predicts the amplitude of the vortex shedding, it exhibits a
transient energy overshoot before settling onto the limit cycle. Figure 6b
shows the reason for this; whereas the DNS is restricted to an approximately
parabolic invariant manifold, the finite relaxation time in the shift mode
dynamics results in the energy of the leading modes continuing its exponential
growth for longer before the energy can be absorbed into the mean flow
deformation.

In contrast, the cubic term in the MMR and invariant manifold models
approximates the relaxation of the shift mode dynamics as instantaneous,
pinning its amplitude to that of the primary mode pair. While the structure
and interpretation of the two models is the same, there are quantitative
differences in their accuracy. In particular, the manifold model slightly
overestimates the energy of the limit cycle, with mode amplitudes larger than
the DNS by $\sim 13\%$, while the multiscale closure accurately predicts the
limit cycle amplitude. This discrepancy can most likely be attributed to the
energy balance procedure used to estimate the subscale damping $\bm{\nu}$ in
the MMR model, which is not used in the invariant manifold model.

Finally, the full flow field can also be reconstructed from the two degree-of-
freedom multiscale model by means of the sparse polynomial approximation to
the nonlinear correlations in the higher-order modes (Fig. 6c). The flow field
predicted by the model gives a good approximation of the DNS solution,
although over long times the phase of the vortex shedding tends to drift from
that of the high-dimensional solution. By systematically averaging over the
fast variables, the multiscale closure method is thus able to reduce the
dimensionality of the model while also improving its fidelity.

###  6.3 Chaotic lid-driven cavity flow

The Stuart-Landau equations derived for the analytic problems in Sec 6.1 and
the cylinder wake in Sec. 6.2 are examples of relatively simple nonlinear
oscillators with either a stable fixed point or periodic limit cycle as an
equilibrium. Here we show that the MMR closure approach can also be applied to
systems with more complex dynamics, as illustrated on the chaotic lid-driven
cavity flow introduced in Sec. 5.2.

Figure 7: Reduced-order models of the lid-driven cavity flow. The Galerkin
system eventually equilibrates at much higher energy levels, while the MMR
closure model produces stable predictions that closely match the power
spectrum of the chaotic flow (left) The reconstructed fields from the MMR
model are also more physically consistent with the DNS than the Galerkin model
(right).

First, we construct a standard POD-Galerkin model with $r=64$ modes,
sufficient to capture $\sim 99.8\%$ of the fluctuation kinetic energy. The
Galerkin system is chaotic and energetically stable. Despite its precise
kinematic resolution, the Galerkin model is inaccurate, overestimating the
kinetic energy by roughly one order of magnitude, as shown by the blue traces
in Fig. 7.

Reducing the model to $r_{0}=21$ degrees of freedom via a multiscale closure
stabilizes it with an equilibrium energy that approximately matches the power
spectrum of the DNS (Fig. 7, orange). The MMR model still retains enough modes
for a $\sim 98\%$ accurate POD reconstruction. The predicted flow fields are
coherent and qualitatively consistent with the DNS fields, even though they do
not match in every detail due to the chaotic nature of the flow.

Varying the dimension $r_{0}$ of the approximate slow variables results in
models that are energetically stable, but not necessarily chaotic. Instead,
many such models eventually settle onto post-transient dynamics that are
periodic or quasi-periodic. The persistence of the dominant frequencies
suggests that even for models with relatively many degrees of freedom
($r_{0}>2$), the generalized Stuart-Landau equations produced by the MMR
closure scheme can be viewed as coupled nonlinear oscillators, especially when
the spatial components of the POD modes suggest that they may be grouped in
pairs with similar temporal frequency content.

###  6.4 Convective instability and vortex pairing in a mixing layer

Advection-dominated flows have historically posed a challenge for reduced-
order models, primarily because the global mode ansatz of POD is a poor
representation of traveling waves. As a final numerical example, we apply the
MMR closure to model the mixing layer introduced in Sec. 5.3, a canonical
example of convective instability in free shear flows.

In one sense, the dynamics of the mixing layer are relatively simple. The flow
close to the inlet is dominated by the shear layer instability and subsequent
vortex roll-up. These vortices are advected downstream at approximately the
midline flow velocity, until the secondary vortex pairing begins at the
subharmonic of the primary instability frequency. However, since this behavior
unfolds as the disturbances are carried downstream, a principal challenge for
low-dimensional models based on global POD modes is to preserve phase
coherence between the harmonics and subharmonics, so that the reconstructed
flow fields consist of approximately discrete vortices undergoing the vortex
pairing.

As discussed in Sec. 5.3, both the POD mode basis and the difficulty of
modeling this flow depend on the streamwise extent of the domain. Since the
subharmonic vortex pairing is energetically dominant downstream, the leading
POD modes computed from a longer domain will be related to vortex pairing,
even though the primary shear layer instability drives the dynamics from a
physically causal perspective. To illustrate this effect, we consider two
domain lengths: one that is dominated by the primary instability and vortex
roll-up, and the full modeling domain with vortex pairing.

\begin{overpic}[width=433.62pt]{figs/mix_L2_results.pdf} \end{overpic} Figure
8:  Models of the mixing layer on a short domain. The phase portraits of the
mode pairs (left) show phase-locking between the fundamental mode pair
$(a_{1},a_{2})$, the subharmonic vortex pairing $(a_{3},a_{4})$, and the first
harmonic $(a_{5},a_{6})$. These phase relationships are preserved by the MMR
closure, resulting in physically consistent flow field estimates (right).

In this case the POD-Galerkin models do not necessarily improve with an
increasing number of modes. For both domains, the best dynamical approximation
is found for models constructed from $r=24$ modes. On the short domain
($L_{x}=150$), this model is energetically stable but eventually approaches an
unphysical fixed point (Fig. 8). On the longer domain ($L_{x}=250$), the model
initially follows the DNS trajectory before eventually becoming unstable (Fig.
9). In either case, increasing the dimension of the POD-Galerkin model
resulted in worse models. We then reduce the $24$-dimensional Galerkin systems
to $r_{0}=16$-dimensional generalized Stuart-Landau models with the MMR
closure.

Figure 8 shows the comparison of the closed model to the DNS and Galerkin
systems for the short domain. The MMR model closely matches the phase
portraits of the DNS, indicating that it preserves the phase relationships
between the various modes, including the early subharmonic vortex pairing
$(a_{3},a_{4})$ and higher harmonics, e.g. $(a_{5},a_{6})$. This remains true
even when integrated to very long times (the final integration time is
$t=2820$ in this case, as with the DNS described in Sec. 5.3). As a result,
the approximated fields are highly coherent, with the same pattern of vortex
roll-up, advection, and pairing as seen in the DNS (Fig. 8, right).

\begin{overpic}[width=390.25534pt]{figs/mix_L3_results.pdf} \end{overpic}
Figure 9:  Mixing layer models on the long domain. While standard POD-Galerkin
models are energetically unstable until at least $r=64$, the multiscale
closure stabilizes the model with only 16 variables (left). The model also
remains coherent on long integration times (see also Fig. 10), producing flow
field predictions that are consistent with the large-scale structure of the
DNS (right).  \begin{overpic}[width=433.62pt]{figs/mix_L3_phase.pdf}
\end{overpic} Figure 10:  Phase portraits of the mixing layer model on the
long domain. The MMR closure (orange) closely matches the DNS trajectories
(black) for the most energetic modes, even though the structure is more
complicated than the Lissajous-type harmonics on the short domain (Fig. 8).
This figure does not show a Galerkin model for comparison because they are
energetically unstable (Fig. 9).

On the longer domain, the multiscale closure also stabilizes the Galerkin
model. Figure 9 compares the true energy content of the POD coefficients to
that predicted by the model, both in the time and frequency domain. Although
the evolution of the energy matches at the dominant frequencies, the MMR model
does not capture the weakly chaotic dynamics exhibited by the DNS as a result
of slight irregularity in the vortex pairing. Instead, after a brief transient
the solution becomes approximately periodic with phase locking between the
modes (Fig. 10). As with the short domain, this synchronization ensures that
the predicted flow fields remain coherent, even when integrated for long times
(Fig. 9, right).

As might be anticipated based on the more complex dynamics on the longer
domain, the fields predicted by the MMR model do not match the DNS quite as
closely as on the shorter domain, particularly closer to the inlet. This is
likely because fine-scale resolution of the upstream region of the flow
requires higher harmonics of the primary instability that do not appear in the
POD modes of the longer domain until much higher-order mode number. Higher-
resolution predictions might therefore be constructed by applying nonlinear
correlations analysis and polynomial regression to selected higher-order
modes, as in the cylinder reconstruction in Sec. 6.2; see also Loiseau et al.
(2018a) and Callaham et al. (2022). Still, the dominant vortices are still
evident in the prediction of the low-dimensional model, even in the downstream
pairing region.

On both domains, the synchronization and long-time phase coherence between the
modes further reinforces the picture outlined in Sec. 6.3 of the closure model
as consisting of coupled nonlinear oscillators. Although the description of
this flow in terms of spatially fixed global modes is not an ideal
representation of the advection-driven dynamics, the proposed MMR closure
method is able to successfully stabilize the POD-Galerkin models and
accurately capture the phase relationships necessary for accurate predictions
of the full flow fields.

##  7 Discussion

In this work, we have developed a multiscale model reduction (MMR) approach to
improve linear-quadratic dynamical systems, derived from POD-Galerkin
projection of the Navier-Stokes equations, with the addition of systematically
computed cubic closure terms. These cubic closure terms are derived through an
adaptation and application of a multiscale stochastic averaging method that
originated in singular perturbation theories of Markov processes (Kurtz, 1973;
Papanicolaou, 1976; Majda et al., 2001; E, 2011; Pavliotis & Stuart, 2012).
Whereas the standard truncation of the Galerkin system disregards the
influence of unresolved variables, the proposed multiscale model reduction
method, accounts for their effect in an average sense through averaging via
the stochastic Koopman operator. In particular, this approach is able to model
nonlinear interactions between resolved and unresolved variables, capturing
key mechanisms such as mean flow deformation and the energy cascade. The
closed model includes cubic terms, taking the form of generalized Stuart-
Landau equations that often act as coupled nonlinear oscillators.

Since the derivation of the MMR closure method relies on an asymptotic
timescale separation between resolved and subscale variables that is generally
absent in fluid flows, we do not attempt to prove the general validity of this
approach when the flow is not close to a bifurcation. Instead, we have
demonstrated by extensive numerical examples that the method dramatically
improves the stability and accuracy of low-dimensional models of unsteady
fluid flow, compared to the standard POD-Galerkin models. After a comparison
to the well-known benchmark problem of vortex shedding in a cylinder wake, we
have shown that the MMR method can reproduce chaotic dynamics in lid-driven
cavity flow. Finally, we developed models of an incompressible mixing layer
and showed that the coupled-oscillator MMR closures not only stabilize the
models, but preserve phase relationships between the modes, which is critical
for physically consistent predictions of advection-driven flow.

These results raise several possibilities for interesting future work. For
example, the mixing layer is an “amplifier” flow characterized by convective
instability and transient energy growth. It is therefore highly sensitive to
the nature of the upstream flow; models developed for a specific inlet forcing
may not generalize. Further consideration of this issue could be the
foundation of a multiscale approach to input/output models in the resolvent
framework (McKeon & Sharma, 2010; Sharma & McKeon, 2013; Zare et al., 2017;
Padovan et al., 2020; Pickering et al., 2021). In this work we neglect the
time-varying boundary condition in the region upstream of the modeling domain,
but an input/output model could more properly treat this explicitly as
forcing, as is done for instance in control-oriented reduced-order models
(Barbagallo et al., 2009; Sipp et al., 2010). Recent work on trajectory-based
optimal oblique projections has also shown promise for capturing dynamically
important, but typically low-energy, modes in an input-output framework for
model reduction of open shear flows (Otto et al., 2022).

As discussed for the chaotic lid-driven cavity flow in Sec. 6.3, we have
observed that the MMR models often tend towards phased-locked periodic or
quasiperiodic solutions rather than chaos, consistent with the idea of a
system of weakly coupled oscillators. However, we expect that in the
development of models of chaotic or turbulent flows it may be helpful to re-
introduce some degree of stochasticity to better match the flow statistics.
Fortunately, this is straightforward in the multiscale modeling framework (E,
2011; Pavliotis & Stuart, 2012). For instance, the homogenization method
(Majda et al., 2001) begins from slightly different scaling assumptions and
results in a closed model in which the fast-scale stochastic forcing appears
as a new diffusion term. A careful investigation of different scale separation
assumptions and stochastic approximations could be valuable in extending this
framework to turbulent flows.

In this work we have only considered linear-quadratic dynamics derived from
POD-Galerkin projection. This has the advantage of being based on the
governing equations, but recent work has shown that data-driven model
discovery can be a powerful alternative (Brunton et al., 2016; Loiseau &
Brunton, 2018; Loiseau et al., 2018a, b; Qian et al., 2020; Rubini et al.,
2020; Deng et al., 2020; Callaham et al., 2022). Data-driven modeling is
especially useful when the governing equations are unknown or projection-based
modeling is not as straightforward as for incompressible fluid flow (Qian et
al., 2020; Loiseau, 2020; Guan et al., 2021). In a sense, data-driven modeling
might circumvent the need for closure modeling by fitting directly to the time
series, including leveraging the intuition that cubic terms will appear in the
effective dynamics (Loiseau & Brunton, 2018). However, the linear-quadratic
system has certain symmetries and conservation properties that can be enforced
in constrained regression to a quadratic model (Schlegel & Noack, 2015;
Loiseau & Brunton, 2018; Kaptanoglu et al., 2021a, b), while extending these
properties to cubic nonlinearity is not straightforward. An interesting future
direction could be exploring a two-step process, first identifying a quadratic
model that satisfies the constraints, and then applying the MMR closure to
reduce it to a lower-dimensional cubic model. Since the empirical models are
often sparse, this could potentially result in models with improved stability
and scaling compared to the projection-based approach. Alternatively, a more
flexible deep-learning model could be used to replace the proposed diagonal
drift approximation, for instance with a modified variational autoencoder
(Kingma & Welling, 2013) that parameterizes the distribution of subscale
variables conditioned on the slow variables.

It will also be interesting to extend this analysis to systems where the
leading-order cubic terms tend to _destabilize_ the system, such as the
subcritical bifurcation in plane Poiseuille flow (Stewartson & Stuart, 1971).
This is an important scenario in the transition to turbulence, particularly
for shear flows in which non-normal energy amplification can activate the
nonlinearity even when the flow is globally stable. Higher-order effective
nonlinearity could potentially be incorporated into the framework via state-
dependent diffusion or an alternative fast drift model, for instance.

Finally, the analytic elimination of small scales in the Fourier-Galerkin
model of Burgers’ equation raises the possibility of a subscale LES closure
model based on multiscale analysis. This might take the form of a multiple-
scale expansion in both space and time or of a filter designed for the
isotropic scales of turbulence, for instance. In either case, the Stuart-
Landau form of the MMR closure model suggests that a spatioteporal analysis
would be akin to a Ginzburg-Landau theory of variational multiscale turbulence
modeling.

Beyond the context of POD-Galerkin reduced-order modeling, the multiscale
model reduction framework represents a general theory of effective cubic
nonlinearity arising in amplitude equations for global modes in quadratically
nonlinear fluid flow. The system of generalized Stuart-Landau equations
produced by the closure model presents a picture of the flow as a set of
coupled nonlinear oscillators describing the evolution of mode pairs. Close to
the threshold of instability the MMR closure is consistent with a weakly
nonlinear analysis, but it does not rely on proximity to a bifurcation. The
multiscale modeling methodology is therefore an alternative to standard
asymptotic expansions that may have implications for a variety of theoretical
and numerical approaches to modeling the large, slow scales of unsteady fluid
flows.

## Acknowledgements

SLB acknowledges funding support from the Army Research Office (ARO
W911NF-19-1-0045) and National Science Foundation AI Institute in Dynamic
Systems (grant number 2112085). JLC acknowledges funding support from the
Department of Defense (DoD) through the National Defense Science & Engineering
Graduate (NDSEG) Fellowship Program.

## Declaration of interests

The authors report no conflict of interest.

##  Appendix A Stability analysis of the mixing layer inlet profile

Figure 11:  An incompressible mixing layer at $\mbox{{Re}}=500$. The
hyperbolic tangent base flow is convectively unstable, amplifying small
perturbations as they are carried downstream. We force the flow at the inlet
with eigenfunctions of inviscid flow equations linearized around the base flow
(left). The central part of the domain is used for further modeling to avoid
boundary effects, while the downstream extent $x\in(250,300)$ is strongly
damped to prevent numerical instability at the outlet.

This appendix describes the eigenfunction forcing applied at the inlet of the
mixing layer configuration described in Sec. 5.3. The flow acts an amplifier
(Huerre & Monkewitz, 1985, 1990; Chomaz, 2005), so the downstream flow is
highly sensitive to the inlet conditions; eigenfunction forcing accentuates
the natural dynamics of the flow.

Since the dominant Kelvin-Helmholtz instability is inviscid (Michalke, 1964),
we derive the forcing from a temporal stability analysis of the inviscid
equations linearized about the parallel hyperbolic tangent inlet profile
$\bm{u}(\bm{x})=(U(y),0)$. Defining a perturbation streamfunction
$\psi(\bm{x},t)$ with real wavenumber $\alpha$ and complex velocity
$c=c_{r}+ic_{i}$:

$\psi(\bm{x},t)=\mbox{Re}\left\\{\hat{\psi}(y)e^{i\alpha(x-ct)}\right\\},$
(74)

the evolution of the perturbation is given by the Rayleigh equation

$\left(U-c\right)\left(\derivative[2]{y}-\alpha^{2}\right)\hat{\psi}-U^{\prime\prime}\hat{\psi}=0.$
(75)

For each wavenumber $\alpha$, this can be written as a generalized eigenvalue
problem

$\mathbfcal{A}(\alpha)\hat{\psi}_{j}(y;\alpha)=c^{(j)}\mathbfcal{B}(\alpha)\hat{\psi}_{j}(y;\alpha)$
(76)

for the $j^{\mathrm{th}}$ velocity $c$ and perturbation streamfunction
$\hat{\psi}$, where we use the convention of sorting by decreasing
$\mbox{Im}\\{c^{(j)}\\}$.

The base profile is odd, so the real phase velocity is the midline value
$c_{r}=\bar{U}$ (Michalke, 1964). With $\bar{U}=1$, perturbations oscillate
with frequency $\omega=\alpha$ and growth rate $\sigma(\alpha)=\alpha
c_{i}(\alpha)$. For further details on linear stability analysis of shear
flows, see for example Huerre & Monkewitz (1985); Schmid & Henningson (2012).

We choose a spectral discretization of Eq. (75) using a Hermite collocation
scheme from a Python implementation of the dmsuite library (Weidemann & Reddy,
2000). We find the least stable wavenumber $\alpha^{*}$ by maximizing the
growth rate $\sigma$ (i.e. minimizing -$\alpha c_{i}^{(1)}$) using an
implementation of the BFGS algorithm available in scipy. Since this departs
from the standard shooting method, we validate it against Michalke (1964), who
uses the profile ${U(y)=0.5(1+\tanh y)}$ based on nondimensionalizing with the
shear layer thickness $\delta_{\omega}/2$ rather than vorticity thickness.
With that profile we find a least-stable wavenumber $\alpha^{*}=0.4449$ and
growth rate $\sigma=0.0949$ compared to the reported results of
$\alpha^{*}=0.4446$ and $\sigma=0.0949$. For the present base profile (64) we
find $\alpha^{*}=0.8912$ and $\sigma=0.2129$.

Once the wavenumber $\alpha^{*}$ and eigenfunction $\hat{\psi}_{1}$
corresponding to the maximum growth rate is identified, the inlet velocity
perturbation can be computed from

$\hat{u}_{1}(y;\alpha^{*})=\hat{\psi}_{1}^{\prime}(y;\alpha^{*}),\hskip
56.9055pt\hat{v}_{1}(y;\alpha^{*})=i\alpha^{*}\hat{\psi}_{1}(y;\alpha^{*}).$
(77)

The optimal perturbations are shown in Fig. 11. Following Colonius et al.
(1997), we also apply forcing at the first three subharmonics of the least-
stable mode computed from solving (75) at $\alpha^{*}/2$, $\alpha^{*}/4$, and
$\alpha^{*}/8$, rescaling each so that the maximum value of $\hat{u}_{1}$ is
$10^{-3}\Delta U$. We perturb the inlet profile by
$\mbox{Re}\\{\hat{\bm{u}}_{1}(y,\alpha_{n})e^{i\omega_{n}t}\\}$ for
$n=1,2,4,8$ and $\alpha_{n}=\omega_{n}=\alpha^{*}/n$ and use a simulation time
step of $\Delta t=0.00705=10^{-3}\times 2\pi/\alpha^{*}$, so that the sampling
rate of the simulation is commensurate with the forcing frequency.

##  Appendix B Variable truncation of the mixing layer model

One of the limitations of projection-based reduced-order modeling for systems
without an explicit scale separation is a lack of principled criteria for
selecting the truncation rank, or dimension of the projection subspace. A
common approach is to use _ad hoc_ criteria such as the number of modes
required to resolve some fraction of the fluctuation kinetic energy (e.g. 90%
or 99%) on average. To some extent this simply replaces the arbitrary choice
of rank truncation with an equally arbitrary choice of approximation accuracy.
This is especially the case when the dynamic approximation accuracy of the
reduced-order model does not monotonically improve with increasing kinematic
approximation accuracy (as can be seen in Fig. 9, for instance).

Unfortunately, this problem is compounded in the secondary reduction step of
the MMR method, since one must select both the rank $r$ of the original POD-
Galerkin system and the rank $r_{0}$ of the final MMR system. Our approach in
this work was to follow the heuristic of simulating the POD-Galerkin system at
various values of $r$ and choosing that which performs best or remains stable
for longest (e.g. $r=24$ for the mixing layer model in Fig. 9). The final MMR
models then typically behave similarly for a range of final ranks $r_{0}$; we
report results for those that best reproduce the dynamics of the fluid flow.

Figure 12:  Alternative truncations of the mixing layer model in Sec. 6.4.
Rather than the “optimal” value of $r=24$, the MMR models shown here are
constructed from Galerkin systems of rank $r=64$. While the POD-Galerkin
models are inaccurate and energetically unstable for any rank truncation,
while the MMR models are less sensitive to the specific values of rank
truncation.

However, in order to show that the method is relatively insensitive to the
specific choices of rank truncation, here we repeat the model reduction for
the mixing layer on the longest domain shown in Fig. 9. These models begin
with $r=64$ rather than the “optimal” POD-Galerkin rank of $r=24$. Although
less accurate than the model reported in Sec. 6.4, the resulting models are
stable for all values $r_{0}$ of the final rank and significantly outperform
any POD-Galerkin model.

Figure 13:  Truncation of the mixing layer in Sec. 6.4 based on the singular
value spectrum. In this case the initial truncation level $r=32$ is selected