# Neural Closure Models for Dynamical Systems

Abhinav Gupta<EMAIL_ADDRESS>Department of Mechanical Engineering,
Massachusetts Institute of Technology, Cambridge, MA 02139  Pierre F.J.
Lermusiaux<EMAIL_ADDRESS>(Corresponding author.) Department of Mechanical
Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139

###### Abstract

Complex dynamical systems are used for predictions in many domains. Because of
computational costs, models are truncated, coarsened, or aggregated. As the
neglected and unresolved terms become important, the utility of model
predictions diminishes. We develop a novel, versatile, and rigorous
methodology to learn non-Markovian closure parameterizations for known-
physics/low-fidelity models using data from high-fidelity simulations. The new
neural closure models augment low-fidelity models with neural delay
differential equations (nDDEs), motivated by the Mori-Zwanzig formulation and
the inherent delays in complex dynamical systems. We demonstrate that neural
closures efficiently account for truncated modes in reduced-order-models,
capture the effects of subgrid-scale processes in coarse models, and augment
the simplification of complex biological and physical-biogeochemical models.
We find that using non-Markovian over Markovian closures improves long-term
prediction accuracy and requires smaller networks. We derive adjoint equations
and network architectures needed to efficiently implement the new discrete and
distributed nDDEs, for any time-integration schemes and allowing nonuniformly-
spaced temporal training data. The performance of discrete over distributed
delays in closure models is explained using information theory, and we find an
optimal amount of past information for a specified architecture. Finally, we
analyze computational complexity and explain the limited additional cost due
to neural closure models.

Keywords— delay differential equations, reduced-order-model, turbulence
closure, ecosystem modeling, data assimilation, machine learning

##  1 Introduction

Most models only resolve spatio-temporal scales, processes, and field
variables to a certain level of accuracy because of the high computational
costs associated with high-fidelity simulations. Such truncation of scales,
processes, or variables often limit the reliability and usefulness of
simulations, especially for scientific, engineering, and societal applications
where longer-term model predictions are needed to guide decisions. There are
many ways to truncate high-fidelity models to low-fidelity models. Examples
abound and three main classes of truncations are: evolving the original
dynamical system in a reduced space, e.g., using reduced-order-models (ROMs)
[90, 43]); coarsening the model resolution to the scales of interest [50, 2];
and reducing the complexity or number of state variables, components, and
parameterizations [52, 14, 91]. In many applications, the neglected and
unresolved terms along with their interactions with the resolved ones can
become important over time, and a variety of modeling techniques have been
developed to represent the missing terms. Techniques that express these
missing terms as functions of modeled state variables and parameters are
referred to as closure models. A main challenge is that no one closure
approach to date is directly applicable to all four main classes of model
truncations. Another is that closure models are only well-defined for either
linear problems or simple cases. Finally, they can easily become ineffective
in the face of nonlinearities.

Due to the explosion of use of a variety of machine learning methods for
solving or simulating dynamical systems, a number of data-driven methods have
been proposed for the closure problem. Most of them attempt to learn a neural
network as the instantaneous map between the low-fidelity solution and the
residual of the high- and low-fidelity solution, or their residual dynamics
[67, 88, 89, 80, 65]. They often use recurrent networks such as long-short
term memory networks (LSTMs), gated recurrent units (GRUs) etc., with
justification based on the Mori-Zwanzig formulation [84, 17, 32] and embedding
theorems by Whitney [92] and Takens [85]. These approaches do not however take
into account accumulation of numerical time-stepping error in the presence of
neural-networks during training. and uniformly-spaced high-fidelity data to be
able to compute the time derivative of the state with high level of accuracy.
Such requirement on the training data can be a luxury in a lot of scenarios.
The requirement of very frequent snapshot data of the system is also true for
methods which achieve model discovery using sparse-regression and provide
interpretable learned models [11, 41, 65]. All of the above issues are
addressed by using neural ordinary differential equations (nODEs; [16]) and
some researchers recently used nODEs for closure modeling. Some directly learn
the ODE system from high-fidelity simulation data without using the available
low-fidelity models [57], which could lead to the requirement of bigger neural
networks. Others combine nODEs with model discovery using sparse-regression
[94] or only learn the values of parameters in existing closure models [68].
Nearly all existing studies primarily only attempt to address the closure for
ROMs. Finally, the existing machine learned closure models are not yet used
for long-term predictions, i.e. forecasting significantly outside of the time-
period to which the training data belonged to.

In the present study, we propose a new neural delay differential equations
(nDDEs) based framework to learn closure parameterizations for low-fidelity
models using data from high-fidelity simulations and to increase the long-term
predictive capabilities of these models. Instead of using ODEs, we learn non-
Markovian closure models using DDEs. We base the theoretical justification for
using DDEs on the Mori-Zwanzig formulation [84, 17, 32] and the presence of
inherent delays in many dynamical systems [64], especially biological systems
[31, 86]. We refer to the new modeling approach as neural closure models. We
demonstrate that our methodology drastically improves the predictive
capability of low-fidelity models for the main classes of model truncations.
Specifically, our neural closure models efficiently account for truncated
modes in ROMs, capture the effects of subgrid-scale processes in coarse
models, and augment the simplification of complex mathematical models. We also
provide adjoint equation derivations and network architectures needed to
efficiently implement nDDEs, for both discrete and distributed delays. In the
case of distributed delays, we propose a novel architecture consisting of two
coupled neural networks, which eliminates the need for using recurrent
architectures for incorporating memory. We find that our nDDE closures
substantially improve nODE closures and outperform classic dynamic closures
such as the Smagorinsky subgrid-scale model. We explain the better performance
of nDDE closures based on information theory and the amount of past-
information being included. Our first two classes of simulation experiments
utilize the advecting shock problem governed by the Burger’s partial
differential equation (PDE), with low-fidelity models derived either by
proper-orthogonal-decomposition Galerkin projection [5] or by reducing the
spatial grid resolution. Our third class of experiments considers marine
biological models of varying complexities [62, 26, 9] and then their physical-
biogeochemical extensions, with low-fidelity models obtained by aggregation of
components and other simplifications of processes and parameterizations.
Finally, we analyze computational complexity and explain the limited
additional computational cost due to the presence of neural closure models.

##  2 Closure Problems

The need for closure modeling in dynamical systems arises for a variety of
reasons. They often involve computational costs considerations, but also
include the lack of data to resolve complex real processes, the limited
understanding of fundamental dynamics, and the inherent nonlinear growth of
uncertainties due to model errors and predictability limits [e.g. 76, 47, 48,
77]. In this section, we examine three main classes of low-fidelity models
that can require closure modeling.

###  2.1 Reduced Order Modeling

Let us consider a nonlinear dynamical system with state variable
$u\in\mathbb{R}^{N}$ and the full-order-model (FOM) dynamics governed by,

$\frac{du(t)}{dt}=Lu(t)+h(u(t)),\quad\text{with}\quad u(0)=u_{0}\,,$ (1)

where $L\in\mathbb{R}^{N\times N}$ is the linear, and
$h(\cdot):\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ the nonlinear, part of the
system. We are mainly interested in dynamical systems whose solution could be
effectively approximated on a manifold of lower dimension,
$\mathcal{V}\in\mathbb{R}^{m}\subset\mathbb{R}^{N}$ [e.g. 27]. Ideally, the
dimension of this manifold is much smaller than that of the system, i.e.
$m<<N$. For the classic Galerkin-based reduced-order modeling, a linear
decomposition of the form,

$u(t)\approx\bar{u}+Va$ (2)

is used, where $\bar{u}\in\mathbb{R}^{N}$ is a reference value, the columns of
$V=[v_{1},...,v_{m}]\in\mathbb{R}^{N\times m}$ a basis of the $m$-dimensional
subspace $\mathcal{V}$, and $a\in\mathbb{R}^{m}$ the vector of coefficients
corresponding to the reduced basis. A popular choice for this basis is the
proper-orthogonal-decomposition (POD) that defines the subspace such that the
manifold $\mathcal{V}$ preserves the variance of the system as much as
possible when projected on $V$ for a given $m$. The reference value
($\bar{u}$) is then commonly chosen as the mean of the system state, in order
to prevent the first reduced coefficient from containing the majority of the
energy of the system and to help stabilize the reduced system [35].

Now, substituting Eq. 2 into Eq. 1, and projecting the result on the
orthonormal modes $V$, we obtain the following set of ordinary differential
equations for the coefficients $a$,

$\frac{da}{dt}=V^{T}LVa+V^{T}h(\bar{u}+Va)+V^{T}L\bar{u}\,,\quad\text{with}\quad
a(0)=V^{T}(u_{0}-\bar{u})\,.$ (3)

This $m$ dimensional system, with $m<<N$, is computationally much cheaper than
the original FOM Eq. 2. This method of dimensionality reduction is commonly
referred to as the POD Galerkin Projection (POD-GP) method. It can suffer from
a number of issues. First, the truncated modes can play an important role in
the dynamical behaviour of the system, and neglecting them can thus lead to a
very different forecast [43]. Second, the error in the reduced state may be
simply too large for truncation, i.e. the POD reduction is not efficient.
Third, if steady POD are employed, they may quickly become irrelevant for the
evolving system state [81, 27, 28]. To address these issues, several methods
try to represent the effect of the truncated modes. The most common approaches
introduce a nonlinear parameterization of the coefficients [e.g. 55] in Eq. 3,
however, they are not yet generally applicable to all classes of closures.

The geometric interpretation of the goal of closure modeling for ROMs is
sketched in Fig. 1. The FOM solution of our dynamical system lies outside the
lower dimension manifold, $\mathcal{V}$. A ROM approximate solution, denoted
by $u^{ROM}$, starts with the projection of the full-order initial condition
onto the manifold, $VV^{T}u(0)$, but quickly diverges from the actual
projection of the full-order solution onto the manifold ($VV^{T}u$), often
leading to a significant source of error. A closure model in this case
basically attempts to keep the updated “closed” solution, $u^{ROM+C}$, as
close as possible to the actual projection of the FOM solution ($VV^{T}u$)
which can be seen as the truth.

Figure 1: Geometric interpretation of the closure for reduced-order-models
(ROMs). $u$ ( ): Solution to the full-order-model (FOM); $VV^{T}u$ ( ):
Projection of $u$ on the subspace $V$; $u^{ROM}$ ( ): Solution to the proper-
orthogonal-decomposition Galerkin-projection (POD-GP) ROM; and, $u^{ROM+C}$ (
): Solution to POD-GP ROM with closure. Adapted from [89].

###  2.2 Subgrid-Scale Processes

A key decision while setting up any numerical simulation is the selection of
spatio-temporal resolution, which is in general limited by the computing power
available. Using a coarse resolution (low-fidelity) model may however lead to
a number of undesired artifacts, such as missing critical scales and processes
for longer-term predictions or numerical diffusion that causes unintended or
unacceptable results [95, 44]. These artifacts become especially important in
the case of ocean models. For example, present-day global observing systems
and global model solutions only resolve open-ocean mesoscale processes
($\mathcal{O}(10-100km)$), but the submesocale (subgrid-scale) processes do
have global consequences, in relation to the mechanisms of energy dissipation
in the general circulation, vertical flux of material concentrations, and
intermediate-scale horizontal dispersal of materials [59, 60, 19]. The
neglected and unresolved scales along with their interactions with the
resolved ones are then at the core of closure parameterizations. Most present
oceanic models consist of a nonlinear system of PDEs, each of a nonlinear
advection type, supplemented by other possible diagnostic nonlinear equations,
and boundary conditions. There is however no unique way of defining such
parameterizations and multiple approaches such as non-dimensional analyses,
physical balance hypotheses, statistical correlation constraints, and other
empirical methods are commonly employed to develop closure models. Similar
statements can be made for atmospheric, Earth system, and climate models [82].
For all these applications, a general approach for subgrid-scale closures
would thus be most useful.

###  2.3 Simplification of Complex Dynamical Systems

Due to incomplete understanding and limited measurements, it is common when
modeling real dynamical systems in nature and engineering that the dynamics
cannot be accurately explained just by using conservation laws and fundamental
process equations. We refer to such systems as complex dynamical systems. The
number of candidate models and equations can then be almost as large as the
number of modelers. The resulting models also vary greatly in terms of their
complexity. More complex models can capture key processes and feedbacks.
Complexity is increased by adding more parameters and parameterizations to the
existing components (state variables) of the dynamical model, but at some
point, it quickly becomes inevitable to add and model new components to
capture the underlying real processes accurately, hence further increasing
model complexity [63, 26, 58]. This is common, for example, in marine
ecosystem models, where simpler models only resolve the broad biogeochemical
classes, while more complex models capture detailed sub-classes [91, 52].
Increasing the number of components however can come at great computational
cost, can increase the overall uncertainty, and can lead to loss of accuracy
or stability due to the nonlinearities. Also, the unknown parameters for
models with more components are calibrated from available data and the
optimization process and parameter estimation quickly become challenging with
the increase in complexity, due to the simultaneous explosion in the number of
unknown parameters [74]. Thus, instead of adding more unknown
parameterizations or increasing the number of components, one might use
adaptive models [46, 45] and in general the present neural closure models with
time delays to incorporate the effects of missing processes in low-complexity
models, enabling them to adapt and emulate the response from high-complexity
models.

##  3 Theory and Methodology

In this section, we develop the theory and methodology for learning data-
assisted closure models for dynamical systems. We first review the Mori-
Zwanzig (MZ) formulation [84, 17, 32] which derives the exact functional form
of the effects of truncated dynamics for common reduced models. Unfortunately,
apart for very simple linear dynamical systems, the use of this formulation is
challenging without making unjustified approximations and simplifications. We
then discuss the presence of delays in complex dynamical systems, and their
impact on modeling [31]. Motivated by both the MZ formulation and the presence
of delays, we finally derive the new nDDEs and neural closure models,
including adjoint equations and network architectures, for both discrete and
distributed delays.

###  3.1 Mori-Zwanzig Formulation and Delays in Complex Dynamical Systems

Without loss of generality, the full nonlinear dynamical system model is
written as,

$\frac{du_{k}(t)}{dt}=R_{k}(u(t),t)\,,\quad\text{with}\quad
u_{k}(0)=u_{0k}\,,\quad k\in\mathfrak{F}\,.$ (4)

The full state vector is
$u=(\\{u_{k}\\}),\;k\in\mathfrak{F}=\mathfrak{R}\cup\mathfrak{U}$, where
$\mathfrak{R}$ is the set corresponding to the resolved variables (e.g. coarse
field or reduced variables), and $\mathfrak{U}$ the set corresponding to the
unresolved variables (e.g. subgrid field or complement variables), which as a
union, $\mathfrak{F}$, form the set for full space of variables. We also
denote $u=\\{\hat{u},\tilde{u}\\}$ where
$\hat{u}=(\\{u_{k}\\}),\;k\in\mathfrak{R}$ and
$\tilde{u}=(\\{u_{k}\\}),\;k\in\mathfrak{U}$. Similarly,
$u_{0}=\\{\hat{u}_{0},\tilde{u}_{0}\\}$, with
$\hat{u}_{0}=(\\{u_{0k}\\}),\;\in\mathfrak{R}$ and
$\tilde{u}_{0}=(\\{u_{0k}\\}),\;k\in\mathfrak{U}$.

The Mori-Zwanzig (MZ) formulation allows rewriting the above nonlinear system
of ODEs as,

$\begin{split}\frac{\partial}{\partial
t}u_{k}(u_{0},t)=\underbrace{R_{k}(\hat{u}(u_{0},t))}_{\text{Markovian}}+\underbrace{F_{k}(u_{0},t)}_{\text{Noise}}+\underbrace{\int_{0}^{t}K_{k}(\hat{u}(u_{0},t-s),s)ds}_{\text{Memory}},\quad
k\in\mathfrak{R}\,,\end{split}$ (5)

where $R_{k}$ is as in the full model dynamics (Eq. 4). Importantly, the above
equation is an exact representation of Eq. 4 for the resolved components. A
derivation is provided in the Supplementary Information (Sec.
LABEL:sec:_MZ_formulation_in_supp_info). Eq. 5 provides useful guidance for
closure modeling. The first term in Eq. 5 is the Markovian term dependent only
on the values of the variables at the present time, while the closure consists
of two terms: the noise term and a memory term that is non-Markovian. We can
further simplify Eq. 5 by applying the $P$ projection and using the fact that
the noise term lives in the null space of $P$ for all times, which could be
easily proved. For ROMs with initial conditions devoid of any unresolved
dynamics, i.e. $\tilde{u}_{0}=0$ and thus $u_{0}=\hat{u}_{0}$, we then retain
the exact dynamics after the projection step, noticing in this case that
$Pu_{k}(u_{0},t)=u_{k}(\hat{u}_{0},t),\forall k\in\mathfrak{R}$,

$\begin{split}\frac{\partial}{\partial
t}u_{k}(\hat{u}_{0},t)=PR_{k}(\hat{u}(\hat{u}_{0},t))+P\int_{0}^{t}K_{k}(\hat{u}(\hat{u}_{0},t-s),s)ds,\quad
k\in\mathfrak{R}\,.\end{split}$ (6)

Hence, for such systems, the closure model would only consider the non-
Markovian memory term. The above derivation of the MZ formulation has been
adapted from [84, 89, 32].

The MZ formulation clearly shows that a non-Markovian closure term requiring
time-lagged state information is theoretically needed to model the unresolved
or missing dynamics. This theoretical basis directly applies to the first two
classes of low-fidelity dynamical systems (Sec. 2), ROMs and coarse resolution
models. For the third category, the simplification of complex dynamical
systems, we emphasize biological and chemical systems. Many are modeled using
ODEs, with one state variable per biological or chemical component. Such ODEs
implicitly assume that information between state variables is exchanged
instantaneously. In reality, however, there are often time-delays for several
reasons. First, changes in populations or reactions have non-negligible time-
scales [e.g. 40, 21, 64]. Such time-scales are introduced in more complex
models by modeling intermediate state variables. Hence, the time response of
lower-complexity models can be comparable to that of high-complexity models
only by explicitly introducing delays [40, 7, 31, 86]. Second, many reactive
systems are modeled assuming smooth concentration fields of state variables
governed by PDEs with fluid flow advection and/or mixing, leading to
advection-diffusion-reaction PDEs [36, 25]. In that case, simplified models
still require time-delays due to the neglected reactive or biogeochemical
dynamics but now also due to truncated modes and/or subgrid-scale processes of
numerical models. For all of these reasons, the need for memory based closure
terms is clearly justified to represent complex dynamical systems.

There are some results for data-assisted / data-driven closure modeling based
on the MZ formulation. Some schemes create a coupled system of stochastic
differential equation using appropriate hidden-variables for approximate
Markovization of the non-Markovian term [39, 8]. Others use a variational
approach to derive nonlinear parameterizations approximating the Markovian
term [15]. Schemes using machine learning to learn non-Markovian residual of
the high- and low-fidelity dynamics limit themselves to specific functional
forms for the residual term, simple Euler time-stepping scheme, and very
frequent and uniformly-spaced training data [65, 89, 88]. They also lack the
rigorous use of the theory for time-delay systems [73].

###  3.2 Neural Delay Differential Equations

The non-Markovian closure terms with time-lagged state information lead us to
delay differential equations (DDEs) [22]. DDEs have been widely used in many
fields such as biology [83, 53], pharmacokinetic-pharmacodynamics [38],
chemistry [78], economics [37], transportation [54], control theory [42],
climate dynamics [29, 6], etc. Next, we summarize the state of the art for
learning and solving differential equations using neural-networks (NNs) and
develop theory and schemes for neural DDEs including adjoint equations for
backpropagation.

The interpretation of residual networks as time integration schemes and flow
maps for dynamical systems has led to pioneering development of neural
ordinary differential equations (nODEs) [16]. A nODE parameterizes an ODE
using a neural-network and solves the initial value problem (IVP) given by,

$\frac{du(t)}{dt}=f_{NN}(u(t),t;\theta)\,,\quad
t\in(0,T]\,,\quad\text{with}\quad u(0)=u_{o}\,,$ (7)

where $f_{NN}$ is the prescribed neural-network and $\theta$ are the weights.
Starting from the initial conditions, the nODE (Eq. 7) is integrated forward
in time using any time-integration scheme, and then gradients are computed
based on a loss function using the adjoint sensitivity method. The gradient
computation boils down to solving a second ODE backwards in time. Using
standard backpropagation for Eq. 7 has however several issues: it would be
very memory expensive as one needs to store the state at every time step; its
computational cost would increase when using higher-order time-integration or
implicit schemes; and, it might become infeasible if the forward time-
integration code does not support automatic differentiation. The adjoint
method, however, provides a backpropagation for nODEs [16] that is memory
efficient and flexible as it treats the time-integration scheme as a ”black-
box”. In our case, we need to incorporate state-delays. Though extending the
nODE framework to incorporate DDEs comes under the ambit of universal
differential equations (UDEs) [71], deeper investigations are warranted.
First, the UDEs are presently implemented using the Julia library
DiffEqFlux.jl [70] which can perform automatic adjoint equation solves, but
other popular open-source languages such as Python and R, and ML-Frameworks
such as TensorFlow [1], PyTorch [66], etc., would require explicit derivation
and coding of the corresponding adjoint equations. Second, we need to study
two different types of DDEs, the discrete and distributed delays, which it
turns out require different architectures. Next, we thus develop the theory
and schemes for efficient implementation of neural delay differential
equations (nDDEs) in any programming language.

####  3.2.1 Discrete Delays

The most popular form of delay differential equations (DDEs) is,

$\begin{split}\frac{du(t)}{dt}=f(u(t),u(t-\tau_{1}),...,u(t-\tau_{K}),t),\quad
t\in(0,T]\,,\quad\text{with}\quad u(t)=h(t),\quad t\leq 0\,,\end{split}$ (8)

where $\tau_{1},...,\tau_{K}$ are $K$ number of discrete-delays (discrete
DDEs). Instead of a single initial value as in the case of ODEs, DDEs require
specification of a history function, $h(t)$. Due to the presence of a given
fixed number delays, we can parameterize the above system by replacing the
time-derivative function with potentially any type of NNs. For example, to use
fully-connected NNs we would concatenate all the delayed states vertically to
form the input vector, or concatenate them horizontally to form an input
matrix for a convolutional NN. However, recurrent NN (RNN) architectures, such
as simple-RNNs, LSTMs, GRUs, etc., are ideal and most efficient for our need
due to the time-series nature of the delayed states. We can assume that the
discrete delays are evenly spaced (this is not a hard requirement as we can
easily extend schemes to irregularly spaced discrete-delays using ODE-RNNs
[79], but for brevity we make this assumption) and use a RNN with weights
$\theta$. Hence, our new discrete-DDE system can be written as,

$\begin{split}\frac{du(t)}{dt}=f_{RNN}(u(t),u(t-\tau_{1}),...,u(t-\tau_{K}),t;\theta),\quad
t\in(0,T]\,,\;\;\text{with}\;\;u(t)=h(t),\quad t\leq 0\,,\end{split}$ (9)

where $f_{RNN}(\bullet;\theta)$ is the recurrent architecture. We refer to
this parameterization of discrete DDEs as discrete-nDDE. The graphical
representation of Eq. 9 in time-discretized form is depicted in Figure 2a. Let
data be available at $M$ times, $T_{1}<...<T_{M}\leq T$. We then optimize the
total loss function given by,
$\mathcal{L}=\int_{0}^{T}\sum_{i=1}^{M}l(u(t))\delta(t-T_{i})dt$, where
$l(\bullet)$ are scalar loss functions such as mean-squared-error (MSE), and
$\delta(t)$ is the Kronecker delta function. To perform this optimization with
any gradient descent algorithm, we need the gradient of the loss function
w.r.t. the weights of the RNN, $\theta$. Using the adjoint sensitivity method
[13] to compute the required gradients, we start by writing the Lagrangian for
the above system,

$\begin{split}L=&\mathcal{L}(u(t))+\int_{0}^{T}\lambda^{T}(t)\left(d_{t}{u(t)}-f_{RNN}(u(t),u(t-\tau_{1}),...,u(t-\tau_{K}),t;\theta)\right)dt\\\
&+\int_{-\tau_{K}}^{0}\mu^{T}(t)(u(t)-h(t))dt\,,\end{split}$ (10)

where $\lambda(t)$ and $\mu(t)$ are the Lagrangian variables. In order to find
the gradients of $L$ w.r.t. $\theta$, we first solve the following adjoint
equation (for brevity we denote,
$\frac{\partial}{\partial(\bullet)}\equiv\partial_{(\bullet)}$ and
$\frac{d}{d(\bullet)}\equiv d_{(\bullet)}$),

$\displaystyle
d_{t}\lambda^{T}(t)=\sum_{i=1}^{M}\partial_{u(t)}l(u(t))\delta(t-T_{i})-\lambda^{T}(t)\partial_{u(t)}f_{RNN}(u(t),u(t-\tau_{1}),...,u(t-\tau_{K}),t;\theta)$
$\displaystyle-\sum_{i=1}^{K}\lambda^{T}(t+\tau_{i})\partial_{u(t)}f_{RNN}\left(u(t+\tau_{i}),u(t-\tau_{1}+\tau_{i}),...,u(t-\tau_{K}+\tau_{i}),t+\tau_{i};\theta\right),\;t\in[0,T)$
$\displaystyle\lambda(t)=0,\qquad t\geq
T\;.\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad$
(11)

Details of the derivation of the above adjoint Eq. 11 are in the accompanying
Supplementary Information. Note that Eq. 11 needs to be solved backward in
time, and one would require access to $u(t),~{}0\leq t\leq T$. In the original
nODE work [16], Eq. 9 is solved backward in time and augmented with the
adjoint Eq. 11, so as to shrink the memory footprint by avoiding the need to
save $u$ at every time-step. Solving Eq. 9 backward can however lead to
catastrophic numerical instabilities as is well known in data assimilation
[93, 75]. Improvements have been proposed, such as the ANODE method [30], but
they are not applicable in case of DDEs. In our present implementation, in
order to access $u(t),~{}0\leq t\leq T$, while solving the adjoint equation,
we create and continuously update an interpolation function using the $u$
obtained at every time-step as we solve Eq. 9 forward in time. To be more
memory efficient, we can, for example, use the method of checkpointing [33],
or the interpolated reverse dynamic method (IRDM) [20]. After solving for
$\lambda$, we can compute the required gradients as,

$d_{\theta}L=-\int_{0}^{T}\lambda^{T}(t)\partial_{\theta}f_{RNN}(u(t),u(t-\tau_{1}),...,u(t-\tau_{K}),t;\theta)dt\;.$
(12)

Finally, using any gradient descent algorithm, we can find the optimal values
of the weights $\theta$.

####  3.2.2 Distributed Delays

In some applications, the delay is distributed over some past time-period
[72],

$\begin{split}\frac{du(t)}{dt}=f\left(u(t),\int_{t-\tau_{2}}^{t-\tau_{1}}g(u(\tau),\tau)d\tau,t\right),\quad
t\in(0,T]\,,\quad\text{with}\quad u(t)=h(t),\quad t\leq 0\,.\end{split}$ (13)

It should be noted that the discrete DDEs can be written as a special case of
distributed DDEs using dirac-delta functions. We can approximate the two
functions $f$ and $g$ using two different neural-networks, and re-write the
above Eqs. 13 as our new coupled discrete DDEs,

$\begin{split}\frac{du(t)}{dt}&=f_{NN}\left(u(t),y(t),t;\theta\right),\quad
t\in(0,T]\\\
\frac{dy(t)}{dt}&=g_{NN}(u(t-\tau_{1}),t-\tau_{1};\phi)-g_{NN}(u(t-\tau_{2}),t-\tau_{2};\phi),\quad
t\in(0,T]\\\ \text{with}\quad u(t)&=h(t),\quad\tau_{2}\leq t\leq
0\,,\quad\text{and}\quad
y(0)=\int_{-\tau_{2}}^{-\tau_{1}}g_{NN}(h(t),t;\phi)dt\,,\end{split}$ (14)

where $f_{NN}(\bullet;\theta)$ and $g_{NN}(\bullet;\phi)$ are the two NNs
parameterized by $\theta$ and $\phi$ respectively. We refer to this
parameterization of distributed DDEs as distributed-nDDE. The graphical
representation of the above system (Eqs. 14) in time-discretized form is
depicted in Figure 2b. Interestingly in the case of distributed-delays, we
obtain a novel architecture consisting of two coupled NNs, which enables us to
incorporate memory without the use of any recurrent networks such as RNN,
LSTMs, GRUs, etc. We can consider $f_{NN}$ as the main network, and $g_{NN}$
as the auxiliary network. Again, we define a scalar loss function given by
$\mathcal{L}=\int_{0}^{T}\sum_{i=1}^{M}l(u(t))\delta(t-T_{i})dt$ for the
available data at $M$ times, $T_{1}<...<T_{M}\leq T$. The Lagrangian for the
above system is,

$\begin{split}L=&\mathcal{L}(u(t))+\int_{0}^{T}\lambda^{T}(t)(d_{t}u(t)-f_{NN}(u(t),y(t),t;\theta))\,dt\\\
&+\int_{0}^{T}\mu^{T}(t)\left({d_{t}y(t)-g_{NN}(u(t-\tau_{1}),t-\tau_{1};\phi)+g_{NN}(u(t-\tau_{2}),t-\tau_{2};\phi)}\right)dt\\\
&+\int_{-\tau_{2}}^{0}\gamma^{T}(t)(u(t)-h(t))dt+\alpha^{T}\left(y(0)-\int_{-\tau_{2}}^{-\tau_{1}}g_{NN}(h(t),t;\phi)dt\right)\,,\end{split}$
(15)

where $\lambda(t),\mu(t),\gamma(t)$, and $\alpha$ are the Lagrangian
variables. In order to find the gradients of $L$ w.r.t. the parameters of the
two NNs, we first solve the following coupled adjoint equations backward in
time,

$\begin{split}d_{t}\lambda^{T}(t)=&\sum_{i=1}^{M}\partial_{u(t)}l(u(t))\delta(t-T_{i})-\lambda^{T}(t)\partial_{u(t)}f_{NN}(u(t),y(t),t;\theta)\\\
&-\mu^{T}(t+\tau_{1})\partial_{u(t)}g_{NN}(u(t),t;\phi)+\mu^{T}(t+\tau_{2})\partial_{u(t)}g_{NN}(u(t),t;\phi)\,,\;t\in[0,T)\\\
d_{t}\mu^{T}(t)=&-\lambda^{T}(t)\partial_{y(t)}f_{NN}(u(t),y(t),t;\theta)\,,\qquad
t\in[0,T)\\\ \lambda^{T}(t)=&~{}0\quad\text{and}\quad\mu^{T}(t)=0,\qquad t\geq
T\,.\end{split}$ (16)

Details of the derivation of the above adjoint Eq. 16 are in the Supplementary
Information. For accessing $u$ values while solving the adjoint equations, we
use the same approach as for our discrete-nDDE (Sec. 3.2.1). After solving for
$\lambda$ and $\mu$, we can compute the required gradients as,

$\begin{split}d_{\theta}L=&-\int_{0}^{T}\lambda^{T}(t)\partial_{\theta}f_{NN}(u(t),y(t),t;\theta)dt\,,\\\
d_{\phi}L=&-\int_{0}^{T}\mu^{T}(t)\left(\partial_{\phi}g_{NN}(u(t-\tau_{1}),t-\tau_{1};\phi)-\partial_{\phi}g_{NN}(u(t-\tau_{2}),t-\tau_{2};\phi)\right)dt\\\
&-\mu^{T}(0)\int_{-\tau_{2}}^{-\tau_{1}}\partial_{\phi}g_{NN}(h(t),t;\phi)dt\,.\end{split}$
(17)

Finally, using any gradient descent algorithm, we can optimize the neural-
networks $f_{NN}$ and $g_{NN}$, and find the optimal values of the weights
$\theta$ and $\phi$.

(a)

(b)

Figure 2: Graphical representation of the time discretized neural delay
differential equations (nDDEs). The blocks labeled RNN and DNN represent any
recurrent or deep neural-network architectures respectively. The block labeled
$\int$ symbolizes any time-integration scheme.

###  3.3 Neural Closure Models

Now that we have the framework for representing delay differential equations
using neural-networks, we can replace the non-Markovian memory term in Eq. 6
using nDDEs to obtain a hybrid closure model which could be trained using data
from high-fidelity simulations or real observations. The modified low-fidelity
dynamical system with the nDDE closures, which approximates the high-fidelity
model would be given by,

$\begin{split}\frac{\partial\hat{u}(t)}{\partial
t}&=\underbrace{PR(\hat{u}(t))}_{\text{Low-
Fidelity}}+\underbrace{f_{RNN}(\hat{u}(t),\hat{u}(t-\tau_{1}),...,\hat{u}(t-\tau_{K}),t;\theta)}_{\text{Neural
Closure}}\\\ \hat{u}(0)&=\hat{u}_{0},\qquad t\leq 0\end{split}$ (18)

using discrete-delays, or by,

$\begin{split}\frac{\partial\hat{u}(t)}{\partial
t}&=\underbrace{PR(\hat{u}(t))}_{\text{Low-
Fidelity}}+\underbrace{f_{NN}\left(\hat{u}(t),\int_{t-\tau_{2}}^{t-\tau_{1}}g_{NN}(\hat{u}(\tau),\tau;\phi)d\tau,t;\theta\right)}_{\text{Neural
Closure}}\\\ \hat{u}(0)&=\hat{u}_{0},\qquad t\leq 0\end{split}$ (19)

using distributed-delays. The initial conditions at $t=0$, $\hat{u}_{0}$, can
be used for $t<0$ as well, as an approximation. Apart from the neural-network
architectures, the amount of delay to be used also becomes a hyperparameter to
tune. These novel _neural closure models_ provides extreme flexibility in
designing the non-Markovian memory term in order to incorporate subject matter
expert insights. At the same time, we can also learn the unknown parts of the
Markovian low-fidelity model using nODEs if the need arises. Next, we will
compare the performance and advantages of using no-delays (nODEs), discrete-
delays (discrete-nDDEs), and distributed-delays (distributed-nDDEs) in closure
terms for various low-fidelity dynamical systems.

##  4 Application Results and Discussion

After presenting the main classes of low-fidelity dynamical models that
require closure (Sec. 2), we derived a novel, versatile, and rigorous
methodology for learning and modeling non-Markovian closure terms using nDDEs
(Sec. 3). The resulting neural closure models have their underpinning in the
Mori-Zwanzig formulation and the presence of inherent delays in models of
complex dynamical systems such as biogeochemical systems. Now, we evaluate the
performance and advantages of these new neural closure models over those of
neural ODEs (Markovian).

We run experiments encompassing each of the classes of low-fidelity models
(Sec. 2). For each experiment, we follow the same training protocol for nODEs
(no-delays) and the two nDDEs, discrete-nDDE (discrete-delays) and
distributed-nDDE (distributed-delays), closure models. The training data are
regularly sampled from high-fidelity simulations in all experiments, but this
is not a requirement. We use performance over the validation period (past the
period for which high-fidelity data snapshots are used for training) to fine
tune various training related hyperparameters. The final evaluation is based
on much longer-term future prediction performance, well past these periods. As
the field of scientific machine learning (SciML; [69]) is relatively new, the
metrics for performance evaluation vary greatly. On the one hand, many
learning studies randomly sample small time-sequences from a given period for
which high-fidelity data are available, and then split them into training,
validation, and test sets [57]. As the training, validation, and test sets
belong to the same time-domain, hence, the learned networks are only evaluated
for their interpolation performance and predicting the unseen data becomes
easy for them. On the other hand, for the few studies where the training and
test (prediction) periods do not overlap, the prediction period is often much
shorter than the training period [96]. In the present work, we consider a more
stringent evaluation. First, our validation period does not overlap the
training period. Second, our future prediction period is equal to or much
longer than the training and validation periods, and has no overlap with
either. Hence, we strictly measure the out-of-sample/generalization
performance of the learned network for its extrapolation capabilities into the
future. Of course, other evaluation metrics are possible and there is indeed a
need for standardization of evaluation procedures in the SciML community. In
the rest of the paper, for all the figure, table, and section references
prefixed with “SI-”, we direct the reader to the Supplementary Information.

###  4.1 Experiments 1: Advecting Shock - Reduced Order Model

For the first experiments, neural closure models learn the closure of proper-
orthogonal-decomposition Galerkin projection (POD-GP) based reduced order
model of the advecting shock problem. The full-order-model (FOM) for this
problem is given by the Burger’s equation,

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

where $\nu$ is the non-dimensional diffusion coefficient. The initial and
boundary conditions are

$\begin{split}u(x,0)&=\frac{x}{1+\sqrt{\frac{1}{t_{0}}}\exp\left(Re\frac{x^{2}}{4}\right)}\,,\quad
u(0,t)=0,\quad\text{and}\quad u(L,t)=0\,,\end{split}$ (21)

where $Re=1/\nu$ and $t_{0}=\exp(Re/8)$. Let the POD of the state variable
$u(x,t)$ be given by, $u(x,t)=\bar{u}(x)+\sum_{i=1}^{m}u_{i}(x)a_{i}(t)$, we
obtain the reduced-order equations as outlined in Section 2.1,

$\begin{split}\frac{da_{k}}{dt}=&-\left\langle\bar{u}\frac{\partial\bar{u}}{\partial
x},u_{k}\right\rangle-a_{i}\left\langle u_{i}\frac{\partial\bar{u}}{\partial
x},u_{k}\right\rangle-a_{j}\left\langle\bar{u}\frac{\partial u_{j}}{\partial
x},u_{k}\right\rangle-a_{i}a_{j}\left\langle u_{i}\frac{\partial
u_{j}}{\partial x},u_{k}\right\rangle\\\
&+\left\langle\nu\frac{\partial^{2}\bar{u}}{\partial
x^{2}},u_{k}\right\rangle+a_{i}\left\langle\nu\frac{\partial^{2}u_{i}}{\partial
x^{2}},u_{k}\right\rangle\,,\quad\text{with}\quad
a_{k}(0)=\langle(u(x,0)-\bar{u}(x)),u_{k}(x)\rangle\,.\end{split}$ (22)

We solve the FOM (Eqs. 20 and 21) for $Re=1000$, $L=1$, and maximum time
$T=4.0$. The singular value decomposition (SVD) of this solution form the POD
modes for the ROM. We only keep the first three modes which capture $60.8\%$
of energy, and evolve the corresponding coefficients using Eq. 22, thus
requiring a closure. The high-fidelity or true coefficients are obtained
solving the FOM (Eq. 20) with initial conditions without the contribution from
the unresolved modes, i.e. $u(x,0)=\bar{u}(x)+\sum_{i=1}^{3}u_{i}(x)a_{i}(0)$,
and projecting the obtained solution onto the first three modes. For
comparison, we also present the true coefficients in Fig. 3, which is what the
ROMs with neural closure are trying to match. For this true data generation,
we solve the FOM using an explicit Runge-Kutta (RK) time-integration of order
(4)5 (dopri5; [34]) with adaptive time-stepping (storing data at time-steps of
$\Delta t=0.01$) and grid spacing of $\Delta x=0.01$, using finite-difference
schemes (upwind for advection and central difference for diffusion).

Our three test periods for the advecting shock ROM (Eq. 22) with three modes
are as follows. For training our neural closure models, we only use the true
coefficient values up to time $t=2.0$. For validation (used only to tune
hyperparameters), we use true coefficient values from $t=2.0$ to $t=4.0$.
Finally for testing, we make a future prediction from $t=4.0$ to final time
$T=6.0$. We compare the three different closures: nODE (no-delays), discrete-
nDDE, and distributed-nDDE with architecture details presented in Table
LABEL:table:_Exp1_and_2_architecture. The architectures are not exactly the
same for the three cases, but they are set-up to be of comparable expressive
power. Mostly, we employ a bigger architecture for the no-delays case in order
to help it compensate the lack of past information. We also ensure that the
networks are neither under-parameterized nor over-parameterized. Along with
the classical hyperparameters such as batch size, number of iterations per
epoch, number of epochs, learning rate schedule, etc., we also have the delay
values ($\tau_{1},...,\tau_{K}$ for discrete-nDDE; and $\tau_{1},~{}\tau_{2}$
for distributed-nDDE) as additional hyperparameters to tune. We chose to use
six discrete delays for the discrete-nDDE in the present experiments. The
values of other hyperparameters are given in Sec. LABEL:SI:_hyperparameters.
For evaluation, at each epoch, we evolve the coefficient of the learned system
$(\\{a^{pred}(T_{i})=\\{a_{k}^{pred}(T_{i})\\}_{k=1}^{3}\\}_{i=1}^{M})$ using
the RK time-integration scheme mentioned earlier, and compare them with the
true coefficients $(\\{a^{true}(T_{i})\\}_{i=1}^{M})$ using the time-averaged
$L_{2}$ error,
$\mathcal{L}=\allowdisplaybreaks\frac{1}{M}\sum_{i=1}^{M}\left(\sqrt{\sum_{k=1}^{3}|a_{k}^{pred}(T_{i})-a_{k}^{true}(T_{i})|^{2}}\right)$,
which is also our loss function for training. The error for the time period
$t=0~{}\text{to}~{}2.0$ forms the training loss, the error for
$t=2.0~{}\text{to}~{}4.0$ the validation loss, and the error for
$t=4.0~{}\text{to}~{}6.0$ the prediction loss.

The performance of the three neural closure models after 200 epochs (the
stochastic gradient descent nearly converges, see Fig.
LABEL:fig:_Exp_1_train_val_loss_all) is evaluated by comparison with the true
coefficients and with the POD-GP coefficients spanning training, validation,
and future prediction periods. Results are shown in Figure 3. The details of
the architectures employed are in Table LABEL:table:_Exp1_and_2_architecture.
We find that using no-delays (nODE), discrete-delays (discrete-nDDE), and
distributed-delays (distributed-nDDE) perform equally well for the training
period, exactly matching the true coefficients. As soon as one enters the
validation period, all the neural closure models starts to slightly diverge,
with the nODE diverging the most by the end of prediction period. Importantly,
both nDDE closures maintain a great improvement over just using the POD-GP
model, and showcase a better performance than the nODE closure, even though
the latter had a deeper architecture with significantly more trainable
parameters. We also find that the performance of the distributed-nDDE closure
is a little better than that of the discrete-nDDE closure for the prediction
period. In a similar set of experiments, Maulik et. al., 2020 (section 3.1,
“Advecting Shock”, [57]) used nODE and LSTM to learn the time evolution of the
first three high-fidelity (true) coefficients without utilizing the known
physics/low-fidelity model. As a result, they required bigger architectures
and more time samples for training data than we do. This confirms our benefits
of learning only the unknown closure model. Due to the highly nonlinear nature
of neural networks, analytical stability analyses are not direct. Nonetheless,
we provide empirical stability results by reporting the evolution of the root-
mean-square-error (RMSE) (Fig. 3). We find that both discrete-nDDE and
distributed-nDDE closures, due to the existence of delays, may have a stronger
dissipative character and thus show much better stability at later times than
the POD-GP and the nODE closure.

Figure 3: Comparison of the true coefficients (solid) with the coefficients
from the POD-GP ROM (dashed-dot) and from the POD-GP ROMs augmented with the
three different learned neural closure models at the end of training (dashed).
For each neural closure, the training period is from $t=0~{}\text{to}~{}2.0$,
the validation period from $t=2.0~{}\text{to}~{}4.0$, and the future
prediction period from $t=4.0~{}\text{to}~{}6.0$. Top-left: neural ODEs with
no-delays (nODE); Top-right: neural DDEs with discrete-delays (Discrete-nDDE);
Bottom-left: neural DDEs with distributed-delays (Distributed-nDDE). Bottom-
right: Evolution of root-mean-squared-error
$(\text{RMSE}(t)=\sqrt{\frac{1}{3}\sum_{k=1}^{3}|a_{k}^{pred}(t)-a_{k}^{true}(t)|^{2}})$
of coefficients from the four different ROMs. These results correspond to the
architectures detailed in Table LABEL:table:_Exp1_and_2_architecture.

One might expect the distributed-delay (distributed-nDDE) to always perform
better than the discrete-nDDE closure because of the presence of the integral
of the state variable over a delay period instead of the state variable at
specific points in the past. The former thus seemingly contains more
information, but there is in fact no guarantee for this being true in all
cases. We can derive an intuition for this from information theory. According
to the data processing inequality [18], let $X$ and $Y$ be two random
variables, then,

$I(g(X);Y)\leq I(X;Y)\,,$ (23)

where $I$ is the mutual information and $g$ is any function which post-
processes $X$. Now, if $X$ is composed of $K$ random variables,
$X=\\{X_{1},...,X_{K}\\}$, and $g(X)=X_{1}+...+X_{K}$, then,

$I(X_{1}+...+X_{K};Y)\leq I(\\{X_{1},...,X_{K}\\};Y)\,.$ (24)

If we consider the effect of the integral of the state variable over the delay
period in the case of distributed-nDDE as a data processing step, this might
actually be decreasing the information content as compared to the discrete-
nDDE closure. We use “might”, even though Eq. 24 is a strong bound, because in
the present experiments we only use six delay values for the discrete-nDDE,
while the integral in the distributed-nDDE is computed using many past state
values, and also the architectures are different. Hence, a direct comparison
using the data processing inequality (Eq. 24) is not possible, but it provides
us with a plausible explanation.

In addition to the results just illustrated, we completed many other
experiments-1 to assess the sensitivity of our framework to various
hyperparameters. In all cases, the time-period corresponding to the training
data should be at least equal to one characteristic time-scale of the
dynamics, otherwise the prediction performance deteriorated, as shown in the
Fig. LABEL:fig:_Exp1_RMSE_with_time_train_time and discussed in Sec.
LABEL:SI:_Sensitivity_to_Network_Size_and_Training_Period_Length. Adding the
neural closure to the low-fidelity model improved its matching with the high-
fidelity data in nearly all cases. Its performance deteriorated with
increasing the length of the time-sequences used to form the batches, and also
with increasing the batch-size (the number of iterations per epoch is a
dependent hyperparameter as mentioned earlier). This also led to an increase
in training time. Depth of the networks affected the performance
significantly, with shallower networks performing poorly than deeper networks
as expected, however, the incremental gain in performance starts to taper off
after certain depths (see Fig.
LABEL:fig:_Exp_1_train_val_loss_all_network_size and Sec.
LABEL:SI:_Sensitivity_to_Network_Size_and_Training_Period_Length). Using an
exponentially decaying learning schedule over a constant learning rate
tremendously improved learning performance and reduced the number of epochs
needed. Further, training times slightly increased when using more delay times
in the case of discrete-nDDE. In general, we found that the training time for
discrete-nDDEs was similar to that for distributed-nDDEs. Such behaviors by
machine learning methods are difficult to anticipate in advance but they
should be mentioned.

Overall, in the experiments-1, we find that using memory-based neural closure
models as we derived from the Mori-Zwanzig formulation is advantageous over
just a Markovian closure. Using the new nDDEs as closure models helps maintain
generalizability of the learned models for longer time-periods, and
significantly reduces the longer-term prediction error of the ROM.

###  4.2 Experiments 2: Advecting Shock - Subgrid-Scale Processes

In the second experiments, we again use the advecting shock problem governed
by the Burger’s equation (Eq. 20), but we now reduce the computational cost of
the FOM by coarsening the spatial resolution, again leading to the need of a
closure model (Sec. 2.2). For the high-fidelity/high-resolution solution, we
employ a fine grid with the number of grid point in the $x$ direction
$N_{x}=100$, while for the low-fidelity/low-resolution solution, we employ a 4
times coarser grid with $N_{x}=25$. A comparison of high- and low- resolution
solutions solved using exactly the same numerical schemes and data stored at
every time-step of $\Delta t=0.01$ is provided in Fig. 4. We observe that by
decreasing the resolution, we introduce numerical diffusion and error in the
location of the shock peak at later times. The goal of the neural closure
models in these experiments is thus to augment the low-resolution model such
that it matches the sub-sampled/interpolated high-resolution solution at the
coarse (low-resolution) grid points.

(a)

(b)

(c)

(d)

Figure 4: Comparison of solutions of Burger’s equation (Eq. 20) for different
grid resolutions. (a): Solution for a high-resolution grid with number of grid
points, $N_{x}=100$; (b): Solution for a low-resolution grid with $N_{x}=25$;
(c): High-resolution solution interpolated onto the low-resolution grid. (d):
Absolute difference between fields in panels (b) and (c). We also provide a
pair of time-averaged errors, specifically: $L_{2}$ error; and RMSE
considering only the grid points where the error is at least 2% of the maximum
velocity value, denoted by RMSE($>$2%).

For training our neural closure models for the low-resolution discretization
with $N_{x}=25$, we use the same training regiment as in Experiments-1 (Sec.
4.1), with architectures details presented in Table
LABEL:table:_Exp1_and_2_architecture. In order to exploit the fact that each
grid point only affects its immediate neighbors over a single time-step, we
use 1-D convolutional layers for these experiments. For the nODE, we again
employ a deeper architecture with more trainable parameters, and for the
discrete-nDDE, six discrete delay values are again used. The values of the
other hyperparameters are in Sec. LABEL:SI:_hyperparameters. the validation
period from $t=1.25~{}\text{to}~{}2.5$, and the future prediction period from
$t=2.5~{}\text{to}~{}5.0$. We have chosen a prediction period of the combined
length of training and validation periods. For time-integration, we use the
Vode scheme [10] with adaptive time-stepping. The true data are generated by
interpolating the high-resolution solution onto the low-resolution grid
$(\\{\\{u^{true}(x_{k},T_{i})\\}_{k=1}^{N_{x}=25}\\}_{i=1}^{M})$, as shown in
Fig. 4c, and we use the time-averaged $L_{2}$ error,
$\mathcal{L}=\frac{1}{M}\sum_{i=1}^{M}\left(\sqrt{\sum_{k=1}^{N_{x}=25}|u^{pred}(x_{k},T_{i})-u^{true}(x_{k},T_{i})|^{2}}\right)$,
as the loss function.

The performance of the three neural closure models after 250 epochs (the
stochastic gradient descent nearly converges, see Fig.
LABEL:fig:_Exp_2_train_val_loss_all) is evaluated by taking the absolute
difference with the high-resolution solution interpolated onto the low-
resolution grid (Fig. 4c) spanning training, validation, and prediction
periods. We further benchmark our performance against the popular Smagorinsky
model [56] used for subgrid-scale turbulence closure in large eddy simulation
(LES). For the Burger’s Eq. 20, it introduces a dynamic turbulent eddy
viscosity ($\nu_{e}$) leading to,

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=\nu\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial}{\partial
x}\left(\nu_{e}\frac{\partial u}{\partial x}\right)\,,$ (25)

where $\nu_{e}=(C_{s}\Delta x)^{2}\big{|}\frac{\partial u}{\partial x}\big{|}$
and $C_{s}$ is the Smagorinsky constant. Results are shown in Fig. 5. The
details of the architectures employed are in Table
LABEL:table:_Exp1_and_2_architecture. As shown by the error fields of the
baseline (Fig. 4d) and closure models (Fig. 5), and by the corresponding pairs
of averaged error numbers (see Figs.), all closures improve the baseline.
However, the nODE and Smagorinsky closures only lead to a 55-60% decrease in
error, while the nDDE closures achieve a 80-90% decrease. Despite the deeper
architecture for the nODE, both the discrete-nDDE and distributed-nDDE (with
smaller architectures) again achieve smaller errors, for the whole period of
$t=0~{}\text{to}~{}5.0$. This means that they have lower numerical diffusion,
thus capturing the targeted subgrid-scale process. As opposed to the findings
of experiments-1 (Fig. 3), in the present experiments-2, the discrete-nDDE
performs slightly better than the distributed-nDDE in the prediction period.

(a)

(b)

(c)

(d)

Figure 5: Solutions of the Burger’s PDE on the low-resolution grid with
different closure models (left-column), and their absolute differences (right-
column) with the high-resolution solution interpolated onto the low-resolution
grid (Fig. 4c). For the trained neural closure models, the training period is
from $t=0~{}\text{to}~{}1.25$, the validation period from
$t=1.25~{}\text{to}~{}2.5$, and the prediction period from
$t=2.5~{}\text{to}~{}5.0$. For each closure, we also provide the pair of time-
averaged errors (see Fig. 4 for description). (a): Smagorinsky LES model with
$C_{s}=1.0$; (b), (c), (d): different neural closure models. These results
correspond to the architectures detailed in Table
LABEL:table:_Exp1_and_2_architecture.

We now study the effect of changing the amount of past information
incorporated in the closure model on the time-averaged $L_{2}$ error. For
this, we fix $\tau_{1}=0$ for the distributed-nDDE closure, and vary the
values of only $\tau_{2}$, keeping the architecture the same (Table
LABEL:table:_Exp1_and_2_architecture). First, for the time-averaged training
loss (not shown), we found no discernible trend and differences were mostly
due the stochastic gradient descent. Second, for the validation period, Fig.
6a shows the statistical summary of the validation losses (time-averaged
$L_{2}$ error) between the last epochs $200$ to $250$, for different delay-
period lengths. In order to ensure statistical soundness of the results,
$10-15$ repeats of the training were done, and the results aggregated for each
delay-period length. Results indicate that the validation loss first decreases
and then increases as we incorporate more-and-more past information, starting
from a very small delay period. For a specified architecture, neither too
little nor too much past information is helpful: there is likely an optimal
amount of information to incorporate. The initial improvement in the
performance of the closure models as the delay period is increased is due to
the increase in information content about the recent past. However, a
particular network architecture of finite size will have a limit on capturing
the increasing information content effectively due to its limited expressive
power, thus leading to a decrease in performance when too much information is
provided. An estimate of the range of delay period lengths to consider can be
obtained from properties of the given dynamical system such as the main time-
scales, e.g. advection and diffusion times-scales in the present system, and
main decorrelation times of state variables. Overall, from Fig. 6a, we can
notice the optimal delay period length to be around $0.075$. Similar trends
between performance of neural closure models and delay period lengths were
also found in Experiments-1 (not shown). Some published studies attempt to
derive analytical expressions for the optimal memory length, making many
approximations in the process [51, 65]. A final option is to learn the delay
amount as a part of the training process itself, however, this requires
modified adjoint equations.

We conducted again a series of experiments-2 to assess the sensitivity to the
various other hyperparameters, and found similar trends (not shown here) as in
the experiments-1. Finally, we noticed that using the dopri5 [34] time-
integration scheme severely impaired the learning ability in the
experiments-2.

(a)

(b)

Figure 6: Variation of distributed-nDDE closure validation loss (time-averaged
$L_{2}$ error) averaged over the last 50 training epoch for Experiments-2 &
3a. All the experiments have $\tau_{1}=0$, and different $\tau_{2}$
(horizontal-axis). Note that $\tau_{2}=0$ corresponds to the nODE closure. We
use boxplots to provide statistical summaries for multiple training repeats
done for each experiment. The box and its whiskers provide a five number
summary: minimum, first quartile (Q1), median (orange solid line), third
quartile (Q3), and maximum, along with outliers (black circles) if any.

Overall, these results demonstrate the superiority of using our new memory-
based closure models in capturing subgrid-scale processes.

###  4.3 Experiments 3a: 0-D Marine biological Models

For our third experiments, we use neural closure models to incorporate the
effects of missing processes and state variables in lower-complexity
biological models, thus targeting the third class of closure modeling (Sec.
2.3). Marine biological models are based on ODEs that describe the food-web
interactions in the ecosystem. They can vary greatly in terms of complexity
[26]. The marine biological models used in our experiments are adapted from
Newberger et. al., 2003 ([62]). They were used to simulate the ecosystem in
the Oregon coastal upwelling zone. They provide hierarchical embedded models
compatible with each other. We employ the three-component NPZ model (nutrients
($N$), phytoplankton ($P$), and zooplankton ($Z$)), and the five-component
NNPZD model (ammonia ($NH_{4}$), nitrate ($NO_{3}$), $P$, $Z$, and detritus
($D$)) in a zero-dimensional setting (0-D; only temporal variation). The low
complexity NPZ model is given by,

$\begin{split}\frac{dN}{dt}&=-G\frac{PN}{N+K_{u}}+\Xi P+\Gamma Z+R_{m}\gamma
Z(1-\exp^{-\Lambda P})\\\ \frac{dP}{dt}&=G\frac{PN}{N+K_{u}}-\Xi
P-R_{m}Z(1-\exp^{-\Lambda P})\\\
\frac{dZ}{dt}&=R_{m}(1-\gamma)Z(1-\exp^{-\Lambda P})-\Gamma Z\\\
\text{with}\quad N(0)&=T_{bio},\quad P(0)=0,\quad\text{and}\quad
Z(0)=0\,,\end{split}$ (26)

with $G$ representing the optical model,

$G=V_{m}\frac{\alpha
I}{(V_{m}^{2}+\alpha^{2}I^{2})^{1/2}},\quad\text{and}\quad
I(z)=I_{0}\exp(k_{w}z)\,,$ (27)

where $z$ is depth and $I(z)$ models the availability of sunlight for photo-
chemical reactions. The parameters in Eqs. 26 and 27 are: $k_{w}$, light
attenuation by sea water; $\alpha$, initial slope of the $P-I$ curve; $I_{0}$,
surface photosynthetically available radiation; $V_{m}$, phytoplankton maximum
uptake rate; $K_{u}$, half-saturation for phytoplankton uptake of nutrients;
$\Xi$, phytoplankton specific mortality rate; $R_{m}$, zooplankton maximum
grazing rate; $\Lambda$, Ivlev constant; $\gamma$, fraction of zooplankton
grazing egested; $\Gamma$, zooplankton specific excretion/mortality rate; and
$T_{bio}$, total biomass concentration. In the NPZ model (Eq. 26), the
nutrient uptake by phytoplankton is governed by a Michaelis-Menten
formulation, which amounts to a linear uptake relationship at low nutrient
concentrations that saturates to a constant at high concentrations. The
grazing of phytoplankton by zooplankton follows a similar behavior: their
growth rate becomes independent of $P$ in case of abundance, but proportional
to available $P$ when resources are scarce, hence zooplankton grazing is
modeled by an Ivlev function. The death rates of both $P$ and $Z$ are linear,
and a portion of zooplankton grazing in the form of excretion goes directly to
nutrients.

In the higher complexity NNPZD model, the nutrients are divided into ammonia
and nitrates, which are the two most important forms of nitrogen in the ocean.
With the intermediate of decomposed organic matter, detritus, the NNPZD model
captures new processes such as: phytoplankton cells preferentially taking up
ammonia over nitrate because the presence of ammonia inhibits the activity of
the enzyme nitrate reductase essential for the uptake kinetics; the pool of
ammonium coming from remineralization of detritus; and part of this ammonium
pool getting oxidized to become a source of nitrate called the process of
nitrification, etc. Overall, the NNPZD model is given by,

$\begin{split}\frac{dNO_{3}}{dt}=&~{}\Omega
NH_{4}-G\left[\frac{NO_{3}}{NO_{3}+K_{u}}\exp^{-\Psi NH_{4}}\right]P\\\
\frac{dNH_{4}}{dt}=&~{}-\Omega NH_{4}+\Phi D+\Gamma
Z-G\left[\frac{NH_{4}}{NH_{4}+K_{u}}\right]P\\\
\frac{dP}{dt}=&~{}G\left[\frac{NO_{3}}{NO_{3}+K_{u}}\exp^{-\Psi
NH_{4}}+\frac{NH_{4}}{NH_{4}+K_{u}}\right]P-\Xi P-R_{m}Z(1-\exp^{-\Lambda
P})\\\ \frac{dZ}{dt}=&~{}R_{m}(1-\gamma)Z(1-\exp^{-\Lambda P})-\Gamma Z\\\
\frac{dD}{dt}=&~{}R_{m}\gamma Z(1-\exp^{-\Lambda P})+\Xi P-\Phi D\\\
\text{with}\quad&NO_{3}(0)=T_{bio}/2,\quad NH_{4}(0)=T_{bio}/2,\quad
P(0)=0,\quad Z(0)=0,\;\;\text{and}\;\;D(0)=0\,,\end{split}$ (28)

where the new parameters are: $\Psi$, $NH_{4}$ inhibition parameter; $\Phi$,
detritus decomposition rate; and $\Omega$, $NH_{4}$ oxidation rate.

Solutions of the above two models are presented in Fig. 7. Different values of
the parameters and initial conditions set these models in different dynamical
regimes. From the responses in time, the present solutions in experiments-3a
are in a stable nonlinear limit-cycle regime. The $N$ class in the NPZ model
is a broader class encompassing $NO_{3}$, $NH_{4}$, and $D$ from the NNPZD
model. Since the NNPZD model resolves many more processes, the concentrations
of $NO_{3}+NH_{4}+D$, $P$, and $Z$ differ significantly from the $N$, $P$, and
$Z$ of the NPZ model. The goal of the neural closure models in these
experiments is thus to augment the low-complexity NPZ model such that it
matches the aggregated states of the high-complexity NNPZD model.

(a)

(b)

(c)

Figure 7:  Solutions of the marine biological models used in Experiments-3a
(concentrations vs. time in days). Parameter values used are (adopted from
[62]): $k_{w}=0.067~{}m^{-1}$, $\alpha=0.025~{}(W~{}m^{-2}~{}d)^{-1}$,
$V_{m}=1.5~{}d^{-1}$, $I_{0}=158.075~{}W~{}m^{-2}$,
$K_{u}=1~{}mmol~{}N~{}m^{-3}$, $\Psi=1.46~{}(mmol~{}N~{}m^{-3})^{-1}$,
$\Xi=0.1~{}d^{-1}$, $R_{m}=1.52~{}d^{-1}$,
$\Lambda=0.06~{}(mmol~{}N~{}m^{-3})^{-1}$, $\gamma=0.3$,
$\Gamma=0.145~{}d^{-1}$, $\Phi=0.175~{}d^{-1}$, $\Omega=0.041~{}d^{-1}$,
$z=-25~{}m$, and $T_{bio}=30~{}mmol~{}N~{}m^{-3}$. (a): Nutrient-
Phytoplankton-Zooplankton (NPZ) model (Eq. 26); (b): Nitrate-Ammonia-
Phytoplankton-Zooplankton-Detritus (NNPZD) model (Eq. 28); (c): Comparison
between $NO_{3}+NH_{4}+D$, $P$, and $Z$ from the NNPZD model (solid) with $N$,
$P$ and $Z$ from the NPZ model (dashed-dot).

For training our neural closure models for the NPZ model, we use the same
training regiment as in Experiments-1 & 2 (Secs. 4.1 & 4.2), with
architectures details presented in Table LABEL:table:_Exp3_and_4_architecture.
For the nODE, we again employ a bigger architecture, and for the discrete-
nDDE, six discrete delay values are again used. The values of other
hyperparameters are given in Sec. LABEL:SI:_hyperparameters. The training
period ranges from $t=0~{}\text{to}~{}30~{}days$, validation period from
$t=30~{}\text{to}~{}60~{}days$; and the prediction period from
$t=60~{}\text{to}~{}330~{}days$. We have chosen a prediction period nine times
longer than the training period. For biological ODE models, there exists
invariant knowledge about the system, such as biological state variables
cannot be negative, and the sum of all the states remains constant with time
(this can be verified by summing the ODEs of NPZ or NNPZD models). We enforce
the constraints as follows. The positivity is enforced as a penalization term
in the loss function. The constant total biomass constraint is embedded in the
architectures of neural closures by introducing a new custom layer named,
BioConstrainLayer. This layer is applied at the end, and expects an input of
size 1. The output of this layer is formed by splitting the input into three
with the proportions, $\beta$, $-1$, and $1-\beta$; where $\beta$ is a
trainable parameter. This ensures that summing the right hand side of the
augmented NPZ model does not leave any new residual due to the neural closure
terms. The stiffness of such biological ODE models also poses a challenge in
maintaining these desired properties [12]. The flexibility of our framework
however allows the use of appropriate time-stepping schemes, such as A-stable
implicit schemes, etc., to overcome stiffness. The true data are generated by
aggregating the variables of the NNPZD model $(N\equiv
NO_{3}+NH_{4}+D,~{}P,~{}\text{and}~{}Z;~{}\\{\\{B^{true}(T_{i})\\}_{B\in\\{N,P,Z\\}}\\}_{i=1}^{M})$.
Finally, we use the dopri5 [34] scheme with adaptive time-stepping and
simulation data were stored at every $\Delta t=0.05$ days for all our time-
integration requirements, along with a $L_{2}$ error loss function,
$\mathcal{L}=\allowdisplaybreaks\frac{1}{M}\sum_{i=1}^{M}\left(\sqrt{\sum_{B\in\\{N,P,Z\\}}|B^{pred}(T_{i})-B^{true}(T_{i})|^{2}}\right)$.

The performance of the three neural closure models augmenting the NPZ is
evaluated after 350 epochs of training (the stochastic gradient descent nearly
converges, as evident from the Fig. LABEL:fig:_Exp_3_train_val_loss_all) by
comparison with the aggregated biology variables from the high-complexity
NNPZD model (Eq. 28) spanning training, validation, and prediction periods.
Results are presented in Figure 8. The details of the architectures employed
are in Table LABEL:table:_Exp3_and_4_architecture. When compared with the
aggregated NNPZD variables (true variables), we find again that despite the
bigger architecture of the nODE, it starts to develop significant errors
around $t=180~{}days$ and quickly gets out-of-phase thereafter (Fig. 8a). The
discrete-nDDE and distributed-nDDE, both with smaller architectures, are
however able to match the true variables for nearly the whole period of
$t=0~{}\text{to}~{}330~{}days$ (Fig. 8b), with only distributed-nDDE starting
to getting out-of-phase after $t=270~{}days$ (Fig. 8c) at the end of the long
prediction period. These results are corroborated by the time evolution of the
RMSE and average cross-correlation for the three variables over the prediction
period (Fig. 8d). From the progression of the time-averaged $L_{2}$ loss
(here, the error between the variables from the closure-model-augmented NPZ
system, and the true variables), the nODE performs either equally well or
better than both discrete-nDDE and distributed-nDDE during training and
validations periods (Fig. LABEL:fig:_Exp_3_train_val_loss_all), however, it is
not able to maintain long-term accuracy. We also notice very large spikes in
the first half of the training regime, which are due to weights of the neural-
networks taking values that lead to negative biology variables. As training
progresses, we however don’t observe this behavior anymore because the
trainable weights starts to converge towards biologically feasible regimes. In
conclusion, using a memory based closure for a low-complexity model can
efficiently help emulate the high-complexity model dynamics.

(a)

(b)

(c)

(d)

Figure 8:  Comparison of the biological variables from the learned NPZ model
augmented with the three neural closure models (dashed), aggregated variables
from the NNPZD model (ground truth; solid), and variables from the NPZ model
(dashed-dot) at the end of training. For each neural closure, the training
period is from $t=0~{}\text{to}~{}30$ days, the validation period is from
$t=30~{}\text{to}~{}60$ days, while prediction period is from
$t=60~{}\text{to}~{}330$ days. (a), (b), (c): different neural closure models;
(d): the left plot shows the evolution of root-mean-squared-error (RMSE), and
the right plot shows the average cross-correlation (only for the prediction
period) w.r.t. the ground truth. These results correspond to the architectures
detailed in Table LABEL:table:_Exp3_and_4_architecture.

As for experiments-2, we conducted a series of experiments-3a, to study the
effect of changing the amount of past information incorporated in the neural
closure models. In Fig. 6b, we show the variation of the average validation
loss (time-averaged $L_{2}$ error) between the last epochs $300$ to $350$, for
different delay-period lengths ($\tau_{1}=0$, and $\tau_{2}$ varying in case
of distributed-nDDE). In order to ensure statistical soundness of the results,
$10-12$ repeats of the training were done, and the results were aggregated for
each delay-period length (excluding the runs which diverged). We again find an
optimal memory length for a specified architecture, however, with more and
more runs failing to converge for longer delay period lengths. For the present
system, estimates of delay period lengths to consider can be obtained from the
time-scales of biological behaviors and adjustments, and from the
decorrelation times of the biological state variables. Taking into account the
limited effectiveness of a network architecture of finite size for capturing
increasing information content, from Fig. 6b, we find an optimal delay period
length to be around $1~{}day$. We also conducted a series of experiments-3a to
study the sensitivity to the various hyperparameters, and found similar trends
(not shown here) as in experiments-1 & 2\. For good performance, we further
found that using a small enough time-step was critical as well as limiting the
number of internal steps in the dopri5 [34] time-integration scheme, while
penalizing negative values in the loss function did not make much of a
difference. Whenever multiple terms are present in the loss function enforcing
different inherent properties of the system, they should be normalized (e.g.
using non-dimensional variables) and given appropriate relative weights.

In general, the ecosystem ODEs are coupled with regional or global ocean
modeling systems, leading to advection-diffusion-reaction PDEs [4]. If highly
complex ecosystem models are employed, a very large number of PDE state
variables need to be solved for, rendering the computations very expensive. A
large number of unknown parameter values as well as uncertain initial
conditions then also need to be estimated, requiring specific methods [e.g.
49]. The available biogeochemical observations are not always sufficient for
calibrating these many unknown parameters and for estimating initial
conditions of high-complexity models. If the corresponding errors are large,
this can lead to integrating models in the wrong dynamical regimes [e.g. 87].
Finally, in some applications, one is only interested in the aggregated state
variables as the output, but cannot use low-complexity models because their
dynamics are too inaccurate for the goals of the applications. Using neural
closure models as shown here, one can increase the accuracy of the low-
complexity models to match the response of high-complexity models (possibly up
to models such as ERSEM [3]) without adding the computational burden of
modeling all the intermediate biological states and processes, while reducing
the effects of other uncertainties listed above. Results of our neural
closures in 1-D PDEs is showcased next.

###  4.4 Experiments 3b: 1-D Marine Biogeochemical Models

For our final set of experiments, we extend the ODE models used in
Experiments-3a (Sec. 4.3) to contain a vertical dimension (thus, 1-D) and
vertical eddy mixing parameterized by the operator, $\partial/\partial
z\left(K_{z}(z,M)\partial/\partial z(\bullet)\right)$, where $K_{z}$ is a
dynamic eddy diffusion coefficient. A mixed layer of varying depth ($M=M(t)$)
is used as a physical input to the ecosystem models. Thus, each biological
state variable $B(z,t)$ is governed by the following non-autonomous PDE,

$\frac{\partial B}{\partial t}=S^{B}+\frac{\partial}{\partial
z}\left(K_{z}(z,M(t))\frac{\partial B}{\partial z}\right)\;,$ (29)

where $S^{B}$ are the corresponding biology source terms, which also makes it
stiff. The dynamic depth-dependent diffusion parameter $K_{z}$ is given by,

$K_{z}(z,M(t))=K_{z_{b}}+\frac{(K_{z_{0}}-K_{z_{b}})(\arctan(-\gamma(M(t)-z))-\arctan(-\gamma(M(t)-D)))}{\arctan(-\gamma
M(t))-\arctan(-\gamma(M(t)-D))}\;,$ (30)

where $K_{z_{b}}$ and $K_{z_{0}}$ are the diffusion at the bottom and surface
respectively, $\gamma$ is the thermocline sharpness, and $D$ is the total
depth. The 1-D model and parameterizations are adapted from Eknes and Evensen,
2002 [23], and Newberger et. al., 2003 [62]. They simulate the seasonal
variability in upwelling, sunlight, and biomass vertical profiles. The dynamic
mixed layer depth, surface photosynthetically-available radiation $I_{0}(t)$,
and biomass fields $B(z,t)$ are shown in Figure 9a. The radiation $I_{0}(t)$
and total biomass concentration, $T_{bio}(z,t)$, affect $S^{B}$ and the
initial conditions.

For these Experiments-3b, we consider $20$ grid points in the vertical and use
the dopri5 [34] scheme with adaptive time-stepping. Data is stored at every
time-step of $\Delta t=0.1~{}days$ for all our time-integration requirements.
Solutions of aggregated states of the high-complexity 1-D NNPZD model (true
data) and their absolute difference with the corresponding low-complexity 1-D
NPZ model states are provided in Figs. 9a & 9b, respectively. For training our
neural closure models for the 1-D NPZ model, we use the same training regiment
as in Experiments-1, 2, & 3a (Secs. 4.1, 4.2, & 4.3). We note that in the 1-D
NPZ model, the local mixing across depths occurs only due to the eddy
diffusion term, and not to the biology source terms. Thus, we employ 1-D
convolutional layers with receptive fields of size 1. We again use the custom
layer, BioConstrainLayer (Sec. 4.3), to ensure that the sum of the biology
source terms of the augmented 1-D NPZ model does not leave any new residual
due to the neural closure terms. Along with this, we define a new custom
layer, called AddExtraChannels, to add additional channels to the input of
this layer. We add one for the depths at different grid points, and the other
for the corresponding values of available sunlight for photo-chemical
reactions ($I(z,t)$). The architectures details for the three closure models
used are presented in Table LABEL:table:_Exp3_and_4_architecture. For the
nODE, we again employ a bigger architecture, and for the discrete-nDDE, only
four discrete delay values are used. The values of other hyperparameters are
given in Sec. LABEL:SI:_hyperparameters. The training period ranges from $t=0$
to $30~{}days$ and validation period from $t=30$ to $60~{}days$, both within
the first season. The prediction period, however, ranges from $t=60$ to
$364~{}days$: it is more than 10 times longer than the training period and
involves the four seasons. Together, the three periods span a full year. For
loss function, we combine the $L_{2}$ errors, considering all the biological
states computed for individual depths, and then averaged over all the depths
and times,
$\mathcal{L}=\allowdisplaybreaks\frac{1}{M}\sum_{i=1}^{M}\left(\frac{1}{N_{z}}\sum_{k=1}^{N_{z}=20}\left(\sqrt{\sum_{B\in\\{N,P,Z\\}}|B^{pred}(z_{k},T_{i})-B^{true}(z_{k},T_{i})|^{2}}\right)\right)$.

The performance of the three neural closure models augmenting the 1-D NPZ
model is evaluated after 200 epochs of training (the stochastic gradient
descent nearly converges, see Fig. LABEL:fig:_Exp_4_train_val_loss_all). The
truth fields are the aggregated biology variables from the high-complexity 1-D
NNPZD model (Eqs. 28 & 29) spanning training, validation, and prediction
periods. Results are presented in Figure 9. We find again that despite the
bigger architecture for the nODE case, it develops spurious oscillations
around $t=250~{}days$. The discrete-nDDE and distributed-nDDE, both with
smaller architectures, however match well with the true variables for nearly
the full year of simulation, about 10 months of which is future prediction.
The distributed-nDDE performs slightly better than its counterpart. In Fig. 9,
we also provide averaged error numbers for the baseline (Fig. 9b) and the
different closure models, all of which improve the baseline. As in
Experiments-3a, we again notice large spikes in the starting of the training
regime, for the same reason as given earlier, and similar trends for
hyperparameter sensitivity. We also found that the Experiments-3b were
affected by the choice of loss function. For example, using $L_{2}$ error
computed for each biological state vector (containing values for all the
depths) and then averaging over the number of biological states and times
deteriorated the quality of learning.

Despite the presence of complex physical processes and relatively large
dimensions as compared to the previous experiments, the nDDEs closures were
found to effectively match the high-complexity model and maintain long-term
accuracy.

(a)

(b)

(c)

(d)

(e)

Figure 9:  Comparison of the 1-D physical-biogeochemical PDE models used in
Experiments-3b with and without closure models. Along with the parameter
values mentioned in Figure 7, we consider: a sinusoidal variation in
$I_{o}(t)$; linear vertical variation in total biomass $T_{bio}(z)$ from
$10~{}mmol~{}N~{}m^{-3}$ at the surface to $30~{}mmol~{}N~{}m^{-3}$ at
$z=100~{}m$; $K_{z_{b}}=0.0864~{}(m^{2}/day)$; $K_{z_{0}}=8.64~{}(m^{2}/day)$;
$\gamma=0.1~{}m^{-1}$; and $D=-100~{}m$, all adapted from [62, 23]. For the
neural closure models, the training period is from $t=0~{}\text{to}~{}30$
days, the validation period from $t=30~{}\text{to}~{}60$ days, and the long
future prediction period from $t=60~{}\text{to}~{}364$ days. (a): Top plots
show the yearly variation of solar radiation and the bottom plots the
aggregated states from the NNPZD model (ground truth) overlaid with the
dynamic mixed layer depth in dashed red lines. In the subsequent plots (b),
(c), (d), and (e), we show the absolute difference of the different neural
closure cases with the ground truth. For each case, we also provide the pair
of time-averaged errors (see Fig. 4 for description). These results correspond
to the architectures given in Table LABEL:table:_Exp3_and_4_architecture.

###  4.5 Computational Complexity

It is crucial to analyze complexity and in particular the cost of adding a
neural closure model to a low-fidelity model. In this section, we analyze the
computational complexity in terms of flop (floating point operations) count
for evaluating the right-hand-side (RHS) of the low-fidelity models, and for
the forward-pass of the neural closure models [61]. We will also comment on
the training costs. The Burger’s PDE considered for Experiments-1&2 (Secs. 4.1
& 4.2) has a nonlinear advection term. Hence, for the POD-GP ROM and the FOM,
the upper flops is of the order of the square of number of resolved modes and
of the spatial grid resolution, respectively. In general, for reaction terms
and biogeochemical systems, the number of nonlinear parameterizations present
are of the order of the number of components in the model. Hence, even for
Experiments-3a (Sec. 4.3), the upper flops is of the order of the square of
the number of biological components. For Experiments-3b (Sec. 4.4), the upper
flops is also affected by the diffusion terms. Let the number of state
variables in the low-fidelity models be $N\in\mathbb{N}$, thus the leading
order computational complexity would be $\mathcal{O}(cN^{2})$, where
$c\in\mathbb{R}^{+}$ is some constant dependent on the numerical schemes used
for spatial discretization, the exact functional form of the RHS, etc.

Now, when neural closure models are added to the low-fidelity models, the time
integration requires a forward pass through the neural-network. This cost
varies with the neural architecture and model type, here either a fully-
connected or convolutional, and discrete-nDDE or distributed-nDDE,
respectively. As observed in our experiments, using delays in the closure
model enables us to use shallower networks, with a depth independent of the
number of state variables ($N$). We also found that the width of the networks
was similar to, or smaller than, $N$. In case of distributed-nDDEs, we
observed that the width of the auxiliary network ($g_{NN}$) could be on an
average nearly half the size of the main network ($f_{NN}$). Let the size of
the hidden state for the RNN in discrete-nDDEs be $N_{h}\in\mathbb{N}$, and
the number of neurons in the hidden layers of the main and auxiliary networks
in the case of distributed-nDDEs be $N_{h}$ and $N_{h}/2$, respectively, with
$N_{h}\lesssim N$. It could be easily shown that the leading order complexity
for a single iteration of RNN would be $\mathcal{O}(N_{h}^{2}+N_{h}N)$, which
is due to the hidden and input state vectors being multiplied by the weight
matrices, while the application of activation function would be
$\mathcal{O}(N_{h})$ only. As the number of discrete-delays in discrete-nDDEs
are independent of $N$ and $\mathcal{O}(1)$, it does not affect the complexity
of the RNN. The complexity of the first hidden layer and/or the output layer
of the deep neural-networks used in discrete-nDDEs and distributed-nDDEs (main
network, $f_{NN}$) will be $\mathcal{O}(N_{h}N)$ and
$\mathcal{O}(N_{h}(N+N_{h}/2))$, respectively, while for the remaining hidden
layers, it will be $\mathcal{O}(N_{h}^{2})$. Focusing on the integral of the
auxiliary network ($g_{NN}$) over the delay period in distributed-nDDEs, if
implemented efficiently, at every time-step, we only need to compute the
integral twice over periods of size $\Delta t$, each adjacent to the ends of
the present delay-period. We can add and subtract these integrals over $\Delta
t$ periods to compute the overall integral in a rolling window sense. Hence,
the contribution to the computational cost by the auxiliary network would be
$\mathcal{O}(N_{h}N/2)$ (for first hidden layer) and
$\mathcal{O}(N_{h}^{2}/4)$ (subsequent hidden layers). Considering
$D\in\mathbb{N}$ as the depth for all of the networks considered, the
complexity for the forward pass through the discrete-nDDE closure is
$\mathcal{O}(DN_{h}^{2}+DN_{h}N)$, while for distributed-nDDE, it is
$\mathcal{O}((3/2)N_{h}N+(7/4)DN_{h}^{2})$. These costs were computed
considering fully-connected layers, however, will only be cheaper in case of
convolutional layers. Thus, the additional computational cost due to the
presence of neural closure models is of similar or lower complexity than the
existing low-fidelity model.

Estimating the computational cost/complexity of training in flops is not
common because apart from time-integrating the forward model and adjoint
equations, there are many other operations such as here: automatic
differentiation through the neural networks; creation and use of interpolation
functions; the integral to compute the final derivatives; the gradient descent
step, etc. The overall cost also depends on the number of epochs needed for
convergence. The present training cost is of course non-negligible, as with
any supervised learning algorithm. However, in applications where one needs to
repeatedly solve a low-fidelity model, investing in a one-time cost of
training a neural closure model can later lead to accuracy close to that of
the high-fidelity model with only a small increase in the computational cost
of the low-fidelity model.

##  5 Conclusion

We developed a novel, versatile, rigorous, and unified methodology to learn
closure parameterizations for low-fidelity models using data from high-
fidelity simulations. The Mori-Zwanzig formulation [84, 17, 32] and the
presence of inherent delays in complex dynamical systems [24], especially
biological systems [31, 86], justify the need for non-Markovian closure
parameterizations. To learn such non-Markovian closures, our new neural
closure models extend neural ordinary differential equations (nODEs; [16]) to
neural delay differential equations (nDDEs). Our nDDEs do not require access
to the high-fidelity model or frequent enough and uniformly-spaced high-
fidelity data to compute the time derivative of the state with high accuracy.
Further, it enables the accounting of errors in the time-evolution of the
states in the presence of neural networks during training. We derive the
adjoint equations and network architectures needed to efficiently implement
the nDDEs, for both discrete and distributed delays, agnostic to the specifics
of the time-integration scheme, and capable of handling stiff systems. For
distributed-delays, we propose a novel architecture consisting of two coupled
deep neural networks, which enables us to incorporate memory without the use
of any recurrent architectures.

Through simulation experiments, we showed that our methodology drastically
improves the long-term predictive capability of low-fidelity models for the
main classes of model truncations. Specifically, our neural closure models
efficiently account for truncated modes in reduced-order-models (ROMs),
capture the effects of subgrid-scale processes in coarse models, and augment
the simplification of complex biological and non-autonomous physical-
biogeochemical models. Our first two classes of simulation experiments utilize
the advecting shock problem governed by the Burger’s PDE, with its low-
fidelity models derived either by proper-orthogonal-decomposition Galerkin
projection or by reducing the spatial grid resolution. Our third class of
experiments considers marine biological ODEs of varying complexities and their
physical-biogeochemical PDE extensions with non-autonomous dynamic
parameterizations. The low-fidelity models are obtained by aggregation of
components and other simplifications of processes and parameterizations. In
each of these classes, results consistently show that using non-Markovian over
Markovian closures improves the accuracy of the learned system while also
requiring smaller network architectures. Our use of the known-physics/low-
fidelity model also helps to reduce the required size of the network
architecture and the number of time samples for the training data. We also
outperform classic dynamic closures such as the Smagorinsky subgrid-scale
model. These results are obtained using stringent evaluations: we compare the
performance of the learned system for the training period (during which high-
fidelity data snapshots are used for training) and validation period (during
which hyperparameter tuning occurs) as often done, but we also compare it for
much longer-term future prediction periods with no overlap with the preceding
two. We even consider a prediction period reaching $10$ times the length of
the training/validation period, thus successfully demonstrating the
extrapolation capabilities of nDDE closures.

In our experiments, we find that just using a few numbers of discrete delays
might perform equally well or better than using a distributed delay which
involves an integral of the state variable over a delay period. We provide a
plausible explanation of this counter-intuitive observation using the data
processing inequality from information theory. We also show that there exists
an optimal amount of past information to incorporate for a specified
architecture and the relevant time-scales present in the dynamical system,
thus indicating that neither too little nor too much past information is
helpful. Finally, a computational complexity analysis using flop (floating
point operation) count proves that the additional computational cost due to
the presence of our neural closure models is of similar or lower complexity
than the existing low-fidelity model.

The present work provides a unified framework to learn non-Markovian closure
parameterization using delay differential equations and neural networks. It
enables the use of the often elusive Mori-Zwanzig formulation [84, 17, 32] in
its full glory without unjustified approximations and simplifications. Our
nDDE closures are not just limited to the shown experiments, but could be
widely extended to other fields such as control theory, robotics,
pharmacokinetic-pharmacodynamics, chemistry, economics, biological regulatory
systems, etc.

Data Accessibility. The codes and data used in this work are available in the
GitHub repository: https://github.com/mit-mseas/neuralClosureModels.git

Author’s Contributions. A.G. conceived the idea of using neural delay
differential equations for closure parameterizations; derived the adjoint
equations; implemented the neural network architectures and the simulation
experiments; interpreted the computational results; and wrote a first draft of
the manuscript. P.F.J.L. supervised the work; conceived the ideas to extend
the methodology to capture the effects of subgrid-scale processes in coarse
models, and to augment the simplification of complex mathematical models;
interpreted the computational results; and edited and wrote significant parts
of the manuscript.

Competing Interests. We declare we have no competing interests.

Funding. We are grateful to the Office of Naval Research for partial support
under grant N00014-20-1-2023 (MURI ML-SCOPE) to the Massachusetts Institute of
Technology. We also thank MathWorks and the Mechanical Engineering Department
at MIT for awarding a competitive 2020-2021 MathWorks Mechanical Engineering
Fellowship for A.G. Some of the computations were made possible due the Google
Cloud Platform research credits, which we gratefully acknowledge.

Acknowledgement. We thank the members of our MSEAS group for their
collaboration and insights, especially Mr. Aaron Charous. We thank our ML-
SCOPE team for many useful discussions, especially Dr. Mickaël Chekroun for
his comments on the final draft of the manuscript. We also thank the anonymous
reviewers for their constructive feedback which helped improve the manuscript.

## References

  * Abadi et al. [2015] M. Abadi et al.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Alfonsi [2009] G. Alfonsi.  Reynolds-averaged navier–stokes equations for turbulence modeling.  _Applied Mechanics Reviews_ , 62(4), 2009.
  * Baretta et al. [1995] J. Baretta, W. Ebenhöh, and P. Ruardij.  The European regional seas ecosystem model, a complex marine ecosystem model.  _Netherlands Journal of Sea Research_ , 33(3-4):233–246, 1995.
  * Beşiktepe et al. [2003] Ş. T. Beşiktepe, P. F. J. Lermusiaux, and A. R. Robinson.  Coupled physical and biogeochemical data-driven simulations of Massachusetts Bay in late summer: Real-time and post-cruise data assimilation.  _J. of Marine Systems_ , 40–41:171–212, 2003.
  * Behzad et al. [2015] F. Behzad et al.  On the sensitivity and accuracy of proper-orthogonal-decomposition-based reduced order models for Burgers equation.  _Computers & Fluids_, 106:19–32, 2015.
  * Bhattacharya et al. [1982] K. Bhattacharya, M. Ghil, and I. Vulis.  Internal variability of an energy-balance model with delayed albedo effects.  _Journal of the Atmospheric Sciences_ , 39(8):1747–1773, 1982.
  * Bocharov and Rihan [2000] G. A. Bocharov and F. A. Rihan.  Numerical modelling in biosciences using delay differential equations.  _Journal of Computational and Applied Mathematics_ , 125(1-2):183–199, 2000.
  * Boers et al. [2017] N. Boers, M. D. Chekroun, et al.  Inverse stochastic–dynamic models for high-resolution greenland ice core records.  _Earth System Dynamics_ , 8(4):1171–1190, 2017\.
  * Borja et al. [2014] A. Borja et al.  Tales from a thousand and one ways to integrate marine ecosystem components when assessing the environmental status.  _Frontiers in Marine Science_ , 1:72, 2014.
  * Brown et al. [1989] P. N. Brown, G. D. Byrne, and A. C. Hindmarsh.  VODE: A variable-coefficient ODE solver.  _SIAM journal on scientific and statistical computing_ , 10(5):1038–1051, 1989.
  * Brunton et al. [2016] S. L. Brunton, J. L. Proctor, and J. N. Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  _PNAS_ , 113(15):3932–3937, 2016.
  * Burchard et al. [2005] H. Burchard, E. Deleersnijder, and A. Meister.  Application of modified patankar schemes to stiff biogeochemical models for the water column.  _Ocean Dynamics_ , 55(3-4):326–337, 2005.
  * Calver and Enright [2017] J. Calver and W. Enright.  Numerical methods for computing sensitivities for ODEs and DDEs.  _Numerical Algorithms_ , 74(4):1101–1117, 2017\.
  * Chassignet and Verron [2012] E. P. Chassignet and J. Verron.  _Ocean modeling and parameterization_ , volume 516.  Springer Science & Business Media, 2012.
  * Chekroun et al. [2019] M. D. Chekroun, H. Liu, and J. C. McWilliams.  Variational approach to closure of nonlinear dynamical systems: Autonomous case.  _Journal of Statistical Physics_ , pages 1–88, 2019.
  * Chen et al. [2018] T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud.  Neural ordinary differential equations.  In _Advances in neural information processing systems_ , pages 6571–6583, 2018.
  * Chorin et al. [2000] A. J. Chorin, O. H. Hald, and R. Kupferman.  Optimal prediction and the Mori–Zwanzig representation of irreversible processes.  _PNAS_ , 97(7):2968–2973, 2000.
  * Cover [1999] T. M. Cover.  _Elements of information theory_.  John Wiley & Sons, 1999.
  * Dauhajre et al. [2019] D. P. Dauhajre, J. C. McWilliams, and L. Renault.  Nearshore lagrangian connectivity: submesoscale influence and resolution sensitivity.  _JGR: Oceans_ , 124(7):5180–5204, 2019.
  * Daulbaev et al. [2020] T. Daulbaev et al.  Interpolation technique to speed up gradients propagation in neural ODEs.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Dell’Anna [2020] L. Dell’Anna.  Solvable delay model for epidemic spreading: the case of covid-19 in italy.  _Scientific Reports_ , 10(1), Sept. 2020.
  * Diekmann et al. [2012] O. Diekmann, S. A. Van Gils, S. M. Lunel, and H.-O. Walther.  _Delay equations: functional-, complex-, and nonlinear analysis_ , volume 110.  Springer Science & Business Media, 2012.
  * Eknes and Evensen [2002] M. Eknes and G. Evensen.  An ensemble kalman filter with a 1-d marine ecosystem model.  _Journal of Marine Systems_ , 36(1-2):75–100, 2002.
  * Erneux [2009] T. Erneux.  _Applied delay differential equations_ , volume 3.  Springer, 2009.
  * Faugeras and Maury [2007] B. Faugeras and O. Maury.  Modeling fish population movements: from an individual-based representation to an advection–diffusion equation.  _J. of Theoretical Bio._ , 247(4):837–848, 2007\.
  * Fennel and Neumann [2014] W. Fennel and T. Neumann.  _Introduction to the modelling of marine ecosystems_.  Elsevier, 2014.
  * Feppon and Lermusiaux [2018a] F. Feppon and P. F. J. Lermusiaux.  A geometric approach to dynamical model-order reduction.  _SIAM Journal on Matrix Analysis and Applications_ , 39(1):510–538, 2018a.
  * Feppon and Lermusiaux [2018b] F. Feppon and P. F. J. Lermusiaux.  Dynamically orthogonal numerical schemes for efficient stochastic advection and Lagrangian transport.  _SIAM Review_ , 60(3):595–625, 2018b.
  * Ghil et al. [2015] M. Ghil, M. D. Chekroun, and G. Stepan.  A collection on ‘climate dynamics: multiple scales and memory effects’, 2015.
  * Gholami et al. [2019] A. Gholami, K. Keutzer, and G. Biros.  Anode: Unconditionally accurate memory-efficient gradients for neural odes.  _arXiv preprint:1902.10298_ , 2019.
  * Glass et al. [2021] D. S. Glass, X. Jin, and I. H. Riedel-Kruse.  Nonlinear delay differential equations and their application to modeling biological network motifs.  _Nature communications_ , 12(1):1–19, 2021.
  * Gouasmi et al. [2017] A. Gouasmi, E. J. Parish, and K. Duraisamy.  A priori estimation of memory effects in reduced-order models of nonlinear systems using the Mori–Zwanzig formalism.  _Proceedings of the Royal Society A_ , 473(2205):20170385, 2017.
  * Griewank [1992] A. Griewank.  Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation.  _Optimization Methods and software_ , 1(1):35–54, 1992.
  * Hairer et al. [1993] E. Hairer, S. P. NÃ¸rsett, and G. Wanner.  _Solving ordinary differential equations I_.  Springer, 1993.
  * Holmes et al. [2012] P. Holmes et al.  _Turbulence, coherent structures, dynamical systems and symmetry_.  CUP, 2012.
  * Hundsdorfer and Verwer [2013] W. Hundsdorfer and J. G. Verwer.  _Numerical solution of time-dependent advection-diffusion-reaction equations_ , volume 33.  Springer Science & Business Media, 2013.
  * Keller [2010] A. A. Keller.  Generalized delay differential equations to economic dynamics and control.  _American-Math_ , 10:278–286, 2010.
  * Koch et al. [2014] G. Koch et al.  Modeling of delays in PKPD: classical approaches and a tutorial for delay differential equations.  _Journal of pharmacokinetics and pharmacodynamics_ , 41(4):291–318, 2014.
  * Kondrashov et al. [2015] D. Kondrashov, M. D. Chekroun, and M. Ghil.  Data-driven non-Markovian closure models.  _Physica D: Nonlinear Phenomena_ , 297:33–55, 2015.
  * Kuang [1993] Y. Kuang.  _Delay differential equations: with applications in population dynamics_.  AP, 1993.
  * Kulkarni et al. [2020] C. S. Kulkarni, A. Gupta, and P. F. Lermusiaux.  Sparse regression and adaptive feature generation for the discovery of dynamical systems.  In _International Conference on Dynamic Data Driven Application Systems_ , pages 208–216. Springer, 2020.
  * Kunisch [1982] K. Kunisch.  Approximation schemes for the linear-quadratic optimal control problem associated with delay equations.  _SIAM J. on Control and Optimization_ , 20(4):506–540, 1982.
  * Kutz et al. [2016] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor.  _Dynamic mode decomposition: data-driven modeling of complex systems_.  SIAM, 2016.
  * Laizet et al. [2015] S. Laizet, J. Nedić, and C. Vassilicos.  Influence of the spatial resolution on fine-scale features in DNS of turbulence generated by a single square grid.  _International Journal of Computational Fluid Dynamics_ , 29(3-5):286–302, 2015.
  * Lermusiaux [2007] P. F. J. Lermusiaux.  Adaptive modeling, adaptive data assimilation and adaptive sampling.  _Physica D: Nonlinear Phenomena_ , 230(1):172–196, 2007.
  * Lermusiaux et al. [2004] P. F. J. Lermusiaux et al.  Adaptive coupled physical and biogeochemical ocean predictions: A conceptual basis.  In _Computational Science - ICCS 2004_ , volume 3038 of _Lecture Notes in Computer Science_ , pages 685–692. Springer Berlin Heidelberg, 2004.  ISBN 978-3-540-22116-6.
  * Lermusiaux et al. [2006] P. F. J. Lermusiaux et al.  Progress and prospects of U.S. data assimilation in ocean research.  _Oceanography_ , 19(1):172–183, 2006.  doi: 10.5670/oceanog.2006.102.
  * Lermusiaux et al. [2006] P. F. J. Lermusiaux et al.  Quantifying uncertainties in ocean predictions.  _Oceanography_ , 19(1):92–105, 2006.  doi: 10.5670/oceanog.2006.93.
  * Lermusiaux et al. [2011] P. F. J. Lermusiaux et al.  Multiscale physical and biological dynamics in the Philippine Archipelago: Predictions and processes.  _Oceanography_ , 24(1):70–89, 2011.
  * Lesieur et al. [2005] M. Lesieur, O. Métais, et al.  _Large-eddy simulations of turbulence_.  Cambridge U. press, 2005.
  * Li and Stinis [2019] J. Li and P. Stinis.  Mori-zwanzig reduced models for uncertainty quantification.  _Journal of Computational Dynamics_ , 6(PNNL-SA-132853), 2019.
  * Los and Blaas [2010] F. Los and M. Blaas.  Complexity, accuracy and practical applicability of different biogeochemical model versions.  _Journal of Marine Systems_ , 81(1-2):44–74, 2010\.
  * MacDonald [2008] N. MacDonald.  _Biological delay systems: linear stability theory_.  CUP, 2008.
  * Matsuya and Kanai [2015] K. Matsuya and M. Kanai.  Exact solution of a delay difference equation modeling traffic flow and their ultra-discrete limit.  _arXiv preprint:1509.07861_ , 2015.
  * Matthies and Meyer [2003] H. G. Matthies and M. Meyer.  Nonlinear galerkin methods for the model reduction of nonlinear dynamical systems.  _Computers & Structures_, 81(12):1277–1286, 2003.
  * Maulik and San [2018] R. Maulik and O. San.  Explicit and implicit les closures for burgers turbulence.  _Journal of Computational and Applied Mathematics_ , 327:12–40, 2018.
  * Maulik et al. [2020] R. Maulik et al.  Time-series learning of latent-space dynamics for reduced-order model closure.  _Physica D: Nonlinear Phenomena_ , 405:132368, 2020.
  * May [2019] R. M. May.  _Stability and complexity in model ecosystems_ , volume 1.  Princeton U. press, 2019.
  * McWilliams [2016] J. C. McWilliams.  Submesoscale currents in the ocean.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 472(2189):20160117, 2016.
  * McWilliams [2017] J. C. McWilliams.  Submesoscale surface fronts and filaments: secondary circulation, buoyancy flux, and frontogenesis.  _Journal of Fluid Mechanics_ , 823:391, 2017.
  * Mizutani and Dreyfus [2001] E. Mizutani and S. E. Dreyfus.  On complexity analysis of supervised MLP-learning for algorithmic comparisons.  In _IJCNN’01. International Joint Conference on Neural Networks. Proceedings (Cat. No. 01CH37222)_ , volume 1, pages 347–352. IEEE, 2001.
  * Newberger et al. [2003] P. A. Newberger, J. S. Allen, and Y. H. Spitz.  Analysis and comparison of three ecosystem models.  _Journal of Geophysical Research: Oceans (1978–2012)_ , 108(C3), 2003.
  * Nowak [2006] M. A. Nowak.  _Evolutionary dynamics: exploring the equations of life_.  Harvard U. press, 2006.
  * Otto et al. [2019] A. Otto, W. Just, and G. Radons.  Nonlinear dynamics of delay systems: An overview.  _Philosophical Transactions of the Royal Society A_ , 377(2153):20180389, 2019.
  * Pan and Duraisamy [2018] S. Pan and K. Duraisamy.  Data-driven discovery of closure models.  _SIAM Journal on Applied Dynamical Systems_ , 17(4):2381–2413, 2018.
  * Paszke et al. [2019] A. Paszke et al.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in neural information processing systems_ , pages 8026–8037, 2019.
  * Pawar et al. [2020] S. Pawar, S. E. Ahmed, O. San, and A. Rasheed.  Data-driven recovery of hidden physics in reduced order modeling of fluid flows.  _Physics of Fluids_ , 32(3):036602, 2020.
  * Portwood et al. [2019] G. D. Portwood et al.  Turbulence forecasting via neural ODE.  _arXiv preprint:1911.05180_ , 2019.
  * [69] C. Rackauckas et al.  SciML Scientific Machine Learning Software.  https://sciml.ai/.
  * Rackauckas et al. [2019] C. Rackauckas et al.  Diffeqflux.jl - A julia library for neural differential equations.  _arXiv preprint:1902.02376_ , 2019.
  * Rackauckas et al. [2020] C. Rackauckas et al.  Universal differential equations for scientific machine learning.  _arXiv preprint:2001.04385_ , 2020.
  * Rasmussen et al. [2003] H. Rasmussen, G. Wake, and J. Donaldson.  Analysis of a class of distributed delay logistic differential equations.  _Mathematical and computer modelling_ , 38(1-2):123–132, 2003.
  * Richard [2003] J.-P. Richard.  Time-delay systems: an overview of some recent advances and open problems.  _automatica_ , 39(10):1667–1694, 2003.
  * Robinson and Lermusiaux [2002] A. R. Robinson and P. F. J. Lermusiaux.  Data assimilation for modeling and predicting coupled physical–biological interactions in the sea.  In _Biological-Physical Interactions in the Sea_ , volume 12 of _The Sea_ , chapter 12, pages 475–536. John Wiley and Sons, New York, 2002\.
  * Robinson et al. [1998] A. R. Robinson, P. F. J. Lermusiaux, and N. Q. Sloan III.  Data assimilation.  In _The Global Coastal Ocean-Processes and Methods_ , volume 10 of _The Sea_ , chapter 20, pages 541–594. John Wiley and Sons, New York, 1998\.
  * Robinson et al. [2002] A. R. Robinson et al.  Predictive skill, predictive capability and predictability in ocean forecasting.  In _Proceedings of ”The OCEANS 2002 MTS/IEEE” conference_ , pages 787–794. Holland Publications, September 2002.
  * Robinson et al. [2017] N. M. Robinson et al.  A systematic review of marine-based species distribution models (SDMs) with recommendations for best practice.  _Frontiers in Marine Science_ , 4:421, 2017.
  * Roussel [1996] M. R. Roussel.  The use of delay differential equations in chemical kinetics.  _The journal of physical chemistry_ , 100(20):8323–8330, 1996.
  * Rubanova et al. [2019] Y. Rubanova, R. T. Chen, and D. Duvenaud.  Latent odes for irregularly-sampled time series.  _arXiv preprint:1907.03907_ , 2019.
  * San and Maulik [2018] O. San and R. Maulik.  Neural network closures for nonlinear model order reduction.  _Advances in Computational Mathematics_ , 44(6):1717–1750, 2018.
  * Sapsis and Lermusiaux [2012] T. P. Sapsis and P. F. J. Lermusiaux.  Dynamical criteria for the evolution of the stochastic dimensionality in flows with uncertainty.  _Physica D: Nonlinear Phenomena_ , 241(1):60–76, 2012.
  * Schneider et al. [2017] T. Schneider et al.  Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations.  _Geophy. Res. L._ , 44(24):12–396, 2017.
  * Smith [2011] H. L. Smith.  _An introduction to delay differential equations with applications to the life sciences_ , volume 57.  Springer New York, 2011.
  * Stinis [2015] P. Stinis.  Renormalized Mori–Zwanzig-reduced models for systems without scale separation.  _Proceedings of the Royal Society A_ , 471(2176):20140446, 2015.
  * Takens [1981] F. Takens.  Detecting strange attractors in turbulence.  In _Dynamical systems and turbulence, Warwick 1980_ , pages 366–381. Springer, 1981.
  * Tokuda et al. [2019] I. T. Tokuda, O. E. Akman, and J. C. Locke.  Reducing the complexity of mathematical models for the plant circadian clock by distributed delays.  _J. of theoretical biology_ , 463:155–166, 2019.
  * Ueckermann and Lermusiaux [2010] M. P. Ueckermann and P. F. J. Lermusiaux.  High order schemes for 2D unsteady biogeochemical ocean models.  _Ocean Dynamics_ , 60(6):1415–1445, Dec. 2010\.
  * Wan et al. [2018] Z. Y. Wan, P. Vlachas, P. Koumoutsakos, and T. Sapsis.  Data-assisted reduced-order modeling of extreme events in complex dynamical systems.  _PloS one_ , 13(5), 2018.
  * Wang et al. [2020] Q. Wang, N. Ripamonti, and J. S. Hesthaven.  Recurrent neural network closure of parametric POD-galerkin reduced-order models based on the Mori-Zwanzig formalism.  _Journal of Computational Physics_ , page 109402, 2020.
  * Wang et al. [2012] Z. Wang, I. Akhtar, J. Borggaard, and T. Iliescu.  Proper orthogonal decomposition closure models for turbulent flows: a numerical comparison.  _CMAME_ , 237:10–26, 2012.
  * Ward et al. [2013] B. A. Ward et al.  When is a biogeochemical model too complex? objective model reduction and selection for North Atlantic time-series sites.  _Progress in Oceanography_ , 116:49–65, 2013.
  * Whitney [1936] H. Whitney.  Differentiable manifolds.  _Annals of Mathematics_ , pages 645–680, 1936.
  * Wunsch [1996] C. Wunsch.  _The ocean circulation inverse problem_.  Cambridge University Press, 1996.
  * Yang et al. [2020] Y. Yang, M. Aziz Bhouri, and P. Perdikaris.  Bayesian differential programming for robust systems identification under uncertainty.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 476(2243):20200290, 2020.
  * Yeung et al. [2018] P. Yeung et al.  Effects of finite spatial and temporal resolution in direct numerical simulations of incompressible isotropic turbulence.  _Physical Review Fluids_ , 3(6):064603, 2018.
  * Yuval et al. [2020] J. Yuval, C. N. Hill, and P. A. O’Gorman.  Use of neural networks for stable, accurate and physically consistent parameterization of subgrid atmospheric processes with good performance at reduced precision.  _arXiv preprint:2010.09947_ , 2020.