textwidth=0.8textheight=0.9 [table]capposition=top

# Physics-informed Neural Networks with Periodic Activation Functions for
Solute Transport in Heterogeneous Porous Media

Salah A Faroughi  Pingki Datta  Seyed Kourosh Mahjour  Shirko Faroughi  Geo-
Intelligence Laboratory, Ingram School of Engineering, Texas State University,
San Marcos, Texas, 78666, USA  Department of Mechanical Engineering, School of
Engineering, Urmia University of Technology, Urmia, Iran

###### Abstract

Solute transport in porous media is relevant to a wide range of applications
in hydrogeology, geothermal energy, underground CO2 storage, and a variety of
chemical engineering systems. Due to the complexity of solute transport in
heterogeneous porous media, traditional solvers require high resolution
meshing and are therefore expensive computationally. This study explores the
application of a mesh-free method based on deep learning to accelerate the
simulation of solute transport. We employ Physics-informed Neural Networks
(PiNN) to solve solute transport problems in homogeneous and heterogeneous
porous media governed by the advection-dispersion equation. Unlike traditional
neural networks that learn from large training datasets, PiNNs only leverage
the strong form mathematical models to simultaneously solve for multiple
dependent or independent field variables (e.g., pressure and solute
concentration fields). In this study, we construct PiNN using a periodic
activation function to better represent the complex physical signals (i.e.,
pressure) and their derivatives (i.e., velocity). Several case studies are
designed with the intention of investigating the proposed PiNN’s capability to
handle different degrees of complexity. A manual hyperparameter tuning method
is used to find the best PiNN architecture for each test case. Point-wise
error and mean square error (MSE) measures are employed to assess the
performance of PiNNs’ predictions against the ground truth solutions obtained
analytically or numerically using the finite element method. Our findings show
that the predictions of PiNN are in good agreement with the ground truth
solutions while reducing computational complexity and cost by, at least, three
orders of magnitude.

###### keywords:

Physics-informed Neural Networks, Solute Transport, Heterogeneous Porous
Media, Advection-Dispersion Equation, Deep Learning, Scientific Computing

##  1 Introduction

Solute transport in porous media is crucial in many environmental and
reservoir engineering applications, including risk and safety assessment of
groundwater contamination [1], secondary recovery processes in petroleum
reservoirs [2], geological storage of radioactive waste [3], hydrogen [4], and
carbon dioxide [5], just to name a few [6]. It is a complicated process due to
the varied temporal and spatial transport characteristics of the flow network
[7]. These features can be divided into flowing areas influenced mostly by
advection and stagnant regions governed mainly by dispersion [8]. The process
is considerably more complicated in heterogeneous porous media, because all
the available pore space does not contribute to flow uniformly; certain
locations are dead ends, or transport is extremely sluggish due to an
inadequate connection to the main flow channels [9]. Despite a conceptual
understanding of the significance of multi-component solute transport
mechanisms in porous media in general, attempts to conduct reliable
experimental investigations are hindered by several barriers, resulting in a
dearth of published information. The major causes are the technical difficulty
of creating flow conditions and pore space connection under-regulated
laboratory or in-situ test circumstances [10]. Hence, a system of field
equations, which are numerically represented by partial differential equations
(PDEs), is practically the only tool to predict the behavior of such non-
linear flows in interconnected systems. In solute transport studies,
traditional solvers, including finite element (FEM), finite difference (FDM),
and finite volume methods (FVM), have been proposed to solve the governing
PDEs. For example, Zhang et al. [11] suggested a 1D model based on FEM for
multi-component solute transport in saturated soil. Bagalkot and Suresh Kumar
[12] presented a 1D numerical evaluation of the FDM for multi-species
radionuclide transport in a single horizontally coupled fracture-matrix
system. Mostaghimi et al. [13] and Maheshwari et al. [14] investigated the
effect of pore structure heterogeneity on reactive transport using a 3D pore
scale model based on FVM. Although most large-scale calculations still use
FEM, FDM, or FVM, they require high resolution meshing to describe the
geometry of the domain being represented. Because the size of an elementary
mesh is normally too large, discrete domain descriptions may also need the
availability of efficient transport equations at the mesh scale and the
consideration of unresolved subgrid-scale effects [15]. Hence, new
methodologies and algorithms, especially mesh-free methods based on deep
learning, are required to accelerate the numerical simulations of solute
transport through porous media.

Deep Learning (DL) methods have the potential to offer mesh-free solvers and
address some of the aforementioned challenges. While most applications
leverage DL (i.e., neural networks) to solve a lack of effective data modeling
processes, explore vast design domains, and identify multidimensional
connections [16, 17], there is a rising interest in using neural networks to
solve PDEs [18, 19, 20, 21]. In general, there are three main neural network
frameworks to accelerate scientific computing while considering underlying
physics: Physics-guided Neural Networks (PgNN), Physics-informed Neural
Networks (PiNN), and Physics-encoded Neural Networks (PeNN). Readers are
referred to a recent review by Faroughi et al. [22], where different neural
network frameworks are compared head-to-head, and their challenges and
limitations are thoroughly discussed. In this study, we adopt PiNN due to its
straightforward mechanism to integrate underlying physics in comparison to
other approaches. By combining a loss function comprised of the residuals of
physics equations, initial conditions and boundary restrictions, PiNN-based
models adhere to physical laws. They employ automated differentiation to
differentiate the outputs of neural networks with respect to their inputs
(i.e., spatiotemporal coordinates and model parameters) [23]. The network can
estimate the solution with high precision by minimizing the loss function
[24]. Therefore, PiNN provides the foundation for a mesh-free solver that
incorporates long-standing advances in mathematical physics into DL [18].

Raissi et al. [25] developed PiNN as a novel computing paradigm for both
forward and inverse modeling in a series of studies [18, 25, 26]. Raissi et
al. [26] developed a PiNN framework, dubbed hidden fluid mechanics (HFM), to
encode the physical laws governing fluid flows, i.e., the Navier-Stokes
equations. They leveraged underlying conservation laws to derive hidden
quantities of objective functions such as velocity and pressure fields from
spatiotemporal visualizations of a passive scalar concentration in arbitrarily
complex domains. Their technique accurately predicted 2D and 3D pressure and
velocity fields in benchmark problems inspired by real-world applications.
Since then, PiNN and its different variants have been aggressively applied to
different fields, including porous media flows. Almajid and Abu-Al-Saud [27]
employed the PiNN framework to solve both the forward and inverse problems of
the Buckley-Leverett PDE equation representing two-phase flow in porous media.
In order to test their implementation, they applied the classic problem of gas
drainage via a porous water-filled medium. Several cases were examined to
demonstrate the significance of the connectivity between observable data and
PiNNs for various parameter spaces. According to their results, PiNNs are
capable of capturing the solution’s broad trend even without observed data,
but their precision and accuracy improve significantly with observed data.
Hanna et al. [28] employed PiNN to simulate 1D and 2D two-phase flow in porous
media based on a new residual-based adaptive algorithm. Applying the PDE
residual to build a probability density function from which additional
collocation points are taken and added to the training set was fundamental to
their work. The approach was individually applied to each PDE in the coupled
system, considering the different collocation points for each PDE. In
addition, the method was applied to enrich the points used to capture the
initial and boundary conditions. They claimed that their approach yielded
superior results with less generalization error than the conventional PiNN
with fixed collocation points. Haghighat et al. [29] introduced a PiNN method
to solve the coupled flow and deformation equations in porous media for
single-phase and multiphase flow. Due to the problem’s dynamic nature, they
reported a dimensionless form of the couplled governing equations for the
optimizer. In addition, they presented a way for sequential training based on
the stress-split algorithms of poromechanics. They demonstrated that
sequential training based on stress-split performs well for a variety of
problems, whereas the conventional strain-split algorithm exhibits instability
comparable to that observed for FEM solvers. He et al. [30] extended a PiNN-
based parameter estimation method to integrate multiphysics measurement. They
studied a subsurface transport problem with sparse conductivity, a hydraulic
head, and solute concentration. In their methodology, they employed the Darcy
and advection-dispersion equations in conjunction with data to train deep
neural networks, reflecting space-dependent conductivity, head, and
concentration fields. They proved that the proposed PiNN method considerably
increased the accuracy of parameter and state estimates for sparse data when
compared to conventional deep neural networks trained with data alone. He and
Tartakovsky [31] also suggested a discretization-free technique based on PiNN
for solving the coupled advection-dispersion equations and Darcy flow equation
with space-dependent hydraulic conductivity. They used PiNN for both one-
dimensional and two-dimensional forward advection-dispersion equations and
compared its performance for various Peclet numbers (Pe) with analytical and
numerical solutions. They found that the PiNN was accurate with errors of less
than 1% and outperformed other discretization-based approaches for large Pe.
In addition, they proved that the PiNN remained accurate for backward
advection-dispersion equations, with relative errors below 5% in the majority
of instances. In another study, Vadyala et al. [32] implemented PiNN with the
use of a machine learning framework such that it can be employed in reduced-
order models to reduce epistemic (model-form) ambiguity associated with the
advection equation. They demonstrated that PiNN provided an accurate and
consistent approximation with the PDEs. Additionally, they showed that PiNN
could transform the physics simulation field by enabling real-time physics
simulation and geometry optimization without costly and time-consuming
simulations on large supercomputers.

Although there has been an increase in the application of PiNN to fluid flow
studies in porous media, the application of PiNN to the advection–dispersion
equation in heterogeneous porous media has received less attention. The main
objective of this study is to develop a PiNN model capable of predicting the
pressure and solute concentration fields simultaneously in heterogeneous
porous media. To this end, we construct PiNN using periodic activation
functions to better represent the complex physical signals and their
derivatives. Using this modification, the derivative of the network inherits
the network’s properties entirely, which enables supervising any derivative of
the network with complex physical signals, such as solute transport in
heterogeneous media. Several 1D and 2D test cases are examined to evaluate how
well the proposed PiNN model can predict the major transport fields.

This paper is structured as follows: In Section 2, the underlying physics for
solute transport in a porous medium is presented. Section 3 discusses the
PiNN’s algorithm with periodic activation functions to predict solute
transport. In Section 4, the PiNN is deployed to several one-dimensional (1D)
and two-dimensional (2D) case studies, and its predictions are examined by
comparing with analytical and/or FEM solutions. Finally, in Section 5, we
summarize the main conclusions of this work.

##  2 Underlying Physics

Solute transport in porous media is governed by three major mechanisms:
advection, diffusion, and dispersion [33, 34, 35]. The mathematical
formulation of the advection-diffusion-dispersion equation is explained
briefly here in terms of implementing the flow and transport factors that
affect the transport process. The mass conservation equation serves as the
basis for the solute transport equation [36, 37] that reads as,

$\nabla.\mathbf{J}=-\frac{\partial}{\partial{t}}(\phi{C}),$ (1)

where $\mathbf{J}$ is the solute mass flux, $\phi$ is the porosity, and $C$ is
the solute concentration. The flux $\mathbf{J}$ contains two mechanisms,

$\mathbf{J}=-\phi{D}\nabla{C}+\mathbf{u}(\phi{C}),$ (2)

where the first term refers to the diffusion flux, and the second term refers
to the advection flux [38, 39]. In Eq. (2), $D$ is the molecular diffusion
coefficient ($D=D_{0}$), and $\mathbf{u}=[u_{x},u_{y},u_{z}]$ is the velocity
vector of pore fluid flow computed using Darcy’s law,

$\mathbf{u}=\mathbf{q}/\phi=-\frac{k(x,y,z)}{\mu\phi}\nabla P,$ (3)

where $\mathbf{q}$ is the Darcy’s flux, $P$ refers to the pressure field,
$k(x,y,z)$ is the permeability, and $\mu$ is the fluid viscosity. To include
the impact of dispersion, the diffusion coefficient in Eq. 2 can be changed to

$D=D_{0}+\alpha{U},$ (4)

where $\alpha$ is known as the dynamic dispersivity, $D$ is now called the
hydrodynamic dispersion coefficient, and $U$ is the magnitude of fluid
velocity defined as $U=\sqrt{{u_{x}}^{2}+{u_{y}}^{2}+{u_{z}}^{2}}$. We allow
$U$ to vary spatially for a heterogeneous porous medium where $\mathbf{u}$
varies due to permeability changes [40, 41]. Substituting Eq. (2) into Eq. (1)
leads to

$\nabla.[\phi{D}\nabla{C}]-\nabla.[\mathbf{u}(\phi{C})]=\frac{\partial({\phi{C}})}{\partial{t}},$
(5)

which can be rewritten as,

$\nabla.[\phi{D}\nabla{C}]-(\nabla.\mathbf{u})[(\phi{C})]-\mathbf{u}.\nabla(\phi{C})=\frac{\partial({\phi{C}})}{\partial{t}},$
(6)

that incorporates diffusion, dispersion, and advection transports and serves
as the foundation of solute transport in porous media [42]. If the solute
transport process has no effect on fluid density (i.e., incompressible flow,
$\nabla.\mathbf{u}=0$), the velocity field can be calculated independently of
the solute concentration using Darcy’s law, and Eq. (6) reduces to,

$\nabla.[\phi{D}\nabla{C}]-\mathbf{u}.\nabla(\phi{C})=\frac{\partial({\phi{C}})}{\partial{t}}.$
(7)

The flow incompressibility condition generates an equation for computing the
pressure field as,

$\nabla.\mathbf{u}=\nabla.(-\frac{k(x,y,z)}{\mu\phi}\nabla P)=0,$ (8)

that reads as,

$\frac{\partial}{\partial{x}}[\zeta(x,y,z)\frac{\partial{P}}{\partial{x}}]+\frac{\partial}{\partial{y}}[\zeta(x,y,z)\frac{\partial{P}}{\partial{y}}]+\frac{\partial}{\partial{z}}[\zeta(x,y,z)\frac{\partial{P}}{\partial{z}}]=0,$
(9)

for three-dimensional (3D) flow problems, where

$\zeta(x,y,z)=\frac{k(x,y,z)}{\mu\phi}.$ (10)

Equations (7) and (9) together form the governing equations for the
concentration and pressure fields in a porous medium with a heterogeneous
permeability distribution [42].

##  3 Methodology

In this study, we leverage Physics-informed Neural Networks (PiNN) [18, 17] to
resolve the solute transport in porous media, which is governed by a strong
mathematical form consisting of the pressure and advection-dispersion
equations as well as the relevant initial and boundary conditions. Unlike
Physics-guided Neural Networks (PgNNs), PiNN is a mesh-free solution learning
method that does not need labeled datasets [22]. In PiNNs, the underlying
physics is included outside of the neural network architecture to constrain
the model during training and ensure that the outputs follow given physical
laws. The most common method to emulate this process is through a weakly
imposed penalty loss that penalizes the network when it does not follow the
physical constraints [25, 43, 44].

Figure 1: A schematic architecture of Physics-informed Neural Networks
(PiNNs). The network digests spatiotemporal coordinates, (x,t), as inputs to
predict a solution set, s, as an approximate to the ground truth solution,
$\hat{\textbf{{s}}}$. The automatic differentiation (AD) is then used to
generates the derivatives of the predicted solution s with respect to inputs.
These derivatives are used to formulate the residuals of the governing
equations in the loss function weighted by different coefficients. $\theta$
and $\lambda$ are the learnable parameters for weights/biases and unknown PDE
parameters, respectively, that can be learned simultaneously while minimizing
the loss function.

A schematic representation of our PiNN architecture is illustrated in Fig. 1.
As illustrated, a PiNN consists of three elements: a neural network (NN), an
automatic differentiation (AD) layer, and a feedback mechanism informed by
physics [45]. A neural network is first built in order to digest the
spatiotemporal features (i.e., x and $t$) as input parameters and approximate
the solution, $s(\mathbf{x},t)$, of a phenomenon described by PDEs with known
boundary values. Assume that
$\mathbb{N}_{L}{(\textbf{x})}:\mathbb{R}^{d_{in}}$ $\rightarrow$
$\mathbb{R}^{d_{out}}$ be a $L$-layer neural network, or a ($L-1$)-hidden
layer neural network, with $N_{l}$ neurons in the $l^{th}$ layer ($N_{0}$ =
$d_{in}$, $N_{L}$ = $d_{out}$) [46]. Each layer in the NN is associated with
its own weight matrix ${\textbf{W}^{l}}\in\mathbb{R}^{N_{l}\times N_{l-1}}$
and bias vector $\mathbf{b}^{l}\in\mathbb{R}^{N_{l}}$. The NN is recursively
defined using an element-wise nonlinear activation function, $\sigma$, as
follows, [47],

$\displaystyle\textnormal{Input\leavevmode\nobreak\
layer:}\leavevmode\nobreak\ \mathbb{N}_{0}(\mathbf{x})$
$\displaystyle=\mathbf{x}\in\mathbb{R}^{d_{in}},$ (11)
$\displaystyle\textnormal{Hidden\leavevmode\nobreak\
layers:}\leavevmode\nobreak\ \mathbb{N}_{l}(\mathbf{x})$
$\displaystyle=\sigma({\textbf{W}^{l}}\mathbb{N}_{l-1}(\textbf{x})+\mathbf{b}^{l})\in\mathbb{R}^{N_{l}},\leavevmode\nobreak\
\leavevmode\nobreak\ for\leavevmode\nobreak\ \leavevmode\nobreak\ 1\leq l\leq
L-1,$ (12) $\displaystyle\textnormal{Output\leavevmode\nobreak\
layer:}\leavevmode\nobreak\ \mathbb{N}_{L}(\mathbf{x})$
$\displaystyle={\textbf{W}^{l}}\mathbb{N}_{L-1}(\textbf{x})+\mathbf{b}^{l}\in\mathbb{R}^{d_{out}},$
(13)

where $\sigma$ is commonly specified as the logistic sigmoid
$\sigma(x)=1/(1-exp(-x))$, hyperbolic tangent
$\sigma(x)=tanh\leavevmode\nobreak\ (x)$, or the rectified linear unit,
$\sigma(x)=max\leavevmode\nobreak\ (x,0)$ [46]. In this work, we propose to
leverage periodic activation functions (e.g., sinusoidal),

$\sigma({\textbf{W}^{l}}\mathbb{N}_{l-1}(\textbf{x})+\mathbf{b}^{l})=\sin({\textbf{W}^{l}}\mathbb{N}_{l-1}(\textbf{x})+\mathbf{b}^{l}),$
(14)

that are ideally suited for representing complex physical signals and their
derivatives [48]. Using periodic activation functions, the derivative of a
network behaves like the network itself (i.e., the derivative of the sine is a
cosine that is a phase-shifted sine). As a result, the derivatives of the
$\mathbb{N}_{L}{(\textbf{x})}:\mathbb{R}^{d_{in}}$ $\rightarrow$
$\mathbb{R}^{d_{out}}$ inherit its properties that enable supervising any
derivative of the network with complex physical signals, including boundary
value problems [48]. In addition, it has been shown that NN with periodic
activation functions converges significantly faster than baseline
architectures [48].

The outputs of the NN are then fed into the automated differentiation (AD)
layer, which is the central property of PiNNs. AD is employed to assess the
derivatives of the network outputs with respect to the network inputs [49, 46,
23, 50]. Consider a strong mathematical form with a PDE specified on the
domain $\Omega$ and parametrzied by $\bm{\lambda}$,

$f\big{(}\textbf{x},t:\Delta s;\bm{\lambda}\big{)}=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \mathbf{x}\in\Omega,$ (15)

and boundary conditions (e.g., Dirichlet, Neumann, Robin boundary condition,
etc.) specified on the boundary of the domain $\partial\Omega$,

$\psi(\mathbf{s}(\mathbf{x},t))=0,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \mathbf{x}\in\partial\Omega,$ (16)

where $\mathbf{s}$ is the unknown solution, and $\Delta$ represents the linear
or nonlinear differential operator (e.g., $\frac{\partial}{\partial{t}}$,
$\frac{\partial}{\partial{x}}$, $\frac{\partial^{2}}{\partial{x^{2}}}$, etc.).
The initial condition can be considered as a kind of Dirichlet boundary
condition on the spatiotemporal domain. AD is used to apply $\Delta$ and
$\psi$ operators on the neural network with respect to inputs, i.e.,
$\mathbf{x}$ and $t$, and generate the required terms in the loss function in
order to optimize the PDE solution, $\mathbf{s}$. Lastly, a feedback mechanism
is constructed to minimize the loss terms through optimizations. The total
loss term is defined as [51, 50],

$\mathcal{L}={w_{i}}\mathcal{L}_{IC}+{w_{b}}\mathcal{L}_{BC}+{w_{d}}\mathcal{L}_{Data}+{w_{p}}\mathcal{L}_{PDE},$
(17)

where $w_{i}$, $w_{b}$, $w_{d}$, and $w_{p}$, are referred to as the weights
for the loss due to the initial conditions, boundary conditions, labeled data,
if any, and PDEs, respectively. The individual loss terms are computed as,

$\displaystyle\mathcal{L}_{IC}$
$\displaystyle=\frac{1}{N_{i}}\sum_{i=1}^{N_{i}}(\hat{s}{|_{\Omega,t_{0}}}-s{|_{\Omega,t_{0}}})^{2},$
(18) $\displaystyle\mathcal{L}_{BC}$
$\displaystyle=\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}(({\partial_{n}}\hat{s}{|_{\partial{\Omega}}}-{\partial_{n}}s{|_{\partial{\Omega}}})-(\hat{s}{|_{\partial{\Omega}}}-s{|_{\partial{\Omega}}}))^{2},$
(19) $\displaystyle\mathcal{L}_{Data}$
$\displaystyle=\frac{1}{N_{d}}\sum_{i=1}^{N_{d}}(\hat{s}{|_{\Omega}}-s{|_{Data}})^{2},$
(20) $\displaystyle\mathcal{L}_{PDE}$
$\displaystyle=\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}(\emph{f}\leavevmode\nobreak\
({s},\leavevmode\nobreak\ {{\partial}_{t}}{s},\leavevmode\nobreak\
{{\partial}_{x}}{s},\leavevmode\nobreak\ \ldots,\leavevmode\nobreak\
\lambda))^{2},$ (21)

where $t_{0}$ is the initial time, $s$ is the neural network prediction,
$\hat{s}$ is the ground truth solution, $N_{i}$, $N_{b}$, $N_{d}$, and $N_{p}$
are, respectively, the number of spatiotemporal points representing the
initial conditions, boundary conditions, labeled data, and collocation points
(i.e., spatiotemporal points within the domain where the neural network
prediction, $s(\textbf{x},t)$ is checked against the constraints of PDEs).
Finally, a gradient-based optimizer (e.g., the Limited-memory
Broyden–Fletcher–Goldfarb–Shanno algorithm, L-BFGS) is employed to minimize
the loss function. L-BFGS is a second-order optimizer that incorporates the
Hessian matrix [52]. It is observed that L-BFGS can find a satisfactory
solution for smooth PDE solutions faster than Adam and with fewer iterations
[53, 54, 55]. At each iteration, the loss is checked against
$\mathcal{L}<\epsilon$, where $\epsilon$ is a specified tolerance. If the
condition is not met, error backpropagation is implemented to update the
learnable parameters ($\theta$ and/or $\lambda$). The entire cycle is repeated
for a given number of iterations until the PiNN model produces learnable
parameters with a loss error less than $\epsilon$ [56]. PiNN training is more
difficult than PgNN training, because PiNNs are composed of sophisticated non-
convex and multi-objective loss functions, which may cause instability during
optimization [57, 58].

We deployed the proposed PiNN to model solute transport in porous media
predicting the pressure and concentration fields, i.e., $\textbf{s}=(P,C)$. We
selected different PiNN architectures for the solute transport case studies in
this paper because of the different PDEs and boundary conditions. However, we
ensured that the selected PiNN had enough layers and width to estimate each
target field. There are two critical normalization steps in the implementation
of PiNNs that must be followed to ensure convergence to the correct solution
[59]. The first is mapping the network’s input and output variables to the
interval $[0,1]\in{R}$. The second step is to scale the pressure and
concentration PDEs such that all terms are in the same order (unity).

##  4 Computational Experiments

In this section, we applied the proposed PiNN model to solve 1D and 2D solute
transport phenomena under different conditions (i.e., a total of four
computational experiments are explored). We validated the PiNN’s prediction
against analytical and/or FEM solutions to examine the potential, capability,
and accuracy of the PiNN model.

###  4.1 Case 1: 1D solute transport with constant velocity

This test case explores solute transport in a 1D domain ($x\in[0,L]$)
representing an isotropic porous medium with a constant velocity field. As the
steady-state velocity field is known, the governing equation, Eq. (7), reduced
to,

$\frac{\partial{C}}{\partial{t}}+{u_{x}}\frac{\partial{C}}{\partial{x}}=\frac{\partial}{\partial{x}}(D_{x}\frac{\partial{C}}{\partial{x}}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
0<x<L,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ 0<t<t_{0},$ (22)

where $L=1$ is the length of the domain, $u_{x}=0.5\leavevmode\nobreak\ m/s$
is the steady state velocity field, and $D_{x}=0.02\leavevmode\nobreak\
m^{2}/s$ refers to the hydrodynamic dispersion coefficient in the $x$
direction. For this case, we assume the following initial and boundary
conditions,

$\begin{cases}C(x,t)=0;\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ t=0,\\\ C(x,t)=C_{0};\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ x=0,\\\
\frac{\partial{}}{\partial{x}}C(x,t)=0;\leavevmode\nobreak\
\leavevmode\nobreak\ x=L.\end{cases}$ (23)

where $C_{0}=1.0\leavevmode\nobreak\ kg/m$ is the injected concentration at
$x=0$ (i.e., Dirichlet boundary condition). This 1D advection–dispersion can
be solved analytically [60]. The analytical solution is,

$C(x,t)=C_{0}\left(1-2\exp(\frac{xu_{x}}{2D_{x}}-\frac{{u_{x}^{2}}t}{4D_{x}})\sum_{i=1}^{\infty}\frac{{\beta_{i}}\sin({\frac{\beta_{i}{x}}{L}})\exp(-\frac{{\beta_{i}^{2}}D_{x}{t}}{L^{2}})}{{\beta_{i}^{2}}+{({\frac{{u_{x}}L}{2D_{x}}})^{2}}+\frac{{u_{x}}L}{2D_{x}}}\right),$
(24)

where $\beta_{i}$ are the roots of,

$\beta\cot{\beta}+\frac{{u_{x}}L}{2D_{x}}=0.$ (25)

Equations (24) and (25) can be approximated using [36],

$C(x,t)=C_{0}\leavevmode\nobreak\ A(x,t),\leavevmode\nobreak\
\leavevmode\nobreak\ 0<t\leq t_{0},$ (26)

where $A(x,t)$ is computed as,

$\displaystyle A(x,t)$
$\displaystyle=\frac{1}{2}\operatorname{erfc}\left(\frac{x-u_{x}{t}}{2\sqrt{{D_{x}}t}}\right)\;+\;\frac{1}{2}\exp(\frac{u_{x}{x}}{D_{x}})\;\operatorname{erfc}\left(\frac{x+{u_{x}}t}{2\sqrt{{D_{x}}t}}\right)$
(27)
$\displaystyle+\frac{1}{2}\left(2+\frac{{u_{x}}(2L-x)}{{D_{x}}}+\frac{{u_{x}^{2}}t}{{D_{x}}}\right)\;\exp(\frac{{u_{x}}L}{D_{x}})\operatorname{erfc}\left(\frac{2L-x+{u_{x}}t}{2\sqrt{{D_{x}}t}}\right)$
$\displaystyle-({\frac{{u_{x}^{2}}t}{\pi{{D_{x}}}}})^{\frac{1}{2}}\exp\left(\frac{{u_{x}}L}{D_{x}}-\frac{1}{4D_{x}t}(2L-x+{u_{x}}t)^{2}\right).$

Figure 2: One-dimensional solute transport with a constant velocity field,
$u_{x}=0.5\leavevmode\nobreak\ m/s$. The upper panel shows the concentration
field predicted by PiNN within the spatiotemporal domain $(0\leq x\leq
1,\leavevmode\nobreak\ 0\leq t\leq 10)$ at an injection rate of $C_{0}=1kg/m$.
The lower panels show the comparison of PiNN’s prediction with the ground
truth obtained analytically using Eq. (26) at three different times. The
comparison demonstrates a good agreement yielding $MSE=1.15\times 10^{-6}$ for
the entire spatiotemporal domain.

A PiNN with the periodic activation function and randomly distributed
collocation points is then used to predict the solute transport in 1D
described by Eq. (22) with the initial and boundary conditions given in Eq.
(23). A manual hyperparameter tuning process was practiced to find the best
architecture based on the MSE accuracy measures. This process yielded a
network with four hidden layers with {32, 32, 16, 16} neurons per layer, and
5,000 spatiotemporal random points within $0\leq x\leq
1and\leavevmode\nobreak\ 0\leq t\leq 10$ as the collocation points to enforce
the PDE loss term. We used identical weights for all the loss terms
($w_{i}=w_{b}=w_{p}=1$), and used 5,000 spatiotemporal random points,
collectively, to represent the boundary and initial conditions informing the
loss terms.

Figure 2 represents the concentration field in the ($\mathbf{x},t$) domain
predicted by PiNN. The PiNN’s prediction is also compared with the ground
truth (analytical solution given by Eq. (27) at different simulation times.
The results of the comparison are shown in Fig. 2, which reveal a good
agreement between PiNN’s prediction and the ground truth resulting in
$MSE=1.15\times 10^{-6}$ for the considered spatiotemporal domain. This
confirms the accuracy of the trained PiNN model for this simplified problem.
In the next test cases, the PiNN is applied to more complex problems
representing the complexity of solute transport in porous media.

###  4.2 Case 2: 2D solute transport with constant velocity

This test case explores solute transport in a 2D rectangular domain
($x\in[0,L],\leavevmode\nobreak\ y\in[0,W]$) representing an isotropic porous
medium with a constant velocity field, $\textbf{u}=(u_{x},u_{y})$. Since the
velocity field is given, the governing equation for solute transport, Eq. (7),
is simplified to,

$\frac{\partial{C}}{\partial{t}}+{u_{x}}\frac{\partial{C}}{\partial{x}}+{u_{y}}\frac{\partial{C}}{\partial{y}}=\frac{\partial}{\partial{x}}(D_{x}\frac{\partial{C}}{\partial{x}})+\frac{\partial}{\partial{y}}(D_{y}\frac{\partial{C}}{\partial{y}}),$
(28)

where $D_{y}=D_{x}=0.02\leavevmode\nobreak\ m^{2}/s$ refers to the
hydrodynamic dispersion coefficient in the $y$ and $x$ directions. As
schematically shown in Fig. 3, the following initial and boundary conditions,

$\begin{cases}C(x,y,t)=0;\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ t=0,\\\ C(x,y,t)=C_{0};\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ x=0\leavevmode\nobreak\
\leavevmode\nobreak\ and\leavevmode\nobreak\ \leavevmode\nobreak\
{y_{1}}\leq{y}\leq{y_{2}},\\\ C(x,y,t)=0;\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ x=0\leavevmode\nobreak\
\leavevmode\nobreak\ and\leavevmode\nobreak\ \leavevmode\nobreak\
y<{y_{1}}\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ y>{y_{2}},\\\
\frac{\partial{}}{\partial{x}}C(x,y,t)=0;\leavevmode\nobreak\
\leavevmode\nobreak\ x=L,\leavevmode\nobreak\ \leavevmode\nobreak\
and\leavevmode\nobreak\ \leavevmode\nobreak\ 0\leq y\leq W,\\\
\frac{\partial{}}{\partial{y}}C(x,y,t)=0;\leavevmode\nobreak\
\leavevmode\nobreak\ y=0\leavevmode\nobreak\ \leavevmode\nobreak\
or\leavevmode\nobreak\ \leavevmode\nobreak\ y=W,\leavevmode\nobreak\
\leavevmode\nobreak\ and\leavevmode\nobreak\ \leavevmode\nobreak\ 0\leq x\leq
L,\end{cases}$ (29)

are considered for this test case. The blue points on the left boundary
bounded by $y_{1}$ and $y_{2}$ represent the positions of point-source solute
injection at the rate of $C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$. The other
three sides of the domain are assigned as zero gradient boundaries, assuming a
free outflow of solute. This 2D advection–dispersion problem defined by Eqs.
(28) and (29) can be solved analytically [60],

$C(x,y,t)=C_{0}\sum_{n=0}^{\infty}{L_{n}}{P_{n}}\cos{({\eta}y)}\left(\exp(\frac{x({u_{x}}-\zeta)}{2{D_{x}}})\operatorname{erfc}(\frac{x-{\zeta}t}{2\sqrt{{D_{x}}t}})+\exp(\frac{x({u_{x}}+\zeta)}{2{D_{x}}})\operatorname{erfc}(\frac{x+{\zeta}t}{2\sqrt{{D_{x}}t}})\right),$
(30)

where $L_{n}$ is defined as,

${L_{n}}=\begin{cases}\frac{1}{2};\leavevmode\nobreak\ \leavevmode\nobreak\
n=0,\\\ 1;\leavevmode\nobreak\ \leavevmode\nobreak\ n>0,\end{cases}$ (31)

and $P_{n}$ is computed as,

${P_{n}}=\begin{cases}({y_{2}}-{y_{1}})/W;\leavevmode\nobreak\
\leavevmode\nobreak\ n=0,\\\
(\sin{(\eta{y_{2}})}-\sin{(\eta{y_{1}})})/(n\pi);\leavevmode\nobreak\
\leavevmode\nobreak\ n>0,\end{cases}$ (32)

where $\eta$ and $\zeta$, respectively, are,

${\eta}=\frac{n{\pi}}{W}\leavevmode\nobreak\ \leavevmode\nobreak\
(n=0,1,2,3,\ldots),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \zeta=\sqrt{{{u_{x}}^{2}}+4{{\eta}^{2}}{D_{x}}{D_{y}}}.$
(33)

We first employed a PiNN with the periodic activation function and randomly
distributed collocation points to predict the solute transport in the 2D
domain ($L=W=1\leavevmode\nobreak\ m$) considering dispersion only with
$D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$ and
$\textbf{u}=(u_{x},u_{y})=(0.0,0.0)\leavevmode\nobreak\ m/s$ in Eq. (28) and
the initial and boundary conditions given by Eq. (29). The solute injection
rate is set to $C_{0}=0.2kg/m^{2}$ between $y_{1}=0.3\leavevmode\nobreak\ m$
and $y_{2}=0.7\leavevmode\nobreak\ m$. To determine the best PiNN
architecture, a manual hyperparameter tuning approach was adopted. This
approach produced a network with four hidden layers and {32,16,16,16} neurons
per layer using 8,000 spatiotemporal random points within $0\leq x,y\leq
1,\leavevmode\nobreak\ 0\leq t\leq 1$ as collocation points to enforce the PDE
loss term. We employed identical weights for the loss terms and utilized 8,000
spatiotemporal random points to represent the boundary and initial conditions
informing the loss terms (sample selected points are shown in Fig. 3).

Figure 3 depicts the comparison of the 2D concentration field predicted by
PiNN and the analytical solution given by Eq. (30). The point-wise error is
also shown, which illustrates a good accordance between the analytic solution
and PiNN’s prediction yielding $MSE=0.67\times 10^{-6}$. It can be observed
that the mismatch is close to zero everywhere in the domain except locations
where extremely high concentration gradients exist (e.g., at points
$(0.0,0.30)$ and $(0.0,0.70)$). This mismatch is due to a difficult-to-
minimize optimization error, i.e., PiNN struggles to converge to ground truth
solutions in high gradient areas [61].

Figure 3: A comparison between PiNN’s prediction and ground truth for solute
transport in a 2D domain representing an isotropic porous medium considering
$D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$,
$\textbf{u}=(u_{x},u_{y})=(0.0,0.0)\leavevmode\nobreak\ m/s$, and solute
injection rate of $C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$ between
$y_{1}=0.3\leavevmode\nobreak\ m$ and $y_{2}=0.7\leavevmode\nobreak\ m$. The
upper panels show the domain with boundary conditions and distributions of
randomly selected points on which different terms in the loss function are
evaluated. The red dots represent the collocation points inside the domain
corresponding to the loss term associated with the 2D solute transport PDE,
and the blue, green, and black dots represent the points on the boundary of
the domain corresponding to the loss terms associated with the boundary
conditions. The lower panel shows the comparison between the PiNN’s prediction
of the concentration field at $t=1.00\leavevmode\nobreak\ s$ and the ground
truth obtained analytically using Eq. (30). The point-wise error shows the
mismatch between the solutions yielding $MSE=0.67\times 10^{-6}$. Figure 4: A
comparison between PiNN’s prediction and ground truth obtained analytically
for the concentration field considering $D_{x}=D_{y}=0.02\leavevmode\nobreak\
m^{2}/s$, $\textbf{u}=(u_{x},u_{y})=(0.5,0.0)\leavevmode\nobreak\ m/s$, and an
injection rate of $C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$ between
$y_{1}=0.3\leavevmode\nobreak\ m$ and $y_{2}=0.7\leavevmode\nobreak\ m$. PiNN
predictions are shown for $t=0.25\leavevmode\nobreak\
s,0.50\leavevmode\nobreak\ s,0.75\leavevmode\nobreak\ s$, and
$1.00\leavevmode\nobreak\ s$. At $t=1.00\leavevmode\nobreak\ s$, the PiNN’s
prediction is compared with the analytical solution obtained using Eq. (30).
The point-wise error shows the mismatch between the solutions yielding
$MSE=1.54\times 10^{-6}$.

A PiNN with the same architecture is then employed to model the
advection–dispersion solute transport in 2D domains with
$D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$ and
$\textbf{u}=(u_{x},u_{y})=(0.5,0.0)\leavevmode\nobreak\ m/s$. The solute is
injected again at a rate of $C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$ between
$y_{1}=0.3\leavevmode\nobreak\ m$ and $y_{2}=0.7\leavevmode\nobreak\ m$. The
PiNN’s predictions for the concentration field at different times
($t=0.25\leavevmode\nobreak\ s,\leavevmode\nobreak\ 0.50\leavevmode\nobreak\
s,\leavevmode\nobreak\ 0.75\leavevmode\nobreak\ s,\leavevmode\nobreak\
1.00\leavevmode\nobreak\ s$) are shown in Fig. 4. The concentration field at
$t=1.00\leavevmode\nobreak\ s$ is also compared with the analytical solution
given by Eq. (28). The point-wise error shows that there is a good match
between the PiNN’s prediction and the analytical solution within the entire
domain yielding $MSE=1.54\times 10^{-6}$.

The same PiNN model is also deployed to predict the solute transport in the 2D
domain discussed above, but with line source injection (i.e.,
$y_{1}=0.0\leavevmode\nobreak\ m$ and $y_{2}=1.0\leavevmode\nobreak\ m$). A
comparison between the PiNN’s prediction and analytical solution for the
concentration field at $t=0.75\leavevmode\nobreak\ s$ is shown in Fig. 5. As
imposed by the boundary conditions, the contours are straight lines parallel
to the $x=0$ boundary, which is predicted correctly by PiNN. The point-wise
error shown in Fig. 5 demonstrates excellent agreement between the PiNN’s
prediction and the ground truth yielding $MSE=0.12\times 10^{-6}$ that clearly
signifies the PiNN model’s capability to solve advection–dispersion equations.

Figure 5:  A comparison between PiNN’s prediction and ground truth obtained
analytically for the concentration field considering
$D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$,
$\textbf{u}=[u_{x},u_{y}]=[0.5,0.0]\leavevmode\nobreak\ m/s$, and an injection
rate of $C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$ between
$y_{1}=0.0\leavevmode\nobreak\ m$ and $y_{2}=1.0\leavevmode\nobreak\ m$ (i.e.,
line-source injection). At $t=0.75\leavevmode\nobreak\ s$, the PiNN’s
prediction is validated against the analytical solution obtained using Eq.
(30). The point-wise error shows the mismatch between the solutions yielding
$MSE=0.12\times 10^{-6}$.

###  4.3 Case 3: 2D solute transport in homogeneous porous media

In this test case, the solute transport in a 2D homogeneous porous medium
($x\in[0,L],\leavevmode\nobreak\ y\in[0,W]$) is investigated. For this
problem, assuming $\zeta(x,y)=1$, the pressure equation defined by Eq. (9)
reduces to,

$\frac{\partial^{2}{P}}{\partial{x}^{2}}+\frac{\partial^{2}{P}}{\partial{y}^{2}}=0,$
(34)

that must be solved to determine the velocity field using Eq. (3). The
governing equation for solute concentration remains the same as Eq. (28). The
following boundary conditions are considered for the pressure field,

$\begin{cases}P(x,y)=1.0;\leavevmode\nobreak\ \leavevmode\nobreak\
y=0\leavevmode\nobreak\ \leavevmode\nobreak\ or\leavevmode\nobreak\
\leavevmode\nobreak\ y=W,\leavevmode\nobreak\ \leavevmode\nobreak\
and\leavevmode\nobreak\ \leavevmode\nobreak\ 0\leq x\leq L,\\\
P(x,y)=0.1;\leavevmode\nobreak\ \leavevmode\nobreak\ x=0\leavevmode\nobreak\
\leavevmode\nobreak\ or\leavevmode\nobreak\ \leavevmode\nobreak\
x=L,\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ 0\leq y\leq W,\end{cases}$ (35)

and for the concentration field,

$\begin{cases}C(x,y,t)=0;\leavevmode\nobreak\ \leavevmode\nobreak\
y=0\leavevmode\nobreak\ \leavevmode\nobreak\ or\leavevmode\nobreak\
\leavevmode\nobreak\ y=W\leavevmode\nobreak\ \leavevmode\nobreak\
and\leavevmode\nobreak\ \leavevmode\nobreak\ x<{x_{1}}\leavevmode\nobreak\
\leavevmode\nobreak\ and\leavevmode\nobreak\ \leavevmode\nobreak\
x>{x_{2}},\\\ C(x,y,t)=C_{0};\leavevmode\nobreak\ \leavevmode\nobreak\
y=0\leavevmode\nobreak\ \leavevmode\nobreak\ or\leavevmode\nobreak\
\leavevmode\nobreak\ y=W\leavevmode\nobreak\ \leavevmode\nobreak\
and\leavevmode\nobreak\ \leavevmode\nobreak\ {x_{1}}\leq{x}\leq{x_{2}},\\\
\frac{\partial{}}{\partial{x}}C(x,y,t)=0;\leavevmode\nobreak\
\leavevmode\nobreak\ x=0\leavevmode\nobreak\ \leavevmode\nobreak\
or\leavevmode\nobreak\ \leavevmode\nobreak\ x=L,\leavevmode\nobreak\
\leavevmode\nobreak\ and\leavevmode\nobreak\ \leavevmode\nobreak\ 0\leq y\leq
W,\end{cases}$ (36)

with $C(x,y,t=0)=0$ as the initial condition. The pressure value at the corner
points that belong to two perpendicular sides is set to $0.55$, as additional
constraints, to keep the symmetry. The ground truth solutions for pressure and
concentration fields are obtained using FEM with a $100\times 100$
quadrilateral elements assuming $L=W=1\leavevmode\nobreak\ m$,
$D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$, and
$C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$.

A PiNN with the periodic activation function and randomly distributed
collocation points is employed to predict the solute transport in the 2D
porous media described above. The PiNN’s inputs are the spatiotemporal
coordinates, $(x,y,t)$, and the outputs are the pressure and concentration
fields, i.e., $\textbf{s}=\\{P,C\\}$ in Fig. 1. The pressure and concentration
fields are decoupled, therefore one may use two PiNNs side-by-side to solve
this problem; however, training a PiNN with multiple outputs was more
desirable in this study to minimize hyperparameter tuning. The hyperparameter
tuning approach was again practiced manually to obtain the best PiNN
architecture based on the MSE accuracy measures. This approach resulted in a
network with five hidden layers and {32, 16, 16, 8, 8} neurons per layer using
12,000 spatiotemporal random points within $0\leq x,y\leq
1,\leavevmode\nobreak\ 0\leq t\leq 1$ as collocation points to enforce the
pressure and concentration PDEs’ loss terms. We employed identical weights for
the loss terms and 12,000 spatiotemporal random points to represent the
boundary and initial conditions informing the loss terms (sample selected
points are shown in Fig. 6).

Figure 6 presents the ground truth (FEM solutions) and the PiNN’s predictions
for the pressure field, velocity field in the $x$ direction, total velocity
field, and flow streamlines for the homogeneous porous medium. The point-wise
error is also shown for all the fields, which shows a good agreement between
PiNN and FEM solutions. Figure 7 shows the same comparison, but for the
concentration field at $t=1.00\leavevmode\nobreak\ s$. The point-wise error
again shows a considerable agreement between the PiNN’s prediction and the FEM
solution yielding $MSE=0.24\times 10^{-6}$. These results imply that a PiNN
can simultaneously solve for the pressure and solute concentration fields in a
2D homogeneous porous medium with a high accuracy. The following section
explores the application of PiNN to solve the problem of solute transport in a
heterogeneous porous medium.

Figure 6: A comparison between PiNN’s predictions and the ground truth for the
pressure field in a 2D homogeneous porous medium. The upper panels show the
domain with pressure and concentration boundary conditions as well as the
distributions of randomly selected points on which different terms in the loss
function are evaluated. The lower panels show the comparison between the
PiNN’s predictions and ground truth solutions obtained using FEM with
$100\times 100$ quadrilateral elements for the pressure field, velocity field
in the $x$ direction, and total velocity field. The flow direction of the pore
fluid is also represented by the streamlines. The point-wise error is also
shown for each field to illustrate the mismatch between the solutions. Figure
7: A comparison between the PiNN’s prediction and ground truth for the
concentration field at $t=1.00\leavevmode\nobreak\ s$ in a 2D homogeneous
porous medium considering $D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$ and
$C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$. The ground truth is obtained using
FEM with $100\times 100$ quadrilateral elements. The point-wise error shows
the mismatch between the solutions yielding $MSE=0.24\times 10^{-6}$.

###  4.4 Case 4: 2D solute transport in heterogeneous porous media

In this test case, the solute transport is analyzed in a 2D rectangular domain
($x\in[0,L],\leavevmode\nobreak\ y\in[0,W]$) representing a heterogeneous
porous medium, i.e., the permeability field $k(x,y)$, and hence, $\zeta(x,y)$
in Eq. (9) vary spatially. The permeability field considered here is shown in
Fig. 9. The initial and boundary conditions for this case are identical to
those reported for Case 3. The pressure constraints at the corner points that
correspond to the two perpendicular sides are not maintained due to
permeability. The ground truth solutions for pressure and concentration fields
are determined using FEM with $100\times 100$ quadrilateral elements
considering $L=W=1\leavevmode\nobreak\ m$, $\alpha=0$,
$D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$, $\mu\phi=0.001$, and
$C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$.

Figure 8: A comparison between the PiNN’s predictions and ground truth
solutions for solute transport in a heterogeneous porous medium with the given
permeability field, $k(x,y)$. The comparison is illustrated for the pressure
field as well as the concentration field at $t=1.00\leavevmode\nobreak\ s$
considering $D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$, $\alpha=0$,
$\mu\phi=0.001$, and $C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$. The PiNN model
here uses the tanh activation function, and the ground truth solutions are
obtained using FEM with $100\times 100$ quadrilateral elements. The point-wise
error for each field demonstrates that the PiNN has difficulty to converge to
the ground truth solution when using the tanh activation function. The
disparity is more pronounced in the concentration field, which is caused by
the rapid variation of the velocity field caused by heterogeneities. Figure 9:
A comparison between the PiNN’s predictions and ground truth solutions for
solute transport in a heterogeneous porous medium with the given permeability
field, $k(x,y)$. The comparison is illustrated for the pressure field as well
as the concentration field at $t=1.00\leavevmode\nobreak\ s$ considering
$D_{x}=D_{y}=0.02\leavevmode\nobreak\ m^{2}/s$, $\alpha=0$, $\mu\phi=0.001$,
and $C_{0}=0.2\leavevmode\nobreak\ kg/m^{2}$. The PiNN model here uses the sin
activation function, and the ground truth solutions are obtained using FEM
with $100\times 100$ quadrilateral elements. The point-wise error for each
field shows that the PiNN, with the Sin activation function, can accurately
predict the complex pressure and concentration fields in heterogeneous porous
media. The comparison yields $MSE=1.19\times 10^{-6}$ and $MSE=0.45\times
10^{-6}$ for the pressure field and concentration fields.

A PiNN model with randomly distributed collocation points is deployed to
simultaneously predict both the pressure and concentration fields. The best
PiNN architecture was determined through a manual hyperparameter tuning
approach. This procedure resulted in a network with five hidden layers and
{32, 32, 16, 16, 16} neurons per layer using 15,000 spatiotemporal random
points within $0\leq x,\leavevmode\nobreak\ y\leq 1,\leavevmode\nobreak\ 0\leq
t\leq 1$ as collocation points to enforce the PDE loss terms. We also utilized
15,000 spatiotemporal random points to represent the boundary and initial
conditions informing the loss terms. To achieve convergence while meeting the
boundary condition requirements satisfactorily, we employed $\omega_{PDE}=0.7$
as the weights for the PDE loss terms.

Figures 8 and 9 compare the PiNN’s predictions with the ground truth (FEM
solutions) obtained for the solute transport in the heterogeneous porous
medium. The comparison is shown for the steady pressure field and the
concentration field at $t=1.00\leavevmode\nobreak\ s$. Figure 8 shows the
PiNN’s prediction using the tanh activation function validated against the FEM
solution with $100\times 100$ quadrilateral elements. The point-wise error for
each field shows that the PiNN, with the tanh activation function, has
difficulty to converge to the ground truth solution. The discrepancy is more
pronounced in the concentration field, which is originated by the rapid
variation of velocity fields due to heterogeneities.

Figure 9 shows the PiNN’s prediction using the sin activation function
validated against FEM solutions with $100\times 100$ quadrilateral elements.
The point-wise error is also shown for both fields, which shows a considerable
accuracy improvement compared to the PiNN with the tanh activation function.
The maximum point-wise error for the pressure field is reduced to $\approx
0.02$ from $\approx 0.49$. Also, the maximum point-wise error for the
concentration field is reduced to $\approx 0.01$ from $\approx 0.13$. The
comparison of the PiNN with the periodic activation function and FEM solutions
leads to $MSE=1.19\times 10^{-6}$ and $MSE=0.45\times 10^{-6}$ for the
pressure field and concentration field, respectively. As illustrated in Fig.
9, the majority of mismatch in this case still occurs near high gradient
areas. These comparative findings show that a PiNN with a periodic activation
function can simultaneously solve for the pressure and solute concentration
fields in 2D heterogeneous media with a higher degree of accuracy compared to
other typical activation function. The results also show that PiNN, in
general, decreases the computational complexity tremendously for solute
transport phenomena in complex porous media.

Table 1: A comparison of the inference time for PiNN with the sin activation function and FEM resolving the solute transport in the 2D homogeneous (Case 3) and heterogeneous (Case 4) porous media. Data is obtained using a 3.00-GHz 48-core Intel Xeon Gold 6248R CPU processor with 128 GB of RAM. |  FEM |  PiNN |  Speed-up Factor
---|---|---|---
|  (s) |  (s) |  (—)
2D Homogeneous |  $183.8$ |  $0.126$ |  1458.7x
2D Heterogeneous |  $363.2$ |  $0.259$ |  1402.3x

Table 1 compares the inference time for PiNN with the sin activation function
and FEM resolving the solute transport in the 2D homogeneous (Case 3) and
heterogeneous (Case 4) porous media. The speed-up factor achieved is around
three orders of magnitude and can considerably escalate when a higher-
resolution mesh (i.e., a higher number of elements) is used by the FEM solver;
a situation that is often faced in cases with higher levels of
heterogeneities. Indeed, in such cases, the number of collection points used
by PiNN also increases, which could negatively impact the cost of the training
but not the inference time. Therefore, the computational efficiency of PiNN is
certainly superior to that of a traditional solver (e.g., FEM) in domains with
complex geometry and heterogeneity. It is also notable that PiNN comes with
several drawbacks. The two most common drawbacks that were also faced in this
study are: (i) PiNN’s learning process and hyperparameter tuning were
conducted manually. The training can face gradient vanishing problems, which
prohibitively slow down the training phase; and (ii) PiNN did not always
converge due to competing non-linear loss terms (e.g., PDE loss terms related
to the pressure and concentration fields in addition to the loss terms related
to initial/boundary conditions for both of the fields). It requires trials and
errors to adjust the loss terms’ weight functions to mitigate the instability.
The lack of a theoretical condition or constraint (e.g., Courant number in
traditional computational fluid dynamics [62]) to assure convergence is an
open research area for investigation. All in all, considering the
effectiveness of PiNN in combining scientific computing and DL as well as
accelerating computation, the addition of PiNNs to the simulations of porous
media flow should be investigated further. For example, a future study could
be the assessment of PiNNs to model non-isothermal reactive flows in highly
heterogeneous and deformable porous media.

##  5 Conclusions

In this study, we presented the application of physics-informed neural
networks (PiNNs) to solve the solute transport problem in homogeneous and
heterogeneous porous media governed by the advection-dispersion equation. We
constructed PiNN using a periodic activation function to accurately predict
the complex pressure and concentration fields as well as reduce the
computational cost. Using this modification, the derivative of the network
inherits the network’s properties entirely, which enables supervising any
derivative of the network with complex physical signals. In order to evaluate
the capabilities of the proposed PiNN, we built a number of case studies in
both 1D and 2D, each of which possessed homogeneous and/or heterogeneous
permeability fields. Using a manual hyperparameter tuning approach, the
optimal PiNN architecture for each test case was determined. The point-wise
error and mean square error (MSE) metrics were used to evaluate the accuracy
of PiNN’s predictions in comparison to the ground truth solutions derived
analytically or numerically using FEM. Our findings demonstrated that the
proposed PiNN can accurately predict the solute transport phenomenon under
different conditions. In addition, the results showed that PiNN’s simultaneous
predictions of the pressure and concentration fields in heterogeneous porous
media are in good agreement with the ground truth solutions while lowering the
computational complexity and expense by three-orders of magnitude.

##  6 Acknowledgements

S.A.F. would like to acknowledge the support by Texas Sate University’s
International Research Accelerator Grant (award no. 9000003039).

##  7 Conflict of Interest

The authors declare no conflict of interests.

## References

  * Cook and Böhlke [2000] Peter G Cook and John-Karl Böhlke.  Determining timescales for groundwater flow and solute transport.  In _Environmental tracers in subsurface hydrology_ , pages 1–30. Springer, 2000.  DOI: https://doi.org/10.1007/978-1-4615-4557-6˙1.
  * Gaus et al. [2008] Irina Gaus, Pascal Audigane, Laurent André, Julie Lions, Nicolas Jacquemet, Pierre Durst, Isabelle Czernichowski-Lauriol, and Mohamed Azaroual.  Geochemical and solute transport modelling for co2 storage, what to expect from it?  _International journal of greenhouse gas control_ , 2(4):605–625, 2008.  DOI: https://doi.org/10.1016/j.ijggc.2008.02.011.
  * Pruess et al. [2002] Karsten Pruess, Steve Yabusaki, Carl Steefel, and Peter Lichtner.  Fluid flow, heat transfer, and solute transport at nuclear waste storage tanks in the hanford vadose zone.  _Vadose Zone Journal_ , 1(1):68–88, 2002.  DOI: https://doi.org/10.2136/vzj2002.6800.
  * Bienert et al. [2006] Gerd P Bienert, Jan K Schjoerring, and Thomas P Jahn.  Membrane transport of hydrogen peroxide.  _Biochimica et Biophysica Acta (BBA)-Biomembranes_ , 1758(8):994–1003, 2006.  DOI: https://doi.org/10.1016/j.bbamem.2006.02.015.
  * Kristensen and Hansen [1999] Erik Kristensen and Kim Hansen.  Transport of carbon dioxide and ammonium in bioturbated (nereis diversicolor) coastal, marine sediments.  _Biogeochemistry_ , 45(2):147–168, 1999.  DOI: https://doi.org/10.1023/A:1006131710076.
  * Li et al. [2020] Xinxin Li, Dianqing Li, Yi Xu, and Xiaobo Feng.  A dfn based 3d numerical approach for modeling coupled groundwater flow and solute transport in fractured rock mass.  _International Journal of Heat and Mass Transfer_ , 149:119179, 2020.  DOI: https://doi.org/10.1016/j.ijheatmasstransfer.2019.119179.
  * Hasan et al. [2020] Sharul Hasan, Vahid Niasar, Nikolaos K Karadimitriou, Jose RA Godinho, Nghia T Vo, Senyou An, Arash Rabbani, and Holger Steeb.  Direct characterization of solute transport in unsaturated porous media using fast x-ray synchrotron microtomography.  _Proceedings of the National Academy of Sciences_ , 117(38):23443–23449, 2020.  DOI: https://doi.org/10.1073/pnas.2011716117.
  * Faraji and Mazaheri [2022] Mojtaba Faraji and Mehdi Mazaheri.  Mathematical model of solute transport in rivers with storage zones using nonlinear dispersion flux approach.  _Hydrological Sciences Journal_ , pages 1–13, 2022.  DOI: https://doi.org/10.1080/02626667.2022.2099280.
  * Yang and Wang [2019] Xiao-Rong Yang and Yan Wang.  Ubiquity of anomalous transport in porous media: Numerical evidence, continuous time random walk modelling, and hydrodynamic interpretation.  _Scientific reports_ , 9(1):1–11, 2019.  DOI: https://doi.org/10.1038/s41598-019-39363-3.
  * Zhao et al. [2011] Zhihong Zhao, Lanru Jing, Ivars Neretnieks, and Luis Moreno.  Numerical modeling of stress effects on solute transport in fractured rocks.  _Computers and Geotechnics_ , 38(2):113–126, 2011\.  DOI: https://doi.org/10.1016/j.compgeo.2010.10.001.
  * Zhang et al. [2018] Zhi-hong Zhang, Jia-pei Zhang, Zhan-ying Ju, and Min Zhu.  A one-dimensional transport model for multi-component solute in saturated soil.  _Water Science and Engineering_ , 11(3):236–242, 2018.  DOI: https://doi.org/10.1016/j.wse.2018.09.007.
  * Bagalkot and Suresh Kumar [2015] Nikhil Bagalkot and G Suresh Kumar.  Effect of nonlinear sorption on multispecies radionuclide transport in a coupled fracture-matrix system with variable fracture aperture: a numerical study.  _ISH Journal of Hydraulic Engineering_ , 21(3):242–254, 2015.  DOI: https://doi.org/10.1080/09715010.2015.1016125.
  * Mostaghimi et al. [2016] Peyman Mostaghimi, Min Liu, and Christoph H Arns.  Numerical simulation of reactive transport on micro-ct images.  _Mathematical Geosciences_ , 48(8):963–983, 2016\.  DOI: https://doi.org/10.1007/s11004-016-9640-3.
  * Maheshwari et al. [2013] P Maheshwari, RR Ratnakar, N Kalia, and V Balakotaiah.  3-d simulation and analysis of reactive dissolution and wormhole formation in carbonate rocks.  _Chemical Engineering Science_ , 90:258–274, 2013.  DOI: https://doi.org/10.1016/j.ces.2012.12.032.
  * Noetinger et al. [2016] Benoit Noetinger, Delphine Roubinet, Anna Russian, Tanguy Le Borgne, Frederick Delay, Marco Dentz, Jean-Raynald De Dreuzy, and Philippe Gouze.  Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale.  _Transport in Porous Media_ , 115(2):345–385, 2016.  DOI: https://doi.org/10.1007/s11242-016-0693-z.
  * Im et al. [2021] Sunyoung Im, Jonggeon Lee, and Maenghyo Cho.  Surrogate modeling of elasto-plastic problems via long short-term memory neural networks and proper orthogonal decomposition.  _Computer Methods in Applied Mechanics and Engineering_ , 385:114030, 2021.  DOI: https://doi.org/10.1016/j.cma.2021.114030.
  * Karniadakis et al. [2021] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  _Nature Reviews Physics_ , 3(6):422–440, 2021\.  DOI: https://doi.org/10.1038/s42254-021-00314-5.
  * Raissi et al. [2019] Maziar Raissi, Paris Perdikaris, and George E 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\.  DOI: https://doi.org/10.1016/j.jcp.2018.10.045.
  * Han et al. [2018] Jiequn Han, Arnulf Jentzen, and Weinan E.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ , 115(34):8505–8510, 2018.  DOI: https://doi.org/10.1073/pnas.1718942115.
  * Berg and Nyström [2018] Jens Berg and Kaj Nyström.  A unified deep artificial neural network approach to partial differential equations in complex geometries.  _Neurocomputing_ , 317:28–41, 2018.  DOI: https://doi.org/10.1016/j.neucom.2018.06.056.
  * Sirignano and Spiliopoulos [2018] Justin Sirignano and Konstantinos Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  _Journal of computational physics_ , 375:1339–1364, 2018\.  DOI: https://doi.org/10.1016/j.jcp.2018.08.029.
  * Faroughi et al. [2022] Salah A Faroughi, Nikhil Pawar, Celio Fernandes, Subasish Das, Nima K Kalantari, and Seyed Kourosh Mahjour.  Physics-guided, physics-informed, and physics-encoded neural networks in scientific computing.  _arXiv preprint arXiv:2211.07377_ , 2022.  DOI: https://doi.org/10.48550/arXiv.2211.07377.
  * Van Merriënboer et al. [2018] Bart Van Merriënboer, Olivier Breuleux, Arnaud Bergeron, and Pascal Lamblin.  Automatic differentiation in ml: Where we are and where we should be going.  _Advances in neural information processing systems_ , 31, 2018.  DOI: https://doi.org/10.48550/arXiv.1810.11530.
  * Baydin et al. [2018] Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.  Automatic differentiation in machine learning: a survey.  _Journal of Marchine Learning Research_ , 18:1–43, 2018\.  DOI: https://doi.org/10.48550/arXiv.1502.05767.
  * Raissi et al. [2017] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations.  _arXiv preprint arXiv:1711.10561_ , 2017.  DOI: https://doi.org/10.48550/arXiv.1711.10561.
  * Raissi et al. [2018] Maziar Raissi, Alireza Yazdani, and George Em Karniadakis.  Hidden fluid mechanics: A navier-stokes informed deep learning framework for assimilating flow visualization data.  _arXiv preprint arXiv:1808.04327_ , 2018.  DOI: https://doi.org/10.48550/arXiv.1808.04327.
  * Almajid and Abu-Al-Saud [2022] Muhammad M Almajid and Moataz O Abu-Al-Saud.  Prediction of porous media fluid flow using physics informed neural networks.  _Journal of Petroleum Science and Engineering_ , 208:109205, 2022.  DOI: https://doi.org/10.1016/j.petrol.2021.109205.
  * Hanna et al. [2022] John M Hanna, Jose V Aguado, Sebastien Comas-Cardona, Ramzi Askri, and Domenico Borzacchiello.  Residual-based adaptivity for two-phase flow simulation in porous media using physics-informed neural networks.  _Computer Methods in Applied Mechanics and Engineering_ , 396:115100, 2022.  DOI: https://doi.org/10.1016/j.cma.2022.115100.
  * Haghighat et al. [2022] Ehsan Haghighat, Danial Amini, and Ruben Juanes.  Physics-informed neural network simulation of multiphase poroelasticity using stress-split sequential training.  _Computer Methods in Applied Mechanics and Engineering_ , 397:115141, 2022.  DOI: https://doi.org/10.1016/j.cma.2022.115141.
  * He et al. [2020] QiZhi He, David Barajas-Solano, Guzel Tartakovsky, and Alexandre M Tartakovsky.  Physics-informed neural networks for multiphysics data assimilation with application to subsurface transport.  _Advances in Water Resources_ , 141:103610, 2020.  DOI: https://doi.org/10.1016/j.advwatres.2020.103610.
  * He and Tartakovsky [2021] QiZhi He and Alexandre M Tartakovsky.  Physics-informed neural network method for forward and backward advection-dispersion equations.  _Water Resources Research_ , 57(7):e2020WR029479, 2021.  DOI: https://doi.org/10.1029/2020WR029479.
  * Vadyala et al. [2022] Shashank Reddy Vadyala, Sai Nethra Betgeri, and Naga Parameshwari Betgeri.  Physics-informed neural network method for solving one-dimensional advection equation using pytorch.  _Array_ , 13:100110, 2022.  DOI: https://doi.org/10.1016/j.array.2021.100110.
  * Zhang et al. [2021] Chunwei Zhang, Kazuki Kaito, Yingxue Hu, Anindityo Patmonoaji, Shintaro Matsushita, and Tetsuya Suekane.  Influence of stagnant zones on solute transport in heterogeneous porous media at the pore scale.  _Physics of Fluids_ , 33(3):036605, 2021.  DOI: https://doi.org/10.1063/5.0038133.
  * Khan et al. [2022] Shahzad Khan, Sharifah E Alhazmi, Fahad M Alotaibi, Massimiliano Ferrara, and Ali Ahmadian.  On the numerical approximation of mobile-immobile advection-dispersion model of fractional order arising from solute transport in porous media.  _Fractal and Fractional_ , 6(8):445, 2022.  DOI: https://doi.org/10.3390/fractalfract6080445.
  * Zhao and Toksoz [1994] Xiaomin Zhao and M Nafi Toksoz.  Solute transport in heterogeneous porous media.  _Massachusetts Institute of Technology. Earth Resources Laboratory_ , 1994.  DOI: https://doi.org/10.1017/S0022112084002858.
  * Van Genuchten and Alves [1982] M Th Van Genuchten and WJ Alves.  Analytical solutions of the one-dimensional convective-dispersive solute transport equation.  _United States. Dept. of Agriculture. Technical bulletin (USA)_ , 1982\.  DOI: https://doi.org/10.22004/ag.econ.157268.
  * Sun et al. [2020] Liwei Sun, Han Qiu, Chuanhao Wu, Jie Niu, and Bill X Hu.  A review of applications of fractional advection–dispersion equations for anomalous solute transport in surface and subsurface water.  _Wiley Interdisciplinary Reviews: Water_ , 7(4):e1448, 2020.  DOI: https://doi.org/10.1002/wat2.1448.
  * Haigh et al. [2021] Michael Haigh, Luolin Sun, James C McWilliams, and Pavel Berloff.  On eddy transport in the ocean. part ii: The advection tensor.  _Ocean Modelling_ , 165:101845, 2021.  DOI: https://doi.org/10.1016/j.ocemod.2021.101845.
  * Lou et al. [2020] Shuai Lou, Shu-sheng Chen, Bo-xi Lin, Jian Yu, and Chao Yan.  Effective high-order energy stable flux reconstruction methods for first-order hyperbolic linear and nonlinear systems.  _Journal of Computational Physics_ , 414:109475, 2020.  DOI: https://doi.org/10.1016/j.jcp.2020.109475.
  * Talon [2022] Laurent Talon.  On the statistical properties of fluid flows with transitional power-law rheology in heterogeneous porous media.  _Journal of Non-Newtonian Fluid Mechanics_ , 304:104789, 2022.  DOI: https://doi.org/10.1016/j.jnnfm.2022.104789.
  * Baioni et al. [2021] Elisa Baioni, Mohaddeseh Mousavi Nezhad, Giovanni Michele Porta, and Alberto Guadagnini.  Modeling solute transport and mixing in heterogeneous porous media under turbulent flow conditions.  _Physics of Fluids_ , 33(10):106604, 2021.  DOI: https://doi.org/10.1063/5.0065734.
  * Berkowitz et al. [2002] Brian Berkowitz, Joseph Klafter, Ralf Metzler, and Harvey Scher.  Physical pictures of transport in heterogeneous media: Advection-dispersion, random-walk, and fractional derivative formulations.  _Water Resources Research_ , 38(10):9–1, 2002\.  DOI: https://doi.org/10.1029/2001WR001030.
  * Jagtap et al. [2020] Ameya D Jagtap, Ehsan Kharazmi, and George Em Karniadakis.  Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems.  _Computer Methods in Applied Mechanics and Engineering_ , 365:113028, 2020.  DOI: https://doi.org/10.1016/j.cma.2020.113028.
  * Jagtap and Karniadakis [2021] Ameya D Jagtap and George E Karniadakis.  Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations.  In _AAAI Spring Symposium: MLPS_ , 2021.  DOI: https://doi.org/10.4208/cicp.OA-2020-0164.
  * Cuomo et al. [2022a] Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maizar Raissi, and Francesco Piccialli.  Scientific machine learning through physics-informed neural networks: Where we are and what’s next.  _arXiv preprint arXiv:2201.05624_ , 2022a.  DOI: https://doi.org/10.48550/arXiv.2201.05624.
  * Depina et al. [2022] Ivan Depina, Saket Jain, Sigurdur Mar Valsson, and Hrvoje Gotovac.  Application of physics-informed neural networks to inverse problems in unsaturated groundwater flow.  _Georisk: Assessment and Management of Risk for Engineered Systems and Geohazards_ , 16(1):21–36, 2022.  DOI: https://doi.org/10.1080/17499518.2021.1971251.
  * Lu et al. [2021] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis.  Deepxde: A deep learning library for solving differential equations.  _SIAM Review_ , 63(1):208–228, 2021.  DOI: https://doi.org/10.1137/19M1274067.
  * Sitzmann et al. [2020] Vincent Sitzmann, Julien N. P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions, 2020\.  URL https://arxiv.org/abs/2006.09661.
  * Rumelhart et al. [1986] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  _nature_ , 323(6088):533–536, 1986.  DOI: https://doi.org/10.1038/323533a0.
  * Guo et al. [2020] Yanan Guo, Xiaoqun Cao, Bainian Liu, and Mei Gao.  Solving partial differential equations using deep learning and physical constraints.  _Applied Sciences_ , 10(17):5917, 2020.  DOI: https://doi.org/10.3390/app10175917.
  * Cai et al. [2022] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks (pinns) for fluid mechanics: A review.  _Acta Mechanica Sinica_ , pages 1–12, 2022.  DOI: https://doi.org/10.1007/s10409-021-01148-1.
  * Liu and Nocedal [1989] Dong C Liu and Jorge Nocedal.  On the limited memory bfgs method for large scale optimization.  _Mathematical programming_ , 45(1):503–528, 1989\.  DOI: https://doi.org/10.1007/BF01589116.
  * Bengio [2000] Yoshua Bengio.  Gradient-based optimization of hyperparameters.  _Neural computation_ , 12(8):1889–1900, 2000\.  DOI: https://doi.org/10.1162/089976600300015187.
  * Kylasa et al. [2019] Sudhir Kylasa, Fred Roosta, Michael W Mahoney, and Ananth Grama.  Gpu accelerated sub-sampled newton’s method for convex classification problems.  In _Proceedings of the 2019 SIAM International Conference on Data Mining_ , pages 702–710. SIAM, 2019.  DOI: https://doi.org/10.1137/1.9781611975673.79.
  * Richardson [2018] Alan Richardson.  Seismic full-waveform inversion using deep learning tools and techniques.  _arXiv preprint arXiv:1801.07232_ , 2018.  DOI: https://doi.org/10.48550/arXiv.1801.07232.
  * Olmo et al. [2022] Alberto Olmo, Ahmed Zamzam, Andrew Glaws, and Ryan King.  Physics-driven convolutional autoencoder approach for cfd data compressions.  _arXiv preprint arXiv:2210.09262_ , 2022.  DOI: https://doi.org/10.48550/arXiv.2210.09262.
  * Cuomo et al. [2022b] Salvatore Cuomo, Fabio Giampaolo, Stefano Izzo, Carlo Nitsch, Francesco Piccialli, and Cristina Trombetti.  A physics-informed learning approach to bernoulli-type free boundary problems.  _Computers & Mathematics with Applications_, 128:34–43, 2022b.  DOI: https://doi.org/10.1016/j.camwa.2022.10.003.
  * Shah et al. [2022] Karan Shah, Patrick Stiller, Nico Hoffmann, and Attila Cangi.  Physics-informed neural networks as solvers for the time-dependent schr$\backslash$” odinger equation.  _arXiv preprint arXiv:2210.12522_ , 2022.  DOI: https://doi.org/10.48550/arXiv.2210.12522.
  * Rasht-Behesht et al. [2022] Majid Rasht-Behesht, Christian Huber, Khemraj Shukla, and George Em Karniadakis.  Physics-informed neural networks (pinns) for wave propagation and full waveform inversions.  _Journal of Geophysical Research: Solid Earth_ , 127(5):e2021JB023120, 2022.  DOI: https://doi.org/10.1029/2021JB023120.
  * Zhou [2009] Jian Guo Zhou.  A lattice boltzmann method for solute transport.  _International journal for numerical methods in fluids_ , 61(8):848–863, 2009.  DOI: https://doi.org/10.1002/fld.1978.
  * Chiu et al. [2022] Pao-Hsiung Chiu, Jian Cheng Wong, Chinchun Ooi, My Ha Dao, and Yew-Soon Ong.  Can-pinn: A fast physics-informed neural network based on coupled-automatic–numerical differentiation method.  _Computer Methods in Applied Mechanics and Engineering_ , 395:114909, 2022.  DOI: https://doi.org/10.1016/j.cma.2022.114909.
  * Atmakidis and Kenig [2010] Theodoros Atmakidis and Eugeny Y Kenig.  A study on the kelvin-helmholtz instability using two different computational fluid dynamics methods.  _The Journal of Computational Multiphase Flows_ , 2(1):33–45, 2010.  DOI: https://doi.org/10.1260/1757-482X.2.1.33.


# Virtual Control Group: Measuring Hidden Performance Metrics

Moshe Tocker  Meta

(June 2022)

##  1 ABSTRACT

Performance metrics measuring in Financial Integrity systems are crucial for
maintaining an efficient and cost effective operation. An important
performance metric is False Positive Rate. This metric cannot be directly
monitored since we don’t know for sure if a user is bad once blocked. We
present a statistical method based on survey theory and causal inference
methods to estimate the false positive rate of the system or a single blocking
policy. We also suggest a new approach of outcome matching that in some cases
including empirical data outperformed other commonly used methods. The
approaches described in this paper can be applied in other Integrity domains
such as Cyber Security.

Keywords— causal inference, financial risk, outcome score matching,
performance metrics, control group, false positive rate

##  2 INTRODUCTION

Measuring performance metrics in Financial Integrity Systems is crucial for
decision making. Integrity systems need to hold a high level of precision and
block bad users to reduce bad activity and comply with regulation and
commercial requirements. On the other hand, risk systems that block many good
users provide a terrible user experience which affects growth and adoption of
the services.

To monitor False Positives one can use a control group or periodically send
transactions for manual review. Control Groups cost is high as they introduce
bad activity into the system. In addition, for Integrity systems in many cases
this approach is not possible due to product or regulatory constraints. Manual
labeling is another approach to estimate false positive rate. This method
however requires manpower and is extremely noisy (fraud labeling is a high-
dimensional problem).

Another shortcoming of the approaches described above is measurement time
delay. In the financial domain the most common indications of fraud are
chargebacks. A chargeback is a return of money to a payer of a transaction,
especially a credit card transaction after a successful dispute. This process
takes time and therefore the indication of fraud takes time to arrive. This
effect introduces a delay in the metrics measurement.

In this paper we present an alternative approach based on survey theory and
causal inference methods that is purely statistical and has many advantages in
terms of operational costs and time to measurement delay. We were able to
predict the fraud prevalence among the blocked transactions with acceptable
errors in several use cases.

We believe that this approach can be useful in many integrity domains such as
Financial Services, Cyber security (Malware detection, Spam detection, …),
Content Integrity and more.

##  3 PRELIMINARIES

Assume we have two functions of covariates vector x: Outcome Y(x) and
Treatment assignment model W(x). In the financial integrity setting we control
the treatment assignment as we define the risk policies based on system risk
signals (in contrast to the normal causal inference settings in which some
treatment covariates are unknown). Units are sent to treatment based on the
binary function W(x). W(x) can model a single fraud policy or the entire risk
policies as a whole. The outcome function Y(x) is known only for units i of
which W($x_{i}$)=0. We seek to estimate the mean value of the outcome of the
treated group $E[Y|W=1]$. This problem can be viewed as estimating population
average of the non respondents from in a survey which has been dealt with by
many authors such as [8] and [5]. In the setting of survey theory surveys are
sent to a sampled population. Some of the people in the sample might not
respond to the survey (a.k.a non-respondents) and thus the survey outcome
would be unknown for them. Estimating the survey results just by the outcome
of the respondents might lead to a biased result. We suggest that estimating
false positive rate in financial risk is an equivalent problem in which the
outcome is fraud or not-fraud and the non-respondents are the blocked
transactions for which the outcome is unknown.

###  3.1 Estimators

The fields of Survey theory and causal inference produced statistical methods
to try and replicate randomized experiments when only observational data is
available. We describe in this section some of the methods of which we have
considered. The method tries to estimate W(x), Y(x) or both based on the
observed covariates and then draw conclusions over the mean of an unobserved
outcome of the treated group.

Many estimation methods are based on the Propensity Score which is the
probability of a unit to be treated and defined by [8].

$\pi(x)=P(W=1|X=x)$ (1)

Propensity Score matching is a well-known approach used widely in practice.
Each unit’s propensity score is estimated. Then units from the treated group
are matched to units from the untreated based on their nearest neighbor in
terms of the propensity score. The matching can be 1:1 also called pair
matching or it can be k-nearest neighbor matching in which every item in the
treated group is matched to k-nearest neighbors with replacement in the
untreated group. Let $\mathcal{I}_{\pi}$ be the set of units in the untreated
group which are the k-nearest neighbors by the propensity score of some unit
in the treated group. Let $y_{i}$ be the true outcome of unit i. The set
$\mathcal{I}_{\pi}$ is then used to estimate the mean outcome of the treated
group:

$\hat{\mu}_{PSM}=\frac{1}{|\mathcal{I}_{\pi}|}\sum_{i\in\mathcal{I}_{\pi}}{y_{i}}$
(2)

Let $\hat{m}(x)$ be a Logistic Regression model fitted on the outcome of the
untreated group. Outcome score matching is done by matching units from the
treated group with units from the untreated group with nearest values of
$\hat{m}(x)$. Let $\mathcal{I}_{y}$ be the set of units in the untreated group
which are the k-nearest neighbors by their predicted outcome of some unit in
the treated group. The mean outcome is then estimating the mean the same way
as in Equation 2.

$\hat{\mu}_{OSM}=\frac{1}{|\mathcal{I}_{y}|}\sum_{i\in\mathcal{I}_{y}}{y_{i}}$
(3)

We used matching on a single continuous covariate i,e the outcome score since
according to [1] for the case when only a single continuous covariate is used
to match, the efficiency loss can be made arbitrarily close to zero by
allowing a sufficiently large number of matches.

We have also considered inverse propensity weighting estimator over the
population of nonrespondents [5]

$\hat{\mu}_{IPW-
NR}=\frac{\sum_{i}{w_{i}{\pi}^{-1}(1-\pi)y_{i}}}{\sum_{i}{w_{i}\pi^{-1}(1-\pi)}}$
(4)

A different group of estimators are regression estimators. These estimators
model the outcome conditioned on the covariates using regression and use this
model to estimate the missing outcomes. Let the outcome model denoted by
$\hat{m}(x)$ and $\mathcal{U}_{t}$ the set of units in the treated group then
the mean predicted outcome estimator is

$\hat{\mu}_{MPO}=\frac{1}{|\mathcal{U}_{t}|}\sum_{i\in\mathcal{U}_{t}}{\hat{m}(x_{i})}$
(5)

Doubly robust methods take into account both the model predictions for $y_{i}$
with inverse-probability weights. These methods are highly efficient when the
y-model is true yet remain asymptotically unbiased when the y-model is
misspecified [6]. We have considered such regression model for which the
Logistic Regression was trained with non-respondents weights
${\pi}^{-1}(1-\pi)$ and refer to to it as $\hat{m}_{w}(x)$ so that the
weighted mean predicted outcome estimator is :

$\hat{\mu}_{WMPO}=\frac{1}{|\mathcal{U}_{t}|}\sum_{i\in\mathcal{U}_{t}}{\hat{m}_{w}(x_{i})}$
(6)

###  3.2 confidence interval Estimation

In contrast to Machine Learning settings causal inference requires reporting
of confidence intervals on top of the prediction. This is highly important so
that the analyst can make a proper decision and know when prediction is of low
quality. We use bootstrapping to estimate the standard error of all our
estimators as recommended by many authors including [2]. In bootstrapping data
variability estimation the dataset is sampled many times with replacement to
produce a sample with the same size. For each sample the estimator is being
calculated. The standard deviation of the estimations is calculated and
reported as the estimation of the standard error of the estimator. From this
standard error (SE) the 95% confidence interval can be simply derived by

$CI95\%=[\hat{\mu}-1.96*SE,\hat{\mu}+1.96*SE]$ (7)

##  4 SIMULATION STUDIES

###  4.1 Treated Group Mean Estimation

A simulation study was done to study the performance of different estimators
for the potential mean outcome of the treated group. Design (1) is based on
similar setting described by [7] and design (2) is similar to one described in
[2]

Treatment Design 1: $W(X)=\mathbbm{1}{[a_{1}X_{1}+a_{2}X_{2}+a_{3}X_{3}>w]}$
Outcome Design 1 :
$Y(X)=\mathbbm{1}{[b_{1}X_{1}+b_{2}X_{2}+b_{3}X_{3}+\epsilon>y]}$ Outcome
Design 2 :
$Y(X)=\mathbbm{1}{[b_{1}X_{1}^{2}+b_{2}X_{2}+b_{3}X_{3}+\epsilon>y]}$
$X_{i}\sim\mathcal{N}(0,1),a_{i}\sim\mathcal{U}(0,1),b_{i}\sim\mathcal{U}(0,1),\epsilon\sim\mathcal{N}(0,\sigma)$
$w\sim\mathcal{U}(-1,1),y\sim\mathcal{U}(-1,1)$

We use $\sigma$ to control the outcome model noise and thus the error rate of
the estimation. For each design we sample 1000 data sets and estimate the mean
outcome using all the estimators described in section 3.1.

The results are summarized in Table 1, Table 2 and Table 3. All results in
this paper are presented in percentage format (actual value times 100) for
reader convenience. The method of Mean predicted Outcome gave best results in
our simulation studies in terms or Root mean square error and Mean Absolute
error. We’ve noticed that an increase in the number of samples from 10K to
100K resulted in an improvement for the misspecified model but in quantities
much lower than $\sqrt{n}$. The 10-NN outcome matching was the second runner
up with low errors suggesting the benefits of matching over more than one
unit. Methods not taking into account the outcome such as Inverse probability
weighting and propensity score matching gave very bad results which makes them
an unreasonable choice for the problem settings.

[6] found that the y-model and -model based estimators are sensitive to
misspecification models as not observing the direct confounders and rather
observing some non-linear function of them. We have observed the same results
for the methods tested. We’ve also noticed that introducing non-linearity by
itself is enough to degrade the mean population estimation even if the
covariates are observed.

Figure 1: Mean Estimation Vs True Mean for 6 different methods for Design 1 with 10K samples Figure 2: Mean Estimation Vs True Mean for 6 different methods for Design 2 with 10K samples Figure 3: Mean Estimation Vs True Mean for 6 different methods for Design 2 with 100K samples Table 1: Performance of the estimators with 10,000 samples over 1000 iterations for design 1 - correctly specified y-model and $\pi$-model  Method | BIAS | RMSE | MAE
---|---|---|---
Mean Predicted Outcome | 0.003 | 3.017 | 1.726
10-NN Outcome Matching | 0.077 | 4.09 | 2.059
IPW Mean Predicted Outcome | 0.111 | 3.879 | 2.268
1-NN Outcome Matching | 0.23 | 6.267 | 2.559
IPW | 0.353 | 13.512 | 10.714
1-NN Propensity Score Matching | 2.088 | 48.309 | 44.327
Table 2:  Performance of the estimators with 10,000 samples over 1000 iterations for design 2 - misspecified y-model and correctly specified $\pi$-model  Method | BIAS | RMSE | MAE
---|---|---|---
Mean Predicted Outcome | -0.06 | 5.473 | 3.626
10-NN Outcome Matching | -0.204 | 6.025 | 3.837
1-NN Outcome Matching | 0.06 | 7.128 | 4.21
IPW Mean Predicted Outcome | 0.163 | 6.567 | 4.321
IPW | 0.192 | 14.462 | 11.218
1-NN Propensity Score Matching | 0.646 | 48.969 | 45.304
Table 3:  Performance of the estimators with 100,000 samples over 1000 iterations for design 2 - misspecified y-model and correctly specified $\pi$-model  Method | BIAS | RMSE | MAE
---|---|---|---
Mean Predicted Outcome | 0.199 | 4.949 | 2.903
10-NN Outcome Matching | 0.359 | 5.192 | 2.95
1-NN Outcome Matching | 0.202 | 5.826 | 3.109
IPW Mean Predicted Outcome | 0.206 | 5.477 | 3.21
IPW | 0.242 | 13.699 | 10.612
1-NN Propensity Score Matching | -1.961 | 49.951 | 46.488

###  4.2 Confidence Interval

We’ve evaluated the confidence interval validity produced by bootstrapping. We
ran the simulation again with 100 samples of design 1. Each sample was of size
10K. Standard errors were estimated using 100 bootstrapped samples for each
iteration. We have calculated the distribution ratio between the error and the
estimated standard error and also the 95% coverage rate presented in Tables 4,
5. For the correctly specified model all confidence intervals were estimated
relatively good besides for the IPW confidence interval which underestimated
the error with only about 24% coverage. For the design 2 with misspecified
outcome model our estimation of the confidence interval degraded with about
75% coverage. The only good estimation of the coverage rate was produced by
the propensity score matching estimator.

Table 4:  Design 1 Coverage Rate of different estimation methods using bootstrapping  Method | Coverage Rate
---|---
Mean Predicted Outcome | 94
1-NN Outcome Score Matching | 96
10-NN Outcome Score Matching | 92
IPW Mean Predicted Outcome | 91
1-NN Propensity Score Matching | 91
IPW | 24
Table 5:  Design 2 Coverage Rate of different estimation methods using bootstrapping  Method | Coverage Rate
---|---
1-NN Propensity Score Matching | 94
1-NN Outcome Score Matching | 76
10-NN Outcome Score Matching | 76
IPW Mean Predicted Outcome | 75
Mean Predicted Outcome | 72
IPW | 23
Figure 4: Design 1 Histograms of ratio between estimation error and the
standard error Figure 5: Design 2 Histograms of ratio between estimation error
and the standard error

##  5 EMPIRICAL RESULTS

###  5.1 Treated Group Mean Estimation

To further investigate the estimator we’ve conducted an analysis over the
public dataset “Credit Card Fraud Detection” from the Kaggle website. The
dataset contains transactions made by credit cards in September 2013 by
European cardholders. This dataset presents transactions that occurred in two
days, with 492 frauds out of 284,807 transactions. The dataset is highly
unbalanced, the positive class (frauds) account for 0.172% of all
transactions. It contains only numeric input variables which are the result of
a PCA transformation. Due to confidentiality issues, the original features
were not provided. Features V1, V2, … V28 are the principal components
obtained with PCA, the only features which have not been transformed with PCA
are ’Time’ and ’Amount’. Feature ’Time’ contains the seconds elapsed between
each transaction and the first transaction in the dataset. The feature
’Amount’ is the transaction Amount, Feature ’Class’ is the response variable
and it takes value 1 in case of fraud and 0 otherwise.

We performed feature selection using the Anova test and selected the top 10
features. To simulate a policy and treated group we sample at each iteration 4
features and create a mock policy using a Gauss Naive Bayes model. The data is
split 50:50 to a train set and test set. The model is trained on the train
set. The treated group size is fixed to 100 units. The treated group is
assigned to the units with top model scores in the test set. Following that we
run all the estimators and evaluate the results which are summarized in Table
6.

For the empirical results 10-NN score matching outperformed all other methods
in terms of RMSE and MAE. IPW remains a bad choice for our domain but
surprisingly 1-NN propensity score matching performance increased compared to
our initial simulation studies.

The empirical studies along with the simulation studies give strong evidence
that these methods can be used in practice to get estimates of fraud
prevalence and catch fraud trends.

Table 6:  Performance of the estimators over 1000 iterations with Kaggle fraud credit card dataset and treated group size of 100  Method | BIAS | RMSE | MAE
---|---|---|---
10-NN Outcome Matching | -1.77 | 3.82 | 3.26
1-NN Outcome Matching | -1.90 | 4.77 | 4.01
Mean Predicted Outcome | -4.56 | 4.31 | 5.10
IPW Mean Predicted Outcome | -5.81 | 7.25 | 7.63
IPW | 7.37 | 6.97 | 8.83
1-NN Propensity Score Matching | 8.79 | 8.02 | 9.81
Figure 6: Mean Estimation Vs True Mean for 6 different methods for Real Fraud
Dataset from Kaggle website

###  5.2 Confidence Interval

We evaluate the performance of bootstrapping in estimating the confidence
interval of our estimators. We sample the Kaggle dataset 100 times for each
estimation method. Each sample contains 900 non-fraud units and 100 fraud
units. We train a Gaussian Naive Bayse model on 10% of the data with 4 random
features. We then predict the model score over the entire set. We add random
uniform noise $\mathcal{U}(0,1)$ to the Bayse Model score to get a higher
variance in the positive ratio within the treated group. We select a random
treated group size n with distribution $\mathcal{U}(50,80)$ and then assign
the top n scores to the treated group. The results are summarized in Table 7
and Figure 7.

Table 7:  Kaggle Fraud Dataset Coverage Rate of different estimation methods using bootstrapping  Method | Coverage Rate
---|---
IPW | 95
1-NN Propensity Score Matching | 98
1-NN Outcome Score Matching | 100
10-NN Outcome Score Matching | 100
IPW Mean Predicted Outcome | 100
Mean Predicted Outcome | 100
Figure 7: Fraud Dataset Histograms of ratio between estimation error and the
standard error

It seems that bootstrapping for the use case of Kaggle fraud dataset with
Naive Bayse policy overestimated the confidence interval and yielded
conservative estimation.

##  6 PRACTICAL CONSIDERATIONS

  1. 1.

Null Ratio - special consideration should be given to the null ratio both
feature level and item level. Since in production many times feature breakage
comes in bursts and items with high null rate $(>10\%)$ will get the wrong
match. We recommend using very reliable features and monitoring for breakage.

  2. 2.

Feature Selection - Features should be selected such that they explain both
the outcome and the treatment. This can be verified using common unitary
feature selection methods such as mutual information. Model should include as
many covariates as possible to ensure correctness of the outcome model.

##  7 DISCUSSION

In this work we’ve presented a method for estimating unobserved performance
metrics with focus on Financial Services. We’ve shown how matching on the
outcome score can predict metrics such as False Positive Rate with acceptable
errors and serve as an additional tool. Moreover it was found to outperform
other methods in a real fraud dataset.

We believe these methods, most of them known for many years in the statistics
community, can be extremely useful in many Integrity domains for which control
groups are not available and manual labeling effort reduction is desirable.

## Acknowledgement

The author would like to thank Shaked Bar for valuable comments.

## References

  * [1] Alberto Abadie and Guido W Imbens.  Large sample properties of matching estimators for average treatment effects.  econometrica, 74(1):235–267, 2006.
  * [2] Peter C Austin and Dylan S Small.  The use of bootstrapping when using propensity-score matching without replacement: a simulation study.  Statistics in medicine, 33(24):4306–4319, 2014.
  * [3] Jianqing Fan and Qiwei Yao.  Efficient estimation of conditional variance functions in stochastic regression.  Biometrika, 85(3):645–660, 1998.
  * [4] Michele Jonsson Funk, Daniel Westreich, Chris Wiesen, Til Stürmer, M Alan Brookhart, and Marie Davidian.  Doubly robust estimation of causal effects.  American journal of epidemiology, 173(7):761–767, 2011.
  * [5] Keisuke Hirano and Guido W Imbens.  Estimation of causal effects using propensity score weighting: An application to data on right heart catheterization.  Health Services and Outcomes research methodology, 2(3):259–278, 2001.
  * [6] Joseph DY Kang and Joseph L Schafer.  Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data.  Statistical science, 22(4):523–539, 2007.
  * [7] Ronnie Pingel.  Estimating the variance of a propensity score matching estimator for the average treatment effect.  Observational Studies, 4(1):71–96, 2018.
  * [8] James M Robins, Andrea Rotnitzky, and Lue Ping Zhao.  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American statistical Association, 89(427):846–866, 1994.
  * [9] Paul R Rosenbaum and Donald B Rubin.  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70(1):41–55, 1983.
  * [10] Elizabeth A Stuart.  Matching methods for causal inference: A review and a look forward.  Statistical science: a review journal of the Institute of Mathematical Statistics, 25(1):1, 2010.
  * [11] Machine Learning Group ULB.  Credit card fraud detection [dataset].  https://www.kaggle.com/mlg-ulb/creditcardfraud, 2018.  Data was retrieved from Kaggle.
  * [12] Stefan Wager and Susan Athey.  Estimation and inference of heterogeneous treatment effects using random forests.  Journal of the American Statistical Association, 113(523):1228–1242, 2018.
  * [13] Yuxiang Xie, Meng Xu, Evan Chow, and Xiaolin Shi.  How to measure your app: A couple of pitfalls and remedies in measuring app performance in online controlled experiments.  In Proceedings of the 14th ACM International Conference on Web Search and Data Mining, pages 949–957, 2021.