# The Phosphine Controversy:
Is it Phosphine? Is there life on Venus?

Priya Hasan

###### Abstract

On 14th September 2020, the Royal Astronomical Society made an official
statement coupled with a webminar on the discovery of phosphine on Venus.
Single-line millimetre-waveband spectral detections of phosphine (with a
signal-to-noise ratio of $\approx$ 15$\sigma$) from the JCMT and ALMA
telescopes indicated a phosphine abundance of 20 ppb111ppb: parts per billion,
1000 times more than that on the Earth. Phosphine is an important biomarker
and immediate speculation in the media about indicators of life being found on
Venus followed. This article presents an analysis of the study and the results
on the observation of the spectral absorption feature of phosphine in the
clouds of Venus, thus implying as a potential biosignature. If phosphine is
produced through biotic, as opposed to abiotic pathways, the discovery could
imply a significant biomass in the Venusian atmosphere. The discovery led to a
major controversy with criticism of the analysis and results and responses to
it. The issue remains unresolved, leading to a fresh interest in the study of
Venus including ground-based observations as well as space-probes that can
answer these questions conclusively.

‘We are alive and we resonate with idea of life elsewhere, but only careful
accumulation and assessment of the evidence can tell us whether a given world
is inhabited.’
Carl Sagan

## Introduction

Phosphine gas is not very pleasant. Pure phosphine is odourless, but technical
grade samples smell like rotting fish, toxic and spontaneously flammable, with
the chemical formula of $\rm{PH_{3}}$.

On 14th September 2020, the Royal Astronomical Society made an official
statement coupled with a webminar on the discovery of phosphine on Venus. At
the heart of this announcement was the paper by [1] titled ‘Phosphine gas in
the cloud decks of Venus’. It was about the detection of phosphine in the
temperate, but hyperacidic clouds of Venus with an abundance of 20 ppb. They
also made a very detailed study of various chemical pathways of generation of
phosphine in the given abundance by lightning, volcanoes or meteorites. They
concluded, that the observed abundance of phosphine could possibly be a result
of biological processes. They also stated that a detection of statistical
significance needed to be backed up by observations to validate the spectral
features of phosphine and of the conditions in the Venusian atmosphere.

Significant work in this context was done by [2] who at made a detailed study
of all the possible chemical processes that can generate these amounts of
phosphine. [3] proposed a Venusian aerial biosphere to explain the presence of
phosphine as a possible hypothesis.

##  1 Biosignatures: Why Phosphine?

Wikipedia defines biosignatures as: A biosignature is any substance – such as
an element, isotope, or molecule – or phenomenon that provides scientific
evidence of past or present life. Measurable attributes of life include its
complex physical or chemical structures and its use of free energy and the
production of biomass and wastes. The obvious way forward is to identify
molecules present on the Earth that could work as biosignatures. An ideal bio-
signature should have its sole source from living organisms be intrinsically
strong and easily identifiable spectroscopically.

The element that comes first to the mind is oxygen. However, as shown by [4],
oxygen is not most suitable, as it creates several false positive mechanisms
on a variety of planet scenarios are possible (see Fig. 1). This implies that
detection as well as non-detection of certain molecules could lead to false
positive biosignatures222A false positive is a set of non-biological processes
that can mimic the detectable feature of a biosignature. False negative
biosignatures occur where life may be present on another planet, but potential
biosignatures are undetectable..

Figure 1: Oxygen as a biosignature: False positive mechanisms for oxygen on a
variety of planet scenarios. The main contributors to a spectrum of the
planet’s atmosphere are shown in large rectangles. The yellow and red circles
represent molecules that would confirm a false positive biosignature and and
molecules crossed confirm false positive biosignatures if not detected.
Cartoon adapted from [4] (Image Credit: Wikipedia).

In the Earth’s atmosphere, methane reacts with oxygen to get converted to
carbon-dioxide and water, very efficiently. However, the atmosphere is found
to have an abundance of 1 molecule per million and this can only be explained
by a constant source of methane production (life) that oxygen cannot keep pace
with in conversion to carbon-dioxide and water. Hence due to an imbalance of
the production and destruction of methane, we still find some of it in our
atmosphere. This is an indicator of life in an oxygen rich environment.

Biogenic methane production is the main contributor to the methane flux coming
from the surface of Earth. Methane has a photochemical sink in the atmosphere,
and hence can be detected only in the presence of biogenic methane production
that is at a rate larger than its oxidation (Fig. 2). [5] stated that the
detection of methane in the atmosphere of another planet can a viable
biosignature, especially if the host star is of $G$ or $K$ type.

Figure 2: Methane as a viable biosignature on the Earth. (Image Credit:
Wikipedia).

A similar argument was used to study methane as a biomarker on Mars. However,
the processes of generation and destruction of methane on Mars are not fully
understood (See Fig. 3). This is still controversial and a mass spectrometer
measuring the isotopic proportions of carbon-12 to carbon-14 in methane is
used to validate the bio-origin of methane [6].

Figure 3: Methane as a biosignature on Mars: Generation and Destruction (Image
Credit: Wikipedia).

Phosphine is similar to ammonia, made of a single atom of phosphorous and
three atoms of hydrogen. On the Earth, we have a phosphine abundance as 1 part
per trillion. In an oxidizing environment, phosphine will immediately react
with oxygen to form phosphorous acid or phosphoric acid.

Phosphine is exclusively associated to life and is not generated by other
natural processes involving the geology or atmosphere of our planet. Hence its
presence in the oxidizing atmosphere, indicates a constant generation, caused
by biological processes exceeding the rate of oxidation, making it a reliable
biomarker.

On the Earth, phosphine is associated with anaerobic ecosystems, and as such
it is a potential biosignature gas in anoxic exoplanets. It has also been
found on the gas giants Jupiter and Saturn. These planets have a hydrogen rich
atmosphere, and the phosphine gets generated in the high temperature and
pressure environment getting mixed well with the gas layers due to convection.
In the case of rocky planets, the surface provides a natural barrier stopping
the phosphine produced at the core from getting mixed with the atmosphere.
Hydrogen is more likely to combine with oxygen to form water, or carbon to
form methane. Hence for phosphorous to combine with Hydrogen in its very low
abundance, is very improbable. Hence the detection of phosphine on a temperate
rocky planet, is a robust indicator of life 333A high temperature could
destroy the phosphine. Hence a process of regeneration of phosphine would be
necessary for detection.. On the Earth, phosphine could be produced directly
by microbial reduction of more oxidized phosphorus species or indirectly by
microbial production of reduced phosphorous compounds, such as hypophosphite
leading to $\rm{PH_{3}}$, implying a bio-origin.

An excess of phosphine would be due to an imbalance in the production and
oxidation processes and can only be produced by life forms on rocky planets,
for example, the Earth. If phosphine is produced through biotic, as opposed to
abiotic pathways, the discovery could imply a significant biomass in the
Venusian atmosphere [7].

There are a few other mechanisms that can produce phosphine. Even if it was
made in the lower layers, it would be destroyed before it reaches layers where
it would be found. Volcanoes, thunder, etc. can produce it, but an appropriate
hydrogen pressure would be required to keep the phosphine. A document of more
than 100 pages was made to study all such mechanisms [2]. There is also strong
ultraviolet (UV) absorbtion in the Venusian atmosphere. Venus was probably
habitable in the past till the green house effect took over to make the planet
inhabitable.

Also, we should always keep in mind that life elsewhere could be very
different from life on Earth from a bio-chemical perspective, specially with
the variety of bio-diversity on the Earth.

##  2 Venus: The Evil Twin

Venus has been called the evil twin of the Earth since its size and mass is
comparable. ‘Heaven and Hell’ were the words used by Carl Sagan to describe
the twins: Earth and Venus. The atmosphere of Venus is composed of 96.5%
$\rm{CO_{2}}$ and 3.5% $\rm{N_{2}}$ and its present surface conditions are
very different from that on Earth.

Life of Venus has been discussed in the past. [8, 9, 10] proposed that the
conditions between the lower and middle atmosphere were conducive to
(terrestrial) biology, while higher altitudes could possibly host
microorganisms. Life in the Venusian cloud layers was also considered by [11]
due to the presence of conducive chemical and physical conditions, like the
presence of sulfur compounds, carbon dioxide ($\rm{CO_{2}}$), and water, and
moderate temperatures ($0-60^{0}$ C) and pressures ($\approx 0.4-2$ atm).

In an hypothesis article [12], proposed that the lower cloud layer of Venus
($47.5-50.5$ km) is an region favorable to microbial life, as it has moderate
temperatures and pressures ($\approx$ $60^{0}$ C and 1 atm), and the presence
of micron-sized sulfuric acid aerosols (Fig. 4).

Figure 4: The potential for microorganisms to survive in Venus’ lower clouds
and contribute to the observed bulk spectra in a hypothesis paper [12] (Image
credit: Limaye et al, doi: 10.1089/ast.2017.1783.)

The formidable Venusian surface has a blistering temperature of 4700 C and a
pressure 92 times higher than that at the surface of Earth. This is similar to
the pressure found in the depth of 900 m of Earth’s oceans. This is an extreme
hostile environment for life, as we know it. The Russian Venera program was a
series of space probes between 1961 and 1984 to study Venus. Ten probes
successfully landed on the surface of the planet, including the two Vega
program and Venera-Halley probes, while thirteen probes successfully entered
the Venusian atmosphere. However, due to the hostile surface conditions, the
probes could only survive for short periods ranging from 23 minutes to two
hours.

The surface of Venus is very hot ($\approx 400^{0}$ C) , but at a height of
45–60 km the temperature is about $30^{0}$ C and the pressure is 1 atm. This
is possible region where the physical conditions are suitable for life, as we
know it. It is believed that there could have been similar conditions (cooler
and wetter) on the Earth and Venus at an earlier time and as conditions on
Venus turned hostile, life could have migrated to this more temperate zone in
the Venusian atmosphere. However, the chemical composition of the atmosphere
is very different from that on the Earth, with clouds or droplets of sulphuric
acid. Hence, life if it exists in these conditions, would be very different
from we know it to be.

##  3 Data: JCMT and ALMA

Phosphine as a probable biomarker hence seemed to be good candidate for
detection on Venus. Phosphine is made up of a single atom of phosphorous and
three atoms of hydrogen. The atoms in the molecule can make transitions
between different rotation energy levels corresponding to certain
characteristic frequencies or wavelengths by absorbing solar radiation. One
such transition from the $J=1-0$ levels is at a frequency of 266944.51 MHz =
1.123 mm. Observations of this spectral line can confirm the presence as well
as the abundance of phosphine and the temperature and pressure of the region.

The James Clerk Maxwell Telescope (JCMT) is a submillimetre-wavelength radio
telescope at Mauna Kea Observatory in Hawaii, USA. The telescope is near the
summit of Mauna Kea at 4,092 m. Its primary mirror has a size of 15 m. It is
the largest single-dish telescope that operates in submillimetre wavelengths
of the electromagnetic spectrum (far-infrared to microwave) (Fig. 5).

Figure 5: Model of James Clerk Maxwell Telescope (JCMT). Photographed at the
Royal Astronomical Society’s National Annual Meeting 2009. (Image Credit:
Wikipedia).

The Fig. 6 shows the spectrum obtained from the JCMT. The spectral
ripples444Spectral ripples are instrumental artefacts in the data that become
apparent when observing an object as bright as Venus. had to be fitted and
removed to obtain the spectrum shown in the left panel with the residual line
present inside velocity ranges of $v=8$ km/s (solid, black) and $v=2$ km/s
(dashed, orange). The data was binned for clarity into histograms on the x
axis; representative $1\sigma$ error bars are also shown. The right panel
shows the adopted mid-range solution with $v=5$ km/s (histogram), overlaid
with our model for 20 ppb abundance by volume. The solid red curve shows this
model after processing with the same spectral fitting as used for the data.

Figure 6: Spectrum of PH3(1-0) with JCMT [1]

The Atacama Large Millimeter/submillimeter Array (ALMA) is an astronomical
interferometer of 66 radio telescopes in the Atacama Desert of northern Chile,
which observe electromagnetic radiation at millimeter and submillimeter
wavelengths (Fig. 7). The array has been constructed at the altitude of 5,000
m at the Chajnantor plateau – near the Llano de Chajnantor Observatory and the
Atacama Pathfinder Experiment. The site is perfect due to its high elevation
and low humidity. The array has much higher sensitivity and higher resolution
than earlier submillimeter telescopes such as the single-dish JCMT or existing
interferometer networks such as the Submillimeter Array or the Institut de
Radio Astronomie Millimétrique (IRAM) Plateau de Bure facility.

Figure 7: Atacama Large Millimeter/submillimeter Array (ALMA)(Image Credit: P.
Horálek/ESO)

Figure 8 shows the PH3 ($1-0$) spectrum of the whole planet, with $1\sigma$
errors of $0.11\times 10^{-4}$ per 1.1 km/s spectral bin in the left panel.
The right panel shows the spectra of the polar (histogram in black), mid-
latitude (in blue) and equatorial (in red) zones on the planet.

Figure 8: Spectrum of PH3 $(1-0)$ with ALMA [1]

A comparison of the production and destruction of phosphine and a better
understanding of the temperature and pressure on Venus is an integral part of
this study. The life time of phosphine in the layers at a height larger than
80 km is about $10^{3}$ s as it is destroyed by UV radiation. Near the base,
due to collisions, the molecules are destroyed in $10^{8}$ s. In mid-regions,
there is no data about the lifetimes. It is expected to be about $10^{3}$
years. [3] proposed a hypothetical lifecycle in the Venusian atmosphere, the
arial biosphere. Life could reside in the droplets, decicated spores that stay
on higher levels of the atmosphere and when they reach the temperate zone,
they get metabollically active, reproduce. They drift upwards and downwards as
they evolve.

Life as we know it, cannot survive in the presence of sulphric acid and the
conditions in the Venusian atmosphere. However, possible life-forms can have
an unknown photochemistry and hence the question is still open.

##  4 The Controversy

There are two big questions in this study. Firstly, does the signal detected
correspond to phosphine and no other molecule? Secondly, if it is phosphine,
is it caused by life?

To answer the first questions, the data has been multiply analysed by various
research groups to verify the presence of phosphine. The next step would be to
reobserve the signal repeatedly to see how it is distributed on Venus, does it
change between day and night, or seasons, or different regions? IF it is truly
related to life, then we would expect some kind of variation. And IF it truly
is life, it won’t produce only one such molecule, there would be many more
such molecules being produced by the complete ecosystem.

The data obtained was very noisy and there were various algorithms were used
to reduce the data and extract the signal without causing an artificial
signal. Several observations were used, 18 months apart, to ensure that the
signal was not spurious and applying different kinds of algorithms to identify
the spectral line in multiple ways in repeated efforts by teams taking in
inputs from experts and referees. It took almost 3 years from the first
detection on JCMT to publication.

The authors of the discovery paper stated clearly that ‘Even if confirmed, we
emphasize that the detection of PH3 is not robust evidence for life, only for
anomalous and unexplained chemistry.’ [1].

Shortly after the announcement, the organizing committee of the International
Astronomical Union (IAU) Commission F3 on Astrobiology released a statement
criticising Greaves’ team for the press coverage of their discovery. ‘It is an
ethical duty for any scientist to communicate with the media and the public
with great scientific rigour and to be careful not to overstate any
interpretation which will be irretrievably picked up by the press,’ adding
that the commission ‘would like to remind the relevant researchers that we
need to understand how the press and the media behave before communicating
with them’.

The IAU statement was condemned by many, including the commission’s own
members, due to which the statement was retracted by the IAU executive. The
IAU Executive stated that the organizing committee of Commission F3 had ‘been
contacted to retract their statement and to contact the scientific team with
an apology’.

Shortly after that, [13] contested the paper by [1] and stated that the
observed PH3 feature with JCMT can be fully explained employing plausible
mesospheric SO2 abundances ( 100 ppbv as per the SO2 profile) and the
identification of PH3 in the ALMA data should be considered invalid due to
severe baseline calibration issues. The team ended its abstract with the
‘suggestion’ that Greaves’ team retract its original paper – seen by some as
unduly aggressive and hence an apology was made by Villanueva’s team. ‘We
agree that the sentence calling for retraction was inappropriate and we
apologise for harm caused to the Greaves et al. team.’

There were a few more criticisms of the paper by some groups like [14, 15,
16].

[17, 18] responded to the paper and reanalysed the data to recover PH3 in
Venus’ atmosphere with ALMA ($\approx 5\sigma$ confidence). They stated that
the ALMA data are reconcilable with the JCMT detection ($\approx$20 ppb) if
there is order-of-magnitude temporal variation and more advanced processing of
the JCMT data is underway to check methods. They concluded that both ALMA and
JCMT were working at the limit of observatory capabilities and hence new
spectra should be obtained as spectral ripples could potentially reducing the
significance of real narrow spectral features.

To add, a recent paper by [19] re-examined data obtained from the Pioneer-
Venus Large Probe Neutral Mass Spectrometer (LNMS) from NASA’s Pioneer Venus
spacecraft to search for evidence of phosphorus and confirmed the detection of
phosphine on Venus.

In the words of Lyman Beecher, ‘No great advance has ever been made in
science, politics, or religion, without controversy’. We shall know.

##  5 The Future

As of now, the phosphine paper is the only set of results that has undergone
peer review. Several papers have been submitted with a reanalysis of the data
and if the findings are still under dispute, new data at different frequencies
are required.

Phosphine is difficult to detect from the ground but NASA’s airborne
Stratospheric Observatory For Infrared Astronomy telescope, which flies at an
altitude of over 13.7 km on board a modified Boeing 747SP, could confirm or
deny the finding.

Observations are also being planned in the Infrared using the NASA Infrared
Telescope Facility (NASA IRTF) which is a 3-meter telescope in Hawaiii. Most
telescopes are designed to look at faint sources. Venus is very bright and
hence methods need to be planned to adapt to this bright source otherwise it
saturates our detectors. JWST can look for such signals in far away planets.
However, Venus is too bright for its detectors.

There will also be new missions for making in-situ measurements of the
Venusian atmosphere to help build more realistic models. The orbiter and
atmospheric balloon mission Shukrayaan-1 by ISRO is under development and
planned under development and planned for 2023. The orbiter and lander
Venera-D by Roscosmos is under development and planned for 2026. NASA has
proposed a secondary payload VAMP for the Venera-D lander.

Small-scaled Venus mission could be launched in the near future to confirm the
presence of phosphine and measure its vertical distribution in the atmosphere.
A whole series of missions could look for signs of life and even life itself.
A golden era of Venus exploration lies ahead.

The hypothesis of life in the clouds of Venus will gain scientific validity
only after we have systematically ruled out all other chemical and geological
processes that can explain the presence of phosphine on Venus and widened our
explorations of this hostile, yet inviting twin.

##  6 Conclusion

The search of extraterrestial life has always been at the heart of human quest
of Life, Universe and Everything. We have reached a stage where we are
exploring planets of other Solar Systems and looking for habitability and
signatures of life on them. Our Solar System is the immediate neighbourhood
where we have started these explorations looking for indicators of water,
habitability and life. The quest is not an easy one, however it is of great
meaning and importance. The exploration of Venus is one such attempt. Many
more such attempts will be planned in the future and this can be a small step
in the right direction in us, the universe, trying to find answers about
itself. Afterall, as Carl Sagan once said, ‘extraordinary claims require
extraordinary evidence’ (ECREE). That’s what we are looking for.

## References

  * [1] J S Greaves, A M Richards, W Bains, et al. , Nat Astron, 2020, https://doi.org/10.1038/s41550-020-1174-4.
  * [2] W Bains, J J Petkowski, S Seager, S Ranjan, C Sousa-Silva, P B Rimmer, Z Zhan, J S Greaves, A M Richards, 2020, arXiv:2009.06499[astro-ph.EP], 2020.
  * [3] S Seager, J J Petkowski, P Gao, W Bains, N Bryan, S Ranjan, J Greaves, Astrobiology, published ahead of print https://www.liebertpub.com/doi/full/10.1089/ast.2020.2244, 2020.
  * [4] V S Meadows, C T Reinhard, G N Arney, M N Parenteau, E W Schwieterman et al., 18 (6): 630–662, 2018.
  * [5] G N Arney, The Astrophysical Journal, 873 (1): L7, 2019.
  * [6] K Zahnle, R S Freedman, D C Catling, Icarus, 212 (2): 493–503, 2011.
  * [7] M Lingam, A Loeb, arXiv:2009.07835v2 [astro-ph.EP], 2020.
  * [8] H. A Morowitz, C Sagan, Nature, 215:1259–1260, 1967.
  * [9] D H Grinspoon, M A Bullock, In: Exploring Venus as a Terrestrial Planet, edited by L W Esposito, E R Stafan, and T E Cravens, American Geophysical Union, pp. 191–206, 2007.
  * [10] Cockell CS, Life on Venus, Planet Space Sci, 47:1487–1501, 1999.
  * [11] D Schulze-Makuch, D H Grinspoon, O Abbas, L N Irwin, M.A.Bullock M.A., 4:11–18, 2004.
  * [12] S Limaye, R Mogul, D J Smith, A H Ansari, G P Słowik, P Vaishampayan, Astrobiology, 1181-1198, 2018.
  * [13] G Villanueva, M Cordiner, P Irwin, I de Pater, et al., arXiv:2010.14305[astro-ph.EP], 2020.
  * [14] I A G Snellen, L Guzman-Ramirez, M R Hogerheijde, A P S Hygate, F F S van der Tak, 644, Astronomy & Astrophysics, L2, 2020.
  * [15] T Encrenaz, T K Greathouse, E Marcq, T Widemann, B Bézard, T Fouchet, R Giles, H Sagawa, J Greaves, C Sousa-Silva, Astronomy and Astrophysics, 643, L5, 2020.
  * [16] M A Thompson, Monthly Notices of the Royal Astronomical Society: Letters, 501, L18-L22, 2020.
  * [17] J S Greaves, A M Richards, W Bains. et al., arXiv:2011.08176 [astro-ph.EP], 2020.
  * [18] C Sousa-Silva, S Seager, S Ranjan, J J Petkowski, Z Zhan, R Hu, W Bains, Astrobiology, 235-268, 2020.
  * [19] R Mogul, S S Limaye, M. J. Way, J A Cordova, arXiv: 2009.12758 [astro-ph.E], 2020.


#  Validation and parameterization of a novel physics-constrained neural
dynamics model applied to turbulent fluid flow

Varun Shankar  Carnegie Mellon University  Gavin D. Portwood  Los Alamos
National Laboratory  Lawrence Livermore National Laboratory  Arvind T. Mohan
Los Alamos National Laboratory  Peetak P. Mitra  University of Massachusetts
Amherst  Palo Alto Research Center  Dilip Krishnamurthy  Carnegie Mellon
University  Christopher Rackauckas  Massachusetts Institute of Technology
Lucas A. Wilson  Dell Technologies  David P. Schmidt  University of
Massachusetts Amherst  Venkatasubramanian Viswanathan  Carnegie Mellon
University<EMAIL_ADDRESS>
###### Abstract

In fluid physics, data-driven models to enhance or accelerate solution methods
are becoming increasingly popular for many application domains, such as
alternatives to turbulence closures, system surrogates, or for new physics
discovery. In the context of reduced order models of high-dimensional time-
dependent fluid systems, machine learning methods grant the benefit of
automated learning from data, but the burden of a model lies on its reduced-
order representation of both the fluid state and physical dynamics. In this
work, we build a physics-constrained, data-driven reduced order model for the
Navier-Stokes equations to approximate spatio-temporal turbulent fluid
dynamics. The model design choices mimic numerical and physical constraints
by, for example, implicitly enforcing the incompressibility constraint and
utilizing continuous Neural Ordinary Differential Equations for tracking the
evolution of the differential equation. We demonstrate this technique on
three-dimensional, moderate Reynolds number turbulent fluid flow. In assessing
the statistical quality and characteristics of the machine-learned model
through rigorous diagnostic tests, we find that our model is capable of
reconstructing the dynamics of the flow over large integral timescales,
favoring accuracy at the larger length scales. More significantly,
comprehensive diagnostics suggest that physically-interpretable model
parameters, corresponding to the representations of the fluid state and
dynamics, have attributable and quantifiable impact on the quality of the
model predictions and computational complexity.

††preprint: APS/123-QED

##  I Introduction

Model reduction is commonly employed by scientists and engineers to obtain
solutions to complex turbulent fluid systems which may otherwise be
prohibitively complex to simulate by their physical governing equations.
Indeed various reduced order modeling approaches have been proposed for
turbulent fluid systems which have the potential to aid in applications such
as simulation-based design optimisation or statistical inverse modeling.
Recently, reduced order modeling of turbulent fluid systems by various deep
learning techniques has shown strong promise in accurately modeling various
aspects of fluid motion while being computationally efficient compared to
their first-principles counterparts. While modeling approaches which
demonstrate favorable accuracy with respect to physical diagnostics are
encouraging, we believe more rigorous and interpretable parameterization of
approximation error is essential for model certification in predictive
science.

The understanding of errors in machine-learned reduced order models is
complicated by the lack of formal approximation theory for many classes of
data driven models, in particular our interest in deep feed-forward and
convolutional neural networks. Combined with the tendency for surrogate models
to handle disparate phenomena with so-called deep “end-to-end” architectures,
it is frequently unclear how to attribute model error in deep neural networks
with respect to model architectures.

However, various efforts have improved the state of model error
interpretability. Firstly, there are various approaches in which physical
constraints can be imposed on ML models [1]. For instance, the development of
loss functions which are informed by governing equations in various
manifestations of physics-informed neural networks (PINNs) have demonstrated
effective surrogate models which approximate physical solutions [2]. Network
architectures which implicitly preserve hard physical constraints have shown
promise in improving model accuracy by enforcing strictly physically-
realizable model outputs in some respects [3, 4].

Beyond enforcing physical constraints, various approaches have been suggested
using convolutional neural networks (CNNs) [5, 3, 6, 7], recurrent neural
networks (RNNs) [1] or long short-term memory (LSTM) [6, 8] architectures.
These models typically encode spatial information from the field using CNNs
and use autoregressive predictions or RNN modules to address the temporal
dependency. Beyond reduced-order modeling, other machine learning
methodologies have been proposed with physics-informed components, such as
methods to relax grid-resolution requirements with learned interpolations or
mesh coarsening corrections for the simulation [9], accelerate intermediate
solver steps [10], and improve predictive capacity with hybrid CFD and ML-
correcting algorithms [11, 7, 12].

However, the adoption of various manifestations of differential equation-based
deep learning models, such as NeuralODEs (NODEs), have effectively
demonstrated the shift of modeling burden to numerically-formal time
integrators [13]. NeuralODEs have previously been used for forecasting limited
turbulence prognostics and system states [14, 15]. More importantly, this
framework establishes a principled way to trade-off accuracy and depth, and
offers a step towards bridging the fields of numerical methods and machine
learning.

Leveraging NeuralODE and PINN design concepts in surrogate models, deep
learning architectures may be developed to mimic formal numerical solution
techniques such that neural network subcomponents are meaningful and
intuitive. Therefore, we utilize intuitive design concepts in the development
of a reduced-order model for three dimensional solutions to the incompressible
Navier-Stokes equations as detailed in section 2 of this manuscript. In
addition to enforcing hard physical constraints, we carefully design our deep
network architecture to explicitly treat two discrete modeling demands in
distinct subnetworks: reduction of physical states and reduction of system
dynamics. By parameterizing the architecture the two subnetworks, we
demonstrate an effective modeling methodology with physically-relevant
diagnostics. More importantly, we investigate the behaviour of model error as
a function of individual subnetwork complexity as a means to systematically
demonstrate trade-offs between modeling error and computational complexity.

##  II Problem statement

We consider incompressible homogeneous isotropic turbulence (HIT) as a
fundamental model flow that features sufficient complexity to evaluate our
proposed modeling methodology introduced above. In addition to exhibiting
strong nonlinearity in the parameter regimes we consider, the multiscale
structure of HIT makes the reduced-order representation of flow states a non-
trivial task. Such modeling challenges are compounded by complexity in the
governing equation, which induce non-local phenomena due to the
incompressibility constraint. The equations of motion we consider are given by
the momentum equations and mass balance:

$\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u}=\nu\nabla^{2}\mathbf{u}-\frac{1}{\rho}\nabla
p+\mathbf{f};\quad\nabla\cdot\mathbf{u}=0$ (1)

where $\mathbf{u}$ is the velocity vector, defined $\\{u,v,w\\}$, $p$ is the
pressure, $\rho$ is the density, $\nu$ is kinematic viscosity and $\mathbf{f}$
is a dynamic large-scale forcing function which enforces statistical
stationary in the turbulent kinetic energy [16], defined

$e\equiv\frac{1}{2}(\langle u^{2}\rangle+\langle v^{2}\rangle+\langle
w^{2}\rangle)\quad.$ (2)

We design our flow system with spatio-temporal homogeneity, by enforcing
statistically stationary kinetic energy and periodic boundary conditions, to
remove correlations with space or time which have the potential to be
trivially memorized by our neural networks.

###  II.1 Modeling objectives

Due to the combined multi-scale and non-linear behaviour of the Navier-Stokes
equations, discretized direct numerical solutions to (1) are subject to
unfavorable scaling of computational complexity with respect to the Reynolds
number, where the number of grid points can be shown to scale with $Re^{9/4}$
[17]. To improve computation tractability, much work has been done in reduced
order modeling of dynamical systems. Proper orthogonal decomposition, often
used with fluid flows, takes a statistical approach, projecting the data onto
a linear subspace comprised of dominant modes [18]. The idea of a latent space
representation is common in many fields, including computer vision.
Convolutional autoencoders (ConvAE) have been used for tasks like image
compression [19, 20] or feature extraction [21], both of which closely
parallel our goal of a rich reduced representation of the turbulent flow
field. Data-driven ROMs are gaining traction in fluid applications as well
[22, 23]. Our approach takes advantage of this technique to reduce
instantaneous velocity snapshots into a much smaller latent space. The
dynamics of the system can then be performed on the latent representation,
easing computation. Using a latent space in conjunction with NeuralODEs has
been previously used successfully for modeling many dynamical systems [24].

The aim of the surrogate model is to track the spatiotemporal dynamics of the
system in order to obtain approximate solutions. Given an initial condition
for the velocity field, the model should estimate the solution to the initial
value problem (IVP) by forecasting the velocity field in time over a finite
horizon. Modeling the dynamics of the spatially discretized velocity field as
a differential equation is one method of achieving temporal evolution. This
construction enables use of the NODE framework, which models the dynamics of a
state variable through a neural network. The network assumes the
responsibility of learning the differential dynamics of the state, but the IVP
can be solved with any number of formalized ODE integrators.

The challenge with combining ML models with formal numerical methods is the
computational cost of backpropagation. As an alternative, adjoint-state
methods has been a well-established numerical method for computing gradients
of some functional. Chen et al. proposed using this method for efficiently
calculating the gradients of a parameterized differential equation by solving
the adjoint problem [13]. This dual form, which is well-defined for many
classes of differential equations, can be solved using traditional numerical
methods as in the forward computation. The resulting parameter sensitivities
can then be used to optimize the cost function, in this case, the
reconstruction loss. In the context of deep-learning, the adjoint approach
allows for a myriad of new architectures that permit differential equation
solves to be easily integrated within the network structure, effectively
instantiating a continuous time, infinite-depth network. Beyond the memory-
efficient adjoint-optimization, NODEs permit arbitrary sampling in time,
flexibility in the ODE integrator, among other benefits. The ability to adapt
neural networks to include differential equations has far-reaching
implications for modeling physical phenomena, which are often strictly
described by ordinary or partial differential equations.

In order to incorporate the spatial information necessary to accurately
approximate the true function, we use a convolutional neural network (CNN) to
learn systems dynamics within the NODE model. This provides two benefits: (1)
the memory costs of the network are drastically reduced as the number of
parameters is much smaller than a fully-connected network. This is due to (2)
the network uses a shared local filter in each layer that limits the message
passing to a constrained receptive field around each output point. In a
finite-layer CNN, the output at a location is determined from input values in
a local neighborhood. The size of this neighborhood, or receptive field, is
influenced by the number of layers in the network and the kernel sizes of the
CNN layers.

The convolutional NODE formulation is a form of inductive, or learning bias
imposed on the model. Although the neural network allows for much flexibility
to be fit to the true data, it is still limited by the formulation, which
enforces that the derivative of the velocity field with respect to time is a
function of spatially local quantities, consistent with our knowledge of the
physical system. Inductive biases are especially advantageous for
generalization beyond the training set, as they encompass the assumptions a
researcher makes about unseen data. In scientific machine learning
applications, physics-based learning biases are a natural way to restrict
models to prevent overfitting and produce outputs that more closely align with
what is known about the system.

##  III Machine learning approach

###  III.1 Neural network architectures

The encoder, denoted by $\mathcal{E}$, is a trainable network that transforms
the input velocity field $\mathbf{u}_{0}$ to a lower dimensional latent space
$\mathbf{z}_{0}$, given by

$\mathcal{E}(\mathbf{u}_{0})=\mathbf{z}_{0}\quad.$ (3)

Compression is achieved with the use of strided convolutional layers in the
encoder to reduce the spatial resolution in the latent space. These strided
layers use the same convolutional spatial filtering technique to transform the
input but coarsen the data by subsampling every other point (for a stride 2
layer). The size of the latent space varies in spatial resolution,
corresponding to the number of strided, or coarsening layers, and latent
channels, greater than the 3 input velocity channels. The encoder/decoder
subnetworks are parameterized by the compression ratio, $z$, defined as the
size of the original space relative to the latent space:

$z=\frac{3\cdot 64^{3}}{C_{z}N_{z}^{3}}\quad,$ (4)

where $N_{z}$ is the number of points in each latent spatial dimension and
$C_{z}$ is the number of latent channels. For reference with respect to
physical quantities, the receptive field size of the encoder outputs are
therefore equivalent to the grid spacing at $z=1$ and the domain size at
$z=64$.

The initial condition in the latent space $\mathbf{z}_{0}$ is used to solve
the ODE parameterized by the network weights. The system is defined by

$\mathrm{\frac{d\mathbf{z}}{dt}=g_{\theta}(\mathbf{z}),\quad\mathbf{z}(0)=\mathbf{z}_{0}\quad},$
(5)

where $g_{\theta}$ represents the dynamical convolutional subnetwork that
learns the continuous dynamics of the latent state. The dynamical subnetwork
consists of 3 convolutional layers and is parameterized by the kernel size of
the layers. The kernel size of the NODE has a direct impact on the receptive
field of the network and thus the sparsity of the ODE system, and in effect
gives an indication of the length scale associated with the dynamics. Larger
kernel sizes can learn higher order approximations of the spatial information,
but this comes at the disadvantage of computational costs scaling as
$O(k^{3})$, where $k$ is the kernel size, for this 3-dimensional application.
Figure 1 illustrates how the cost of a training iteration scales with kernel
size and therefore, the number of network parameters. The points indicate
models with $z=6$ and odd kernel sizes from 3 to 11. The scaling is consistent
with using higher order numerical schemes.

Figure 1: Training time in seconds per sample per epoch for 5 kernel sizes (3,
5, 7, 9, 11), shown on a log scale.

Although $\mathbf{z}(t)$ is a continuous function, and evaluated as such, the
solution is saved at discrete time points corresponding to the training and
test set. The decoder, $\mathcal{D}$, is another trainable network with
mirrored architecture to that of the encoder, which decodes the latent
sequence back to the original velocity space as

$\mathcal{D}(\\{\mathbf{z}_{0},\mathbf{z}_{1},...,\mathbf{z}_{n}\\})=\\{\mathbf{\hat{u}}_{0},\mathbf{\hat{u}}_{1},...,\mathbf{\hat{u}}_{n}\\}\quad.$
(6)

The network architecture is visualized in figure 2.

Figure 2: Schematic of the overall model architecture. (Encoder): initial
conditions are given and encoded into a latent space with a CNN. (Hidden
dimensions): augmented channels are concatenated to increase dimensionality of
the latent space. (NeuralODE): the latent dynamics are forecasted through the
NeuralODE approximator. (Project): augmented channels are removed by
projecting down into the original latent space. (Decoder): the sequence is
decoded via another CNN. (Spectral projection): the divergence-free condition
is enforced through the last layer.

In order to build a more detailed understanding of the model, we run a series
of sensitivity studies to determine the effect of the two subnetwork
parameterizations on model accuracy. In our experiments, we found that certain
hyperparameters such as the size of the hidden layers did not have a
significant impact on the results. However, the two hyperparameters with
strong connections to traditional numerical methods, the compression ratio and
the kernel size of the NeuralODE, were found to be a good indicator of model
performance at varying scales.

In our a-posteriori analysis, we test 13 different models: 4 compression
ratios (6, 12, 24, and 48) and 3 kernel sizes (3, 5, and 7), shown in table 1,
to determine their impact. We also examined a “deep” model, which used a
compression ratio of 6 and a 9-layer CNN with kernel sizes of 3 in the
NeuralODE, as opposed to the 3-layer networks in the other models. We will
reference these models with the key-names $\mathbf{k}\alpha\mathbf{z}\beta$,
where $\mathbf{k}$ and $\mathbf{z}$ represent the kernel size and compression
ratio respectively, and $\alpha$ is the kernel size and $\beta$ is the
compression ratio. The “deep” model is indicated with $deep$. The $deep$ model
investigates whether the larger depth of the model, which increases the
overall NeuralODE receptive field, can accommodate for the smaller kernels in
each layer. The receptive field of the 9-layer $deep$ network is equivalent to
the 3-layer $k7z6$ network.

| k=3 | k=5 | k=7 |
---|---|---|---|---
z=48 | k3z48 | k5z48 | k7z48 | -
z=24 | k3z24 | k5z24 | k7z24 | -
z=12 | k3z12 | k5z12 | k7z12 | -
z=6 | k3z6 | k5z6 | k7z6 | deep
Table 1: Outline of the models tested in our studies. $k$ and $z$ denote the
kernel size and compression ratio respectively. $deep$ denotes a 9-layer
network with kernel sizes of 3.

###  III.2 Physically-inspired model constraints

We incorporate a few architectural elements to improve prediction and enforce
physical constraints. As the problem includes periodic boundary conditions,
this constraint is imposed by padding inputs to all convolutions circularly.
We also employ the augmented NeuralODE approach to the latent dynamics, which
is capable of learning a richer set of dynamics by representing the ODE on a
higher dimensional space than the latent space input [25]. This is
accomplished by concatenating channels of zeros to the output of the encoder,
with the predicted results projected back into the original latent space. We
use 2 augmented channels in our model, selected through parameter sweeps.

Lastly, we apply a spectral projection operation to the final model output to
enforce the divergence-free field condition required for a constant density
fluid [26]. The divergence-free velocity field comes from the mass balance in
(1). In this layer formulation, the input $\mathbf{u}^{\prime}$ is transformed
into Fourier space

$\mathbf{U}^{\prime}=\mathcal{F}(\mathbf{u}^{\prime})\quad,$ (7)

where the divergence operator can be applied linearly. With this linear
constraint, we can construct a quadratic optimization problem to minimize the
squared difference between the input and output in Fourier space:

$\displaystyle\min_{\mathbf{\hat{U}}}$
$\displaystyle\frac{1}{2}(\mathbf{U}^{\prime}-\mathbf{\hat{U}})^{T}(\mathbf{U}^{\prime}-\mathbf{\hat{U}})=$
(8) $\displaystyle\min_{\mathbf{\hat{U}}}$
$\displaystyle\frac{1}{2}\mathbf{\hat{U}}^{T}\mathbf{I}\mathbf{\hat{U}}-\mathbf{U}^{\prime
T}\mathbf{\hat{U}}$ s.t.
$\displaystyle\mathbf{A}\mathbf{\hat{U}}=\mathbf{0}\quad,$

where $\mathbf{A}$ is the linear divergence operator. The solution can be
obtained directly:

$\mathbf{\hat{U}}=\mathbf{U}^{\prime}-\frac{\mathbf{k}\cdot\mathbf{U}^{\prime}}{\mathbf{k}\cdot\mathbf{k}}\mathbf{k}\quad,$
(9)

where $\mathbf{k}$ is a Fourier wavenumber vector. The final output comes from
taking the inverse Fourier transform of the optimum,

$\mathbf{\hat{u}}=\mathcal{F}^{-1}(\mathbf{\hat{U}})\quad.$ (10)

Each operation in the layer is fully differentiable and thus backpropagation
through the layers with automatic differentiation is accomplished by the
software without manual intervention.

###  III.3 Data and Training

We obtain solutions to (1) for a triply-periodic domain with direct numerical
simulation (DNS) using standard Fourier psuedo-spectral method [c.f. 27]. The
system is forced at low-wavenumbers to keep the total energy in the system
constant [16]. The resulting data are statistically stationary with
approximately stationary energy spectra. Our training, validation and test
datasets span 100 integral time scales $\mathrm{\tau}\equiv e/\epsilon$, where
$\epsilon$ is the turbulent kinetic energy dissipation rate. The turbulent
Reynolds number, $\mathrm{Re_{T}\equiv e^{2}/(\nu\epsilon)}$ is approximately
380. The solutions are discretized with 64 collocation points in each
direction, for a total of $64^{3}$ points and de-aliased with a standard two-
thirds Fourier truncation filter. We collect velocity fields every $\tau/100$,
which are used to train the neural network.

While the algorithm accepts a single instantaneous velocity field as input,
during training, the trajectory is regressed with sampled temporal DNS
snapshots from the forecasting window. The dynamics are learned from gradient-
based optimization of the predicted sequence sampled at the same temporal rate
as the DNS. From the NODE formulation, the prediction horizon can be adjusted
on-the-fly, during training or validation. We use the normalized mean squared
error (MSE) as the cost function for optimization. For this problem, we
normalize the loss with the average energy in the flow. The loss can be
written as:

$\mathrm{L(\mathbf{u},\mathbf{\hat{u}})=\frac{\sum_{t=0}^{n}\langle\lVert\mathbf{u}_{t}-\mathbf{\hat{u}}_{t}\rVert^{2}\rangle}{e}\quad},$
(11)

where $\mathbf{u}_{t}$ and $\mathbf{\hat{u}}_{t}$ are the true and predicted
velocity fields at time step $\mathrm{t}$, and $\mathrm{e}$ is given by eq.
(2) and averaged over time.

All models were trained with the ADAM optimizer using a decaying learning rate
from $10^{-3}$ \- $10^{-5}$. Additionally, models were trained by
progressively increasing the prediction horizon from $\mathrm{\tau^{*}=0.25}$
to $\mathrm{\tau^{*}=1.0}$, where $\mathrm{\tau^{*}=t/\tau}$. Each model was
trained for a total of 3500 epochs, 2000 at $\mathrm{\tau^{*}=0.25}$, 1000 at
$\mathrm{\tau^{*}=0.5}$, 200 at $\mathrm{\tau^{*}=0.75}$, and 300 at
$\mathrm{\tau^{*}=1.0}$. Based on a-posteriori analysis of the turbulence
prognostics described in the subsequent section, this was determined to be
sufficient for convergence. The training set spanned $10\mathrm{\tau}$,
validation another $10\mathrm{\tau}$, and the test data encompassed the rest
of the $80\mathrm{\tau}$ in the dataset. Our model is implemented in Julia
using the _DifferentialEquations.jl_ and _Flux.jl_ packages [28]. We use the
Tsitouras 5/4 Runge-Kutta method [29] for explicit time integration of the ODE
dynamics. Model training was accomplished with NVIDIA V100 GPUs with a memory
of 32GBs. GPU memory was a critical resource due to the large size of the
dataset.

##  IV Results

We analyze the predictive capabilities of our model by examining a series of
turbulent statistical metrics and comparing them to those of the DNS dataset.
While standard statistical measures like local RMSE are often employed in
machine learning applications to judge performance of models, and indeed we
use the local MSE in our loss function, these values are limited in their
physical insight and ultimately only offer a holistic qualitative measure of
accuracy compared to the true data. In this section, we analyze data obtained
from the DNS solutions at the suite of models described in the previous
section, evaluated for one large eddy turnover time from an initial condition
unseen during training.

Figure 3a depicts the turbulent kinetic energy in the flow, given by eq. 2, as
a function of $\tau^{*}$ for each of the 13 models in table 1. We normalize
each of the trajectories by the kinetic energy of the DNS, which deviates by
less than one percent in this time window. At $\tau^{*}=0$, we observe each
model to systematically underpredict the kinetic energy density at the initial
condition. We expect that our encoder subnetwork, which reduces the
representation of the velocity field state as embodied by the compression
ratio $z$, to resolve a subset of turbulent fluctuations. Indeed, at
$\tau^{*}=0$ we observe the kinetic energy resolved by the velocity field
reduction to decrease as the $z$ increases. Such phenomenology can be
described by the encoder’s exposure to larger and more energetic scales as
$z$, and therefore the receptive field size, increases.

Figure 3: We show three statistics that vary with time. The prediction window
is normalized to the integral timescale of the flow. (a) turbulent kinetic
energy in the flow. The predicted flow loses some energy but remains within
reasonable bounds to the DNS data throughout. (b) velocity autocorrelation
function on a log-log plot. (c) the MSE of the prediction normalized to the
average energy in the flow.

Such phenomenology is observed in figure 4, which provides a qualitative
visual comparison of the output of the network with the DNS ground truth
midway through the prediction window at $\tau^{*}=0.5$. Generally, large-scale
structures appear to be preserved, while smaller scale fluctuations are
filtered out by the model.

Figure 4: Qualitative visual comparison of the ground truth DNS data and the
output of the $k7z6$ network. Shown here is a 3D snapshot of the velocity
magnitude at $\tau*=0.5$. Large scale structures appear to be preserved while
small scale fluctuations are filtered out.

At later times in figure 3a, with the exception of $k7z48$, which appears to
be an outlier, the models tend to slightly underpredict the flow energy in the
initial phase before diverging more significantly after $\tau^{*}=0.5$.
Solutions obtained from many models go unstable by $\tau^{*}=1$. Models
$k7z6$, $k5z6$, and $deep$ are the few that disregard this trend, which appear
to have stable, perhaps overdissipative, solutions by $\tau^{*}=1$. Clearly,
$k7z6$ tracks the DNS kinetic energy the closest.

Next, kinetic energy powerspectra normalized by the true DNS spectra are shown
in figure 5 for the model predictions. In panel 5a, we show spectra from
fields of the neural network model at the initial condition. The wavenumber
$x-$axis is non-dimensionalized by $L_{E}$, given by:

$L_{E}=\frac{e^{3/2}}{\epsilon}\quad.$ (12)

By comparing the DNS and neural network initial condition, denoted
$\tau^{\ast}=0$, we first note the impact of the dimensionality reduction in
removing smallest scale velocity fluctuations, whereby we observe a rapid
attenuation of the predicted power spectra near $||k||L_{E}=20$. Again,
$k7z48$ appears to be an outlier. We note that at this initial condition, the
NeuralODE subnetwork has zero impact on the power spectra such that sources of
error are constrained to the encoder/decoder subnetworks. With respect to the
DNS solutions, we observe the neural network predictions to be accurate at
small wavenumbers approximately less than 20, consistent with alternative
reduced-order solution techniques such as large-eddy simulation.

Fig 5b shows model spectra at $\tau^{\ast}=1$, where the DNS powerspectra are
almost identical to those of the initial condition, owing to the forcing
procedure applied to the system. Here, we see the effect of the temporal
evolution of the models, which now differentiate themselves significantly at
all wavenumbers. In general, the $z=6$ models are superior in accuracy, in the
low wavenumber region below 20 and even in the highest wavenumber region,
greater than 100. We remark that the wavenumber regime where small scale
fluctuations begin to attenuate at $\tau=0$ in the neural network model is
strongly correlated with the wavenumber regime where deviations develop at
$\tau=1$, indicating that the NeuralODE subnetwork is able to accurately model
dynamics for the flow scales which are resolved by the encoder/decoder
subnetworks.

Figure 5: Both plots show the energy of the predicted flow as a function of
non-dimensionalized wavenumber, normalized to the energy of the true DNS data.
(a) A snapshot of the energy spectra for all models corresponding to the
initial time step $\tau^{*}=0$. Only the encoder/decoder networks impact the
solution. (b) The energy spectra at the final time step $\tau^{*}=1$. For
models with large receptive fields and latent spaces, low wavenumbers show
better agreement with DNS, which are sufficient for many flow of engineering
interest. Models with high compression in the latent and dynamics
representations produce less realistic results with greater instability.

We show the temporal velocity autocorrelations of the models, given as:

$\rho_{ii}(\tau^{*})=\frac{\langle\mathbf{u}_{i}(0)\mathbf{u}_{i}(\tau^{*})\rangle}{\langle\mathbf{u}_{i}(0)^{2}\rangle}\quad,$
(13)

again normalized to the true DNS data, in figure 3b. The models begin to
diverge at $\tau^{*}=0.1$ with the three best performing networks, $k7z6$,
$k5z6$, and $deep$, once again separating themselves from the rest. We also
see the greater accuracy of models with larger NODE kernels by comparing
autocorrelations within a family of models with the same compression ratio.
$k7z48$ is an exception here as well. Figure 3c presents the TKE-normalized
mean-squared error over time. The 4 $z=6$ models exhibit the lowest error of
all tested.

Real turbulent flows exhibit certain universal small-scale structures [30],
such as a preference for vorticity alignment with the intermediate strain-rate
eigenvalue. The flow topology can be described using the $\mathrm{Q-R}$ plane,
which represent the second and third velocity gradient tensor invariants
respectively. Isolines of probability in the $\mathrm{Q-R}$ plane, expressing
intimate features of the turbulent flow geometry, have a nontrivial shape
documented in the literature. Different parts of the $\mathrm{Q-R}$ plane are
associated with different structures of the flow. Thus, the lower right corner
(negative $Q$ and $R$), which has higher probability than other regions,
corresponds to a pancake type of structure (two expanding directions, one
contracting) with the direction of rotation (vorticity) aligned with the
second eigenvector of the stress. We additionally perform coarse
$\mathrm{Q-R}$ tests [31], which provide a description of the flow topology at
varying length scales. The scale is described with the radius of the coarse-
graining filter given in terms of the Kolmogorov length scale $\eta$. The
tear-drop shape of the probability isolines becomes more prominent with
decrease of the coarse-graining scale.

Figures 6 and 7 compare the various network prediction $\mathrm{Q-R}$
distributions with the ground truth DNS distributions. Figure 6 demonstrates
the effect of NODE kernel size using the outputs from the $k7z6$, $k5z6$, and
$k3z6$ models. The 4 panels provide insight into the PDFs at different coarse-
graining scales and temporal snapshots. The rows correspond to the initial and
final snapshots, $\tau^{*}=0$ and $\tau^{*}=1$, while the columns represent
length scales $r=0\eta$, where no filter is used, and $r=5.33\eta$, $1/4$ of
the length scale $L$ of the domain. As expected, at $\tau^{*}=0$, where the
NODE effect is null, all three models perform almost identically. It is clear
that at both length scales the flow morphology is well reconstructed by the
encoder/decoder subsystem. At $\tau^{*}=1$, we start to see the differences
between models. With no coarse-graining applied, there is a trend towards more
symmetric PDFs as the kernel size is reduced. A symmetric distribution
suggests the presence of random noise in the prediction that strays from the
true statistical trends found in real turbulence. At $r=5.33\eta$, the
teardrop shape of the DNS data is far less pronounced and there is greater
variation between models. However, $k7z6$ best estimates the underlying true
distribution.

Figure 7 demonstrates the effect of compression ratio using the outputs from
the $k7z6$, $k7z12$, and $k7z24$ models. At $\tau^{*}=0$, perhaps
counterintuitively, we see that the morphology is well reconstructed at both
length scales despite increases in compression ratio. All models perform
similarly, and it appears the compression ratios tested have little bearing on
the predicted morphology when no dynamics are involved. In contrast, at
$\tau^{*}=1$ and $r=0\eta$, the trend towards more symmetric distributions as
the compression ratio is increased is observed again. $k7z6$ also has the best
large-scale structure conformity at $r=5.33\eta$ when compared to the other
models with larger $z$.

Figure 6: $\mathrm{Q-R}$ planes of the flow shown at two length scales and two
temporal snapshots for models with varying $k$. $z=6$ for all models shown.
Rows correspond to different time steps, $\tau^{*}=0$ and $1.0$, and columns
indicate the coarse-grained scales given in terms of the Kolmogorov microscale
$\eta$. Figure 7: $\mathrm{Q-R}$ planes of the flow shown at two length scales
and two temporal snapshots for models with varying $z$. $k=7$ for all models
shown. Rows correspond to different time steps, $\tau^{*}=0$ and $1.0$, and
columns indicate the coarse-grained scales given in terms of the Kolmogorov
microscale $\eta$.

Additionally, figure 8 shows the magnitude of divergence in the velocity field
over one integral timescale for two models. The impact of the spectral
projection step is significant in reducing the output divergence to magnitudes
comparable to the true field. The formulation of the layer as a constrained
optimization problem minimizing the deviations between $\mathbf{U}^{\prime}$
and $\mathbf{\hat{U}}$, from eq. 8, implies inputs $\mathbf{u}^{\prime}$ and
outputs $\mathbf{\hat{u}}$ are highly comparable. The average correlation
coefficient between $\mathbf{u}^{\prime}$ and $\mathbf{\hat{u}}$ is 0.9986.
This lends additional credence to the idea that the latent NeuralODE is
learning the majority of the turbulent dynamics, and this step acts as a
simple projection to correct the velocity field. The strong correlation also
means training dynamics are not significantly affected by using the projection
step, and MSE and other evaluation metrics are comparable between the models.

Figure 8: The model uses a divergence-constraining layer to induce physical
solution fields. The magnitude of the divergence is plotted as a function of
time for three cases: the true divergence, the prediction from the model
without any constraint layer, and the prediction from the same model with the
divergence constraint. The difference between the predictions is significant,
and it is clear the constrained model performs close to the expected output.

##  V Conclusions

Encouraged by broader adoption of deep neural network architectures for
reduced order modeling of physics problems, we explore a model design paradigm
which attempts to reduce ambiguity in the parameterization of model error. We
have approached this objective by considering an idealized problem setting of
the incompressible Navier Stokes equations with spatio-temporal homogeneity, a
system with strong nonlinearity and solution complexity. In addition to
implementing hard physical constraints, such as the implicit preservation of
incompressibility, we have attempted to explicitly decouple individual
modeling demands into discrete and interpretable subnetworks. Firstly
utilizing convolutional encoders and decoders to reduce velocity field
complexity, and also utilizing a NeuralODE architecture to model the dynamics
imposed by the Navier Stokes equations.

We assess the methodology by analyzing forecasts from this scheme with varying
model parameters and comparing them to ground truth DNS data. While
qualitative comparisons to labeled data may be sufficient for non-scientific
applications, we are interested in a rigorous study of the turbulent physics
observed in the model’s outputs. Statistical metrics like energy spectra,
$\mathrm{Q-R}$ plane PDFs, and autocorrelations are used to illustrate the
quality of predictions.

Across the board, we see this surrogate approach’s inclination to prefer
preserving the majority energy-containing large eddies in the flow while
filtering the smaller length scale fluctuations. To an extent, decreasing the
latent space size still allows for conservation of energy at the large scales,
but we see increasingly larger deviations in the temporal correlations and
small-scale performance. In particular, there is a shift towards more
symmetric distributions in the $\mathrm{Q-R}$ plane, indicating that the
network is predicting noise instead of correlated statistics. The kernel size
of layers in the dynamical CNN also imparts an effect on physical
interpretations of the results. The smallest kernel network tested showed non-
stationarity and energy loss even at the larger scales from the spectra, as
well as divergence in the flow morphology at the small-, and to a lesser
extent the medium- and large-scales. At the limit of our computational
capabilities, model $k7z6$, we find excellent agreement with the ground truth.
Some of the energy at the smaller length scales is lost, also reflected in the
TKE, but because the smaller eddies contain less energy, over 90% of the total
energy in the flow is conserved. The divergence-free velocity fields from the
forecasts imply mass in the system is conserved as well.

Building on previous research in scientific ML and turbulence modeling with
ML, this work further demonstrates the effectiveness of data-driven models to
be useful surrogates for approximating fluid dynamics. Data-driven
convolutional NeuralODEs are multiple orders of magnitude faster as surrogates
than traditional simulation methods and still capture important details of the
flow, which supports use-cases where rapid predictions are essential, such as
in design or parameter-space exploration applications. The addition of
physics-grounded learning biases in these models is especially pragmatic for
constraining the model and producing more physical outputs. We have also taken
steps towards parameterizing the quality of such convolutional NeuralODE
architectures, particularly for fluid dynamics applications, paving an avenue
for practitioners to approach real-world engineering or other flow problems
with data-driven methods. As the field of machine learning and fluid dynamics
rapidly develops, others have been pushing the bounds of data-driven physics-
informed models, even incorporating traditional solver machinery into their
architectures. We see the space of possible models as a spectrum, where there
is a tradeoff in computational cost and accuracy as one incorporates
increasingly more physics into the model. It will be vital moving forward for
practitioners to strike the right balance for their specific application, and
therefore in depth studies of models along this spectrum will be crucial to
the community.

##  VI Acknowledgements

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. DGE 1745016 awarded to VS. The
authors from CMU acknowledge the support from the Technologies for Safe and
Efficient Transportation University Transportation Center, and Mobility21, A
United States Department of Transportation National University Transportation
Center. This work was supported in part by Oracle Cloud credits and related
resources provided by the Oracle for Research program. The authors further
acknowledge the support of NVIDIA via The ICEnet Consortium, for providing
necessary GPU compute. A.T.M and G.D.P. have been supported by the LDRD
(Laboratory Directed Research and Development) program at LANL under project
20190059DR.

This work has been co-authored by employee(s) of Triad National Security, LLC
which operates Los Alamos National Laboratory (LANL) under Contract No.
89233218CNA000001 with the U.S. Department of Energy/National Nuclear Security
Administration.

## References

  * Brunton  _et al._ [2020] S. L. Brunton, B. R. Noack, and P. Koumoutsakos, Machine learning for fluid mechanics, Annual Review of Fluid Mechanics 52, 477 (2020).
  * Raissi  _et al._ [2019] M. Raissi, P. Perdikaris, and G. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378, 686–707 (2019).
  * Wang  _et al._ [2020] R. Wang, K. Kashinath, M. Mustafa, A. Albert, and R. Yu, Towards physics-informed deep learning for turbulent flow prediction, Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining 10.1145/3394486.3403198 (2020).
  * Mohan  _et al._ [2020] A. T. Mohan, N. Lubbers, D. Livescu, and M. Chertkov, Embedding hard physical constraints in neural network coarse-graining of 3D turbulence, arXiv preprint arXiv:2002.00021 (2020).
  * Thuerey  _et al._ [2020] N. Thuerey, K. Weißenow, L. Prantl, and X. Hu, Deep learning methods for Reynolds-Averaged Navier–Stokes Simulations of airfoil flows, AIAA Journal 58, 25–36 (2020).
  * Mohan  _et al._ [2019] A. Mohan, D. Daniel, M. Chertkov, and D. Livescu, Compressed convolutional LSTM: An efficient deep learning framework to model high fidelity 3D turbulence, arXiv: Fluid Dynamics (2019).
  * de Avila Belbute-Peres  _et al._ [2020] F. de Avila Belbute-Peres, T. D. Economon, and J. Z. Kolter, Combining differentiable PDE solvers and graph neural networks for fluid flow prediction (2020), arXiv:2007.04439 [cs.LG] .
  * Vlachas  _et al._ [2018] P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, and P. Koumoutsakos, Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474, 20170844 (2018).
  * Kochkov  _et al._ [2021] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer, Machine learning–accelerated computational fluid dynamics, Proceedings of the National Academy of Sciences 118, e2101784118 (2021).
  * Tompson  _et al._ [2017] J. Tompson, K. Schlachter, P. Sprechmann, and K. Perlin, Accelerating Eulerian fluid simulation with convolutional networks (PMLR, International Convention Centre, Sydney, Australia, 2017) pp. 3424–3433.
  * Portwood  _et al._ [2021] G. D. Portwood, B. T. Nadiga, J. A. Saenz, and D. Livescu, Interpreting neural network models of residual scalar flux, Journal of Fluid Mechanics 907, A23 (2021).
  * [12] O. Watt-Meyer, N. D. Brenowitz, S. K. Clark, B. Henn, A. Kwa, J. McGibbon, W. A. Perkins, and C. S. Bretherton, Correcting weather and climate models by machine learning nudged historical simulations, Geophysical Research Letters , e2021GL092555.
  * Chen  _et al._ [2019] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, Neural ordinary differential equations, arXiv:1806.07366 [cs, stat] (2019), arXiv: 1806.07366.
  * Portwood  _et al._ [2019] G. D. Portwood, P. P. Mitra, M. D. Ribeiro, T. M. Nguyen, B. T. Nadiga, J. A. Saenz, M. Chertkov, A. Garg, A. Anandkumar, A. Dengel, and et al., Turbulence forecasting via neural ODE, arXiv:1911.05180 [physics] (2019), arXiv: 1911.05180.
  * Mohan  _et al._ [2021] A. T. Mohan, K. Nagarajan, and D. Livescu, Learning stable galerkin models of turbulence with differentiable programming, arXiv preprint arXiv:2107.07559 (2021).
  * Overholt and Pope [1998] M. R. Overholt and S. B. Pope, A deterministic forcing scheme for direct numerical simulations of turbulence, Comput. Fluids 27, 11 (1998).
  * Pope [2000] S. B. Pope,  _Turbulent Flows_ (Cambridge University Press, 2000).
  * Weiss [2019] J. Weiss, A tutorial on the proper orthogonal decomposition, in  _AIAA Aviation 2019 Forum_ (American Institute of Aeronautics and Astronautics, 2019).
  * Toderici  _et al._ [2016] G. Toderici, S. M. O’Malley, S. J. Hwang, D. Vincent, D. Minnen, S. Baluja, M. Covell, and R. Sukthankar, Variable rate image compression with recurrent neural networks, in  _International Conference on Learning Representations_ (2016).
  * Rippel and Bourdev [2017] O. Rippel and L. Bourdev, Real-time adaptive image compression, in  _Proceedings of the 34th International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 70, edited by D. Precup and Y. W. Teh (PMLR, 2017) pp. 2922–2930.
  * Masci  _et al._ [2011] J. Masci, U. Meier, D. Cireşan, and J. Schmidhuber, Stacked convolutional auto-encoders for hierarchical feature extraction, in  _Artificial Neural Networks and Machine Learning – ICANN 2011_ , edited by T. Honkela, W. Duch, M. Girolami, and S. Kaski (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 52–59.
  * Erichson  _et al._ [2019] N. Erichson, M. Muehlebach, and M. Mahoney, Physics-informed autoencoders for Lyapunov-stable fluid flow prediction, ArXiv abs/1905.10866 (2019).
  * Gonzalez and Balajewicz [2018] F. Gonzalez and M. Balajewicz, Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems, ArXiv abs/1808.01346 (2018).
  * Rubanova  _et al._ [2019] Y. Rubanova, R. T. Q. Chen, and D. K. Duvenaud, Latent ordinary differential equations for irregularly-sampled time series, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 5320–5330.
  * Dupont  _et al._ [2019] E. Dupont, A. Doucet, and Y. W. Teh, Augmented neural odes, in  _Advances in Neural Information Processing Systems_ (2019) pp. 3140–3150.
  * Jiang  _et al._ [2020] C. M. Jiang, K. Kashinath, Prabhat, and P. Marcus, Enforcing physical constraints in neural neural networks through differentiable PDE layer (2020).
  * Orszag and Patterson Jr [1972] S. A. Orszag and G. Patterson Jr, Numerical simulation of three-dimensional homogeneous isotropic turbulence, Physical Review Letters 28, 76 (1972).
  * Rackauckas and Nie [2017] C. Rackauckas and Q. Nie, Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia, Journal of Open Research Software 5 (2017).
  * Tsitouras [2011] C. Tsitouras, Runge–kutta pairs of order 5(4) satisfying only the first column simplifying assumption, Computers & Mathematics with Applications 62, 770 (2011).
  * Elsinga and Marusic [2010] G. Elsinga and I. Marusic, Universal aspects of small-scale motions in turbulence, Journal of Fluid Mechanics 662, 514 (2010).
  * Chertkov  _et al._ [1999] M. Chertkov, A. Pumir, and B. I. Shraiman, Lagrangian tetrad dynamics and the phenomenology of turbulence, Physics of Fluids 11, 2394 (1999).