# A learning-based multiscale model for reactive flow in porous media

Mina Karimi and Kaushik Bhattacharya  California Institute of Technology,
Pasadena CA 91125

###### Abstract

We study solute-laden flow through permeable geological formations with a
focus on advection-dominated transport and volume reactions. As the fluid
flows through the permeable medium, it reacts with the medium, thereby
changing the morphology and properties of the medium; this in turn, affects
the flow conditions and chemistry. These phenomena occur at various lengths
and time scales, and makes the problem extremely complex. Multiscale modeling
addresses this complexity by dividing the problem into those at individual
scales, and systematically passing information from one scale to another.
However, accurate implementation of these multiscale methods are still
prohibitively expensive. We present a methodology to overcome this challenge
that is computationally efficient and quantitatively accurate. We introduce a
surrogate for the solution operator of the lower scale problem in the form of
a recurrent neural operator, train it using one-time off-line data generated
by repeated solutions of the lower scale problem, and then use this surrogate
in application-scale calculations. The result is the accuracy of concurrent
multiscale methods, at a cost comparable to those of classical models. We
study various examples, and show the efficacy of this method in understanding
the evolution of the morphology, properties and flow conditions over time in
geological formations.

##  1 Introduction

The transport of water and solutes through permeable geological formations
couple various phenomena [19, 5]. As it flows through the permeable medium,
solute-laden water reacts with the medium changing the morphology and
properties of the medium; this in turn affects the flow and chemistry.
Permeable medium supports a menagerie of microbial life, and these are
affected by and affect the flow and the chemistry. Chemical interactions,
biological activity and flow may cause mechanical failure which in turn
changes morphology thereby affecting the flow. Further, these varied phenomena
occur, and manifest themselves, at various length and time scales (see for
example [26] and the citations there): at the molecular scale where chemistry
happens, at the pore scale where morphology, biology, mechanical properties
and flow interact, the core scale where averaged features emerge, and finally
at the geological scale relevant to aquifer and reservoir hydrology.
Furthermore, these geological formations are heterogenous at various scales
ranging from the individual grains, to pores, to formations, fissures, caverns
and cracks to geological strata and faults. Finally, they are highly
anisotropic at various scales ranging from the clay particles to faults to
geological formations.

This enormous complexity makes modeling extremely difficult. There are well-
developed models for individual processes at particular scales, and in some
cases conceptual frameworks to link various scales. Still, there is a critical
gap in accurate methods of passing information from one scale to another,
especially when it concerns multiple phenomena and history (or time)-dependent
phenomena. In this paper, we address this coupling between the pore, core and
geological formation scales with a focus on advection-dominated transport and
volume reactions.

We begin with a brief review of the existing approaches, and then describe our
contributions.

Geological scale models. Continuum models for reactive flow at the geological
scale focus on predicting the overall flow and transport. These models usually
employ empirical relations to relate the evolving transport parameters to
changes in porosity, reactivity, and specific area [5]. To do so, many studies
estimate the change in porosity [3], diffusivity [33], permeability [3, 43,
29], reaction rate, and specific area [35] due to chemical reactions. However,
these approaches face some limitations. These models are derived assuming
simple geometries (e.g. spherical grains) and therefore can not provide
detailed insight into the pore structure and specific area, and their
evolution. Further, they have a limited range of accuracy in the presence of
precipitation and dissolution. In particular, they often fail when porosity
approaches critical values or percolation.

Core scale models. Core-scale simulations focus on predicting the evolution of
the pore structure and average fluid transport properties. Using conventional
finite element methods for these simulations is computationally expensive due
to the need to update the geometry [5]. Several alternative numerical
strategies have been developed to overcome this limitation, such as pixel-
based approaches [12, 13], level set methods [18, 39], smoothed particle
hydrodynamics [37], and adaptive discontinuous Galerkin finite element methods
[36].

Multiscale models. These models employ upscaling approaches to bridge the gap
between pore and geological scales. Volume averaging methods [6] estimate the
effective transport properties and derive macroscopic equations. Initially
explored for porous media dispersion by [8, 31], these methods were later
extended to spatially periodic porous media for reactive fluid transport by
Paine [30]. However, such an approach necessitates incorporating ad hoc
assumptions for closure relations.

Early studies by Auriault and Adler [4] and Rubinstein and Mauri [32]
investigate multiscale expansions to upscale Taylor dispersion [38] in
periodic porous media, primarily focusing on diffusion-dominated, non-reactive
transport. Later, this methodology was extended to reactive flow with surface
reactions by Mauri [25]. However, this approach failed to derive effective
properties for advection-dominated flows. Other researchers [24, 2, 10]
introduced a two-scale expansion approach with drift to address advection-
dominated flows. They employed a moving coordinate system with an effective
drift velocity to investigate advection’s impact on effective equations. This
method was further applied by Allaire et al., [1] to study reactive transport
with surface reactions while considering a linear adsorption/desorption rate
at the solid interface.

Machine learning. In recent years, machine learning (ML) techniques have been
applied to expedite reactive transport simulations and estimate effective
porous medium properties. Liu et al [22] have used an ML approach to predict
the effective reaction rate using features of pore structures. Artificial
neural networks have been used to accelerate the geological scale reactive
flow simulations [9, 34], and deep learning methods have been used to estimate
the permeability considering non-reactive transport in the porous medium [42,
40]. These consider a snap-shot in time and do not address evolution over long
periods. Wang and Battiato [41] have developed a deep-learning multiscale
model to predict the clogging caused by solute precipitation in a microcrack
network. Lu et al. [23] have used a neural network model to predict the
evolution of the uranium distribution coefficient in the subsurface due to
thermal and chemical processes.

In the context of mechanical properties, recurrent neural networks (RNNs)
address the history-dependent behaviors [27], Long Short-Term Memory (LSTM)
effectively remembers information over long time intervals [11], gated
recurrent unit (GRU), a simplified variant of LSTM captures similar temporal
relationships in data [28, 44]. However, LSTM-based approaches require
millions of variables to be trained from data. Recently, a recurrent neural
operator (RNO) has been introduced inspired by internal variable theory,
offering an efficient approach for multiscale modeling [21, 20].

##### Our Contributions.

In this work, we develop a methodology for investigating flow through
underground geological formations over long periods of time. We specifically
focus on an advection dominated transport in a porous medium and reactions at
the fluid-solid interface characterized by nonlinear reaction kinetics. We
first use a two-scale expansion method with effective drift velocity [1] to
obtain the governing equations are the core and geological formulations. This
is described in Section 2.

We then introduce a recurrent neural operator (RNO), and use it to learn the
solution operator of the core scale problem. To elaborate, we use data
generated by repeated solutions of the core scale problem to train the RNO to
learn the map between geological scale variables over time (e.g., velocity and
solute concentration histories) and the effective transport properties (e.g.,
permeability, diffusivity, advection velocity, porosity, and specific area).
This is described in Section 3.

Finally, we use the trained RNO surrogate model as a surrogate in geological
scale computations to investigate the long time evolution of these formations
in the presence of volume reactions induced by flow and transport. We
demonstrate the accuracy of the approach, as well as its ability to reveal
non-trivial interactions between the core and geological scale in Section 4.

We conclude in Section 5 with a discussion of implications, promises and open
issues.

##  2 Two scale model

###  2.1 Governing equation and non-dimensionalization

We consider a porous geological formation $\tilde{\Omega}$ composed of a
porous region $\tilde{\Omega}_{p}$ and solid region $\tilde{\Omega}_{s}$ with
$\tilde{\Omega}_{p}\cap\tilde{\Omega}_{s}=\phi,\
\tilde{\Omega}^{c}=\tilde{\Omega}_{p}^{c}\cup\tilde{\Omega}_{s}^{c}$ where the
superscript ‘c’ denotes the closure. We have an incompressible fluid flow in
the pore governed by the steady Stokes equation

$\displaystyle\begin{cases}-\tilde{\nabla}\tilde{p}+\tilde{\nu}\Delta\tilde{{\mathbf{v}}}+\tilde{{\mathbf{f}}}={\mathbf{0}}\quad&\text{in}\
\tilde{\Omega}_{p},\\\ \tilde{\nabla}\cdot\tilde{{\mathbf{v}}}=0&\text{in}\
\tilde{\Omega}_{p},\\\ \tilde{{\mathbf{v}}}={\mathbf{0}}&\text{in}\
\tilde{\Omega}_{s}\end{cases}$ (1)

where $\tilde{{\mathbf{v}}}$ is the particle velocity, $\tilde{p}$ is
pressure, $\tilde{\nu}$ is the viscosity and $\tilde{{\mathbf{f}}}$ is the
body force. We assume that the fluid carries a solute that is transported in
the fluid through a combination of diffusion and advection and reacts with the
surface of the solids, resulting either in deposition or dissolution

$\displaystyle\begin{cases}\tilde{c}_{\tilde{t}}+\tilde{{\mathbf{v}}}\cdot\tilde{\nabla}\tilde{c}-D\tilde{\Delta}\tilde{c}=0\quad&\text{in}\
\tilde{\Omega}_{p},\\\
-D\tilde{\nabla}\tilde{c}\cdot{\mathbf{n}}=-(\tilde{c}-\tilde{m})\tilde{v}_{n}=\tilde{q}_{c}(\tilde{c})&\text{on}\
\partial\tilde{\Omega}_{p}\end{cases}$ (2)

where $\tilde{c}$ is the concentration of the solute, $D$ is the diffusion
constant, ${\mathbf{n}}$ the unit normal vector to the surface
$\partial\tilde{\Omega}_{p}$ oriented outward with respect to
$\tilde{\Omega}_{p}$, $\tilde{m}$ is the concentration of the solute in the
solids, $\tilde{v}_{n}$ is the normal speed of the solid/fluid interface due
to dissolution or deposition and $\tilde{q}_{c}(c)$ is the reaction rate that
depends on the solute concentration. Above, the second term of the bulk
equation describes advection, while the third describes diffusion. On the
interface, the first equation describes the mass balance between the flux of
the solute from the fluid to the interface and the growth of the interface,
while the second relates the growth of the interface to the interfacial
reaction rate. These equations have to be supplemented with appropriate
boundary conditions.

Note that we have assumed that the fluid flow is steady while the solute
transport is time-dependent. We assume that the reaction rate at the interface
and, consequently, the rate of reconstruction of the porous medium is slow
compared to the time-scale of the fluid flow. So, the fluid flow reaches a
steady state at each time as the medium reconstructs.

Now, the characteristic length of the geological formation $L$ is very large
compared to the characteristic length of the core or representative volume
$\ell$, $L>>\ell$. This makes the system (1, 2) difficult to solve: we have to
resolve the flow and transport with a resolution small compared to $\ell$, but
on a domain of size $L$. Therefore, we resort to a two-scale asymptotic
expansion under the assumption that the ratio of length-scales
$\epsilon=\ell/L$ is small, $\epsilon<<1$.

In order to do so, we change to non-dimensional units by scaling length with
the characteristic length $L$ of the geological formation, the velocity with
the characteristic velocity $V$ and the pressure with characteristic pressure
$\Pi$. It follows that time is rescaled by the characteristic time $T=L/V$. We
expect slow flows with small $V$ over long distances $L$, which means that the
characteristic time $T$ is large and consistent with steady state. Now, recall
that the characteristic length of the pores is small, and therefore, in order
to have non-trivial flow, we need the viscosity to be extremely small.
Therefore, we assume that the characteristic viscosity is $\Pi
T/\epsilon^{2}$. The non-dimensional flow equations are given by

$\displaystyle\begin{cases}-{\nabla}{p^{\epsilon}}+\epsilon^{2}\nu\Delta{{\mathbf{v}}^{\epsilon}}+{{\mathbf{f}}}={\mathbf{0}}\quad&\text{in}\
{\Omega}^{\epsilon}_{p},\\\
{\nabla}\cdot{{\mathbf{v}}^{\epsilon}}=0&\text{in}\
{\Omega}^{\epsilon}_{p},\\\ {{\mathbf{v}}^{\epsilon}}={\mathbf{0}}&\text{in}\
{\Omega}^{\epsilon}_{s}.\end{cases}$ (3)

Above, we have introduced $\epsilon$ as a superscript in the non-dimensional
variables to signify that the porosity and therefore variations in these
quantities are at a scale $\epsilon$.

We now turn to solute transport. We non-dimensionalize the concentration with
a characteristic concentration $C$ and introduce two non-dimensional numbers,
the Péclet number Pe and the Damköhler number Da:

$\text{Pe}^{\epsilon}=\frac{LV}{D},\quad\text{Da}^{\epsilon}=\frac{LK}{D}$ (4)

where $K$ is a characteristic reaction rate. The non-dimensional equations of
solute transport are

$\begin{cases}c_{t}^{\epsilon}+\text{Pe}^{\epsilon}{\mathbf{v}}\cdot\nabla
c^{\epsilon}-\Delta c^{\epsilon}=0\quad&\text{in}\ \Omega^{\epsilon}_{p}\\\
-\nabla
c^{\epsilon}\cdot{\mathbf{n}}=-v_{n}^{\epsilon}(c^{\epsilon}-m)=\text{Da}^{\epsilon}q_{c}^{\epsilon}&\text{on}\
\partial\Omega^{\epsilon}_{p}\end{cases}$ (5)

Now, we are interested in situations where we have significant advection at
the pore scale, and very slow reactions at the interface. So, we assume that

$\text{Pe}^{\epsilon}=\frac{\widehat{\text{Pe}}}{\epsilon},\quad\text{Da}^{\epsilon}=\epsilon\widehat{\text{Da}},\quad
v_{n}^{\epsilon}=\epsilon\hat{v}_{n}$ (6)

where $\widehat{\text{Pe}},\widehat{\text{Da}},\hat{v}_{n}$ are all $O(1)$.

###  2.2 Two-scale model

We assume that the porous media is almost periodic, i.e., it is periodic on
the scale $\epsilon$ but can change over long distances compared to
$\epsilon$. To be precise, we assume
$\Omega^{\epsilon}_{p}(x)=\Omega_{p}(x,x/\epsilon)$ where $\Omega_{p}$ is
periodic with period 1 in the second variable; so porous medium is periodic on
$\epsilon Y$ where $Y$ unit cube or unit cell in the vicinity of the point
${\mathbf{x}}$. Further, $Y=Y_{p}({\mathbf{x}})\cup Y_{s}({\mathbf{x}})$ where
$Y_{p}({\mathbf{x}})$ is the pore in the unit cell in the vicinity of the
point ${\bf x}$ in the geological formation. We show that under this
assumption, we can approximate the solution of the system (3,5) by solving a
geological scale problem where the constitutive behavior is determined by
solving a core scale problem. We first describe the two problems, and then
justify this derivation in the following subsections.

####  2.2.1 Geological scale model

We can find the overall flow and solute transport at the geological scale by
solving the following system on $\Omega$:

$\displaystyle\begin{cases}\nabla\cdot{\mathbf{v}}_{0}=0,\quad{\mathbf{v}}_{0}=\frac{1}{\nu}{\mathbf{K}}^{*}({\mathbf{f}}-\nabla
p_{0}),\\\ \lambda{c}_{0t}+\lambda\text{Pe}\overline{\mathbf{v}}\cdot\nabla
c_{0}-\nabla\cdot({\mathbf{D}}^{*}\nabla
c_{0})=-\gamma\widehat{\text{Da}}q_{c}\end{cases}$ (7)

for the overall velocity ${\mathbf{v}}_{0}$, pressure $p_{0}$, concentration
$c_{0}$ at the geological scale subject to boundary conditions. Above, the
parameters ${\mathbf{K}}^{*}={\mathbf{K}}^{*}(t,{\mathbf{x}})$ is the
permeability tensor, ${\mathbf{D}}^{*}={\mathbf{D}}^{*}(t,{\mathbf{x}})$ is
the effective diffusivity tensor, $\lambda=\lambda(t,{\mathbf{x}})$ is the
pore volume fraction,
$\overline{\mathbf{v}}=\overline{\mathbf{v}}(t,{\mathbf{x}})$ is an effective
advection velocity, and $\gamma=\gamma(t,{\mathbf{x}})$ is a local surface
area per unit volume. Note that these parameters are all functions of time,
and they are specified through the core scale problem below.

Note that this geological scale problem is solved at the entire geological
$\Omega$ and the coefficients vary only on the scale of the geological
formation. All the information about the pores has been subsumed into
parameters that only vary at the geological scale.

####  2.2.2 Core scale or unit cell model

Given a porous unit cell at the macroscopic point ${\mathbf{x}}$ at time $t$,
i.e., given $Y_{p}$, the unit cell problem is to solve

$\displaystyle\begin{cases}-\nabla_{y}q^{i}+\nu\Delta_{y}{\mathbf{u}}^{i}+{\mathbf{e}}^{i}={\mathbf{0}},\
\ \nabla_{y}\cdot{\mathbf{u}}^{i}=0\quad&\text{in}\ Y_{p}\\\
{\mathbf{u}}^{i}={\mathbf{0}}&\text{in}\ Y_{s}\\\ \widehat{\text{Pe}}\
{\mathbf{v}}\cdot\left(\nabla_{y}\chi^{j}\right)-\Delta_{y}\chi^{j}=\
({\mathbf{v}}^{*}-\widehat{\text{Pe}}\
{\mathbf{v}})\cdot{\mathbf{e}}^{j}\quad&\text{in}\ Y_{p}\\\
-\left(\nabla_{y}\chi^{j}\right)\cdot{{\mathbf{n}}}={\mathbf{e}}^{j}\cdot{{\mathbf{n}}}&\text{on}\
\partial Y_{p}\end{cases}$ (8)

periodic velocity fluctuation ${\mathbf{u}}^{i}$, pressure fluctuation
$q^{i}$, chemical fluctuation $\chi^{j}$ when the overall flow is in the
direction ${\mathbf{e}}^{i}$ and overall solute transport in the direction
${\mathbf{e}}^{j}$ for $i,j=1,\dots,d$ (dimension $d$). Above, we use
$\nabla_{y},\Delta_{y}$ to signify that these are derivatives with respect to
the spatial variable in the unit cell.

We can then use it to find the parameters

$\displaystyle{\mathbf{K}}^{*}_{ij}$
$\displaystyle=\int_{Y_{p}}\nabla_{y}{\mathbf{u}}^{i}\cdot\nabla_{y}{\mathbf{u}}^{j}\
dy,$ (9) $\displaystyle{\mathbf{D}}^{*}_{ij}$
$\displaystyle=\int_{Y_{p}}{\mathbf{e}}^{i}\cdot{\mathbf{e}}^{j}\
dy+\widehat{\text{Pe}}\int_{Y_{p}}(\overline{\mathbf{v}}_{i}-{\mathbf{v}}_{i})\chi^{j}\
dy+\int_{Y_{p}}\nabla_{y}\chi^{j}\cdot{\mathbf{e}}^{i}\ dy,$ (10)
$\displaystyle\overline{\mathbf{v}}$
$\displaystyle=\frac{1}{|Y_{p}|}\int_{Y_{p}}{\mathbf{v}}\
dy,\quad\lambda=|Y_{p}|,\quad\gamma=|\partial Y_{p}|.$ (11)

Finally, the microstructure evolves according to the equation

$\hat{v}_{n}=\frac{\widehat{\text{Da}}}{m-\tilde{c}_{0}}q_{c}(\tilde{c}_{0})\approx\frac{\widehat{\text{Da}}}{m}q_{c}(\tilde{c}_{0}).$
(12)

The last approximation uses the fact that $\tilde{c}_{0}<<m$. This evolution
happens on the geological time scale, and therefore the parameters above
change on the geological time scale.

####  2.2.3 Summary

We summarize the resulting multiscale formulation. We solve (7) on the
geological scale. This requires us to obtain the parameters
${\mathbf{K}}^{*}(t,{\mathbf{x}}),{\mathbf{D}}^{*}(t,{\mathbf{x}}),\bar{\mathbf{v}}(t,{\mathbf{x}}),\gamma(t,{\mathbf{x}}),\lambda(t,{\mathbf{x}})$.
To do so, at each point ${\mathbf{x}}$ at the geological scale, we provide a
history of the flow ${\mathbf{v}}_{0}$ and solute concentration $c_{0}$ to a
unit cell, solve (8) and obtain the parameters from (9), (10) and (11) as a
function of time while evolving the microstructure according to (12). This is
illustrated in Figure 1.

While this formulation separates the original problem into two separate
problems, it is still computationally demanding: we have to solve a core scale
problem at at every quadrature point of the geological formation and at every
instant of time. The direct implementation of this framework is often referred
to as concurrent multiscale approach. We propose an alternate approach in
Sections 3 and 4 below.

Figure 1: Schematic figure of the two-scale framework.

###  2.3 Details of the asymptotic expansion

This sub-section outlines the derivation of the two-scale formulation above.
It closely follows Allaire et al. [1].

####  2.3.1 Porous media flow

We look for a solution to the system (3) with the ansatz

$\displaystyle{\mathbf{v}}^{\epsilon}({\mathbf{x}})=\sum_{i=0}^{+\infty}\epsilon^{i}{\mathbf{v}}_{i}\left({\mathbf{x}},\frac{{\mathbf{x}}}{\epsilon}\right)\,\quad
p^{\epsilon}({\mathbf{x}})=\sum_{i=0}^{+\infty}\epsilon^{i}p_{i}\left({\mathbf{x}},\frac{{\mathbf{x}}}{\epsilon}\right).$
(13)

By collecting terms and then using the Fredholm alternative, Lévy [17] showed
that ${\mathbf{v}}_{0}={\mathbf{v}}_{0}({\mathbf{x}}),\ p_{0}=p({\mathbf{x}})$
(i.e, they are independent of the fast variable), and that these satisfy (7)1
with ${\mathbf{K}}^{*}$ given by (9) where ${\mathbf{u}}^{i}$ satisfies
(8)1,2. Further, up to leading order the velocity and the pressure are

$\displaystyle{\mathbf{v}}^{\epsilon}({\mathbf{x}})={\mathbf{v}}_{0}({\mathbf{x}})+\epsilon\sum_{i=1}^{d}({\mathbf{v}}_{0})_{i}({\mathbf{x}}){\mathbf{u}}^{i}\left(\frac{{\mathbf{x}}}{\epsilon}\right),\quad
p^{\epsilon}({\mathbf{x}})=p_{0}({\mathbf{x}})+\epsilon\sum_{i=1}^{d}q^{i}\left(\frac{{\mathbf{x}}}{\epsilon}\right)\left({\mathbf{f}}_{i}({\mathbf{x}})-\frac{\partial
p_{0}}{\partial x_{i}}({\mathbf{x}})\right).$ (14)

####  2.3.2 Solute transport

We look for a solution to the system (5) with the ansatz

$\displaystyle
c^{\epsilon}(t,{\mathbf{x}})=\sum^{+\infty}_{i=0}\epsilon^{i}\tilde{c}_{i}\left(t,{\mathbf{x}}-\frac{{\mathbf{v}}^{*}}{\epsilon}t,\frac{{\mathbf{x}}}{\epsilon}\right)$
(15)

where ${\mathbf{v}}^{*}$ is a drift velocity to be determined. We have chosen
a macroscopic coordinate
${\mathbf{x}}^{\prime}={\mathbf{x}}-\frac{{\mathbf{v}}^{*}}{\epsilon}t$ that
is not stationary in the geological formation in order to account for the
advection. We substitute the ansatz into (5), apply the chain rule,

$\displaystyle c_{t}=\frac{\partial{\tilde{c}}}{\partial
t}-\frac{{\mathbf{v}}^{*}}{\epsilon}\cdot\nabla^{\prime}\tilde{c},\quad\nabla
c=\nabla^{\prime}\tilde{c}+\frac{1}{\epsilon}\nabla_{y}\tilde{c}$ (16)

and collect terms with different powers of $\epsilon$.

The leading order (smallest power) is $\epsilon^{-2}$, and we have the
following equations at that order:

$\begin{cases}\hat{\text{Pe}}{\mathbf{v}}\cdot\nabla_{y}\tilde{c}_{0}-\Delta_{y}\tilde{c}_{0}=0\quad&\text{in}\
Y_{p},\\\ -\nabla_{y}\tilde{c}_{0}\cdot{{\mathbf{n}}}=0&\text{on}\ \partial
Y_{p}\end{cases}$ (17)

which yields to $\tilde{c}_{0}=\tilde{c}_{0}({\mathbf{x}}^{\prime})$.

At the next order, $\epsilon^{-1}$, we have

$\begin{cases}\widehat{\text{Pe}}\
{\mathbf{v}}\cdot\nabla_{y}\tilde{c}_{1}-\Delta_{y}\tilde{c}_{1}=({\mathbf{v}}^{*}-\widehat{\text{Pe}}\
{\mathbf{v}})\cdot\nabla^{\prime}\tilde{c}_{0}\quad&\text{in}\ Y_{p},\\\
-\left(\nabla_{y}\tilde{c}_{1}\right)\cdot{{\mathbf{n}}}=\nabla^{\prime}\tilde{c}_{0}\cdot{{\mathbf{n}}}&\text{on}\
\partial Y_{p}.\end{cases}$ (18)

We seek to solve this equation for $\tilde{c}_{1}$. This is possible if and
only if the following compatibility condition (Fredholm alternative) is
satisfied.

$\int_{Y_{p}}({\mathbf{v}}^{*}-\widehat{\text{Pe}})\cdot\nabla^{\prime}\tilde{c}_{0}\
dy-\int_{\partial
Y_{p}}\left(\nabla^{\prime}\tilde{c}_{0}\cdot{{\mathbf{n}}}\right)dS=0\quad\Rightarrow\quad{\mathbf{v}}^{*}=\frac{\widehat{\text{Pe}}}{|Y_{p}|}\int_{Y_{p}}{\mathbf{v}}\
dy.$ (19)

Recalling the expansion (14)1 for ${\mathbf{v}}$, we can re-write the drift
velocity as follow

${\mathbf{v}}_{i}^{*}=\frac{\widehat{\text{Pe}}}{|Y_{p}|}\sum_{i=1}^{d}({\mathbf{v}}_{0})_{i}\int_{Y_{p}}{\mathbf{u}}^{i}\
dy$ (20)

Returning to (18), we notice that the solution $\tilde{c}_{1}$ will depend
linearly on $\nabla^{\prime}c_{0}$. So, set

$\tilde{c}_{1}(t,{\mathbf{x}}^{\prime},{\mathbf{y}})=\sum^{d}_{i=1}\frac{\partial\tilde{c}_{0}}{\partial
x^{\prime}_{i}}(t,{\mathbf{x}}^{\prime})\chi^{i}({\mathbf{y}}).$ (21)

It follows that $\chi^{i}$ satisfy (8)3,4.

Now, turning to order $\epsilon$, we have

$\begin{cases}\widehat{\text{Pe}}\
{\mathbf{v}}\cdot\nabla_{y}\tilde{c}_{2}-\Delta_{y}c_{2}\\\
\quad\quad=-\frac{\partial\tilde{c}_{0}}{\partial
t}+({\mathbf{v}}^{*}-\widehat{\text{Pe}}\
{\mathbf{v}})\cdot\nabla^{\prime}\tilde{c}_{1}+\nabla^{\prime}\cdot\left(\nabla^{\prime}\tilde{c}_{0}+\nabla_{y}\tilde{c}_{1}\right)+\nabla_{y}\cdot(\nabla^{\prime}\tilde{c}_{1})\quad&\text{in}\
Y_{p},\\\
-\nabla_{y}\tilde{c}_{2}\cdot{{\mathbf{n}}}-\nabla^{\prime}\tilde{c}_{1}\cdot{\mathbf{n}}=\hat{v}_{n}(\tilde{c}_{0}-m)=\widehat{\text{Da}}\
q_{c}(\tilde{c}_{0})&\text{on}\ \partial Y_{p}.\end{cases}$ (22)

This equation has a solution for $\tilde{c}_{2}$ if and only if the following
compatibility condition (Fredholm alternative) is satisfied

$\int_{Y_{p}}\left(-\frac{\partial\tilde{c}_{0}}{\partial
t}+({\mathbf{v}}^{*}-\widehat{\text{Pe}}\
{\mathbf{v}})\cdot\nabla^{\prime}\tilde{c}_{1}+\nabla^{\prime}\cdot\left(\nabla^{\prime}\tilde{c}_{0}+\nabla_{y}\tilde{c}_{1}\right)\right)\
dy=\int_{\partial{Y_{p}}}\widehat{\text{Da}}\ q_{c}(\tilde{c}_{0})\ dS.$ (23)

Substituting for $\tilde{c}_{1}$, we obtain the homogenized equation

$-\lambda\frac{\partial\tilde{c}_{0}}{\partial
t}=-\nabla^{\prime}\cdot({\mathbf{D}}^{*}\nabla^{\prime}\tilde{c}_{0})+\gamma\widehat{\text{Da}}\
q_{c}(\tilde{c}_{0})$ (24)

where $\lambda=|Y_{p}|$, and $\gamma=|\partial Y_{p}|$, and the effective
diffusion tensor ${\mathbf{D}}^{*}$ is given by (10). We then transform back
to stationary coordinates by setting.

$c_{0}(t,{\mathbf{x}})=\tilde{c}_{0}\left(t,{\mathbf{x}}-\frac{{\mathbf{v}}^{*}}{\epsilon}t\right)$
(25)

to obtain (7)2.

Finally, the boundary condition at $O(\varepsilon)$ gives (12).

##  3 Learning the core scale behavior

The two-scale formulation above requires us to solve the core scale problem at
each time step and each quadrature point in the geological formation. This is
prohibitively expensive. So, we seek to “learn” the solution operator of the
core scale problem. Specifically, for a given initial microstructure $T^{0}$,
we view the core scale problem as a map from the velocity and concentration
history to the current permeability, diffusivity, drift velocity, specific
area, and pore volume fraction.

$\Phi:{\mathcal{I}}[0,t]\to{\mathcal{O}}(t),\quad{\mathcal{I}}(\tau)=\\{{\mathbf{v}}_{0}({\tau}),c_{0}(\tau)\\},\
{\mathcal{O}}(\tau)=\\{{\mathbf{K}}^{*}(\tau),{\mathbf{D}}^{*}(\tau),\bar{\mathbf{v}}(\tau),\gamma(\tau),\lambda(\tau)\\}.$
(26)

where the input
${\mathcal{I}}[0,t]=\\{{\mathbf{v}}_{0}(\tau),c_{0}(\tau):\tau\in[0,t]\\}$ is
specified over the time interval $[0,t]$ and the output
${\mathcal{O}}(t)=\\{{\mathbf{K}}^{*}(t),{\mathbf{D}}^{*}(t),\bar{\mathbf{v}}(t),\gamma(t),\lambda(t)\\}$
at time $t$. We seek an approximation in the form of a parametrized map

$\Psi:{\mathcal{I}}[0,t]\times{\mathbb{R}}^{p}\to{\mathcal{O}}$ (27)

and train it using data $\\{{\mathcal{I}}^{n},{\mathcal{O}}^{n}\\}_{n=1}^{N}$
that is generated using numerical simulation of $\Phi$. In other words, we
shall postulate a form for $\Psi$ and find the parameters $\Theta^{*}$ that
minimizes a loss function for data generated by repeated solutions of the
core-scale problem.

###  3.1 Recurrent neural operator

There are two issues we have to address in postulating an approximation
$\Psi$. First, the map $\Phi$ (and hence $\Psi$) has as its input a function
defined on an interval of time. Thus our approximation has to be an operator.
One idea is to discretize the functions in time and then find a neural network
approximation with the discretized function as input. Unfortunately, this
approximation will depend on the discretization (time step), and hence can
only be used at that discretization. However, it is natural in a multi-scale
setting to use different discretization for the core scale problem (generating
data) and the geological scale problem. Further, one may use an adaptive
discretization in the geological scale problem. For these and other reasons,
we want approximation to be independent of discretization. Second, the output
at time $t$ depends on the history of the input. So, we want our map to be
history dependent.

Following experience and practice in continuum physics, we postulate that the
history can be encapsulated in $k$ state or internal variables
$\\{\xi_{\alpha}\\}_{\alpha=1}^{k}$ that evolve in time. Then, following
recent work in the multi-scale modeling of mechanical properties of materials
[7, 21], we postulate $\Psi$ to be a recurrent neural operator :

$\Psi:\begin{cases}{\mathcal{O}}(t)=f\left({\mathcal{I}}(t),\\{\xi_{\alpha}(t)\\}_{\alpha=1}^{k};\Theta\right),\\\
\dot{\xi}_{i}(t)=g_{i}\left({\mathcal{I}}(t),\\{\xi_{\alpha}(t)\\}_{\alpha=1}^{k};\Theta\right),\quad
i=1,\dots,k\end{cases}$ (28)

where $f,g_{i}$ are feed-forward deep neural networks parametrized by $\Theta$
(weights and biases). The architecture (28) is formulated to be in continuous
time. To implement it with time discretization, we use a backward Euler
discretization:

$\begin{cases}{\mathcal{O}}^{n}=f\left({\mathcal{I}}^{n},\\{\xi^{n}_{\alpha}\\}_{\alpha=1}^{k}\right)\\\
\xi^{n}_{i}=\xi^{n-1}_{i}+(\Delta
t_{n})g_{i}\left({\mathcal{I}}^{n},\\{\xi^{n-1}_{\alpha}\\}_{\alpha=1}^{k}\right),\quad
i=1,\dots,k.\end{cases}$ (29)

Note that $f,g_{i}$ and the internal variables, and therefore the
approximation, is independent of the discretization.

The number $k$ of internal variables has to be chosen a priori, but the actual
internal variables are identified from the data as a part of the learning. As
noted, they encapsulate the history dependence. They do not necessarily have
any intrinsic physical meaning. Indeed, note that the form of the architecture
(28) is invariant under the reparametrization $\xi^{\prime}=\Xi(\xi)$ for any
diffeomorphism $\Xi$. In some special examples, it is possible to choose a
parametrization so that the internal variables are interpretable [21];
however, this is not always the case. We refer the reader to [21] for a
discussion of these and other aspects of this architecture.

###  3.2 Data and training

We generate the data by solving the core scale problem (8) over some interval
$[0,T]$ to yield our data in the form
$\\{{\mathcal{I}}^{n}[0,t],{\mathcal{O}}^{n}[0,T]\\}_{n=1}^{N}$. To do so, we
have to sample the inputs $\\{{\mathcal{I}}^{n}[0,t]\\}_{n=1}^{N}$ in such a
manner that is rich enough to represent the actual trajectories encountered in
the geological scale model. Broadly, we anticipate trajectories of velocity
and concentration that can vary over time, and also change slope as some
region gets clogged or fully dissolved. So we use the following strategy. We
take our interval $[0,T]$ and divide them into $M$ sub-intervals
$\\{[t^{m-1},t^{m}]\\}_{m=1}^{M}$ with $t^{m}\leq t^{m+1},\ t^{0}=0,\ t^{M}=T$
where $\\{t^{m}\\}_{m=1}^{M-1}$ are chosen from a uniform distribution (and
relabelled to be increasing). We then set each component ${\mathcal{I}}_{i}$
of the input at times $\\{t^{m}\\}_{m=1}^{M-1}$

$({\mathcal{I}}_{i})(t^{m})=(\mathcal{I}_{i})^{m-1}+\nu^{m}{\mathcal{I}}_{i}^{\max}\sqrt{t},~{}~{}\text{with}~{}~{}~{}i=1~{}~{}\text{for}~{}~{}c_{0},~{}~{}i=1,\cdots,d~{}~{}\text{for}~{}~{}{\mathbf{v}}_{0}$
(30)

where $\nu_{n}\in\\{-1,1\\}$. We then obtain ${\mathcal{I}}_{i}[0,T]$ via a
cubic spline interpolation. We refer the reader to [21, 20] for a discussion.
We clarify that $\\{t^{m}\\}$ is distinct from the time steps used for
generating the data. We consider $T=1$, $M=5$ and use 200 time-steps to
calculate the solution.

We emphasize that the data provided to the RNO is at the geological scale:
inputs (velocity and concentration history at a point), and output
(permeability, diffusivity, advection velocity, specific area, and pore volume
fraction). There is no information about the pore scale. The internal
variables are inferred from this data as a part of the training process.

After generating the data, we proceed to train the RNO. The training process
involves sequentially feeding the input data into the network, computing
outputs at each time step, and comparing them to the target values during the
forward propagation. We then use the backpropagation through time to calculate
gradients over the entire sequence, and optimize the parameters.

###  3.3 Results

Figure 2: Initial core-scale microstructure

We now demonstrate the ability of the RNO to approximate the core-scale
problem. We consider an initial microstructure shown in Figure 2. We generate
our data by solving the core-scale problem described in Section 2.2.2 using a
numerical algorithm described elsewhere [14].

We consider a fully connected $4$-layer neural network, with each layer
consisting of $200$ nodes. We use the nonlinear activation function scaled
exponential linear unit (SELU) [16], and optimize the parameters using the
ADAM optimization algorithm [15]. We consider the following loss function

${\mathcal{L}}=\frac{1}{D_{\text{train}}}\sum_{d=1}^{D_{\text{train}}}\frac{\int^{T}_{0}\left||\overline{\mathcal{O}}^{\text{truth}}_{d}-\overline{\mathcal{O}}^{\text{approx}}_{d}\right|^{2}dt}{\int^{T}_{0}\left||\overline{\mathcal{O}}^{\text{truth}}_{d}\right|^{2}dt},$
(31)

where $d$ indexes the trajectory in the training dataset and
$\overline{\mathcal{O}}$ is the normalized output. To compute this, we use
min-max normalization on each physical component of the output,

$(\overline{{\mathbf{D}}}^{*}_{ij})_{d}=\frac{({\mathbf{D}}^{*}_{ij})_{d}-\min{p,q,r}({\mathbf{D}}^{*}_{pq})_{r}}{\max{p,q,r}({\mathbf{D}}^{*}_{pq})_{r}-\min{p,q,r}({\mathbf{D}}^{*}_{pq})_{r}}$
(32)

and so forth. We define and train the RNO with $\log{\mathbf{K}}^{*}$ instead
of ${\mathbf{K}}^{*}$ to properly emphasize the almost clogged regime.

Recall that we have to fix the number of internal variables $k$ before
training. We do so, repeating the training for $k=0,1,\dots,4$. Thus we have
five fully trained RNOs with differing numbers of internal variables.

Figure 3:  Training and testing the RNO. (a) Training error vs. training
epochs, (b) Average test loss of the trained RNOs vs. the number of hidden
variables, (c) Average test loss of an RNO with one internal variable vs.
training set size. (d) Average normalized test error vs. the number of
internal variables.

The results are shown in Figure 3 in terms of the loss (31). Figure 3(a) shows
how the training loss changes with the number of epochs with varying number of
internal variables. We see that an RNO with no internal variables (i.e., no
history dependence) is unable to reduce the training error beyond a certain
point. However, RNOs with one or more internal variable can be trained to a
high degree of accuracy. Figure 3(b) shows the test loss (the same loss as
(31), but computed for the test data set) as a function of the number of
internal variables for the trained RNO. We see that the trained RNO with no
internal variable provides a very poor approximation, but RNOs with one or
more internal variable provides excellent approximation. We repeat the
training for various sizes of training data and the average test loss is shown
in Figure 3(c) for the case of a single internal variable. We see that the
size (800) of our training data set is adequate, and the average test error is
small in each component.

Figure 3(d) shows the average normalized test error of the various physical
outputs. For each physical output, we define the normalized test error as

$\text{average normalized test error in
${\mathbf{D}}^{*}$}=\left(\frac{1}{D_{\text{test}}}\sum_{d=1}^{D_{\text{test}}}\frac{\int^{T}_{0}\left|({\mathbf{D}}^{*})^{\text{truth}}_{d}-({\mathbf{D}}^{*})^{\text{approx}}_{d}\right|^{2}dt}{\int^{T}_{0}\left|({\mathbf{D}}^{*})^{\text{truth}}_{d}\right|^{2}dt}\right)^{1/2}$
(33)

and so forth. We observe that the average normalized test error about 1% for
each of the physical quantities except for the effective advection velocity
$\bar{{\mathbf{v}}}$ where the error is about 6% for a trained RNO with one or
more internal variables. The test includes cases where the effective advection
velocity is zero up to machine precision and these lead to large apparent
errors.

Figure 4 elaborates on the results by focussing on a typical trajectory chosen
arbitrarily from the test data set. Figures 4(a,b) show the input, while
Figures 4(c-g) compare the ground truth and RNO predictions for the outputs:
permeability (${\mathbf{K}}^{*}$), diffusivity (${\mathbf{D}}^{*}$), advection
velocity ($\bar{{\mathbf{v}}}$), specific area ($\gamma$), and pore volume
fraction ($\lambda$). We see excellent agreement.

Finally, we examine the time discretization independence of the trained RNO.
The RNO is initially trained with $\Delta t$. Figure 5 shows the predictions
of the RNO for permeability and diffusivity components evaluated with
different time steps – $0.25\Delta t$, $0.5\Delta t$, $\Delta t$, and $2\Delta
t$ – all for the same input trajectory. We see that the results are
independent of the time step.

In summary, we conclude that an RNO with one internal variable is able to
provide an excellent approximation to the solution operator of the core scale
model.

Figure 4: Input test trajectories, (a) concentration, (b) velocity
trajectories, comparison of estimated values of the RNO with one hidden
variable with ground truth, (c) permeability, (d) diffusivity, (e) advection
velocity, (f) specific area and (g) pore volume fraction.  Figure 5:
Comparison of estimated values of the RNO, permeability components (a)-(c),
and diffusivity components (d)-(f) considering various time steps.

##  4 Multiscale simulation

We now consider a geological scale simulation, but one that constantly updates
the properties, from core scale calculations according to the framework
described in Section 2. However, instead of solving the core scale problem at
each point at each instant, we use the trained neural approximation described
in Section 3 as a surrogate for the core scale problem at each point at each
instant. Of particular interest is to understand how the formation and its
properties, as well as the flow, change over long periods of time, and how
such changes are magnified by heterogeneities in the formation. We implement
the geological scale problem with the Python finite elements library, FEniCS,
and an unstructured mesh with triangular elements.

###  4.1 Reactive flow with uniform initial properties

Figure 6: Geometry of geological domain.

We first consider an example with uniform initial material properties. The
geometry of the geological formation and microstructure are shown in Figure 6:
a solution with high concentration of solute is injected at high pressure
through the well on the left and removed from a well at the right. The
properties and boundary conditions are as follows.

$\displaystyle\text{Pe}=1000,\ \ \text{Da}=0.001,\ \ c^{*}=0.5,$
$\displaystyle\text{left well:}~{}p_{0}=10^{5},~{}~{}c_{0}=0.6,$
$\displaystyle\text{right well:}~{}p_{0}=10^{4},~{}~{}c_{0}=0.4,$
$\displaystyle\text{rest of the boundary: zero flux}.$

Figure 7: Variation of concentration (a)-(c), and pressure (d)-(f), and pore
volume fraction (g)-(i) profiles in the geological formation due to chemical
reactions.

Figures 7 shows three snapshots of the pressure and concentration profile at
the geological scale, at $t_{1}$, $t_{2}$, and $t_{3}$ associated with time
steps $1$, $110$, and $220$, respectively. The pressure and concentration
change gradually away from the wells in the early times. However, at later
times, the pressure and concentration change rapidly in the vicinity of the
wells. Recall that the equations (7) at the geological scale are steady state
equations, i.e., do not involve time explicitly. Therefore, the evolution with
time is related to the change in properties related to the reconstruction of
the porous medium at the core scale.

The concentration prescribed at the inlet exceeds the equilibrium
concentration, resulting in the gradual precipitation of solutes on the
solid/fluid interface in this region. Conversely, the concentration prescribed
at the outlet lies below the equilibrium concentration, leading to a
dissolution of solid structure close to the outlet. As illustrated in Figure
7, precipitation in the inlet increases pore volume fraction that decreases
permeability and diffusivity values while increasing the pressure gradient
over time. This reduction in permeability and diffusivity reduces flow and
chemical transport in the medium, eventually leading to clogging.

Figure 8: Comparison of the (a) permeability, (b) diffusivity, (c) advection
velocity, (d) pore volume fraction, and (e) specific area components obtained
from the RNO with ground truth, and (f) microstructure geometry at $230$th
time step, at three random points.

We explore this further at the three points marked in Figure 6: Points 1 and 3
are in the vicinity of the left well (inlet), and right well (outlet),
respectively; Point 2 is located between the inlet and outlet. The changes in
properties and pore structure with time are shown in Figure 8 as the solid
curves. The left column of the figure shows the changes at Point 1 that is
close to the inlet. The deposition on the surface of solid leads to a decrease
in permeability, diffusivity, advection velocity and pore volume fraction. The
decrease is steady initially, but at some time the permeability, diffusivity
and advection velocity effectively goes to zero: this time coincides with when
the pores are blocked by deposition. The porosity does not go to zero, but
there is no change in its value beyond this time. The specific surface area
initially increases, but saturates earlier than full blocking. The deposition
leads to an increase of surface area, but this is eventually balanced by a
decrease as the pores are blocked and opposite sides of a pore begin to touch.

We see the opposite trends on the right column corresponding to Point 3 that
is close to the outlet. We see that the dissolution leads to an increase of
permeability, diffusivity, advection velocity and pore volume fraction, and
decrease of specific surface area. The middle column shows the results for
point 2 located between the inlet and outlet. This only experiences slight
changes in microstructure and effective properties, including permeability,
diffusivity, pore volume fraction, and specific area. The reduction in the
advection velocity at point 2 primarily results from the decrease in the flow
reaching this point due to the clogging of the inlet.

Figure 8 provides an a posteriori evaluation of the error in using the RNO
surrogate. We take the velocity and concentration experienced by the 3 points
during the macroscale calculations. We compare the output of our RNO surrogate
(solid lines) and the results of core scale calculations (dashed lines) for
these histories. We observe excellent overall agreement. This tells us that
our RNO surrogate performs well not only for the histories in the distribution
used to train the model, but also real histories experienced in actual
calculations.

Figure 9: Variation of flow flux ($q$, solid line) in inlet and outlet, and
chemical flux difference ($\Delta J$, dashed line) between inlet and outlet,
considering various values of Peclet and Damkohler numbers.

We now turn to the overall flow through the geometric formulation and the
amount of material deposited in the formation as a function of time. This is
shown in the center of Figure 9 for the parameters chosen above. We see that
the flow decreases and the amount of material deposited increases steadily
till they saturate. This is consistent with the observations above.

Figure 9 also shows a parameter study for various Péclet and Damköhler numbers
that reveal an interesting interplay between the pore and geological scale
phenomena. Recall that the core scale problem does not depend on the exact
values of these non-dimensional quantities (as long as they satisfy the
scaling). It follows that the RNO is independent of them, and we do not have
to retrain it for each pair of numbers. Instead, these non-dimensional
quantities only appear in the effective geological scale equations. So, we
repeat the geological scale calculation as described above for three distinct
values (differing by a factor 2) of both these non-dimensional quantities. We
see that the flow decreases and the amount of material deposited increases
steadily till they saturate in each case. The change in flow and amount of
material deposited increases significantly with increasing Péclet number, but
only slightly with increasing Damköhler number.

Recall that we have very high Péclet number and a small Damköhler number. This
means that we are in the reaction controlled regime. Therefore, one could
expect that the overall deposition rate would be more sensitive to the
Damköhler number, and less sensitive to the Péclet number. However, we see the
opposite trend in Figure 9. This is because the morphological changes in the
formation lead to clogging, and therefore the flow is transport limited at the
macroscopic scale. The inset in each graph in Figure 9 shows a snapshot of the
concentration at a fixed time: we see that the Péclet number induces a greater
effect than the Damköhler number. In other words, even when the core scale is
reaction limited, the formation can become transport limited.

Table 1: Comparison of computation cost (wall-clock time in seconds) Method | Calculation | Cost
---|---|---
Classical constitutive relations | Geological scale calculation | 800 (CPU))
Proposed method | Geological scale calculation |  $900$ (CPU)
| Off-line data (sequential) |  $3.6\times 10^{6}$ (CPU)
| Off-line data (parallel) |  $4\times 10^{4}$ (GPU)
| Off-line training |  $5\times 10^{3}$ (GPU)
Concurrent multiscale (estimated) | Geological scale calculation |  $1.64\times 10^{7}$ (CPU))

Having established the accuracy and its efficacy in parameter study, we study
the computational cost of the proposed approach and compare it with both the
cost of the classical empirical constitutive model and the cost of concurrent
multiscale models. The results are shown in Table 1 for the base simulation
above. These calculations were conducted on a single core of an Intel Skylake
CPU with a clock speed of 2.1 GHz and an NVIDIA P100 GPU with 3584 cores and
1.3 MHz clock speed. The classical empirical model is described in [5]. We
find that the computational cost of solving the macroscopic problem using the
trained RNO is comparable to the cost of classical constitutive relations.
These are significantly (by factor $10^{5}$) smaller than the estimated cost
of the concurrent multiscale approach (we estimate this by using the cost of
the unit cell problem and multiplying it by the product of the time steps and
spatial grid). The proposed approach has a one-time off-line cost of
generating the data and fitting the RNO. This is also smaller (by a factor
$~{}10^{2}$) than solving a single concurrent multiscale calculation. This
off-line cost can further be reduced by parallelization.

In summary, we conclude that our approach is able to provide the accuracy of a
concurrent multiscale model at the computational cost comparable to that of a
classical constitutive relation.

###  4.2 Reactive flow with non-uniform initial properties

We now turn to examples with geological formation characterized by initially
non-uniform properties as shown in Figure 10. This heterogeneity is defined by
the initial value for the internal variable in the geological scale problem.

Figure 10: Domain geometry and location of blocks.

####  4.2.1 High initial permeability and diffusivity inclusions

We consider a $4\times 8~{}\text{m}$ domain with embedded blocks of higher
initial permeability and diffusivity in the geometry illustrated in Figure 10.
A horizontal flow is introduced from the left boundary and withdrawn from the
right. The parameters are the same as in Section 4.1, and boundary conditions
are

$\displaystyle\text{left boundary:}~{}p_{0}=10^{5},~{}~{}c_{0}=0.4,$
$\displaystyle\text{right boundary:}~{}p_{0}=10^{4},~{}~{}c_{0}=0.5,$
$\displaystyle\text{rest of the boundary: zero flux}.$

Figures 11 (a), and (b) show the temporal evolution of the concentration and
pressure fields over time due to changes in the morphology. The left (inlet)
concentration is lower than the equilibrium concentration, and this leads to
dissolution at that end. The right (outlet) concentration is at the
equilibrium concentration, and so we do not see significant morphological
changes there. Note that the flow and concentration are not uniform across the
(vertical) cross-section even though the boundary conditions are. This is a
result of the inclusions. The flow preferentially enters the regions with high
permeability and diffusivity, and seeks to connect these regions together. The
greater flow leads to greater chemical reaction and further increases
permeability; this is clear from the evolution of the pore volume fraction
shown in Figure 11 (c)). This further aids the channeling from one inclusion
to another.

The differential change in permeability, diffusivity, and advection velocity
is emphasized in Figure 12 that shows the evolution of these quantities at two
points (marked 1 and 2 in Figure 10). The solid lines correspond to Point 1
located within the high permeability block, while the dashed lines represent
the change of effective properties at Point 2 located outside the high
permeability block. As expected, the rate of change in effective properties at
Point 1 is significantly greater compared to Point 2.

Figure 11: Variation of concentration (a), pressure (b), and pore volume
fraction (c) fields, considering non-uniform initial material properties.
Figure 12: Change in effective properties, (a) permeability, (b) diffusivity,
and (c) advection velocity, at points 1 (solid lines), and point 2 (dashed
lines).

####  4.2.2 Low initial permeability and diffusivity inclusions

We now consider the complementary situation with blocks of lower initial
permeability and diffusivity. The parameters are the same as in Section 4.1,
and boundary conditions are the same as Section 4.2.1. Figure 13 (a)-(b) shows
that the flow moves around the block following regions with higher initial
permeability and diffusivity. Consequently, there is more reaction and
dissolution in the matrix compared to the blocks, Figure 13(c). This further
channels the flow through between the blocks. Figure 14 shows the change in
material properties at Points 3 and 4 (as shown in figure 10). We observe a
substantial change in effective properties at Point 4 outside the blocks and
relatively little change at Point 3 inside.

Figure 13: Variation of concentration (a), pressure (b), and pore volume
fraction (c) fields, considering non-uniform initial material properties.
Figure 14: Change in effective properties, (a) permeability, (b) diffusivity,
and (c) advection velocity, at points 3 (dashed lines), and point 4 (solid
lines).

##  5 Conclusion

The transport of water and solutes through permeable geological formations is
a complex multiscale phenomenon. In this work, we have proposed a framework
that harnesses the power of continuum modeling and machine learning to address
this complexity with accuracy at reasonable computational cost. We have
demonstrated the framework on flow through permeable geological formations
with advection-dominated transport and volume reactions.

We begin with an asymptotic analysis of the governing equations that enable us
to split the problem into two interconnected problems: one at the scale of the
geological formation and another at the scale of a core. A key idea in this
analysis is the invocation of a drifting coordinate system to capture the
core. We then introduce a recurrent neural operator (RNO) to approximate the
solution operator of the core-scale problem. This consists of two feed-forward
neural networks and internal variables. The neural networks are trained and
the internal variables are identified from data that is generated by
repeatedly solving the core-scale problem. The key features of this neural
architecture are that it is consistent with common formulations of continuum
physics, that it is relatively simple and that it is independent of the time-
discretization. We demonstrate that it is able to accurately capture the
behavior of the core scale over long periods of time including the
morphological changes in the pores and the resulting change in effective
permeability, diffusivity, advective velocity, pore volume fraction and
specific area. Finally, we solve the problem of transport and morphological
changes at the scale of the geological formation by using the trained RNO as a
surrogate for the small scale problem. We thus obtain the accuracy of a
concurrent multiscale simulations at a cost that is comparable to classical
constitutive relations. We demonstrate the ability of this approach to learn
subtle features of the interaction between the scales including the change
from reaction limited to transport limited regime due to clogging and the
positive feedback of channeling in heterogeneous situations.

We now emphasize a few notable aspects of the proposed approach. First, our
RNO neural approximation is able to capture morphological and property changes
over long periods of time. It is formulated in a time-continuous manner and
discretized as necessary for training and use. It follows that the
approximation is independent of the discretization. This is important because
it is common to use different time discretization for the core scale problem
used to generate data and the application at the geological scale. It also
enables the use of data generated at different discretizations. Second, the
two scale formulation makes the core scale problem to be independent of the
physical characteristics of the flow and reaction rate, specifically the
Péclet and the Damköhler numbers (as long as we are in the advection dominated
regime with slow reactions). This means that we need to generate data and
train the RNO only once, and can use the method for different situations as
these quantities change. Third, we not only demonstrate accuracy over the
distribution used to train the RNO, but also a posteriori over the actual
histories encountered in the geological scale calculations. Finally, the
approach is highly transferable. In this work, we used examples in two
dimensions to demonstrate the framework with modest computational cost.
However, the approach holds in three dimensions. Similarly, one can
incorporate other phenomena including, for example, multi-phase flows, more
complex chemistry, poroelasticity and phase change (melting), as long as we
can use scale separation. Importantly, one can extend this to more than two
scales as long as they interact pairwise.

We close with a discussion of a few limitations and avenues for future work.
First, the approach requires us to train the RNO with a starting pore
morphology. This is not an issue as long as the core scale is chosen to be
large enough to be statistically representative of the underlying medium.
However, this adds to the computational cost. So, one may consider training
the RNO on an ensemble of cores. Second, while we demonstrate a posteriori
accuracy, it would be useful to have a systematic approach to quantifying the
overall uncertainties. Such a quantification may also enable the use of active
learning where we initially train the RNO over synthetic samples as we do in
this work, but then progressively add more data based on examples we
encounter. Third, in this work, we only use geological scale information
obtained by averaging the results of core scale simulations to train the RNO.
However, we have access to core scale information. It is possible that this
data has insights that may lead to a more robust efficient training procedure.
Fourth, we have exclusively used data generated numerically to train the RNO.
it would be interesting to use a combination of experimental and computational
data. Finally, it remains to use the framework established in this work to
actual geological problems. All of this remains a topic of current and future
research.

## Data availability

The data is available at CaltechData through
https://doi.org/10.22002/yd0c5-q5s36

## Acknowledgments

We are delighted to acknowledge numerous helpful discussions with Professor
Xiaojing (Ruby) Fu. We gratefully acknowledge the financial support of the
Resnick Sustainability Institute at the California Institute of Technology. KB
also acknowledges the support of the Army Research Office through grant number
W911NF-22-1-0269. The simulations reported here were conducted on the Resnick
High Performance Computing Cluster at the California Institute of Technology.

## References

  * [1] G. Allaire, R. Brizzi, A. Mikelić, and A. Piatnitski.  Two-scale expansion with drift approach to the taylor dispersion for reactive transport through porous media.  Chemical Engineering Science, 65:2292–2300, 2010.
  * [2] G. Allaire and A.-L. Raphael.  Homogenization of a convection–diffusion model with reaction in a porous medium.  Comptes Rendus Mathematique, 344:523–528, 2007.
  * [3] C. A. J. Appelo and D. Postma.  Geochemistry, groundwater and pollution.  CRC Press, 2004.
  * [4] J.-L. Auriault and P. Adler.  Taylor dispersion in porous media: analysis by multiple scale expansions.  Advances in Water Resources, 18:217–226, 1995.
  * [5] Y. Baqer and X. Chen.  A review on reactive transport model and porosity evolution in the porous media.  Environmental Science and Pollution Research, 29:47873–47901, 2022.
  * [6] J. Bear and Y. Bachmat.  Introduction to modeling of transport phenomena in porous media.  Springer Science & Business Media, 2012.
  * [7] K. Bhattacharya, B. Liu, A. Stuart, and M. Trautner.  Learning markovian homogenized models in viscoelasticity.  Multiscale Modeling & Simulation, 21:641–679, 2023.
  * [8] R. Carbonell and S. Whitaker.  Dispersion in pulsed systems. II: Theoretical developments for passive dispersion in porous media.  Chemical Engineering Science, 38:1795–1802, 1983.
  * [9] E. Demirer, E. Coene, A. Iraola, A. Nardi, E. Abarca, A. Idiart, G. de Paola, and N. Rodríguez-Morillas.  Improving the performance of reactive transport simulations using artificial neural networks.  Transport in Porous Media, 149:1–27, 2022.
  * [10] P. Donato and A. Piatnitskí.  Averaging of nonstationary parabolic operators with large lower order terms.  In Multi Scale Problems and Asymptotic Analysis, pages 153–165. Gakkotosho, 2006.
  * [11] F. Ghavamian and A. Simone.  Accelerating multiscale finite element simulations of history-dependent materials using a recurrent neural network.  Computer Methods in Applied Mechanics and Engineering, 357:112594, 2019.
  * [12] Q. Kang, D. Zhang, and S. Chen.  Simulation of dissolution and precipitation in porous media.  Journal of Geophysical Research: Solid Earth, 108:2505, 2003.
  * [13] Q. Kang, D. Zhang, S. Chen, and X. He.  Lattice Boltzmann simulation of chemical dissolution in porous media.  Physical Review E, 65:036318, 2002.
  * [14] M. Karimi and K. Bhattacharya.  A fast-Fourier transform method for reactive flow in porous media.  In preparation, 2023.
  * [15] D. Kingma and J. Ba.  ADAM: A method for stochastic optimization.  In International conference on learning representations (ICLR), volume 5, page 6, 2015.
  * [16] G. Klambauer, T. Unterthiner, A. Mayr, and S. Hochreiter.  Self-normalizing neural networks.  Advances in Neural Information Processing Systems, 30, 2017.
  * [17] T. Lévy.  Fluids in porous media: Darcy’s law.  In Homogenization Techniques for Composite Media, pages 75–84. Springer, 1985.
  * [18] X. Li, H. Huang, and P. Meakin.  Level set simulation of coupled advection-diffusion and pore structure evolution due to mineral precipitation in porous media.  Water Resources Research, 44, 2008.
  * [19] P. C. Lichtner, C. I. Steefel, and E. H. Oelkers.  Reactive transport in porous media, volume 34.  Walter de Gruyter GmbH & Co KG, 2018.
  * [20] B. Liu, N. Kovachki, Z. Li, K. Azizzadenesheli, A. Anandkumar, A. M. Stuart, and K. Bhattacharya.  A learning-based multiscale method and its application to inelastic impact problems.  Journal of the Mechanics and Physics of Solids, 158:104668, 2022.
  * [21] B. Liu, E. Ocegueda, M. Trautner, A. M. Stuart, and K. Bhattacharya.  Learning macroscopic internal variables and history dependence from microscopic models.  Journal of the Mechanics and Physics of Solids, 178:105329, 2023.
  * [22] M. Liu, B. Kwon, and P. K. Kang.  Machine learning to predict effective reaction rates in 3d porous media from pore structural features.  Scientific Reports, 12:5486, 2022.
  * [23] H. Lu, D. Ermakova, H. M. Wainwright, L. Zheng, and D. M. Tartakovsky.  Data-informed emulators for multi-physics simulations.  Journal of Machine Learning for Modeling and Computing, 2(2), 2021.
  * [24] E. Marušić-Paloka and A. L. Piatnitski.  Homogenization of a nonlinear convection-diffusion equation with rapidly oscillating coefficients and strong convection.  Journal of the London Mathematical Society, 72:391–409, 2005.
  * [25] R. Mauri.  Dispersion, convection, and reaction in porous media.  Physics of Fluids A: Fluid Dynamics, 3:743–756, 1991.
  * [26] G. Medici, L. Smeraglia, A. Torabi, and C. Botter.  Review of modeling approaches to groundwater flow in deformed carbonate aquifers.  Groundwater, 59:334–351, 2021.
  * [27] L. R. Medsker and L. Jain.  Recurrent neural networks.  Design and Applications, 5(64-67):2, 2001.
  * [28] M. Mozaffar, R. Bostanabad, W. Chen, K. Ehmann, J. Cao, and M. Bessa.  Deep learning predicts path-dependent plasticity.  Proceedings of the National Academy of Sciences, 116:26414–26420, 2019.
  * [29] A. P. Oron and B. Berkowitz.  Flow in rock fractures: The local cubic law assumption reexamined.  Water Resources Research, 34:2811–2825, 1998.
  * [30] M. Paine, R. Carbonell, and S. Whitaker.  Dispersion in pulsed systems—i: Heterogenous reaction and reversible adsorption in capillary tubes.  Chemical Engineering Science, 38:1781–1793, 1983.
  * [31] M. Quintard and S. Whitaker.  Transport in ordered and disordered porous media: volume-averaged equations, closure problems, and comparison with experiment.  Chemical Engineering Science, 48:2537–2564, 1993.
  * [32] J. Rubinstein and R. Mauri.  Dispersion and convection in periodic porous media.  SIAM Journal on Applied Mathematics, 46:1018–1023, 1986.
  * [33] N. Seigneur, K. U. Mayer, and C. I. Steefel.  Reactive transport in evolving porous media.  Reviews in Mineralogy and Geochemistry, 85:197–238, 2019.
  * [34] R. Sprocati and M. Rolle.  Integrating process-based reactive transport modeling and machine learning for electrokinetic remediation of contaminated groundwater.  Water Resources Research, 57:e2021WR029959, 2021.
  * [35] C. I. Steefel and P. C. Lichtner.  Multicomponent reactive transport in discrete fractures: I. controls on reaction front geometry.  Journal of Hydrology, 209:186–199, 1998.
  * [36] S. Sun and M. F. Wheeler.  A dynamic, adaptive, locally conservative, and nonconforming solution strategy for transport phenomena in chemical engineering.  Chemical Engineering Communications, 193:1527–1545, 2006.
  * [37] A. M. Tartakovsky, N. Trask, K. Pan, B. Jones, W. Pan, and J. R. Williams.  Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media.  Computational Geosciences, 20:807–834, 2016.
  * [38] G. I. Taylor.  Dispersion of soluble matter in solvent flowing slowly through a tube.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 219:186–203, 1953.
  * [39] C. Varloteaux, M. T. Vu, S. Békri, and P. M. Adler.  Reactive transport in porous media: pore-network model approach compared to pore-scale model.  Physical Review E, 87:023010, 2013.
  * [40] K. Wang, Y. Chen, M. Mehana, N. Lubbers, K. C. Bennett, Q. Kang, H. S. Viswanathan, and T. C. Germann.  A physics-informed and hierarchically regularized data-driven model for predicting fluid flow through porous media.  Journal of Computational Physics, 443:110526, 2021.
  * [41] Z. Wang and I. Battiato.  Upscaling reactive transport and clogging in shale microcracks by deep learning.  Water Resources Research, 57(4):e2020WR029125, 2021.
  * [42] Z. Wang, L. Chen, H. Wei, Z. Dai, Q. Kang, and W.-Q. Tao.  Pore-scale study of mineral dissolution in heterogeneous structures and deep learning prediction of permeability.  Physics of Fluids, 34, 2022.
  * [43] P. A. Witherspoon, J. S. Wang, K. Iwai, and J. E. Gale.  Validity of cubic law for fluid flow in a deformable rock fracture.  Water Resources Research, 16:1016–1024, 1980.
  * [44] L. Wu, N. G. Kilingar, L. Noels, et al.  A recurrent neural network-accelerated multi-scale model for elasto-plastic heterogeneous materials subjected to random cyclic and non-proportional loading paths.  Computer Methods in Applied Mechanics and Engineering, 369:113234, 2020.


without the temporal notions we study. The ML approaches aim to learn
mitigation strategies trained on a fixed input graph which limits their
transferability to new graph instances. Recently, Juul et al. [16] repeat the
comparison of the spread of fake and true news conducted by [46] on two
subsampled datasets with exactly the same size distribution. Interestingly,
they find that under these conditions, the propagation characteristics become
indistinguishable, which has important consequences for misinformation
detection. We note that our work focuses on the general case where the
propagation differences observed by Vosoughi et al. [46] do hold. In
particular, the authors of [16] confirm that false and true information
propagate at different rates when the naturally occurring cascade size
distributions remain unaltered.

##  9 Conclusion and Future Work

In this paper, we address a major shortcoming of existing MM propagation
models by introducing the TCIC model, which captures important temporal
aspects of fake news diffusion and formulate a time-sensitive variant of the
MM problem. We prove our mitigation objective is non-submodular and develop
submodular upper and lower bounding functions to sandwich the objective and
provide data-dependent approximation guarantees. Finally, we propose a reverse
sampling framework that provides $(1-1/e-\epsilon)$-approximate solutions to
our bounding functions and present an anytime version of our approach. Using
experiments over five datasets, we demonstrate that our _NAMM_ algorithm
outperforms various baselines including those that are oblivious to time-
critical aspects.

The new techniques developed here are applicable to other problem settings.
E.g., in the filter bubble problem, multiple conflicting opinions propagate in
a social network and the goal is to ensure balanced exposure. Here, the set of
users adopting any particular opinion is stochastic and thus our _NAMM_
algorithm can be applied to help obtain a more balanced exposure for the
users. Further, revisit the “classic” competitive IM problem, where
realistically the diffusion rates of companies/brands may differ owing to
differences in reputation and marketing strategies. In this case, _NAMM_ can
be applied for the competitive IM problem from a follower perspective. Notice
that previous techniques for competitive IM do not apply to this setting.

For simplicity of exposition, we considered a fixed set of fake seeds.
However, our proposed solution retains its guarantees even when the fake seeds
are not known exactly but are chosen from a distribution. It is interesting to
study the scenario where the fake seed set is dynamically growing. This is
closely tied with adaptive influence maximization, which we leave for future
work.

## References

  * [1] G. M. Beal and J. M. Bohlen.  The diffusion process.  Technical report, 1956.
  * [2] S. Bhagat, A. Goyal, and L. V. Lakshmanan.  Maximizing product adoption in social networks.  In Proceedings of the fifth ACM international conference on Web search and data mining, pages 603–612, 2012.
  * [3] S. Bharathi, D. Kempe, and M. Salek.  Competitive influence maximization in social networks.  In WINE’07, pages 306–311, Berlin, Heidelberg, 2007. Springer-Verlag.
  * [4] A. A. Bian, J. M. Buhmann, A. Krause, and S. Tschiatschek.  Guarantees for greedy maximization of non-submodular functions with applications.  In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 498–507. JMLR. org, 2017.
  * [5] C. Borgs, M. Brautbar, J. Chayes, and B. Lucier.  Maximizing social influence in nearly optimal time.  In Proceedings of the twenty-fifth annual ACM-SIAM symposium on Discrete algorithms, pages 946–957. SIAM, 2014.
  * [6] C. Budak, D. Agrawal, and A. El Abbadi.  Limiting the spread of misinformation in social networks.  In WWW’11, 2011.
  * [7] W. Chen, L. V. S. Lakshmanan, and C. Castillo.  Information and Influence Propagation in Social Networks.  Synthesis Lectures on Data Management. Morgan & Claypool Publishers, 2013.
  * [8] W. Chen, W. Lu, and N. Zhang.  Time-critical influence maximization in social networks with time-delayed diffusion process.  In AAAI, volume 2012, pages 1–5, 2012.
  * [9] D. R. Cox and E. J. Snell.  A general definition of residuals.  Journal of the Royal Statistical Society: Series B (Methodological), 30(2):248–265, 1968.
  * [10] L. Fan, Z. Lu, W. Wu, B. Thuraisingham, H. Ma, and Y. Bi.  Least cost rumor blocking in social networks.  In Distributed Computing Systems (ICDCS), 2013 IEEE 33rd International Conference on, pages 540–549. IEEE, 2013.
  * [11] M. Farajtabar, J. Yang, X. Ye, H. Xu, R. Trivedi, E. Khalil, S. Li, L. Song, and H. Zha.  Fake news mitigation via point process based intervention.  In International Conference on Machine Learning, pages 1097–1106. PMLR, 2017.
  * [12] M. Gabielkov, A. Ramachandran, A. Chaintreau, and A. Legout.  Social clicks: What and who gets read on twitter?  In Proceedings of the 2016 ACM SIGMETRICS international conference on measurement and modeling of computer science, pages 179–192, 2016\.
  * [13] X. He, G. Song, W. Chen, and Q. Jiang.  Influence blocking maximization in social networks under the competitive linear threshold model.  In Proceedings of the 2012 SIAM International Conference on Data Mining, pages 463–474. SIAM, 2012.
  * [14] K. Huang, S. Wang, G. Bevilacqua, X. Xiao, and L. V. Lakshmanan.  Revisiting the stop-and-stare algorithms for influence maximization.  Proceedings of the VLDB Endowment, 10(9):913–924, 2017.
  * [15] D. Hughes-Hallett, A. M. Gleason, and W. G. McCallum.  Calculus: Single and multivariable.  John Wiley & Sons, 2020.
  * [16] J. L. Juul and J. Ugander.  Comparing information diffusion mechanisms by matching on cascade size.  Proceedings of the National Academy of Sciences, 118(46), 2021.
  * [17] S. Kalish.  A new product adoption model with price, advertising, and uncertainty.  Management science, 31(12):1569–1585, 1985.
  * [18] D. Kempe, J. Kleinberg, and E. Tardos.  Maximizing the spread of influence through a social network.  In KDD’03, 2003.
  * [19] E. B. Khalil, B. Dilkina, and L. Song.  Scalable diffusion-aware optimization of network topology.  In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 1226–1235, 2014.
  * [20] A. Khan.  Towards time-discounted influence maximization.  In Proceedings of the 25th ACM International on Conference on Information and Knowledge Management, pages 1873–1876, 2016.
  * [21] J. Kim, B. Tabibian, A. Oh, B. Schölkopf, and M. Gomez-Rodriguez.  Leveraging the crowd to detect and reduce the spread of fake news and misinformation.  In Proceedings of the eleventh ACM international conference on web search and data mining, pages 324–332, 2018.
  * [22] Y. Li, J. Fan, Y. Wang, and K.-L. Tan.  Influence maximization on social graphs: A survey.  IEEE Transactions on Knowledge and Data Engineering, 30(10):1852–1872, 2018.
  * [23] Y. Li, D. Zhang, and K.-L. Tan.  Real-time targeted influence maximization for online advertisements.  Proceedings of the VLDB Endowment, 2015.
  * [24] Y. Lin and J. C. Lui.  Analyzing competitive influence maximization problems with partial information: An approximation algorithmic framework.  Performance Evaluation, 91:187–204, 2015.
  * [25] B. Liu, G. Cong, D. Xu, and Y. Zeng.  Time constrained influence maximization in social networks.  In 2012 IEEE 12th international conference on data mining, pages 439–448. IEEE, 2012.
  * [26] W. Lu, W. Chen, and L. V. Lakshmanan.  From competition to complementarity: comparative influence diffusion and maximization.  Proceedings of the VLDB Endowment, 9(2):60–71, 2015.
  * [27] A. Mitchell, G. Stocking, and K. E. Matsa.  Long-form reading shows signs of life in our mobile news world.  Pew Research Center, 5, 2016.
  * [28] H. T. Nguyen, T. P. Nguyen, T. N. Vu, and T. N. Dinh.  Outward influence and cascade size estimation in billion-scale networks.  Proceedings of the ACM on Measurement and Analysis of Computing Systems, 1(1):1–30, 2017.
  * [29] H. T. Nguyen, M. T. Thai, and T. N. Dinh.  Stop-and-stare: Optimal sampling algorithms for viral marketing in billion-scale networks.  In Proceedings of the 2016 International Conference on Management of Data, pages 695–710, 2016.
  * [30] N. P. Nguyen, G. Yan, M. T. Thai, and S. Eidenbenz.  Containment of misinformation spread in online social networks.  In Proceedings of the 4th Annual ACM Web Science Conference, WebSci ’12, pages 213–222, New York, NY, USA, 2012. ACM.
  * [31] C. V. Pham, Q. V. Phu, and H. X. Hoang.  Targeted misinformation blocking on online social networks.  In Asian Conference on Intelligent Information and Database Systems, pages 107–116. Springer, 2018.
  * [32] C. V. Pham, Q. V. Phu, H. X. Hoang, J. Pei, and M. T. Thai.  Minimum budget for misinformation blocking in online social networks.  Journal of Combinatorial Optimization, 38(4):1101–1127, 2019.
  * [33] B. A. Prakash, L. Adamic, T. Iwashyna, H. Tong, and C. Faloutsos.  Fractional immunization in networks.  In Proceedings of the 2013 SIAM International Conference on Data Mining, pages 659–667. SIAM, 2013.
  * [34] A. Saxena, W. Hsu, M. L. Lee, H. Leong Chieu, L. Ng, and L. N. Teow.  Mitigating misinformation in online social network with top-k debunkers and evolving user opinions.  In Companion Proceedings of the Web Conference 2020, pages 363–370, 2020.
  * [35] M. Simpson, V. Srinivasan, and A. Thomo.  Clearing contamination in large networks.  IEEE Transactions on Knowledge and Data Engineering, 28(6):1435–1448, 2016.
  * [36] M. Simpson, V. Srinivasan, and A. Thomo.  Reverse prevention sampling for misinformation mitigation in social networks.  In 23rd International Conference on Database Theory (ICDT 2020). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2020.
  * [37] C. Song, W. Hsu, and M. L. Lee.  Temporal influence blocking: Minimizing the effect of misinformation in social networks.  In 2017 IEEE 33rd International Conference on Data Engineering (ICDE), pages 847–858. IEEE, 2017.
  * [38] J. Tang, X. Tang, X. Xiao, and J. Yuan.  Online processing algorithms for influence maximization.  In Proceedings of the 2018 International Conference on Management of Data, pages 991–1005, 2018.
  * [39] Y. Tang, Y. Shi, and X. Xiao.  Influence maximization in near-linear time: A martingale approach.  In Proceedings of the 2015 ACM SIGMOD International Conference on Management of Data, pages 1539–1554. ACM, 2015.
  * [40] Y. Tang, X. Xiao, and Y. Shi.  Influence maximization: Near-optimal time complexity meets practical efficiency.  In Proceedings of the 2014 ACM SIGMOD International Conference on Management of Data, SIGMOD ’14, pages 75–86, New York, NY, USA, 2014. ACM.
  * [41] A. Tong, D.-Z. Du, and W. Wu.  On misinformation containment in online social networks.  In Advances in neural information processing systems, pages 341–351, 2018.
  * [42] G. Tong.  Stratlearner: Learning a strategy for misinformation prevention in social networks.  In Advances in Neural Information Processing Systems, volume 33, pages 15546–15555, 2020.
  * [43] G. Tong, W. Wu, L. Guo, D. Li, C. Liu, B. Liu, and D.-Z. Du.  An efficient randomized algorithm for rumor blocking in online social networks.  IEEE Transactions on Network Science and Engineering, 7(2):845–854, 2017.
  * [44] G. A. Tong and D.-Z. Du.  Beyond uniform reverse sampling: A hybrid sampling technique for misinformation prevention.  In IEEE INFOCOM 2019-IEEE Conference on Computer Communications, pages 1711–1719. IEEE, 2019.
  * [45] N. Vo and K. Lee.  The rise of guardians: Fact-checking url recommendation to combat fake news.  In The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval, pages 275–284, 2018.
  * [46] S. Vosoughi, D. Roy, and S. Aral.  The spread of true and false news online.  Science, 359(6380):1146–1151, 2018.
  * [47] M. M. Vu and H. X. Hoang.  Minimizing the spread of misinformation on online social networks with time and budget constraint.  In 2017 9th International Conference on Knowledge and Systems Engineering (KSE), pages 160–165. IEEE, 2017.
  * [48] C. Wang, W. Chen, and Y. Wang.  Scalable influence maximization for independent cascade model in large-scale social networks.  Data Mining and Knowledge Discovery, 25(3):545–576, 2012.
  * [49] H. Wu, Z. Zhang, Y. Fang, S. Zhang, Z. Jiang, J. Huang, and P. Li.  Containment of rumor spread by selecting immune nodes in social networks.  Mathematical Biosciences and Engineering, 18(3):2614–2631, 2021\.
  * [50] Y. Zhang and B. A. Prakash.  Data-aware vaccine allocation over large networks.  ACM Transactions on Knowledge Discovery from Data (TKDD), 10(2):1–32, 2015.