# Reduced order modeling for flow and transport problems with Barlow Twins
self-supervised learning

Teeratorn Kadeethum  Sandia National Laboratories, New Mexico, USA  Cornell
University, New York, USA  Francesco Ballarin  Catholic University of the
Sacred Heart, Brescia, Italy  Daniel O’Malley  Los Alamos National Laboratory,
New Mexico, USA  Youngsoo Choi  Lawrence Livermore National Laboratory,
California, USA  Nikolaos Bouklas  Cornell University, New York, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>Hongkyu Yoon  Sandia National
Laboratories, New Mexico, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a unified data-driven reduced order model (ROM) that bridges the
performance gap between linear and nonlinear manifold approaches. Deep
learning ROM (DL-ROM) using deep-convolutional autoencoders (DC-AE) has been
shown to capture nonlinear solution manifolds but fails to perform adequately
when linear subspace approaches such as proper orthogonal decomposition (POD)
would be optimal. Besides, most DL-ROM models rely on convolutional layers,
which might limit its application to only a structured mesh. The proposed
framework in this study relies on the combination of an autoencoder (AE) and
Barlow Twins (BT) self-supervised learning, where BT maximizes the information
content of the embedding with the latent space through a joint embedding
architecture. Through a series of benchmark problems of natural convection in
porous media, BT-AE performs better than the previous DL-ROM framework by
providing comparable results to POD-based approaches for problems where the
solution lies within a linear subspace as well as DL-ROM autoencoder-based
techniques where the solution lies on a nonlinear manifold; consequently,
bridges the gap between linear and nonlinear reduced manifolds. We illustrate
that a proficient construction of the latent space is key to achieving these
results, enabling us to map these latent spaces using regression models. The
proposed framework achieves a relative error of 2% on average and 12% in the
worst-case scenario (i.e., the training data is small, but the parameter space
is large.). We also show that our framework provides a speed-up of $7\times
10^{6}$ times, in the best case, and $7\times 10^{3}$ times on average
compared to a finite element solver. Furthermore, this BT-AE framework can
operate on unstructured meshes, which provides flexibility in its application
to standard numerical solvers, on-site measurements, experimental data, or a
combination of these sources.

## Introduction

A reduced order model (ROM) is devised to provide an acceptable accuracy while
utilizing a much lower computational cost compared to the full order model
(FOM) [1]. In recent years, a non-intrusive or data-driven ROM approach has
grasped attention because (1) it has a straightforward implementation (i.e.,
does not require any modifications of FOM), (2) it easily lends itself to
different kinds of physical problems, and (3) it allows for more stable and
much faster prediction than intrusive ROM for nonlinear problems [2, 3, 4, 5,
6, 7]. Traditionally, proper orthogonal decomposition (POD) is used as a data
compression tool (i.e., linear subspace approach), which is the optimal way to
construct the linear reduced manifolds. However, POD-based solutions on a
linear subspace are often restrictive for highly nonlinear problems where
reduced spaces lie in nonlinear manifolds. More recently, nonlinear
compression using autoencoder-based deep learning (DL) architectures or
nonlinear manifold approach [5, 6, 8, 9] has been suggested to reconstruct
these nonlinear manifolds, resulting in generic and more refined predictive
capabilities than linear subspace approaches for nonlinear problems. Recent
extensive comparisons, however, show a performance deficit for DL-ROM
approaches in some cases [6].

Kadeethum et al. [6] illustrate that there are two essential issues for DL-
ROM. First, the nonlinear approach outperforms its linear counterpart in
specific settings (e.g., boundary conditions and domain geometry), but the
opposite can occur in other settings. This is because POD provides the optimal
data compression in a linear subspace for the problems with fast-decaying
Kolmogorov’s n-width that measures the degree of accuracy by n-dimensional
linear subspaces [10, 11, 12, 13]. Therefore, the DL-ROM approach could not
exceed the level of POD accuracy for problems that naturally lie within linear
manifolds. However, for problems with slowly decaying Kolmogorov’s width, the
nonlinear manifold approach outperforms the linear subspace one. Even though
the authors hypothesize that a visual comparison between principal component
analysis (PCA) and t-Distributed Stochastic Neighbor Embedding (t-SNE) could
indicate which method will perform better before employing any specific
compression strategy, there is no unified model that could be used across
problem settings without an extensive case-based hyperparameter search.
Second, although the nonlinear approach excels in very complex settings, it
relies on convolutional operators, hindering its application for unstructured
meshes and limiting DL-ROM approaches to less practical problems. Hence, these
limitations in DL-ROM methods need to be resolved and tested with varying
degrees of complex problems.

Convection in porous media is an important process in various applications in
natural and engineered environments (e.g., biomedical engineering, multiphase
flow in the subsurface, seawater intrusion, geothermal energy, and storage of
nuclear and radioactive waste) [14, 15, 16, 17]. As the media temperature and
composition (fluid concentration) are altered, the dynamics of fluid density
and viscosity variations could drive the flow field through flow instabilities
[18, 19]. The gravity-driven flow problem is usually characterized by Rayleigh
number ($\mathrm{Ra}$) in which if the $\mathrm{Ra}$ is low, the flow field is
laminar, while if the $\mathrm{Ra}$ is high, the flow turns into a turbulent
regime. In cases where the driving force is strong enough (very high
$\mathrm{Ra}$), the flow might also exhibit fingering behavior [20].

Numerical simulation of gravity-driven flow in porous media has been a subject
of extensive research. Notable examples of full order model (FOM) include: (1)
TOUGH software suite, which includes multi-dimensional numerical models for
simulating the coupled thermo-hydro-mechanical-chemical (THMC) processes in
porous and fractured media [21, 22], (2) SIERRA Mechanics, which has
simulation capabilities for coupling thermal, fluid, aerodynamics, solid
mechanics and structural dynamics [23], (3) PyLith, a finite-element code for
modeling dynamic and quasi-static simulations of coupled multiphysics
processes [24], (4) OpenGeoSys project, which is developed mainly based on the
finite element method using object-oriented programming THMC processes in
porous media [25], (5) IC-FERST, a reservoir simulator based on control-volume
finite element methods and dynamic unstructured mesh optimization [26], (6)
DYNAFLOW™, a nonlinear transient finite element analysis platform [27], (7)
DARSim, multiscale multiphysics finite volume based simulator [28], (8) the
CSMP, an object-oriented application program interface, for the simulation of
complex geological processes, e.g. THMC, and their interactions [29], and (9)
PorePy, an open-source modeling platform for multiphysics processes in
fractured porous media [30]. In this study, we utilize the FOM developed in
the previous works, a locally conservative mixed finite element framework for
coupled hydro-mechanical–chemical processes in heterogeneous porous media [31,
32] in which interior penalty enriched Galerkin and mixed finite element are
employed. This FOM, however, is computationally expensive for two reasons. The
first one is the problem of interest is highly nonlinear; hence, it takes more
nonlinear iterations to converge. The second reason is to satisfy the
Courant–Friedrichs–Lewy (CFL) condition, the FOM needs to march through many
intermediate time-steps to reach the time-steps of interest [33, 34, 35].

Kadeethum et al. [6] propose a data-driven reduced order model (ROM) that can
reduce computation cost while maintaining an acceptable accuracy for natural
convection in porous media problems. The model is applicable to parameterized
problems [36, 37, 38, 1, 39, 40, 41, 13], depending on a set of parameters
($\bm{\mu}$) which could correspond to physical properties, geometric
characteristics, or boundary conditions. This model sequentially follows (1)
the offline and (2) online stages [1, 42]. The offline stage begins with
initializing a set of input parameters, which we call a training set. Then the
FOM is solved corresponding to each member in the training set (in the
following, we will refer to the corresponding solutions as snapshots). Either
linear, relying on POD, [4, 43] or nonlinear compression, depending on deep
convolutional autoencoder (DL-AE or DL-ROM) [5, 6, 8], is then used to
compress FOM snapshots to produce basis functions that span reduced spaces of
very low dimensionality, but still guarantee accurate reproduction of the
snapshots [44, 45]. The ROM can then be solved during the online stage for any
new value of $\bm{\mu}$ by seeking an approximated solution in the reduced
space.

In this work, we propose a unified data-driven ROM using a combination of
Barlow Twins (BT) self-supervised learning and an autoencoder (BT-AE) that
bridges the performance gap between linear and nonlinear manifold approaches.
In particular, we use BT self-supervised learning to maximize the information
content of the embedding with the latent space through a joint embedding
architecture [46]. With four different example cases that span the degree of
complexity to cover both linear and nonlinear problems, a comparison of the
proposed BT-AE framework with both linear (POD) and nonlinear (DL-AE) ROM
approaches is conducted to demonstrate the performance of the unified data-
driven ROM framework that (1) excels in all test cases (whether the solution
can be captured in a linear or nonlinear manifold) and (2) operates on either
structured or unstructured meshes. Importantly, this model is fully data-
driven; it could be trained by data produced by FOM, on-site measurement,
experimental data, or a combination of them. This characteristic can provide
flexibility across the spectrum in more complex problems. Since it is not
limited by the Courant–Friedrichs–Lewy condition for conventional FOMs, it
could deliver quantities of interest at any given time contrary to the FOM
[6].

## Results

### Data generation

We present a summary of all geometries and boundary conditions we use in
Figure 1. In short, Examples 1, 2, and 3 represent cases where $\bm{\mu}$ is a
scalar quantity, namely $\mathrm{Ra}$, while Example 4 illustrates a case
where $\bm{\mu}$ is a four-dimensional vector, composed of $\mathrm{Ra_{1}}$,
$\mathrm{Ra_{2}}$, $\mathrm{Ra_{3}}$, and $\mathrm{Ra_{4}}$. The information
of each example is presented in Table 1. We note that
$\mathrm{M}_{\mathrm{validation}}$ and $\mathrm{M}_{\mathrm{test}}$ represent
the number of the validation and testing sets with varying Rayleigh number
($\mathrm{Ra}$), respectively (Table 1). Due to time dependence, the total
number of training, validation, and test samples is the product of
$\mathrm{M}$ and $N^{t}$ with varying $N^{t}$ ranges. Specifically the
validation samples, $\mathrm{M}_{\mathrm{validation}}N^{t}$, is determined by
$\mathrm{M}_{\mathrm{validation}}N^{t}=0.1\mathrm{M}N^{t}$ by randomly
sampling 10% of the sum of training/validation sets ($\mathrm{M}N^{t}$).

Figure 1: Domain and boundary conditions for (a) Example 1 (heating from the
left boundary), (b) Example 2 (Elder problem), (c) Example 3 (unit cell of
micromodel), and (d) Example 4 (modified Hydrocoin with four subdomains). The
red line indicates the region of the boundary where the temperature is
elevated. Table 1: Summary of main information for each example.

| Example 1 | Example 2 | Example 3 | Example 4 | remark
---|---|---|---|---|---
$\mathrm{M}$ | 40 | 40 | 40 | 81 | training set for the parameter space $\bm{\mu}$
$\mathrm{M}_{\mathrm{validation}}N^{t}$ | 10% of $\mathrm{M}N^{t}$ | 10% of $\mathrm{M}N^{t}$ | 10% of $\mathrm{M}N^{t}$ | 10% of $\mathrm{M}N^{t}$ | validation set - randomly select from $\mathrm{M}N^{t}$
$\mathrm{M}_{\mathrm{test}}$ | 10 | 10 | 10 | 10 | test set for the parameter space $\bm{\mu}_{\mathrm{test}}$
$\mathrm{M}N^{t}$ | 16802 | 36110 | 44354 | 90175 | total training/validation data
$\mathrm{M}_{\mathrm{test}}N^{t}$ | 3260 | 8951 | 11238 | 11432 | total testing data
$N^{t}$ range | [226, 477] | [790, 1010] | [951, 1265] | [907, 1280] | for training, validation, and test sets
$N_{h}^{T}$ | 7110 | 9600 | 17064 | 11382 | degree of freedom (DOF): $T_{h}$
$t$ range | [0.0, 0.1] | [0.0, 0.1] | [0.0, 0.1] | [0.0, 0.1] | $\left[t^{0},t^{N}\right]$
$\bm{\mu}$ | $\mathrm{Ra}\in[40,80]$ | $\mathrm{Ra}\in[350,450]$ | $\mathrm{Ra}\in[350,450]$ | $\mathrm{Ra_{1}}\in[350,450]$ | only Example 4 has four parameters
$\mathrm{Ra_{2}}\in[350,450]$
$\mathrm{Ra_{3}}\in[350,450]$
$\mathrm{Ra_{4}}\in[350,450]$

It is noted that the final total training samples are $0.9\mathrm{M}N^{t}$
because we allocate 10% of the training samples for the validation set. The
total of testing data is $\mathrm{M_{test}}N^{t}$. We want to emphasize that
our $N^{t}$ is not constant, but it is a function of $\bm{\mu}$. To elaborate,
the higher $\mathrm{Ra}$ value will result in the higher $N^{t}$ to satisfy
CFL condition.

The summary of each model, including the subspace dimension and compression
method, is presented in Table 2. The detailed description of POD, AE, and DC-
AE models is provided in Kadeethum et al. [6], and our newly developed BT-AE
models are described in the Methodology section. In short, for POD models, we
use proper orthogonal decomposition as a compression tool. The AE models use
an autoencoder as a compression method. We employ a deep convolutional
autoencoder to compress our training snapshots ($\mathrm{M}N^{t}$) for DC-AE
models. The BT-AE models utilize a combination of an autoencoder and Barlow
Twins self-supervised learning in their compression procedure. For the POD
models, linear compression, subspace dimension refers to the number of reduced
basis or $\mathrm{N}$ as well as the number of intermediate reduced basis or
$\mathrm{N_{int}}$. We assume $\mathrm{N}=\mathrm{N_{int}}$ for all models for
simplicity. The subspace dimension is the number of latent space
($\mathrm{Q}$) for the nonlinear compression, AE, DC-AE, and BT-AE models.

Table 2: Summary of naming for each model.

model name | compression | subspace dimension | compression techniques
---|---|---|---
POD 16 RB | linear | 16 | proper orthogonal decomposition
POD 50 RB | linear | 50 | proper orthogonal decomposition
POD 100 RB | linear | 100 | proper orthogonal decomposition
POD 500 RB | linear | 500 | proper orthogonal decomposition
AE 4 Q | nonlinear | 4 | autoencoder
DC-AE 4 Q | nonlinear | 4 | deep convolutional autoencoder
BT-AE 4 Q | nonlinear | 4 | Barlow twins + autoencoder
AE 16 Q | nonlinear | 16 | autoencoder
DC-AE 16 Q | nonlinear | 16 | deep convolutional autoencoder
BT-AE 16 Q | nonlinear | 16 | Barlow twins + autoencoder
AE 256 Q | nonlinear | 256 | autoencoder
DC-AE 256 Q | nonlinear | 256 | deep convolutional autoencoder
BT-AE 256 Q | nonlinear | 256 | Barlow twins + autoencoder

Details of POD, AE, DC-AE models are provided in Kadeethum et al. [6].

### Comparisons of BT-AE with POD, AE, and DC-AE models in simple domains

We first compare the BT-AE model accuracy (for different numbers of
$\mathrm{Q}$) with the models developed by Kadeethum et al. [43, 6] (i.e.,
POD, AE, and DC-AE models) in relatively simple model domains. Example 1
illustrates a case where a linear manifold is optimal, while Example 2
presents a case where a nonlinear manifold is optimal. The results of POD, AE,
and DC-AE models presented in Kadeethum et al. [6] demonstrated that the POD-
based and DL-ROM approaches are more suitable for the linear and nonlinear
manifold problems, respectively, and they are used in this manuscript to
evaluate the performance of BT-AE models.

#### Example 1: Heating from the left boundary

The geometry and boundary conditions are shown in Figure 1a, and we adopt this
example from Zhang et al. and Kadeethum et al. [47, 6]. This example
represents a case where its fluid flow is driven by buoyancy as the fluid is
heated on the left side of the domain. The fluid then flows upwards and
rotates to the right side of the domain. We set $\bm{\mu}=(\mathrm{Ra})$, and
its admissible range of variation is $[40.0,80.0]$, see Table 1. For the
training set, we use $\mathrm{M}=40$, which lead to, in total,
$\mathrm{M}N^{t}=16802$ training data points.

We present the test case results of the BT-AE model (BT-AE 16 Q) as
supplimental information (SI-Animation-Example 1). The difference between
solutions produced by the FOM and ROM (DIFF) is calculated by

$\operatorname{DIFF}_{\varphi}(t^{k},\bm{\mu}_{\mathrm{test}}^{(i)})=\left|\varphi_{h}(\cdot;t^{k},\bm{\mu}_{\mathrm{test}}^{(i)})-\widehat{\varphi}_{h}(\cdot;t^{k},\bm{\mu}_{\mathrm{test}}^{(i)})\right|$
(1)

where $\varphi_{h}$ is a finite-dimensional approximation of the set of
primary variables corresponding to velocity, pressure, and temperature fields.
$\widehat{\varphi}_{h}$ is an approximation of $\varphi_{h}$ produced by the
ROM. Thus, $\varphi_{h}(\cdot;t^{k},\bm{\mu}_{\mathrm{test}}^{(i)})$ and
$\widehat{\varphi}_{h}(\cdot;t^{k},\bm{\mu}_{\mathrm{test}}^{(i)})$ represent
$\varphi_{h}$ and $\widehat{\varphi}_{h}$ at all space coordinates (i.e.,
evaluations at each DOF) at time $t^{k}$ with input parameter
$\bm{\mu}_{\mathrm{test}}^{(i)}$, respectively. Note that we only present the
results of the temperature field. Hence, $\varphi_{h}$ and
$\widehat{\varphi}_{h}$ represent $T_{h}$ and $\widehat{T}_{h}$, respectively.
From SI-Animation-Example 1, we observe that BT-AE 16 Q provides a reasonable
approximation of the temperature field.

Figure 2: Example 1 - results: (a) mean squared error (MSE) of each model
(please refer to Table 2), and blue texts represent a mean value of the box
plots - here we show that BT-AE 16 gives performance similar to POD-based
approaches, but AE and DC-AE models don’t, (b) data compression loss for
validation set (Equation (18)), (c) mapping using ANN loss for validation set
(Equation (19)), (d) latent space plot of DC-AE 16 Q model, and (e) latent
space plot of BT-AE 16 Q model. Latent space plots are constructed using
t-Distributed Stochastic Neighbor Embedding (t-SNE). Different colors
represent each value of $\mathrm{Ra}$ value. We calculate the t-SNE plots
using Scikit-Learn package using its default setting and perplexity of 15.

The results of Example 1 is presented in Figure 2. In Figure 2a, The
performance of the different models (Table 2) is evaluated with the mean
square error ($\mathrm{MSE}_{\varphi}(:,\bm{\mu}_{\mathrm{test}}^{(i)})$) of
the test cases defined as follows

${\mathrm{MSE}_{\varphi}(:,\bm{\mu}_{\mathrm{test}}^{(i)}):=\frac{1}{N^{t}}{\sum_{k=0}^{N^{t}}\left|\varphi_{h}(\cdot;t^{k},\bm{\mu}_{\mathrm{test}}^{(i)})-\widehat{\varphi}_{h}(\cdot;t^{k},\bm{\mu}_{\mathrm{test}}^{(i)})\right|_{\varphi}^{2}}}$
(2)

where $\mathrm{MSE}_{\varphi}(:,\bm{\mu}_{\mathrm{test}}^{(i)})$ represents
the MSE values of all $t$ for each $\bm{\mu}_{\mathrm{test}}^{(i)}$. The
$\mathrm{MSE}$ results show that BT-AE models perform better than AE and DC-AE
models. Besides, BT-AE 16 Q delivers similar $\mathrm{MSE}$ results to those
of the POD models. In contrast to the findings presented in Kadeethum et al.
[6] where the linear compression (POD) outperforms nonlinear compression (AE
and DC-AE), BT-AE models in this study could perform similar to the POD
models. To be accurate, BT-AE models still underperform, but errors are
comparable.

We then investigate how the performance of BT-AE models compares to DC-AE.
First, we examine the data compression loss of the validation set (see
Equation (18)) which is presented in Figure 2b. From this Figure, the data
compression losses of BT-AE models are slightly better than those of the DC-AE
models. Subsequently, we illustrate the mapping using ANN loss of the
validation set, see Equation (19), in Figure 2c. From Figure 2c, we observe
that the mapping losses of the BT-AE models are six orders of magnitude less
than those of the DC-AE models. This behavior shows that the BT-AE’s latent
spaces are easier to be mapped (i.e., ANN loss of the validation set for the
BT-AE mapping is much lower than that of the DC-AE.). This speculation is
explained by Figures 2d-e, using a t-Distributed Stochastic Neighbor Embedding
(t-SNE) plot. From Figure 2d, one could see that all latent variables of DC-AE
16 Q blend (i.e., you cannot differentiate among cases with different
$\mathrm{Ra}$ values.). The latent variables of the BT-AE 16 Q model, on the
other hand, shown in Figure 2e, behave in a much better structure (i.e., we
can differentiate among cases with different $\mathrm{Ra}$ values.).

#### Example 2: Elder problem

The Elder problem [48] is a significantly more complicated and ill-posed
problem [48, 49]. High $\mathrm{Ra}$ numbers considered in this case may cause
the flow instability to be fingering behavior. The domain and boundary
conditions are presented in Figure 1b [50, 47, 6]. In short, the model domain
is heated from the half of the bottom boundary (Figure 1b), and the flow is
driven upwards by buoyancy force. We set $\bm{\mu}=(\mathrm{Ra})$, and its
admissible range as $[350.0,400.0]$ (Table 1). Compared to Example 1, this
higher range of $\mathrm{Ra}$ values affects the minimum and maximum $N^{t}$
as its range increases to $[790,1010]$.

The results of Example 2 are presented in Figure 3. From Figure 3a, we observe
that all the models using nonlinear compression (AE, DC-AE, and BT-AE) perform
better than the linear compression (POD). Furthermore, the BT-AE model
accuracy is comparable to that of the DC-AE models. However, the BT-AE model
results seem to be insensitive to the number of $\mathrm{Q}$, while the DC-AE
model results are affected by the number of $\mathrm{Q}$ (i.e., the DC-AE 16 Q
and DC-AE 256 Q are more accurate than the DC-AE 4 Q.). We also present the
results of the test cases for the BT-AE 16 Q model in the supplemental
animation (SI-Animation-Example 2). From these results, we observe that the
BT-AE 16 Q model delivers a reasonable approximation of the solution $T_{h}$
(i.e., $\widehat{T}_{h}$).

Figure 3: Example 2 - results: (a) mean squared error (MSE) of each model
(please refer to Table 2), and blue texts represent a mean value of the box
plots - here we show that BT-AE models provide performance similar to DC-AE
models, but POD-based approaches and AE models don’t, (b) data compression
loss for validation set (Equation (18)), (c) mapping using ANN loss for
validation set (Equation (19)), (d) latent space plot of DC-AE 16 Q model, and
(e) latent space plot of BT-AE 16 Q model. Latent space plots are constructed
using t-Distributed Stochastic Neighbor Embedding (t-SNE). Different colors
represent each value of $\mathrm{Ra}$ value. We calculate the t-SNE plots
using Scikit-Learn package using its default setting and perplexity of 15.

We present the data compression loss of the validation set (Equation (18)) in
Figure 3b. In contrast to the ones shown in Figure 2b, the DC-AE models have a
slightly lower loss than that of the BT-AE models. We then investigate the ANN
mapping loss (see Equation (19)) of the validation set in Figure 3c. Similar
to those presented in Figure 2c, the BT-AE models have much lower mapping
losses compared to those of the DC-AE models. Among the BT-AE models, BT-AE
256 Q has the highest value of ANN mapping loss, which is expected since it
has the highest output dimension (i.e., we are mapping $t$ and $\bm{\mu}$ to
$\bm{z}^{\mathrm{Q}}$). Again, we observe a much better structure of the BT-AE
16 Q latent space than the one from DC-AE 16 Q (see Figures 3d-e). To
elaborate, the latent variables of the DC-AE 16 Q are overlapped to hinder us
from differentiating among each case (different $\mathrm{Ra}$ values). The
latent variables of the BT-AE, on the contrary, are structured in a way that
one can clearly observe different parts that represent different $\mathrm{Ra}$
values as shown in Figures 3e.

### Model performance of BT-AE models on complex geometries

From Examples 1 and 2, we have observed that the BT-AE models could provide
good results while operating on unstructured data. In this section, more
challenging geometries which require an unstructured mesh for the FOM are
evaluated with BT-AE models only since other methods are not suitable for
unstructured mesh problems.

#### Example 3: Unit cell of micromodel

Example 3 uses a unit cell of micromodel where a central part of honeycomb
shape and four corners are impermeable for flow. Still, the heat can conduct
through these five subdomains as presented in Figure 1c. Over the past decade,
the micromodel has been used to study multiple coupled processes, including
flow, reactive transport, bioreaction, and flow instability [51, 52, 53, 54,
20, 19]. The flow is initiated from an influx at the bottom of the domain.
This geometry is more complex than those utilized in Examples 1 and 2 (see
Figures 1a-b). The higher temperature at the bottom surface (shown in red)
alters a fluid density at the bottom, and subsequently, a buoyancy force
drives the flow upwards from the bottom (shown in red) to the top of the
domain. Five subdomains contain very low flow conductivity, but they can
conduct heat. Again, we set $\bm{\mu}=(\mathrm{Ra})$ and its range as
$[350.0,400.0]$ (Table 1). The range of $\mathrm{Ra}$ can also cause flow
instability. We use $\mathrm{M}=40$,
$\mathrm{M}_{\mathrm{validation}}N^{t}=10$% of $\mathrm{M}N^{t}$, and
$\mathrm{M}_{\mathrm{test}}=10$. We have in total $\mathrm{M}N^{t}=44354$
training data points.

The summary of the Example 3 results is shown in Figure 4. For all test cases
the MSE values over time in Figures 4a-c are in the range of $\approx 1\times
10^{-5}$. The MSE values tend to decrease over time until the temperature
field becomes a steady state. Besides, BT-AE models with different
$\mathrm{Q}$ values provide approximately similar results (in line with our
findings from Examples 1 and 2). The behavior infers that utilizing only a
small number of latent spaces; the model can achieve the same level of
accuracy as the one with a large number of latent spaces. This behavior is
very beneficial because the mapping between parameter space and latent space
becomes more manageable. We also present the results of the test cases for the
BT-AE 16 Q model in the supplemental animation (SI-Animation-Example 3).
Overall, BT-AE 16 Q delivers a reasonable approximation of the $T_{h}$ (i.e.,
DIFF results are low, and the relative error lies within 2 %).

Figure 4: Example 3 results: the moving average (a window size of 50) of mean
squared error (MSE) of (a) BT-AE 4 Q, (b) BT-AE 16 Q, (c) BT-AE 256 Q (please
refer to Table 2). (d) Data compression loss for validation set (Equation
(18)), (e) Barlow Twins loss for validation set (Equation (14)), (f) mapping
using ANN loss for validation set (Equation (19)), and (g) latent space plot
of BT-AE 16 Q model. Latent space plots are constructed using t-Distributed
Stochastic Neighbor Embedding (t-SNE). Different colors represent each value
of $\mathrm{Ra}$ value. We calculate the t-SNE plots using Scikit-Learn
package using its default setting and perplexity of 15.

The data compression loss (Equation (18)) is in the range of $\approx 1\times
10^{-5}$ to $1\times 10^{-6}$ (Figure 4d) which is similar to that of Example
1 (Figure 2b), but slightly lower than that of the Example 2 (Figure 3b). The
data compression loss seems to be invariant to $\mathrm{Q}$ values. We also
present the Barlow Twins loss (Equation (14)) in Figure 4e. We observe that
the Barlow Twins loss increases with increasing the $\mathrm{Q}$ value as in
Zbontar et al. [46]. This can be explained that as the $\mathrm{Q}$ value
grows larger, the cross-correlation matrix
$\mathbf{C}^{T}\left(t,\bm{\mu}\right)$ becomes bigger, resulting in more
members in Equations (15) and (16). As stated by Zbontar et al. [46], the
absolute value of Equations (15) and (16) is not as important as their trend.
To elaborate this, in Figure 4e, all models (different $\mathrm{Q}$ values)
reach their saturated points around 40 epochs, meaning that the minimization
of Equations (15) and (16) is completed.

The mapping of the latent space using ANN loss (Equation (19)) is presented in
Figure 4f. Similar to Examples 1 and 2 (Figures 2c and 3c), the mapping loss
is in range of $\approx 1\times 10^{-5}$ to $1\times 10^{-7}$. The higher
$\mathrm{Q}$ values, the mapping loss grows larger because there are more
outputs to map. We present the latent space structure in Figure 4g (only for
BT-AE 16 Q). Following the results shown in Figures 2e and 3e, the latent
structure of the BT-AE model has a good structure since we can differentiate
among different $\mathrm{Ra}$ values. This behavior stems from the fact that
the BT losses maximize the information content of the embedding with the
latent space through a joint embedding architecture.

#### Example 4: Modified Hydrocoin with four subdomains

Example 4 uses the hydrocoin problem [55, 56] with the domain geometry shown
in Figure 1d. In this example, the domain is subdivided into four subdomains
with different $\mathrm{Ra}$ values (i.e.,
$\bm{\mu}=\left(\mathrm{Ra_{1}},\mathrm{Ra_{2}},\mathrm{Ra_{3}},\mathrm{Ra_{4}}\right)$).
The range of $\mathrm{Ra}$ values is $[350.0,400.0]$. Similar to the previous
examples, this $\mathrm{Ra}$ range causes fingering behavior as shown in the
supplemental animation (SI-Animation-Example4). We use $\mathrm{M}=81$,
$\mathrm{M}_{\mathrm{validation}}N^{t}=10$% of $\mathrm{M}N^{t}$, and
$\mathrm{M}_{\mathrm{test}}=10$. We have in total $\mathrm{M}N^{t}=90175$
training data points. We note that as we use $\mathrm{M}=81=3^{4}$ equally
spaced samples, for each parameter $\mathrm{Ra}_{i}$, $i=1,2,3,4$, we only
have three values. As an example, for $\mathrm{Ra_{1}}$ we only sample
$\mathrm{Ra_{1}}=\left(350,400,450\right)$ for the training set. The same goes
for $\mathrm{Ra_{2}},\mathrm{Ra_{3}},$ and $\mathrm{Ra_{4}}$. As a result,
training with relatively sparse samples of each parameter $\mathrm{Ra}_{i}$
makes it very challenging to obtain an accurate data-driven framework in
general [1, 4].

Even though this setting is very challenging, we still observe that the BT-AE
16 Q delivers a reasonable approximation of the $T_{h}$ as seen in the
supplemental animation (SI-Animation-Example4). The summary of the Example 4
results is shown in Figure 5. We present the MSE values as a function of time
in Figure 5a-c. We can observe that the MSE values for all test cases are in
the range of $\approx 1\times 10^{-1}$ to $1\times 10^{-5}$, which are
significantly higher than those of Examples 1, 2, and 3. Moreover, the MSE
values generally increase as we approach steady-state solutions, unlike the
behaviors shown in Example 3. Again, BT-AE models with different $\mathrm{Q}$
provide approximately similar results (in line with our finding from Examples
1, 2, and 3).

Figure 5: Example 4 results: the moving average (a window size of 50) of mean
squared error (MSE) of (a) BT-AE 4 Q, (b) BT-AE 16 Q, (c) BT-AE 256 Q (please
refer to Table 2). (d) Data compression loss for validation set (Equation
(18)), (e) Barlow Twins loss for validation set (Equation (14)), (f) mapping
using ANN loss for validation set (Equation (19)), and latent space plot of
BT-AE 16 Q model for (g) $Ra_{1}$ and (h) $Ra_{4}$. Latent space plots are
constructed using t-Distributed Stochastic Neighbor Embedding (t-SNE).
Different colors represent each value of $\mathrm{Ra}$ value. We calculate the
t-SNE plots using Scikit-Learn package using its default setting and
perplexity of 15.

The data compression loss (Equation (18)) is in the range of $\approx 1\times
10^{-2}$ to $1\times 10^{-4}$ (Figure 5d), which is significantly higher than
that of Examples 1, 2, and 3. This behavior illustrates that this example is
the most challenging case for the BT-AE models. The data compression loss is
the lowest for $\mathrm{Q}=256$ and the highest for $\mathrm{Q}=4$, but the
difference is not critical. As shown in the Barlow Twins loss (Equation (14))
in Figure 5e, the higher values of $\mathrm{Q}$ the larger Barlow Twins loss
is (as we discussed in the previous example.).

The mapping of the latent space using ANN loss (Equation (19)) is presented in
Figure 5f. The mapping loss is in the range of $\approx 1\times 10^{-4}$ to
$1\times 10^{-5}$, which is significantly higher than those of Examples 1, 2,
and 3 (see Figures 2c, 3c, and 4f). This behavior also contributes to the
higher MSE values of the BT-AE models. We present the latent space structure
(only for BT-AE 16 Q) in Figure 5g-h for $\mathrm{Ra_{1}}$ and
$\mathrm{Ra_{4}}$, respectively. Since Example 4 has different $\mathrm{Ra}$
values in four subdomains, the differentiation of the latent space of
individual $\mathrm{Ra}$ does not provide good solutions as each latent space
of each subdomain might also be interconnected.

## Discussion

Recent developments in ML-based data-driven reduced order modeling (DL-ROM or
DC-AE in this study) [5, 6] have shown promising results in capturing
parametrized solutions of systems of nonlinear equations. These models,
however, rely on convolutional operators, which hinders the applicability of
these models to complex geometries where an unstructured mesh is required for
FOMs, as in Examples 3 and 4. Though we could utilize an autoencoder without
convolutional layers, the model could not achieve the same level of accuracy
as DL-ROM [6]. Kadeethum et al. [6] also illustrate that in a specific setting
(simple geometry and boundary conditions), a linear compression approach using
POD can outperform the DL-ROM model (Example 1). We have demonstrated that the
autoencoder model through Barlow Twins self-supervised learning (BT-AE) could
achieve the same accuracy as DL-ROM (Example 2 where POD models perform much
worse than DL-ROM) by regularizing the latent space or nonlinear manifolds.
Besides, it also yields optimal results in the case where the linear
compression model outperforms the DL-ROM (Example 1). It means that the BT-AE
model excels in all the test cases (Examples 1 and 2) while it still can
operate on an unstructured mesh. This behavior has a significant advantage in
scientific computing since most realistic problems require unstructured mesh
representations. Besides, the BT-AE’s performance is insensitive to the number
of latent spaces, suggesting that with only a small number of latent spaces,
the model can achieve the same level of accuracy as the one with a large
number of latent spaces. This behavior is very beneficial because the mapping
between parameter space and latent space becomes more manageable.

The computational time used to develop our ROM can be broken down into three
primary parts: (1) generation of training data through FOM (the second step in
Figure 6), (2) training BT-AE (the third step in Figure 6), and (3) mapping of
$t$ and $\bm{\mu}$ to reduced subspace (the fourth step in Figure 6). Each FOM
model (corresponding to each set of $\bm{\mu}$ or $\mathrm{Ra}$ in this work)
takes, on average, about two hours on AMD Ryzen Threadripper 3970X (4
threads). We note that our FOM utilizes the adaptive time-stepping; hence,
each $\bm{\mu}^{(i)}$ may require a substantially different computational
time. To elaborate, cases that have higher $\mathrm{Ra}$ usually have a
smaller time-step ($N^{t}$ becomes larger), and subsequently, they require
more time to complete.

The wall time used to train BT-AE is approximately 0.4 hours using a single
Quadro RTX 6000. It is noted that this computational cost is much cheaper than
that of the DC-AE model, taking around four to six hours [6]. This is because
DC-AE relies on convolutional layers, dropout, and batch normalization, which
require much higher computational resources. The BT-AE, on the other hand,
utilizes only a plain autoencoder. The BT-AE model is also cheaper than the
POD model. However, we note that this may not be a fair comparison as we
perform POD and BT-AE using different machines (i.e., our POD framework only
works on CPU, but our BT-AE is trained using GPU). Please refer to Kadeethum
et al. [6] for detailed wall time comparisons among POD and DC-AE models. The
mapping of $t$ and $\bm{\mu}$ to reduced subspace through artificial neural
networks (ANN) takes around half an hour to one hour using a single Quadro RTX
6000. As mentioned in the Methodology section, we do not terminate the
training of both BT-AE and mapping of $t$ and $\bm{\mu}$ to reduced subspace
through ANN early, but rather use the model with the best validation loss
through the final epochs. For example, we train for 50 epochs, but the model
that offers the best validation loss might be the model at 20 epochs. However,
the training time we report here is for 50 epochs. Thus, our training time
provided here is considered conservative.

Even though the ROM training time is not trivial, it could provide a fast
prediction during the online phase. Using AMD Ryzen Threadripper 3970X, the
ROM takes approximately several milliseconds for a query of a pair of $t^{k}$
and $\bm{\mu}^{(i)}$. We also note that, as discussed previously, our ROM is
needed to be trained on GPU for the problems at hand. Still, it could utilize
CPU during an online time since we do not have to deal with back-propagation
or optimization during the prediction time. On the contrary, one FOM
simulation (for each $\bm{\mu}^{(i)}$ for all $t$ $\in$
$0=:t^{0}<t^{1}<\cdots<t^{N}:=\tau$) takes about two hours. So, assuming that
we query all $t$ similar to those of the FOM, ROM takes only a matter of
several seconds. In practice, however, we might not need to evaluate all
timestamps in $0=:t^{0}<t^{1}<\cdots<t^{N}:=\tau$ because the quantities of
interest at the specific time may be more important. Since ROM is not bound by
the CFL condition and can predict the quantities of interest at any specific
time without intermediate computation, we could simply perform one query -
$t^{N}$ and $\bm{\mu}^{(i)}$, resulting in saving computational time
significantly. Our ROM could provide a speed-up of $7\times 10^{6}$ at any
specific time step for Example 2, and a speed-up of $7\times 10^{3}$ to
$7\times 10^{6}$ for all examples considered in this work.

Our model is developed upon the data-driven paradigm, which is applicable to
any FOM. Besides, it could be trained using data produced by FOM, on-site
measurements, experimental data, or a combination among them. This
characteristic provides flexibility, which intrusive approaches could not
provide. The data-driven model, though, is usually hungry for training
samples. We have illustrated that as the dimensionality of our parameter space
grows, the model requires more training samples, or it will suffer by losing
its accuracy significantly as in Example 4 compared to accurate prediction in
Example 3. We speculate that an adaptive sampling technique [57, 58, 59],
incorporating physical information [60, 61], or including multimodal
unsupervised training [62] might provide a resolution to this issue in the
future work. Another gap in data-driven machine learning ROM is that a
posteriori error is exceptionally challenging to quantify. An error estimator
developed by Xiao [63] for linear manifolds could be adapted and extended to
the nonlinear manifold paradigm. Additionally, epistemic uncertainty could
also be quantified by adopting the ensemble technique proposed by Jacquier et
al. [64].

## Methodology

A graphical summary of our procedure is presented in Figure 6: the
computations are divided into an offline phase for the ROM construction, which
we will show through four consecutive main steps and (single-step) online
stage for the ROM evaluation.

Figure 6: The summary of procedures taken to establish the proposed BT-AE.

The first step of the offline stage represents an initialization of a training
set ($\bm{\mu}$), validation set ($\bm{\mu}_{\mathrm{validation}}$), and test
set ($\bm{\mu}_{\mathrm{test}}$) of parameters used to train, validate, and
test the framework, of cardinality $\mathrm{M}$,
$\mathrm{M}_{\mathrm{validation}}$, $\mathrm{M}_{\mathrm{test}}$. For the rest
of sections we will discuss only $\bm{\mu}$. The same analogy goes for
$\bm{\mu}_{\mathrm{validation}}$ and $\bm{\mu}_{\mathrm{test}}$. Let
$\mathbb{P}\subset\mathbb{R}^{P}$, $P\in\mathbb{N}$, be a compact set
representing the range of variation of the parameters $\bm{\mu}\in\mathbb{P}$.
For the sake of notation we denote by $\mu_{p}$, $p=1,\ldots,P$, the $p$-th
component of $\bm{\mu}$. To explore the parametric dependence of the
phenomena, we define a discrete training set of $\mathrm{M}$ parameter
instances. Each parameter instance in the training set will be indicated with
the notation $\bm{\mu}^{(i)}$, for $i=1,\ldots,\mathrm{M}$. Thus, the $p$-th
component of the $i$-th parameter instance in the training set is denoted by
$\mu_{p}^{(i)}$ in the following. The choice of the value of $\mathrm{M}$, as
well as the sampling procedure from the range $\mathbb{P}$, is typically user-
and problem-dependent. In this work, we use an equispaced distribution for the
training set as done in [43, 6].

In the second step, we query the FOM, based on the finite element solver
proposed and made publicly available in Kadeethum et al. [32, 6], for each
parameter $\bm{\mu}$ in the training set. In short, we are interested in
gravity driven flow in porous media, and here we briefly describe all the
equations used in this study: (1) mass balance and (2) heat advection-
diffusion equations. Let $\Omega\subset\mathbb{R}^{d}$ ($d\in\\{1,2,3\\}$)
denote the computational domain and $\partial\Omega$ denote the boundary.
$\bm{X}^{*}$ are spatial coordinates in $\Omega$ (e.g.,
$\bm{X}^{*}=[x^{*},y^{*}]$ when $d=2$, which we will focus on throughout this
study). The time domain is denoted by $\mathbb{T}=\left(0,\tau\right]$ with
$\tau>0$ (i.e., $\tau$ is the final time). Primary variables used in this
paper are $\bm{u}^{*}(\cdot,t^{*}):\Omega\times\mathbb{T}\to\mathbb{R}^{d}$,
which is a vector-valued Darcy velocity ($\mathrm{m}\mathrm{/}\mathrm{s}$),
$p^{*}(\cdot,t^{*}):\Omega\times\mathbb{T}\to\mathbb{R}^{d}$, which is a
scalar-valued fluid pressure ($\mathrm{P}\mathrm{a}$), and
$T^{*}(\cdot,t^{*}):\Omega\times\mathbb{T}\to\mathbb{R}^{d}$, which is a
scalar-valued fluid temperature ($\mathrm{C}$). Time is denoted as $t^{*}$
($\mathrm{s}$).

Following Joseph [65], the Boussinesq approximation to the mass balance
equations results in the density difference only appearing in the buoyancy
term. The mass balance equation reads

$\bm{u}^{*}+\bm{\kappa}\left(\nabla
p^{*}+{\mathbf{y}}\left(\rho-\rho_{0}\right)g\right)=0,$ (3)

and

$\nabla\cdot\bm{u}^{*}=0$ (4)

where $\bm{\kappa}=\bm{k}/\mu_{f}$ is the porous medium conductivity, $\bm{k}$
is the matrix permeability tensor, $\mu_{f}$ is the fluid viscosity,
$\mathbf{y}$ is a unit vector in the direction of gravitational force, $g$ is
the constant acceleration due to gravity, $\rho$ and $\rho_{0}$ are the fluid
density at current and initial states, respectively. We assume that $\rho$ is
a linear function of $T^{*}$ [66, 47]

$\rho=\rho_{0}\left(1-\alpha\left(T^{*}-T_{0}^{*}\right)\right),$ (5)

where $\alpha$ is the thermal expansion coefficient, and $T_{0}^{*}$ is the
reference fluid temperature. We note that Equation (5) is the simplest
approximation, and one may easily adapt the proposed method when employing a
more complex relationship provided in [67]. The heat advection-diffusion
equation defined as

$\gamma\frac{\partial T^{*}}{\partial t^{*}}+\bm{u}^{*}\cdot\nabla
T^{*}-K\nabla^{2}T^{*}-{f_{c}}^{*}=0.$ (6)

Here, $\gamma$ is the ratio between the porous medium heat capacity and the
fluid heat capacity, $K$ is the effective thermal conductivity, and
${f_{c}}^{*}$ is a sink/source. We follow Nield and Bejan [18] and define
dimensionless variables as follows

$\bm{X}:=\frac{1}{H}\bm{X}^{*},\quad t:=\frac{\kappa}{\mu\gamma
H^{2}}t^{*},\quad
p:=\frac{\kappa}{K}p^{*},\quad\bm{u}:=\frac{H}{K}\bm{u}^{*},\quad
T:=\frac{T^{*}-T_{0}^{*}}{\Delta T^{*}},\quad f_{c}:=\frac{t^{*}}{\Delta
T^{*}}{f_{c}}^{*},$ (7)

where $H$ is the dimensional layer depth, and $\Delta T^{*}$ is the
temperature difference between two boundary layers. From these dimensionless
variables, we could rewrite our Equations (3) and (4) as

$\begin{split}\bm{u}+\nabla p-\mathbf{y}\operatorname{Ra}T=0,&\text{ \> in \>
}\Omega\times\mathbb{T},\\\ \nabla\cdot\bm{u}=0,&\text{ \> in \>
}\Omega\times\mathbb{T},\\\ p=p_{D}&\text{ \> on \>
}\partial\Omega_{p}\times\mathbb{T},\\\ \bm{u}\cdot\mathbf{n}=q_{D}&\text{ \>
on \>}\partial\Omega_{q}\times\mathbb{T},\\\ p=p_{0}&\text{ \> in \>
}\Omega\text{ at }t=0,\end{split}$ (8)

where $\partial\Omega_{p}$ and $\partial\Omega_{q}$ are the pressure and flux
boundaries (i.e., Dirichlet and Neumann boundary conditions), respectively.
Here, $\mathrm{Ra}$ is the Rayleigh number

$\mathrm{Ra}:=\frac{g\alpha\kappa\Delta T^{*}H}{K}.$ (9)

We then write Equation (6) in dimensionless form as follows

$\begin{split}\frac{\partial T}{\partial t}+\bm{u}\cdot\nabla
T-\nabla^{2}T-f_{c}=0,&\mbox{ \> in \> }\Omega\times(0,\mathbb{T}],\\\
T=T_{D}&\mbox{ \> on \> }\partial\Omega_{{T}}\times(0,\mathbb{T}],\\\
(-\bm{u}T+\nabla T)\cdot\mathbf{n}={T_{\rm in}}\bm{u}\cdot\mathbf{n}&\mbox{ \>
on \> }\partial\Omega_{{\rm in}}\times(0,\mathbb{T}],\\\ \nabla
T\cdot\mathbf{n}=0&\mbox{ \> on \> }\partial\Omega_{{\rm
out}}\times(0,\mathbb{T}],\\\ T={T_{0}}&\text{ \> in \> }\Omega\text{ at
}t=0,\end{split}$ (10)

where $\partial\Omega_{{T}}$ is temperature boundary (Dirichlet boundary
condition), $\partial\Omega_{\rm in}$ and $\partial\Omega_{\rm out}$ denote
inflow and outflow boundaries, respectively, which are defined as

$\partial\Omega_{\rm
in}:=\\{\bm{X}\in\partial\Omega:\bm{u}\cdot\mathbf{n}<0\\}\quad\mbox{ and
}\quad\partial\Omega_{\rm
out}:=\\{\bm{X}\in\partial\Omega:\bm{u}\cdot\mathbf{n}\geq 0\\}.$ (11)

The detail of discretization could be found in Kadeethum et al. [32, 6], and
the FOM source codes are provided in Kadeethum et al. [32]. After the second
step, we have $\mathrm{M}$ snapshots of FOM results associated with the
different parametric configurations in $\bm{\mu}$. Since the problem
formulation is time-dependent, the output of the FOM solver for each parameter
instance $\bm{\mu}^{(i)}$ collects the time series representing the time
evolution of the primary variables for each time-step $t$. Thus, each snapshot
contains approximations of the primary variables ($\bm{u}_{h}$, $p_{h}$, and
$T_{h}$) at each time-step of the partition of the time domain $\mathbb{T}$.
Therefore, based on the training set cardinality $\mathrm{M}$ and the number
$N^{t}$ of time-steps, we have a total of $N^{t}\mathrm{M}$ training data to
be employed in the subsequent steps. We note that as our finite element solver
utilizes an adaptive time-stepping [32, 6], each snapshot may have a different
number of time-steps $N^{t}$, i.e. $N^{t}=N^{t}(\bm{\mu})$.

The third step aims to compress the information provided by the training
snapshots provided by the second step. Kadeethum et al. [6] provide detailed
derivations and comparisons between linear and nonlinear compression.
Especially the convolutional layers, in their classical form, could not deal
with an unstructured data structure (unstructured mesh), which is very common
in scientific computing or, more specifically, finite element analysis. Hence,
our goal is to develop a nonlinear compression that (1) consistently
outperforms (or at least delivers similar accuracy) the linear compression and
(2) is compatible with an unstructured data structure.

To achieve this goal, we propose a nonlinear compression utilizing feedforward
layers in combination with self-supervised learning (SSL) of Barlow Twins (BT)
( Figure 6). The BT for redundancy reduction is proposed by Zbontar et al.
[46]. It operates on a joint embedding of noisy images by producing two
distorted images from an original one through a series of random cropping,
resizing, horizontal flipping, color jittering, converting to grayscale,
Gaussian blurring, and solarization. Since we do not operate on structured
data (image) but unstructured data produced by finite element solver, we only
employ random noise and Gaussian blur operations to produce our noisy data
set, see Figure 6.

Let ${z}_{1}^{\bm{u}},\cdots,{z}_{\mathrm{Q}}^{\bm{u}}$,
${z}_{1}^{p},\cdots,{z}_{\mathrm{Q}}^{p}$, and
${z}_{1}^{T},\cdots,{z}_{\mathrm{Q}}^{T}$ be the nonlinear manifolds of the
$\bm{u}_{h}$, $p_{h}$, and $T_{h}$, respectively. For the sake of compactness,
we will only discuss primary variable $T_{h}$. The same procedure holds for
$\bm{u}_{h}$ and $p_{h}$. Our goal is to achieve
$\mathrm{Q}\ll\mathrm{M}N^{t}$ where $\mathrm{M}N^{t}$ is the total training
data, which implies that our nonlinear manifolds could represent our training
data using much lower dimension. We employ a vanilla AE (using only
feedforward layers) that is regularized by Barlow Twins SSL to obtain
$\bm{z}^{T}=\left[{z}_{1}^{T},\cdots,{z}_{\mathrm{Q}}^{T}\right]$. We do not
use any batch normalization or dropout. The summary of the training process is
presented in Algorithm 1. We will provide the detailed implementation in
https://github.com/sandialabs.

Algorithm 1 Training autoencoder (AE) with Barlow Twins (BT) regularization

# Integrate BT into AE architecture to regularize AE’s latent spaces
$\bm{z}^{T}$
###
# Training data $T_{h}$ and distorted data $T_{h,A}$, $T_{h,B}$ are input of
encoder
# latent spaces $\bm{z}^{T}$, $\bm{z}_{A}^{T}$, and $\bm{z}_{B}^{T}$ are
output encoder
###
# latent space $\bm{z}^{T}$ is output of decoder
# Approximation of $T_{h}$, i.e., $\widehat{T_{h}}$ is output of decoder
###
# latent spaces $\bm{z}_{A}^{T}$ and $\bm{z}_{B}^{T}$ are input of projector
# cross-correlation matrix $\mathbf{C}^{T}$ is output of projector
###

1:Initialize (or load pre-trained models) encoder, decoder, and projector
$\triangleright$ size of latent space $\mathrm{Q}$ has to be specified.

2:Initialize (or load pre-trained optimizers) two optimizers for AE and BT.

3:Load training set $\bm{\mu}$ $\triangleright$ the total training data is
$\mathrm{M}N^{t}$

4:Randomly select 10% of $\mathrm{M}N^{t}$ for validation set
$\bm{\mu}_{\mathrm{validation}}$ $\triangleright$ the total training data
becomes 90% of $\mathrm{M}N^{t}$

5:Add random noise $\triangleright$ see Equation (12)

6:Add Gaussian blur $\triangleright$ see Equation (13)

7:From step 5 and 6, we obtain $T_{h,A}\left(t,\bm{\mu}\right)$ and
$T_{h,B}\left(t,\bm{\mu}\right)$ from $T_{h}\left(t,\bm{\mu}\right)$

8:for each epoch do

9:  _Outer loop: training BT_ $\triangleright$ Batch size
$\mathbf{B}_{\mathrm{outer}}=512$

10: for each $\mathrm{B}_{\mathrm{outer},i}$ $\in$
$\mathbf{B}_{\mathrm{outer}}$ do

11:
$\bm{z}^{T}_{h,A}\left(t,\bm{\mu}\right)=\mathrm{encoder}\left({T}_{h,A}\left(t,\bm{\mu}\right)\right)$

12:
$\bm{z}^{T}_{h,B}\left(t,\bm{\mu}\right)=\mathrm{encoder}\left({T}_{h,B}\left(t,\bm{\mu}\right)\right)$

13:
$\mathbf{C}^{T}\left(t,\bm{\mu}\right)=\mathrm{projector}\left(\bm{z}^{T}_{h,A}\left(t,\bm{\mu}\right),\bm{z}^{T}_{h,B}\left(t,\bm{\mu}\right)\right)$

14: Calculate BT loss $\mathcal{L}_{\mathrm{BT}}^{T}$ $\triangleright$ see
Equation (14)

15: Back-propagation of BT loss w.r.t. each
$\mathrm{encoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ and
$\mathrm{projector}\left({\mathbf{W}},{\mathbf{b}}\right)$

16: Update $\mathrm{encoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ and
$\mathrm{projector}\left({\mathbf{W}},{\mathbf{b}}\right)$ using BT optimizer

17: Update learning rate $\eta_{c}$ of BT optimizer $\triangleright$ see
Equation (17)

18:  _Inner loop: training AE_ $\triangleright$ Batch size
$\mathbf{B}_{\mathrm{inner}}=32$

19: for each $\mathrm{B}_{\mathrm{inner},i}$ $\in$
$\mathbf{B}_{\mathrm{inner}}$ do

20:
$\bm{z}_{h}^{T}\left(t,\bm{\mu}\right)=\mathrm{encoder}\left({T}_{h}\left(t,\bm{\mu}\right)\right)$

21:
$\widehat{{T}_{h}}\left(t,\bm{\mu}\right)=\mathrm{decoder}\left(\bm{z}_{h}^{T}\left(t,\bm{\mu}\right)\right)$

22: Calculate AE loss ${\mathcal{L}_{\mathrm{AE}}^{T}}$ (data compression
loss) $\triangleright$ see Equation (18)

23: Back-propagation of AE loss w.r.t. each
$\mathrm{encoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ and
$\mathrm{decoder}\left({\mathbf{W}},{\mathbf{b}}\right)$

24: Update $\mathrm{encoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ and
$\mathrm{decoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ using AE optimizer

25: Update learning rate $\eta_{c}$ of AE optimizer $\triangleright$ see
Equation (17)

26: end for

27: end for

28:end for

It is noted that, we here only discuss primary variable $T_{h}$ for the sake
of compactness. The same procedures are hold for $\bm{u}_{h}$ and $p_{h}$.
Moreover, this algorithm only reflects the third step in Figure 6.

In short, during the training phase, our BT-AE model is composed of one
encoder, one decoder, and one projector. The training entails two sub-tasks;
the first is BT (encoder and projector), which takes place in the outer loop.
The second sub-task is responsible for the training of AE (encoder and
decoder), which takes place inside the inner loop. The main reasons for this
procedure are two folds. The first reason is Zbontar et al. [46] states that
the BT works better with large batch sizes. The AE, however, generally
requires a small batch size [68, 69]. Our previous numerical experiments based
on DC-AE [6] also align with this statement. Consequently, we set our batch
size of the outer loop as $\mathbf{B}_{\mathrm{outer}}=512$, and the batch
size of the inner loop as $\mathbf{B}_{\mathrm{inner}}=32$

Prior to the training, we distort our training set (i.e., creating
$T_{h,A}\left(t,\bm{\mu}\right)$ and $T_{h,B}\left(t,\bm{\mu}\right)$ from
$T\left(t,\bm{\mu}\right)$) through a series of two operations. First, _add
random noise_ is added as follows

$\widetilde{T_{h,A}}\left(t,\bm{\mu}\right),\widetilde{T_{h,B}}\left(t,\bm{\mu}\right)=T\left(t,\bm{\mu}\right)+\epsilon\operatorname{SD}\left(T\left(t,\bm{\mu}\right)\right)\mathcal{G}\left(0,1\right)$
(12)

where
$\widetilde{T_{h,A}}\left(t,\bm{\mu}\right),\widetilde{T_{h,B}}\left(t,\bm{\mu}\right)$
are intermediate distorted input data. The constant $\epsilon$, which is set
to $0.1$, determines the noise level as it is multiplied with the standard
deviation of the input field. $\mathcal{G}\left(0,1\right)$ is a random value
which is sampled from the standard normal distribution with mean and standard
deviation of zero and one, respectively.

Subsequently, we pass
$\widetilde{T_{h,A}}\left(t,\bm{\mu}\right),\widetilde{T_{h,B}}\left(t,\bm{\mu}\right)$
through Gaussian blur operation, which reads

${T_{h,A}}\left(t,\bm{\mu}\right),{T_{h,B}}\left(t,\bm{\mu}\right)=\frac{1}{\sqrt{2\pi\operatorname{SD}\left(\widetilde{T_{h,A}}\left(t,\bm{\mu}\right),\widetilde{T_{h,B}}\left(t,\bm{\mu}\right)\right)^{2}}}\mathrm{exp}\left({-\frac{{\widetilde{T_{h,A}}\left(t,\bm{\mu}\right),\widetilde{T_{h,B}}\left(t,\bm{\mu}\right)}^{2}}{2\operatorname{SD}\left(\widetilde{T_{h,A}}\left(t,\bm{\mu}\right),\widetilde{T_{h,B}}\left(t,\bm{\mu}\right)\right)^{2}}}\right)$
(13)

to obtain $T_{h,A}\left(t,\bm{\mu}\right)$ and
$T_{h,B}\left(t,\bm{\mu}\right)$.

We use a number of the epoch of 50, see Algorithm 1. The outer loop works as
follows: training BT begins with passing $T_{h,A}\left(t,\bm{\mu}\right)$ and
$T_{h,B}\left(t,\bm{\mu}\right)$ to the encoder (it is noted we have only one
encoder) resulting in $\bm{z}_{A}^{T}\left(t,\bm{\mu}\right)$ and
$\bm{z}_{B}^{T}\left(t,\bm{\mu}\right)$. We then use
$\bm{z}_{A}^{T}\left(t,\bm{\mu}\right)$ and
$\bm{z}_{B}^{T}\left(t,\bm{\mu}\right)$ as input to the projector resulting in
cross-correlation matrix $\mathbf{C}^{T}\left(t,\bm{\mu}\right)$.
$\mathbf{C}^{T}\left(t,\bm{\mu}\right)$ is a square matrix with the
dimensionality of the projector’s output, and its values range between -1,
perfect anti-correlation, and 1, perfect correlation.

The Barlow Twins loss $\mathcal{L}_{\mathrm{BT}}^{T}$ (BT loss) is then
calculated using

$\mathcal{L}_{\mathrm{BT}}^{T}:=\mathcal{L}_{\mathrm{I}}^{T}+\mathcal{L}_{\mathrm{RR}}^{T}$
(14)

where

$\mathcal{L}_{\mathrm{I}}^{T}:=\sum_{i}\left(1-\mathbf{C}_{ii}^{T}\left(t,\bm{\mu}\right)\right)^{2},$
(15)

and

$\mathcal{L}_{\mathrm{RR}}^{T}:=\lambda\sum_{i}\sum_{j\neq
i}{\mathbf{C}_{ij}^{T}\left(t,\bm{\mu}\right)}^{2},$ (16)

where $\mathbf{C}_{ii}^{T}\left(t,\bm{\mu}\right)$ denotes the $i$-th diagonal
entry of $\mathbf{C}^{T}\left(t,\bm{\mu}\right)$, $\lambda$ is a positive
constant, which is set to $5\times 10^{-3}$ as recommended by Zbontar et al.
(2021) [46], and $\mathbf{C}_{ij}^{T}$ are off-diagonal entries of
$\mathbf{C}^{T}$. In short, we train our BT part by trying to force
$\mathcal{L}_{\mathrm{I}}^{T}$ to 1, but $\mathcal{L}_{\mathrm{RR}}^{T}$ to 0
resulting in teaching the encoder and projector learn how to get rid off noise
from the distorted data, $T_{h,A}\left(t,\bm{\mu}\right)$ and
$T_{h,B}\left(t,\bm{\mu}\right)$, and construct a representation that
conserves as much $T\left(t,\bm{\mu}\right)$ information as possible.

Here, we follow the training procedures used by Kadeethum et al. [70, 6]. We
use the ADAM algorithm[71] to adjust learnable parameters of
encoder($\mathrm{W}$ and $\mathrm{b}$) and projector($\mathrm{W}$ and
$\mathrm{b}$). The learning rate ($\eta$) is calculated as [72]

$\eta_{c}=\eta_{\min}+\frac{1}{2}\left(\eta_{\max}-\eta_{\min}\right)\left(1+\cos\left(\frac{\mathrm{step_{c}}}{\mathrm{step_{f}}}\pi\right)\right)$
(17)

where $\eta_{c}$ is a learning rate at step $\mathrm{step_{c}}$, $\eta_{\min}$
is the minimum learning rate, which is set as $1\times 10^{-16}$,
$\eta_{\max}$ is the maximum or initial learning rate, which is selected as
$1\times 10^{-4}$, $\mathrm{step_{c}}$ is the current step, and
$\mathrm{step_{f}}$ is the final step. We note that each step refers to each
time we perform back-propagation, including updating both encoder and
projector’s parameters.

The inner loop is as follows: training AE starts with obtaining
$\bm{z}^{T}\left(t,\bm{\mu}\right)$ by passing $T_{h}\left(t,\bm{\mu}\right)$
to the encoder. We then use $\bm{z}^{T}\left(t,\bm{\mu}\right)$ to reconstruct
$\widehat{T}_{h}\left(t,\bm{\mu}\right)$ through the decoder. Subsequently, we
calculate our data compression loss or AE loss
($\mathcal{L}_{\mathrm{AE}}^{T}$) using

${\mathcal{L}_{\mathrm{AE}}^{T}}:={\mathrm{MSE}^{T}}=\frac{1}{\mathrm{M}N^{t}}\sum_{i=1}^{\mathrm{M}}\sum_{k=0}^{N^{t}}\left|\widehat{T}_{h}\left(t^{k},\bm{\mu}^{(i)}\right)-T_{h}\left(t^{k},\bm{\mu}^{(i)}\right)\right|^{2}.$
(18)

Similar to the training of BT, we use ADAM to adjust learnable parameters of
encoder($\mathrm{W}$ and $\mathrm{b}$) and decoder($\mathrm{W}$ and
$\mathrm{b}$) according the gradient of Equation (18). The $\eta_{c}$ is
adjusted by Equation (17). In contrast to the training of BT, we use
$\eta_{\min}=1\times 10^{-16}$, and $\eta_{\max}=1\times 10^{-5}$.

Following the training of the BT-AE, we now establish the manifold
$\bm{z}^{T}\left(t,\bm{\mu}\right),\quad\forall
t\in\mathbb{T}\>\text{and}\>\forall\bm{\mu}\in\mathbb{P}$ during the fourth
step shown in Figure 6. The data available for this task are the pairs
$(t,\bm{\mu})$ and $\bm{z}^{T}\left(t,\bm{\mu}\right)$ in the training set. We
achieve this through the training of artificial neural networks (ANN).
Following Kadeethum et al. [43, 6], our ANN has five hidden layers, and each
layer has seven neurons. We use tanh as our activation function. Here, we use
a mean squared error (${\mathrm{MSE}^{\bm{z}^{T}}}$) as the metric of our
network loss function, defined as follows

${\mathcal{L}_{\mathrm{ANN}}^{T}}:={\mathrm{MSE}^{\bm{z}^{T}}}=\frac{1}{\mathrm{M}N^{t}}\sum_{i=1}^{\mathrm{M}}\sum_{k=0}^{N^{t}}\left|\widehat{\bm{z}}^{T}\left(t^{k},\bm{\mu}^{(i)}\right)-{\bm{z}^{T}}\left(t^{k},\bm{\mu}^{(i)}\right)\right|^{2}.$
(19)

To minimize Equation (19), we use the ADAM algorithm to adjust each neuron
$\mathrm{W}$ and $\mathrm{b}$, a batch size of 32, a learning rate of 0.001, a
number of epoch of 10,000, and we normalize both our input ($t,\bm{\mu}$) and
output ($\bm{z}^{T}$) to $[0,1]$. To prevent our networks from overfitting
behavior, we follow early stopping and generalized cross-validation criteria
[4, 73, 74]. Note that instead of literally stopping our training cycle, we
only save the set of trained weight and bias to be used in the online phase
when the current validation loss is lower than the lowest validation from all
the previous training cycle.

During the online phase (the fifth step shown in Figure 6), we utilize the
trained ANN and the trained decoder to approximate
$\widehat{{T}}_{h}\left(\cdot;t,\bm{\mu}\right)$ for each inquiry (i.e., a
pair of $(t,\bm{\mu})$ ) through

$\widehat{\bm{z}}^{T}\left(\cdot;t,\bm{\mu}\right)=\operatorname{ANN}\left(t,\bm{\mu}\right),$
(20)

and, subsequently,

$\widehat{{T}}_{h}\left(\cdot;t,\bm{\mu}\right)=\operatorname{decoder}\left(\widehat{\bm{z}}^{T}\left(\cdot;t,\bm{\mu}\right)\right).$
(21)

We note that, for the prediction phase, our ROM could be evaluated using any
timestamps, including one that does not exist in the training phase (i.e., any
$t$ that lies within $[t^{0},\tau$]) because our ROM treats the time domain
$\mathbb{T}$ continuously. Besides, in contrast with the FOM, the ROM is not
bound by the CFL condition and can predict the quantities of interest at any
specific time without intermediate computation. Hence, our proposed framework
can reduce the computational time significantly.

## References

  * [1] Hesthaven, J., Rozza, G. & Stamm, B.  _Certified reduced basis methods for parametrized partial differential equations_ (Springer, 2016).
  * [2] Xiao, D., Fang, F., Pain, C. & Hu, G.  Non-intrusive reduced-order modelling of the Navier–Stokes equations based on rbf interpolation.  _International Journal for Numerical Methods in Fluids_ 79, 580–595 (2015).
  * [3] Xiao, D. _et al._ Non-intrusive reduced order modelling of the Navier–Stokes equations.  _Computer Methods in Applied Mechanics and Engineering_ 293, 522–541 (2015).
  * [4] Hesthaven, J. & Ubbiali, S.  Non-intrusive reduced order modeling of nonlinear problems using neural networks.  _Journal of Computational Physics_ 363, 55–78 (2018).
  * [5] Fresca, S., Dede, L. & Manzoni, A.  A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized PDEs.  _Journal of Scientific Computing_ 87, 1–36 (2021).
  * [6] Kadeethum, T. _et al._ Non-intrusive reduced order modeling of natural convection in porous media using convolutional autoencoders: comparison with linear subspace techniques.  _Advances in Water Resources_ 104098 (2022).
  * [7] Ahmed, S., San, O., Rasheed, A. & Iliescu, T.  Nonlinear proper orthogonal decomposition for convection-dominated flows.  _Physics of Fluids_ 33, 121702 (2021).
  * [8] Kim, Y., Choi, Y., Widemann, D. & Zohdi, T.  A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder.  _Journal of Computational Physics_ 110841 (2021).
  * [9] Kim, Y., Choi, Y., Widemann, D. & Zohdi, T.  Efficient nonlinear manifold reduced order model.  _arXiv preprint arXiv:2011.07727_ (2020).
  * [10] Chatterjee, A.  An introduction to the proper orthogonal decomposition.  _Current science_ 808–817 (2000).
  * [11] Willcox, K. & Peraire, J.  Balanced model reduction via the proper orthogonal decomposition.  _AIAA journal_ 40, 2323–2330 (2002).
  * [12] Choi, Y., Coombs, D. & Anderson, R.  SNS: a solution-based nonlinear subspace method for time-dependent model order reduction.  _SIAM Journal on Scientific Computing_ 42, A1116–A1146 (2020).
  * [13] Kim, Y., Wang, K. & Choi, Y.  Efficient space–time reduced order model for linear dynamical systems in python using less than 120 lines of code.  _Mathematics_ 9, 1690 (2021).
  * [14] Taron, J. & Elsworth, D.  Thermal-hydrologic-mechanical-chemical processes in the evolution of engineered geothermal reservoirs.  _International Journal of Rock Mechanics and Mining Sciences_ 46, 855–864 (2009).
  * [15] Nick, H., Raoof, A., Centler, F., Thullner, M. & Regnier, P.  Reactive dispersive contaminant transport in coastal aquifers: numerical simulation of a reactive henry problem.  _Journal of contaminant hydrology_ 145, 90–104 (2013).
  * [16] Zheng, C. & Bennett, G.  _Applied contaminant transport modeling_ , vol. 2 (Wiley-Interscience New York, 2002).
  * [17] Rutqvist, J. _et al._ A numerical study of THM effects on the near-field safety of a hypothetical nuclear waste repository—BMT1 of the DECOVALEX III project. part 3: Effects of THM coupling in sparsely fractured rocks.  _International Journal of Rock Mechanics and Mining Sciences_ 42, 745–755 (2005).
  * [18] Nield, D. & Bejan, A.  _Convection in porous media_ , vol. 3 (Springer, 2006).
  * [19] Park, S. W., Lee, J., Yoon, H. & Shin, S.  Microfluidic investigation of salinity-induced oil recovery in porous media during chemical flooding.  _Energy & Fuels_ 35, 4885–4892 (2021).
  * [20] Davison, S. M., Yoon, H. & Martinez, M. J.  Pore scale analysis of the impact of mixing-induced reaction dependent viscosity variations.  _Advances in water resources_ 38, 70–80 (2012).
  * [21] Pruess, K.  _TOUGH user’s guide_ (1987).
  * [22] Rutqvist, J.  An overview of TOUGH-based geomechanics models.  _Computers & Geosciences_ 108, 56–63 (2017).
  * [23] Bean, J., Sanchez, M. & Arguello, J.  Sierra mechanics, an emerging massively parallel hpc capability, for use in coupled thmc analyses of hlw repositories in clay/shale.  _5th International meeting Book of abstracts_ (2012).
  * [24] Aagaard, B., Williams, C. & Knepley, M.  PyLith: A finite-element code for modeling quasi-static and dynamic crustal deformation.  _Eos Trans. AGU_ 89 (2008).
  * [25] Kolditz, O. _et al._ OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media.  _Environmental Earth Sciences_ 67, 589–599 (2012).
  * [26] Obeysekara, A. _et al._ Modelling stress-dependent single and multi-phase flows in fractured porous media based on an immersed-body method with mesh adaptivity.  _Computers and Geotechnics_ 103, 229–241 (2018).
  * [27] Prévost, J. H.  Dynaflow.  _Princeton University, Princeton, NJ_ 8544 (1983).
  * [28] HosseiniMehr, M., Vuik, C. & Hajibeygi, H.  Adaptive dynamic multilevel simulation of fractured geothermal reservoirs.  _Journal of Computational Physics: X_ 100061 (2020).
  * [29] Matthai, S. _et al._ Numerical simulation of multi-phase fluid flow in structurally complex reservoirs.  _Geological Society, London, Special Publications_ 292, 405–429 (2007).
  * [30] Keilegavlen, E. _et al._ Porepy: An open-source software for simulation of multiphysics processes in fractured porous media.  _arXiv preprint arXiv:1908.09869_ (2019).
  * [31] Kadeethum, T., Lee, S. & Nick, H.  Finite element solvers for biot’s poroelasticity equations in porous media.  _Mathematical Geosciences_ 52, 977–1015 (2020).
  * [32] Kadeethum, T., Lee, S., Ballarin, F., Choo, J. & Nick, H.  A locally conservative mixed finite element framework for coupled hydro-mechanical-chemical processes in heterogeneous porous media.  _Computers & Geosciences_ 104774 (2021).
  * [33] Diersch, H.  Finite element modelling of recirculating density-driven saltwater intrusion processes in groundwater.  _Advances in Water Resources_ 11, 25–43 (1988).
  * [34] Frolkovič, P. & De Schepper, H.  Numerical modelling of convection dominated transport coupled with density driven flow in porous media.  _Advances in Water Resources_ 24, 63–72 (2000).
  * [35] Kolditz, O., Ratke, R., Diersch, H. & Zielke, W.  Coupled groundwater flow and transport: 1. verification of variable density flow and transport models.  _Advances in water resources_ 21, 27–46 (1998).
  * [36] Carlberg, K., Choi, Y. & Sargsyan, S.  Conservative model reduction for finite-volume models.  _Journal of Computational Physics_ 371, 280–314 (2018).
  * [37] Ballarin, F., D’amario, A., Perotto, S. & Rozza, G.  A POD-selective inverse distance weighting method for fast parametrized shape morphing.  _International Journal for Numerical Methods in Engineering_ 117, 860–884 (2019).
  * [38] Venturi, L., Ballarin, F. & Rozza, G.  A weighted POD method for elliptic PDEs with random inputs.  _Journal of Scientific Computing_ 81, 136–153 (2019).
  * [39] Choi, Y. & Carlberg, K.  Space–time least-squares petrov–galerkin projection for nonlinear model reduction.  _SIAM Journal on Scientific Computing_ 41, A26–A58 (2019).
  * [40] Copeland, D., Cheung, S., Huynh, K. & Choi, Y.  Reduced order models for lagrangian hydrodynamics.  _Computer Methods in Applied Mechanics and Engineering_ 388, 114259 (2022).
  * [41] Hoang, C., Choi, Y. & Carlberg, K.  Domain-decomposition least-squares petrov–galerkin (dd-lspg) nonlinear model reduction.  _Computer Methods in Applied Mechanics and Engineering_ 384, 113997 (2021).
  * [42] Choi, Y., Brown, P., Arrighi, W., Anderson, R. & Huynh, K.  Space–time reduced order model for large-scale linear dynamical systems with application to Boltzmann transport problems.  _Journal of Computational Physics_ 424, 109845 (2021).
  * [43] Kadeethum, T., Ballarin, F. & Bouklas, N.  Data-driven reduced order modeling of poroelasticity of heterogeneous media based on a discontinuous galerkin approximation.  _GEM-International Journal on Geomathematics_ 12, 1–45 (2021).
  * [44] DeCaria, V., Iliescu, T., Layton, W., McLaughlin, M. & Schneier, M.  An artificial compression reduced order model.  _SIAM Journal on Numerical Analysis_ 58, 565–589 (2020).
  * [45] Cleary, J. & Witten, I.  Data compression using adaptive coding and partial string matching.  _IEEE transactions on Communications_ 32, 396–402 (1984).
  * [46] Zbontar, J., Jing, L., Misra, I., LeCun, Y. & Deny, S.  Barlow twins: Self-supervised learning via redundancy reduction.  _arXiv preprint arXiv:2103.03230_ (2021).
  * [47] Zhang, C., Zarrouk, S. & Archer, R.  A mixed finite element solver for natural convection in porous media using automated solution techniques.  _Computers & Geosciences_ 96, 181–192 (2016).
  * [48] Elder, J.  Transient convection in a porous medium.  _Journal of Fluid Mechanics_ 27, 609–623 (1967).
  * [49] Simpson, M. & Clement, T.  Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models.  _Advances in Water Resources_ 26, 17–31 (2003).
  * [50] Diersch, H. & Kolditz, O.  Variable-density flow and transport in porous media: approaches and challenges.  _Advances in water resources_ 25, 899–944 (2002).
  * [51] Yoon, H., Valocchi, A., Werth, C. & Dewers, T.  Pore-scale simulation of mixing-induced calcium carbonate precipitation and dissolution in a microfluidic pore network.  _Water Resources Research_ 48 (2012).
  * [52] Yoon, H., Kang, Q. & Valocchi, A.  Lattice boltzmann-based approaches for pore-scale reactive transport.  _Reviews in Mineralogy and Geochemistry_ 80, 393–431 (2015).
  * [53] Yoon, H., Chojnicki, K. & Martinez, M.  Pore-scale analysis of calcium carbonate precipitation and dissolution kinetics in a microfluidic device.  _Environmental Science & Technology_ 53, 14233–14242 (2019).
  * [54] Yoon, H. _et al._ Adaptation of delftia acidovorans for degradation of 2, 4-dichlorophenoxyacetate in a microfluidic porous medium.  _Biodegradation_ 25, 595–604 (2014).
  * [55] Inspectorate, S. N. P.  The international hydrocoin project—background and results.  _Organization for Economic Co-operation and Development, Paris_ (1987).
  * [56] Flemisch, B. _et al._ Benchmarks for single-phase flow in fractured porous media.  _Advances in Water Resources_ 111, 239–258 (2018).
  * [57] Paul-Dubois-Taine, A. & Amsallem, D.  An adaptive and efficient greedy procedure for the optimal training of parametric reduced-order models.  _International Journal for Numerical Methods in Engineering_ 102, 1262–1292 (2015).
  * [58] Vasile, M. _et al._ Adaptive sampling strategies for non-intrusive pod-based surrogates.  _Engineering computations_ (2013).
  * [59] Choi, Y., Boncoraglio, G., Anderson, S., Amsallem, D. & Farhat, C.  Gradient-based constrained optimization using a database of linear reduced-order models.  _Journal of Computational Physics_ 423, 109787 (2020).
  * [60] Raissi, M., Perdikaris, P. & Karniadakis, G.  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).
  * [61] Kadeethum, T., Jørgensen, T. M. & Nick, H. M.  Physics-informed neural networks for solving nonlinear diffusivity and Biot’s equations.  _PloS one_ 15, e0232683 (2020).
  * [62] Huang, X., Liu, M., Belongie, S. & Kautz, J.  Multimodal unsupervised image-to-image translation.  In _Proceedings of the European conference on computer vision (ECCV)_ , 172–189 (2018).
  * [63] Xiao, D.  Error estimation of the parametric non-intrusive reduced order model using machine learning.  _Computer Methods in Applied Mechanics and Engineering_ 355, 513–534 (2019).
  * [64] Jacquier, P., Abdedou, A., Delmas, V. & Soulaïmani, A.  Non-intrusive reduced-order modeling using uncertainty-aware deep neural networks and proper orthogonal decomposition: Application to flood modeling.  _Journal of Computational Physics_ 424, 109854 (2021).
  * [65] Joseph, D.  _Stability of fluid motions I_ , vol. 27 (Springer Science & Business Media, 2013).
  * [66] Chen, Z., Huan, G. & Ma, Y.  _Computational methods for multiphase flows in porous media_ , vol. 2 (Siam, 2006).
  * [67] Lake, L., Johns, R., Rossen, B. & Pope, G.  _Fundamentals of enhanced oil recovery_ (Society of Petroleum Engineers Richardson, TX, 2014).
  * [68] Wang, H., Ren, K. & Song, J.  A closer look at batch size in mini-batch training of deep auto-encoders.  In _2017 3rd IEEE International Conference on Computer and Communications (ICCC)_ , 2756–2761 (IEEE, 2017).
  * [69] Karras, T., Aila, T., Laine, S. & Lehtinen, J.  Progressive growing of gans for improved quality, stability, and variation.  _arXiv preprint arXiv:1710.10196_ (2017).
  * [70] Kadeethum, T. _et al._ A framework for data-driven solution and parameter estimation of pdes using conditional generative adversarial networks.  _Nature Computational Science_ 1, 819–829, DOI: https://doi.org/10.1038/s43588-021-00171-3 (2021).
  * [71] Kingma, D. & Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * [72] Loshchilov, I. & Hutter, F.  Sgdr: Stochastic gradient descent with warm restarts.  _arXiv preprint arXiv:1608.03983_ (2016).
  * [73] Prechelt, L.  Early stopping-but when?  In _Neural Networks: Tricks of the trade_ , 55–69 (Springer, 1998).
  * [74] Prechelt, L.  Automatic early stopping using cross validation: quantifying the criteria.  _Neural Networks_ 11, 761–767 (1998).

## Acknowledgements

TK and HY were supported by the Laboratory Directed Research and Development
program at Sandia National Laboratories and US Department of Energy Office of
Fossil Energy and Carbon Management, Science-Informed Machine Learning to
Accelerate Real Time Decisions-Carbon Storage (SMART-CS) initiative. FB thanks
the project “Numerical modeling of flows in porous media” funded by the
Catholic University of the Sacred Heart, and the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie Actions,
grant agreement 872442 (ARIA). DO acknowledges support from Los Alamos
National Laboratory’s Laboratory Directed Research and Development Early
Career Award (20200575ECR). YC acknowledges LDRD funds (21-FS-042) from
Lawrence Livermore National Laboratory (LLNL-JRNL-831095). Lawrence Livermore
National Laboratory is operated by Lawrence Livermore National Security, LLC,
for the U.S. Department of Energy, National Nuclear Security Administration
under Contract DE-AC52-07NA27344. NB acknowledges startup support from the
Sibley School of Mechanical and Aerospace Engineering, Cornell University.
Sandia National Laboratories is a multimission laboratory managed and operated
by National Technology and Engineering Solutions of Sandia LLC, a wholly owned
subsidiary of Honeywell International Inc. for the U.S. Department of Energy’s
National Nuclear Security Administration under contract DE-NA0003525. This
paper describes objective technical results and analysis. Any subjective views
or opinions that might be expressed in the paper do not necessarily represent
the views of the U.S. Department of Energy or the United States Government.

## Author contributions statement

T. Kadeethum: Conceptualization, Formal analysis, Software, Validation,
Writing - original draft, Writing - review & editing. F. Ballarin:
Conceptualization, Formal analysis, Supervision, Validation, Writing - review
& editing. D. O’Malley: Conceptualization, Formal analysis, Supervision,
Validation, Writing - review & editing. Y. Choi: Conceptualization, Formal
analysis, Supervision, Validation, Writing - review & editing. N. Bouklas:
Conceptualization, Formal analysis, Funding acquisition, Supervision, Writing
- review & editing. H. Yoon: Conceptualization, Formal analysis, Funding
acquisition, Supervision, Writing - review & editing.

## Competing interests

The authors declare no competing interests

## Code and data availability

Our model scripts and all data generated or analyzed during this study will be
available publicly through the Sandia National Laboratories software portal —
a hub for GitHub-hosted open source projects (https://github.com/sandialabs).