# Estimating permeability of 3D micro-CT images by physics-informed CNNs based
on DNS

Stephan Gärttner1,*, Faruk O. Alpak2, Andreas Meier1,
Nadja Ray1, Florian Frank1,*
1 Friedrich-Alexander-Universität Erlangen-Nürnberg, Department Mathematik,
Cauerstraße 11, 91058 Erlangen, Germany 2 Shell Technology Center, 3333
Highway 6 South, Houston, TX 77082, USA

###### Abstract

In recent years, convolutional neural networks (CNNs) have experienced an
increasing interest for their ability to perform fast approximation of
effective hydrodynamic parameters in porous media research and applications.
This paper presents a novel methodology for permeability prediction from
micro-CT scans of geological rock samples. The training data set for CNNs
dedicated to permeability prediction consists of permeability labels that are
typically generated by classical lattice Boltzmann methods (LBM) that simulate
the flow through the pore space of the segmented image data. We instead
perform direct numerical simulation (DNS) by solving the stationary Stokes
equation in an efficient and distributed-parallel manner. As such, we
circumvent the convergence issues of LBM that frequently are observed on
complex pore geometries, and therefore, improve on the generality and accuracy
of our training data set.

Using the DNS-computed permeabilities, a physics-informed CNN (PhyCNN) is
trained by additionally providing a tailored characteristic quantity of the
pore space. More precisely, by exploiting the connection to flow problems on a
graph representation of the pore space, additional information about confined
structures is provided to the network in terms of the maximum flow value,
which is the key innovative component of our workflow. As a result,
unprecedented prediction accuracy and robustness are observed for a variety of
sandstone samples from archetypal rock formations.

††footnotetext: * Corresponding authors. _E-mail address:_
<EMAIL_ADDRESS><EMAIL_ADDRESS>
_Keywords:_ digital rock, neural networks, deep learning, permeability, porous
media.

_MSC classification:_ 05C21, 68T07, 76D07, 76M10, 76S05.

##  1 Introduction

Artificial neural networks can accelerate or replace classical methods for
estimating various hydrological and petrophysical properties of artificial and
natural rock [1]. In [2], deep neural networks were successfully used to
approximate tomography operators to reconstruct wave velocity models from
seismic data. Likewise, convolutional neural networks (CNNs) have proven their
ability to provide fast predictions of scalar permeability values directly
from images of the pore space of geological specimens [3, 4]. This paper aims
at contributing to the last-mentioned application of machine learning
techniques by presenting a finite-element-based forward simulation approach as
well as introducing a novelly considered characteristic quantity for physics-
informed neural network (PhyCNN) models.

Specimens of natural rock are typically obtained by microcomputed tomography
(µCT) scanning [5, 6] or indirectly from 2D colored images [7]. As an
important characteristic quantity of porous media, permeability measures the
resistance of a fluid to travel through a considered pore space. However, in
general, the precision and generality of the estimates performed by neural
networks heavily depend on the quality of the underlying training data set as
shown in [8]. In terms of permeability prediction, the labeling process of
geological specimens is connected to the computation of 3D stationary flow
fields of a single-phase fluid within the pore space. For permeability
prediction under regular geometries, a broader set of methods is available as
described and compared in [9].

In the literature, various methods are available to compute the stationary
flow on complex geometries such as pore spaces of porous media, a thorough
comparison of which is found in [10]. Most commonly, lattice Boltzmann methods
(LBM) are used to tackle the flow problem on a discrete modeling basis, cf.
[11]. In this approach, a transient interacting many-particle system is driven
to equilibrium state, heavily exploiting the inherent parallelism of the
underlying mathematical structure. For a detailed description of LBM
fundamentals and numerics, we refer to [12]. Yet, since a global equilibrium
has to be reached, complex geometries containing thin channels may cause LBM
to converge extremely slowly or even diverge [13, 4]. This behavior is
frequently observed with common pore-scale flow LBM implementations that rely
on the conventional single-relaxation-time Bhatnagar–Gross–Krook (BGK) scheme
[14]. However, the novel multiple-relaxation-time (MRT) scheme addresses this
problem to a great extent (e. g., [15]) at the cost of a moderate additional
computational overhead. We further note that classical bounce-back rules used
for implementing boundary conditions tend to develop boundary layers that
might locally dominate fluid behavior in scenarios with thin channels, cf.
[16]. Severe limitations arising from simplistic implementations of LBM
restrict the training data sets to ones derived from sphere-packs instead of
real CT data in many publications, see for instance [13]. On the other hand,
industrial-grade state-of-the-art LBM solvers are typically proprietary
without public access to the code base for academic research purposes (e. g.,
[17, 15]). Moreover, a comprehensive permeability computation benchmarking
study [18] demonstrated that not all of the simulation methods deliver a good
compromise of accuracy against computational performance indicating that there
exists a clear need for accurate and computationally efficient alternative
methods for permeability computation.

In practice, LBM simulations are often aborted after a maximal number of
iterations in case one or more convergence criteria (typically linked to the
relative changes in the velocity and/or computed permeability) cannot be met
[4]. As such, training sets for neural networks can be artificially filtered
by the numerical characteristics of the forward simulation, potentially
resulting in biased data. To improve the generality and quality of our PhyCNN
training sets, we base our machine learning data set on the stationary Stokes
equation by performing direct numerical simulations (DNS) on the pore geometry
[19]. More precisely, our forward simulation is based on a distributed-
parallel Stokes solver utilizing the finite element library MFEM [20]. As
studied in [9] for simple cylindrical obstacles, such DNS approaches (in this
case FEM) deliver results comparable to LBM. However, our implementation
successfully alleviates the drawback of an impractically large number of
iterations to obtain the desired accuracy on complex geometries. As such, our
approach allows overcoming prior restrictions in setting up representative
training sets including also confined and complex structures. Moreover, no
artificial data augmentation schemes such as pore space dilation are needed in
our approach to increase the number of training data samples or enhance the
porosity range covered.

Likewise, novel strategies are employed to our neural networks. More
precisely, we exploit the concept of PhyCNNs, where the CNN is provided with
additional (physics-related) input quantities to improve the reliability and
the accuracy of its predictions [3, 21]. As shown in [22], carefully chosen
specific quantities derived from the pore space such as connectivity indices
can deliver reasonable approximation quality for permeability estimation. As
illustrated in [2], also training a network solely on previously extracted
features from the raw data may lead to satisfactory prediction quality.

The outstanding performance of our methodology is achieved by considering the
_maximum flow_ value, a graph-network derived quantity being highly correlated
to the target permeability value and simple to compute. As such, we solve
maximum flow problems on a graph representation of the pore space based on
[23]. Thus, we approximate the Stokes flow through the pore space by an
abstract flow through a graph. By using the scalar quantity of maximum flow as
a second input to our neural network, we additionally provide our CNN with
information that reflects possible thin channel-like structures. As discussed
in [24], graph representations have been shown highly capable of
characterizing the pore-space connectivity in fractured rock and allow for a
convenient way to deduce topological quantities of interest. We demonstrate
that involving the maximum flow value in our newly designed PhyCNN in
combination with the DNS-based forward simulation approach, our methodology
delivers superior prediction accuracy and robustness compared to what is found
in the literature.

The paper is organized as follows: In Section 2, we describe the sampling and
preprocessing procedure of sandstone specimens including the forward
simulation. Subsequently, Section 3 is dedicated to the network architecture
used in our study. Finally, we validate the training performance of our PhyCNN
on different types of sandstone in Section 4.

##  2 Methodology and data preparation

In this section, we describe the workflow and methodology by which we acquire
the data set necessary to train and validate a CNN using a supervised learning
approach. The preprocessing includes the selection and preparation of a set of
pore-space geometries in form of voxel sets and the labeling with their
computed permeability value $k_{\mathrm{cmp}}$. Our complete workflow is
presented in Figure 1. We note that all steps except for the data labeling
procedure (green box) are implemented in Matlab 2021a [25].

Data acquisition • Extract $100^{3}$-voxel subsamples. • Remove disconnected
pore space. • Eliminate impermeable subsamples. Data labeling Perform forward
simulation to determine permeability labels $k_{\mathrm{cmp}}$. Network
training • Compute physics input. • Setup training and validation data set. •
Optimize network parameters via SGD. Network validation against • unseen data
samples of same sandstone type, • data samples of different sandstone type, •
artificially distorted data. Figure 1: Overall workflow in flow chart
representation.

###  2.1 Sampling and preprocessing

The training procedure of our PhyCNN is based on a segmented X-ray µCT scan of
a Bentheimer sandstone sample, see Figure 2, with experimentally measured
porosity $\phi_{\mathrm{exp}}=22.64\%$ and permeability
$k_{\mathrm{exp}}=$386\text{\,}\mathrm{m}\mathrm{D}$$, provided by [26, 27].
Further characteristic quantities for this sandstone as well as two other
sandstone types used below for validation purposes are listed in Table 1.
Bentheimer sandstone is known to exhibit a broad range in pore volume
distribution and high pore connectivity in comparison to other types of
sandstone [28]. As such, this sample is expected to contain a representative
collection of geometrical and topological properties of the pore space in
natural rocks. The data set used in this paper is a
$$1000${\times}$1000${\times}$1000$$ binary voxel image, in which each voxel
either belongs to the pore space (“fluid voxels”) or the rock matrix (“solid
voxels”). The voxel edge length is $2.25\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ yielding an overall cube side length of $2.25\text{\,}\mathrm{mm}$.

Figure 2: 10002 2D slices and 10003 3D pore-space µCT image of all samples considered in this study, illustrating characteristic pore features for Bentheimer, Berea, and Castlegate sandstone. type | $\phi_{\mathrm{exp}}$ | $\phi_{\mathrm{cmp}}$ |  $k_{\mathrm{exp}}$ [mD] | MCD [$\mathrm{\SIUnitSymbolMicro m}$] |  $A_{\mathrm{cmp}}$ [$\mathrm{mm}$]
---|---|---|---|---|---
Bentheimer | 22.64% | 26.72% | 386 | 30.0 | 355
Berea | 18.96% | 21.67% | 121 | 22.3 | 284
Castlegate | 26.54% | 24.67% | 269 | 24.7 | 335
Table 1: Characteristics quantities of sandstone samples used in this paper:
experimentally determined porosity, $\phi_{\mathrm{exp}}$ ($\pm 0.5\%$),
porosity computed from the µCT scan, $\phi_{\mathrm{cmp}}$, experimentally
determined permeability, $k_{\mathrm{exp}}$ ($\pm 10\%$), mean capillary
diameter, MCD, and interior surface area, $A_{\mathrm{cmp}}$, computed via the
Matlab function isosurface. First four quantities are provided in [27].

In the first step, we extract subsamples of $100{\times}100{\times}100$ voxels
from the Bentheimer sample. Henceforth, we use the term ‘subsample’ to refer
to segmented µCT-scan pieces of this specific size. For subsample extraction,
we make use of the sliding frame technique. This approach is commonly used to
further exploit a given data set beyond the partition into disjoint subsets
[29]. More precisely, we sweep a $100{\times}100{\times}100$ voxel frame along
the coordinate axes of the original $1000$-voxel cube and displace it by steps
of $50$ voxels (half a subsample size) resulting in a total of $19^{3}=$6859$$
subsamples. Although data are sampled redundantly, the set of obtained
subsamples can be regarded independently when training neural networks [4].
Furthermore, by rotating the original $100{\times}100{\times}100$ voxel
subsamples by $90^{\circ}$ around the $y$ and $z$ axis, the number of
extractable data is increased further by a factor of three. As such,
artificial data augmentation techniques like shrinkage or inflation of the
pore space from the µCT image are not required here for obtaining a sufficient
and representative amount of training data. Even though we can produce
$3{\cdot}19^{3}{=}$20\,577$$ subsamples by the method above, we select only
the first $10\,000$ ones, since this number is sufficient to train our PhyCNN.
In particular, this number of available data samples exceeds that of similar
studies, cf. [4, 13].

Second, fluid voxels belonging to disconnected pore space with respect to the
$x$ direction possibly occurring within the subsamples are turned into solid
voxels. As disconnected pores do not contributed to Stokes flow being driven
from the inhomogeneous boundary conditions eq. 1c placed on opposing sides of
the subsample, this procedure maintains permeability properties while
facilitating their calculation. To this end, a simple graph walking algorithm
is exploited to identify connected subdomains of the pore space. Starting from
a random fluid node, neighboring nodes within the fluid domain are
successively added until the scheme converges. A thorough description of this
algorithm is found in [30].

By projecting the voxels of the connected pore space onto the $x$ axis, we
conclude whether each subsample has a nontrivial permeability. Subsamples with
zero permeability are excluded from the later workflow (the maximum flow value
is exactly zero if and only if the permeability is zero—therefore no training
on such data samples is necessary, cf. Section 3.2.2). By encoding the
simplified pore geometries in 1-bit raw format, memory consumption is
$125\text{\,}\mathrm{k}\mathrm{B}$ per subsample. Accordingly, our whole
library of data samples allocates $1.25\text{\,}\mathrm{G}\mathrm{B}$ of disk
space and is therefore still manageable on standard personal computers.

We note that a subsample of size $225\text{\,}\mathrm{\SIUnitSymbolMicro m}$
is too small to be a representative elementary volume (REV) for most sandstone
types, cf. [27, 28]. As such, computed effective properties of one subsample
cannot be expected to represent the whole segmented porous medium out of which
it was extracted. On the other hand, the obtained training data set for our
PhyCNN is therefore expected to be highly diverse, i. e. to contain highly
permeable as well as narrow and confined pore geometries. Consequently, this
setup is well suited to underline the robustness and prediction quality of our
proposed methodology. In the conclusions Section 5, we will point out how the
results presented in this paper could be used to approach REV-scale samples.

###  2.2 Forward simulation

To apply a supervised learning approach as outlined in Section 3, each
subsample within the data set derived by the methods of Section 2.1 needs to
be labeled with a computed permeability $k_{\mathrm{cmp}}$ that we use as the
reference value.

To this end, in Section 2.2.1, we perform flow simulations on the pore space.
More precisely, for each of the $10\,000$ subsamples, a stationary flow field
along the $x$ direction is computed by solving the Stokes equations on the
union of fluid voxels $\Omega$ for the fluid velocity ${\boldsymbol{u}}$ and
pressure $p$. The discretization uses arbitrary order (stable) Taylor–Hood or
reduced-order stabilized Taylor–Hood mixed finite elements.

In Section 2.2.2, the (absolute, scalar) permeability is computed by averaging
the pressure gradient and velocity field across the subsample in $x$
direction, cf. Section 2.1. Apparently, as the data set contains $y$ and
$z$-rotated versions of each subsample, this relates to the determination of
the permeability with respect to all three main axes. By accurate bookkeeping,
the diagonal permeability tensor of each initial subsample is retrievable.

Finally, in Section 2.2.3, we justify our choice of discretization parameters,
i. e., mesh refinement level and finite element spaces used to produce the
permeability values $k_{\mathrm{cmp}}$ to train and validate our PhyCNN.

|  Subsample 0 |  Subsample 9213
---|---|---
pressure |  |
|  $p_{h}\in[$-1.5100$,$1.4100$]$ |  $p_{h}\in[$-1.5100$,$1.8100$]$
velocity magnitude |  |
|  $|{\boldsymbol{u}}_{h}|\in[0,$7.510-3$]$ |  $|{\boldsymbol{u}}_{h}|\in[0,$1.910-3$]$
Figure 3: Examples of geometries used for the training process, with pressure
fields  _(top)_ and velocity magnitudes  _(bottom)_. Subsample 0  _(left)_
exhibits moderate pressure gradients due to the wide and highly conductive
channels in the pore space. Contrarily, in Subsample 9213  _(right)_ Stokes
pressure is dominated by thin structures leading to an ill-conditioned
problem, i. e., small changes in the diameter of narrow pores have drastic
impact on the permeability. Thin pore throats are indicated by black arrows in
the pressure plot.

####  2.2.1 Computation of the flow field

We consider the stationary Stokes equation for a Newtonian fluid in the
nondimensionalized form,

$\displaystyle-\frac{1}{\mathrm{Re}}\upDelta{\boldsymbol{u}}+{\boldsymbol{\nabla}}p$
$\displaystyle={\boldsymbol{0}}\quad\text{in}\leavevmode\nobreak\ \Omega,$
(1a) $\displaystyle{\boldsymbol{\nabla}}\cdot{\boldsymbol{u}}$
$\displaystyle=0\quad\text{in}\leavevmode\nobreak\ \Omega,$ (1b) where
$\Omega\subset(0,1)^{3}$ is a domain that consists of the union of fluid
voxels of a considered $100{\times}100{\times}100$ voxel subsample (i. e. the
subsample is inscribed into the unit cube). In (1),
${\boldsymbol{u}}={\boldsymbol{u}}(x,y,z)$ denotes the (dimensionless) fluid
velocity, $p=p(x,y,z)$ the (dimensionless) pressure and
$\mathrm{Re}\coloneqq\rho\,U_{\mathrm{c}}L_{\mathrm{c}}/\mu$ the Reynolds
number of the system with characteristic length $L_{\mathrm{c}}$
$\mathrm{[}\mathrm{m}\mathrm{]}$, characteristic velocity $U_{\mathrm{c}}$
$\mathrm{[}\mathrm{m}\,\mathrm{s}^{-1}\mathrm{]}$, fluid density $\rho$
$\mathrm{[}\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}\mathrm{]}$, and fluid
viscosity $\mu$ $\mathrm{[}\mathrm{Pa}\,\mathrm{s}\mathrm{]}$. Note that the
permeability $k_{\mathrm{cmp}}\leavevmode\nobreak\
$\mathrm{[}\mathrm{m}^{2}\mathrm{]}$$ is invariant with respect to
$\mathrm{Re}$ (see below). The data set is constructed in such a way that
there are connected fluid voxels (i. e. sharing a common face) reaching from
the $x{=}0$ plane to the $x{=}1$ plane, cf. Figure 3, since impermeable
subsamples were excluded from the workflow in the ‘data acquisition’ step, cf.
Figure 1 and Section 2.1. The flow field is driven by a pressure gradient in
$x$ direction induced by the boundary condition
$\left(\frac{1}{\mathrm{Re}}{\boldsymbol{\nabla}}{\boldsymbol{u}}-p\,\mathbfsf{I}\right){\boldsymbol{n}}={\boldsymbol{e}}_{x}\quad\text{on}\leavevmode\nobreak\
\Gamma_{\mathrm{N}},$ (1c) where
$\Gamma_{\mathrm{N}}\coloneqq\big{\\{}(x,y,z)\in\partial\Omega\,|\,x\in\\{0,1\\}\big{\\}}$
and ${\boldsymbol{e}}_{x}$ denoting the unit vector in $x$ direction. On the
remaining boundary
$\Gamma_{\mathrm{D}}\coloneqq\partial\Omega\setminus\Gamma_{\mathrm{N}}$, no-
slip boundary conditions are prescribed,
${\boldsymbol{u}}={\boldsymbol{0}}\quad\text{on}\leavevmode\nobreak\
\Gamma_{\mathrm{D}}.$ (1d)

The weak formulation of (1) is discretized by generalized Taylor–Hood pairs of
spaces, $\mathds{Q}_{\ell+1}^{3}/\mathds{Q}_{\ell}$, where $\mathds{Q}_{\ell}$
denotes the local space of polynomials of degree at most $\ell$ in each
variable $x,y,z$, cf. [31, 32]. For $\ell=0$, the pressure space
$\mathds{Q}_{0}=\mathds{P}_{0}$ is discontinuous and consists of elementwise
constants. The respective “reduced Taylor–Hood pair”
$\mathds{Q}_{1}^{3}/\mathds{P}_{0}$ requires stabilization (see below) due to
a lack of discrete inf-sup-stability [31]. Let
${\boldsymbol{\phi}}_{i}={\boldsymbol{\phi}}_{i}(x,y,z):\Omega\rightarrow\mathds{R}^{3}$,
$i=1,\ldots,n$ and $\psi_{i}=\psi_{i}(x,y,z):\Omega\rightarrow\mathds{R}$,
$i=1,\ldots,m$ denote the basis functions for the global discrete spaces for
velocity and pressure, respectively. The discrete velocity
${\boldsymbol{u}}_{h}={\boldsymbol{u}}_{h}(x,y,z)$ and discrete pressure
$p_{h}=p_{h}(x,y,z)$ then have the representation

${\boldsymbol{u}}_{h}=\sum_{i=1}^{n}[{\boldsymbol{x}}^{\boldsymbol{u}}]_{i}\,{\boldsymbol{\phi}}_{i},\qquad
p_{h}=\sum_{i=1}^{m}[{\boldsymbol{x}}^{p}]_{i}\,\psi_{i}$ (2)

with degree-of-freedom vectors
${\boldsymbol{x}}^{\boldsymbol{u}}\in\mathds{R}^{n}$,
${\boldsymbol{x}}^{p}\in\mathds{R}^{m}$, which are unique solutions of the
linear system

$\displaystyle\begin{bmatrix}\mathbfsf{A}&\phantom{-}\mathbfsf{B}^{\mathrm{T}}\\\
\mathbfsf{B}&-\mathbfsf{C}\end{bmatrix}\begin{bmatrix}{\boldsymbol{x}}^{\boldsymbol{u}}\\\
{\boldsymbol{x}}^{p}\end{bmatrix}=\begin{bmatrix}{\boldsymbol{b}}^{\boldsymbol{u}}\\\
{\boldsymbol{0}}\end{bmatrix}\qquad\Longleftrightarrow:\qquad\mathcal{A}\,{\boldsymbol{x}}={\boldsymbol{b}}$
(3)

with right-hand side ${\boldsymbol{b}}^{\boldsymbol{u}}\in\mathds{R}^{n}$,
$[{\boldsymbol{b}}^{\boldsymbol{u}}]_{i}\coloneqq\int_{\Gamma_{\mathrm{N}}}{\boldsymbol{e}}_{x}\cdot{\boldsymbol{\phi}}_{i}$.
The sparse blocks in $\mathcal{A}$ are the vector-Laplacian matrix
$\mathbfsf{A}\in\mathds{R}^{n,n}$ and the divergence matrix
$\mathbfsf{B}\in\mathds{R}^{m,n}$,

$[\mathbfsf{A}]_{i,j}\coloneqq\int_{\Omega}{\boldsymbol{\nabla}}{\boldsymbol{\phi}}_{i}:{\boldsymbol{\nabla}}{\boldsymbol{\phi}}_{j},\qquad[\mathbfsf{B}]_{k,j}\coloneqq-\int_{\Omega}\psi_{k}{\boldsymbol{\nabla}}\cdot{\boldsymbol{\phi}}_{j},$

and $\mathbfsf{C}$ is a stabilization matrix that is required only for lowest
order $\ell=0$ to guarantee the full rank of $\mathcal{A}$. We choose
$\mathbfsf{C}$ as in [31], (3.84), in which case, $\mathbfsf{C}$ can be
interpreted as a pressure-Laplacian discretized by cell-centered finite
volumes [33]. For $\ell>0$, $\mathbfsf{C}$ is set to zero. In either case,
$\mathcal{A}$ is symmetric and indefinite with $n$ positive and $m$ negative
eigenvalues.

In order to solve the saddle point system (3) efficiently, a preconditioned
MINRES method is applied, as it is the best choice of Krylov subspace methods
for symmetric indefinite systems [34]. The chosen precondition operator for
$\mathcal{A}$ in (3) is the symmetric and positively definite block-diagonal
matrix

$\mathcal{P}\coloneqq\mathrm{diag}\,(\mathbfsf{A},\mathbfsf{W})$

with $\mathbfsf{W}$ being the pressure-mass matrix
$[\mathbfsf{W}_{k,l}]\coloneqq\int_{\Omega}\psi_{k}\,\psi_{l}$. The action of
the inverse $\mathcal{P}^{-1}$ in each Krylov iteration is approximated block-
wise by one V-cycle of an algebraic multigrid method (Boomer AMG from the
Hypre library [35]). Since $\mathbfsf{W}$ is spectrally equivalent to the
(negative) Schur complement
$\mathbfsf{B}\mathbfsf{A}^{-1}\mathbfsf{B}^{\mathrm{T}}$ of $\mathcal{A}$ (for
$\ell>0$), and due to the utilization of AMG, the number of Krylov iterations
required to reach a given relative tolerance is bounded independently of the
mesh size [36, 37, 38] (however, it highly depends on the geometry of the
domain). MINRES with preconditioning as described above belongs to the state-
of-the art Stokes solvers in the high-performance computing context [39]. In
Section 2.2.3, we discuss appropriate choices regarding mesh size and
discretization order.

Figure 3 illustrates Stokes velocity and pressure fields for two exemplary
subsamples exhibiting qualitatively highly different pore-space geometries
such as wide pore throats and narrow channels. We will emphasize the impact of
highly and merely permeable rock samples on the behavior of DNS-based
permeability computations and PhyCNN predictions throughout the paper.

####  2.2.2 Permeability estimation

We deduce the permeability value $k_{\mathrm{cmp}}$ of interest from the
previously calculated Stokes velocity ${\boldsymbol{u}}$ and pressure $p$ (we
suppress the discretization index $h$ in this section). In the Stokes
equations (1), the inflow and outflow boundaries are defined as

$\Gamma_{\mathrm{in}}\coloneqq\big{\\{}(x,y,z)\in\partial\Omega\,|\,x=0\big{\\}}\,,\qquad\Gamma_{\mathrm{out}}\coloneqq\
\big{\\{}(x,y,z)\in\partial\Omega\,|\,x=1\big{\\}}$

and thus disjoint subsets of $\Gamma_{\mathrm{N}}$. From the solution
$({\boldsymbol{u}},p)$ of (1), we compute the approximated permeability
$k_{\mathrm{cmp}}$ $\mathrm{[}\mathrm{m}^{2}\mathrm{]}$ by the classical Darcy
law, which reads in nondimensionalized form,

$Q=-\mathrm{Da}\,\mathrm{Re}\,(P_{\mathrm{out}}-P_{\mathrm{in}}),$ (4a) with
(dimensionless) volume flow rate $Q$, and (dimensionless) area-averaged inflow
and outflow pressures, $P_{\mathrm{in}}$, $P_{\mathrm{out}}$, given by [40]
$Q\coloneqq\int\limits_{\Gamma_{\mathrm{out}}}{\boldsymbol{u}}\cdot{\boldsymbol{n}}\,,\qquad
P_{\mathrm{in}}\coloneqq\frac{1}{|\Gamma_{\mathrm{in}}|}\int\limits_{\Gamma_{\mathrm{in}}}p\,,\qquad
P_{\mathrm{out}}\coloneqq\frac{1}{|\Gamma_{\mathrm{out}}|}\int\limits_{\Gamma_{\mathrm{out}}}p\,.$
From (4a), $k_{\mathrm{cmp}}$ is derived from the Darcy number
$\mathrm{Da}\coloneqq\frac{k_{\mathrm{cmp}}}{L_{\mathrm{c}}^{2}}\,,$ (4b)

where $L_{\mathrm{c}}$ is the characteristic length, in our case, the edge
length of a $100^{3}$-voxel subsample, i. e.,
$L_{\mathrm{c}}=$225\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. We choose
$\mathrm{Re}$ equal to one, since $\mathrm{Re}$ does not influence the
permeability value $k_{\mathrm{cmp}}$ (a rescaling of ${\boldsymbol{u}}$ by
$\mathrm{Re}^{-1}$ implies a rescaling of $Q$ by $\mathrm{Re}^{-1}$ and
therefore, $\mathrm{Re}$ cancels out in (4a)).

Formula (4a) determines an approximation of the _scalar_ permeability in $x$
direction by area averaging. If this approach was applied to all three
principal directions, it yielded a diagonal permeability tensor. In [41], a
volume-averaged approach is proposed that is capable of determining the full
permeability tensor. Application to $x$ direction only, yields a column
vector, whose entries are in general non-trivial. Since we want to train our
PhyCNN with one scalar permeability value only, we decided to utilize the
area-averaging approach in this study.

####  2.2.3 Choice of discretization parameters

In this section, we investigate the influence of mesh refinement and choice of
polynomial order on the solution quality by comparing the computed
permeability $k_{\mathrm{cmp}}$ obtained from (4a) to the analytical value
$k_{\mathrm{ana}}$. In order to have analytical results available, we restrict
our considerations to viscous flow through rectangular channels
$\Omega_{a,b}$. This is supposed to constitute a sufficient benchmark, since
the approximation quality of the computed permeability $k_{\mathrm{cmp}}$ is
dominated by the discretization error in narrow pore throats due to high local
gradients.

For the Stokes equations (1), consider a rectangular channel

$\Omega_{a,b}\coloneqq(0,1)\times\left(\frac{1}{2}-\frac{a}{2},\frac{1}{2}+\frac{a}{2}\right)\times\left(\frac{1}{2}-\frac{b}{2},\frac{1}{2}+\frac{b}{2}\right)\subset\mathds{R}^{3}$

of width $a\in\left(0,\frac{1}{2}\right)$ and height
$b\in\left(0,\frac{1}{2}\right)$. An analytical expression for its
permeability $k_{\mathrm{ana}}$ is [42]

$\displaystyle k_{\mathrm{ana}}$
$\displaystyle\coloneqq\frac{K_{\infty}\cdot\min(a,b)^{3}\cdot\max(a,b)}{12},\qquad\text{where}$
$\displaystyle K_{j}$ $\displaystyle\coloneqq
1-\sum_{n=1}^{j}\frac{1}{(2n-1)^{5}}\cdot\frac{192}{\pi^{5}}\cdot\frac{\min(a,b)}{\max(a,b)}\tanh\left((2n-1)\frac{\pi}{2}\frac{\max(a,b)}{\min(a,b)}\right),$

denoting the limit of $K_{j}$ for $j\to\infty$ by $K_{\infty}$. By application
of the triangular inequality, $0<\min(a,b)\leq\max(a,b)$, and
$|\tanh(x)|<1\;\forall x\in\mathds{R}$, we obtain the following approximation
error bound:

$\displaystyle|K_{\infty}-K_{j}|\leq\sum\limits_{n=j+1}^{\infty}\frac{1}{(2n-1)^{5}}\cdot\frac{192}{\pi^{5}}.$

Piecewise application of Jensen’s inequality to the convex function
$(2x-1)^{-5}$ finally yields the estimate:

$\displaystyle|K_{\infty}-K_{j}|\leq\frac{192}{\pi^{5}}\int\limits_{j+0.5}^{\infty}\frac{1}{(2x-1)^{5}}\;\mathrm{d}x=\frac{192}{128}\cdot\frac{1}{\pi^{5}j^{4}}\approx
0.0049\,j^{-4}.$

As such, we obtain at least six significant digits using $K_{10}$ for the
calculation of the analytical reference permeability $k_{\mathrm{ana}}$.

order $\ell$ | ref. level | $k_{\mathrm{cmp}}$ | $\dfrac{|k_{\mathrm{cmp}}-k_{\mathrm{ana}}|}{k_{\mathrm{ana}}}$ |  $m$ (DOF ${\boldsymbol{u}}_{h}$) |  $n$ (DOF $p_{h}$)
---|---|---|---|---|---
0 | 0 | $8.4410-8$ | $8.8210-2$ | $8\,484$ | $1\,800$
0 | 1 | $9.0610-8$ | $2.2210-2$ | $54\,873$ | $14\,400$
0 | 2 | $9.2210-8$ | $4.6210-3$ | $390\,975$ | $115\,200$
1 | 0 | $9.2510-8$ | $1.2010-3$ | $54\,873$ | $2\,828$
1 | 1 | $9.2610-8$ | $4.0610-4$ | $390\,975$ | $18\,291$
Table 2: Computed permeabilities $k_{\mathrm{cmp}}$ of a channel with
$6{\times}3$ voxel rectangular cross section for different polynomial orders
$\ell$ of the pressure space (cf. Section 2.2.1) and mesh refinement levels.
As in (2), $m$ and $n$ are the numbers of DOF for velocity
${\boldsymbol{u}}_{h}$ and pressure $p_{h}$, respectively. The relative
tolerance of MINRES is set to $1.010-6$, which has shown to be sufficient for
three significant digits in the permeability value using a small scouting test
set.

Table 2 lists the relative error for the computed permeability
$k_{\mathrm{cmp}}$ obtained from (4a) on $\Omega_{0.06,0.03}$ using different
polynomial orders $\ell$ and global mesh refinement levels. As expected, both
refinement and higher-order discrete spaces consistently reduced the error
with respect to the analytical solution. The best cost-to-approximation-
quality ratio is achieved by $\mathds{Q}_{2}^{3}/\mathds{Q}_{1}$ elements
($\ell=1$) on the original grid. In particular, this choice poses a
significant improvement in accuracy over the use of lowest-order spaces
$\mathds{Q}_{1}^{3}/\mathds{Q}_{0}$ ($\ell=0$) with two global mesh refinement
levels.

The advantage of higher order discretizations presented above for narrow
rectangular channels seems to generalize to our actual geological samples,
which include confined structures. For Subsample 0, cf. Figure 3, where the
main flow channels are wide with respect to the voxel resolution,
$k_{\mathrm{cmp}}$ is hardly affected by increasing the polynomial order $p$
from zero to one ($4\%$ deviation). On the other hand, Subsample 9213 (Figure
3) exhibiting narrow and branchy structures experienced a significant relative
increase in computed permeability of $43\%$. Therefore, we approve the
increased computational complexity of lowest order classical Taylor–Hood
elements ($\ell=1$) for permeability labeling in our forward simulations, cf.
Section 2.2.1.

This completes the methodology description for data acquisition and yields a
training set suitable to train our neural network for permeability prediction.
We conclusively note that the solver converged for every subsample in the
database particularly maintaining the generality of our training set.

###  2.3 Evaluation metrics

In our study, we use the following well-known metrics to characterize and
compare different measures of a data set statistically.

First, we define the standard deviation $\sigma$ by

$\displaystyle\sigma=\sqrt{\frac{1}{N-1}\sum\limits_{i=1}^{N}\left(t_{i}-\bar{t}\right)^{2}}$

for a number $N$ of real-valued data $t_{i}$ and arithmetic mean value
$\bar{t}$, measuring the expected deviation from $\bar{t}$. In order to
quantify the correlation between different measures of a data set, we further
introduce the coefficient of determination $R^{2}$

$\displaystyle
R^{2}=1-\frac{\sum\limits_{i=1}^{N}(t_{i}-y_{i})^{2}}{\sum\limits_{i=1}^{N}(t_{i}-\bar{t})^{2}}$

for targets $t_{i}$ and corresponding predictions $y_{i}$, encoding the share
of variance in the targets that is covered by the predictions.

Finally, we define the mean-squared error MSE by

$\displaystyle\text{MSE}=\sum\limits_{i=1}^{N}\frac{(t_{i}-y_{i})^{2}}{N}.$

##  3 Machine learning model

In this section, we motivate and present our suggested neural network
architecture for permeability prediction. Therefore, we provide a brief
introduction to the class of convolutional neural networks in Section 3.1 as
the basis for our specific network architecture. For further reading, we refer
to [43]. In Section 3.2, this model is augmented by an additional physical
input quantity, which is derived in Section 3.2.2 in detail. For simplicity,
we refer to the $100^{3}$-voxel geological subsamples derived in Section 2.1
as (data) samples in the context of supervised machine learning.

###  3.1 Convolutional neural networks

In the following, we discuss the building blocks required to set up an
artificial neural network tailored to the task of permeability predictions
from 3D binary images of pore-space geometries. In a broader sense, we start
with the more general concept of feed-forward neural networks.

As indicated by the name, feed-forward neural networks (FFNs) process their
input data by successively propagating the input information through a finite
number of layers. Each of these layers $i\in\\{0,\dots,L\\}$ consists of
$N_{i}$ neurons, where the indices $0$ and $L$ correspond to the network’s
input and output layer, respectively. Typically, the action of a layer $1\leq
i\leq L$ on input data $x_{i-1}\in\mathds{R}^{N_{i-1}}$ is given as an affine-
linear transformation, followed by application of a nonlinearity (activation
function) $\sigma_{i}$:

$\displaystyle x_{i}=\sigma_{i}\left(W_{i}x_{i-1}+b_{i}\right).$

We refer to $W_{i}\in\mathds{R}^{N_{i},N_{i-1}}$ as the weight matrices and to
$b_{i}\in\mathds{R}^{N_{i}}$ as the biases. Their combined entries are the
learnable parameters of the network.

Convolutional neural networks (CNNs) constitute a subclass of FFNs exhibiting
a dedicated internal structure specialized to the interpretation of image
data. Furthermore, CNNs naturally incorporate translational invariance with
respect to the absolute spatial position of features within the input data. In
the following, we give a short introduction to this kind of networks. For a
broader overview of common network designs and layer types, we refer to the
respective literature, e. g. [43]. A list of the different layer types used in
our study is given in Table 3.

Typically, CNNs consist of several convolutional (conv) lower layers (filter
layers) including pooling operations, followed by dense neuron layers. As
indicated by the name, convolutional layers perform the discrete analog of a
convolution, where the convolution kernel is a learnable parameter of the
network. As such, they possess a highly structured weight matrix $W$ and are
therefore superior to fully connected layers in terms of computational
complexity, especially for highly localized kernels (sparsity). This
restriction is justified in our application scenario, since most (low-level)
geometry features are based on highly localized correlations among neighboring
voxels.

MaxPooling (MP) layers are used to reduce the resolution of the output image.
Presenting a natural bottleneck, the most significant pieces of information
are selected and passed to the subsequent neuron layers. Application of
several different filters to the same input within each layer of the network
extracts increasingly complex characteristic features of the pore geometry.
Finally, this information is interpreted by the subsequent dense neural
layers.

Convergence of the training procedure is improved by exploiting additional
batch normalization (BN) layers, calibrating the mean and spread of the
provided inputs. Throughout the network, we apply LeakyReLU nonlinearities in
an elementwise manner. For a parameter $\alpha\in(0,1)$, the
LeakyReLU($\alpha$) function is defined by

$\displaystyle\texttt{LeakyReLU}(\alpha;x)\coloneqq\max(0,x)-\alpha\max(0,-x).$

Due to the slope of this nonlinearity being bounded away from zero, we
circumvent the dying neuron problem occurring in deep neural networks with
ReLU nonlinearities [44], which arise formally by setting $\alpha=0$.
Furthermore, backpropagation during the learning phase is facilitated. Our
implementation is conducted using the Deep Learning Toolbox in Matlab R2021a
[25].

###  3.2 Physics-informed convolutional neural networks

In the following, we use the building blocks for general CNNs as introduced in
Section 3.1 and derive a network architecture specifically tailored to our
application in permeability prediction. Therefore, we first give an
introduction to PhyCNNs and thereafter discuss the choice of meaningful
additional input quantities to the network.

####  3.2.1 General architecture

As discussed above, standard CNNs are required to extract all important
features and information solely from the input image data to predict the
correct output. Due to the complex interplay of weights and biases within the
layers of a CNN, it is not possible to input a priori knowledge of important
features that one would wish the network to take into account. Nevertheless,
it is possible to provide the network with additional input derived from the
image data in a preprocessing step, leading to the field of PhyCNNs. An
exemplary schematic representation of such a network structure is provided in
Figure 5. As the additional physical input is typically not of image type (and
hence does not exhibit spatially correlated features), it omits the
convolutional layers and is directly coupled to the upper dense layers. As
such, PhyCNNs predict by considering both the extracted geometry features as
well as the provided physical input quantities. Depending on the quality and
quantity of the extracted features as well as the specific application, these
may even provide sufficient information to the network to accurately perform
predictions, cf. [2], underlining the potential of physics-informed
strategies.

This approach has been successfully applied to the prediction of
permeabilities from pore-space geometries: In [3], porosity and surface area
have been used to increase predictions performance for 2D image data. The
specific choice of these parameters is related to the Kozeny–Carman formula
estimating permeability from porosity and surface area. As indicated in [45],
the proposed relation deteriorates in quality for increasingly complex
geometries. However, [21] found porosity to be the most prominent input factor
for their CNN among a set of various derived quantities such as coordinate
number or mean pore-size. Furthermore, characteristics obtained by pore
network modeling have shown low general accuracy [21].

As indicated by our data set statistics in Figure 4, porosity and surface area
are only weakly correlated with the computed permeability $k_{\mathrm{cmp}}$
for our complex 3D setting. For a given surface area $A_{\mathrm{cmp}}$, we
observe data samples covering four orders of magnitude in their respective
permeability values, about half the range is observed for fixed porosities
$\phi_{\mathrm{cmp}}$. Hence, drastic increases in accuracy cannot be expected
by regarding these quantities as additional physical input in our case.
Instead, we construct another more specific physical quantity that can be
computed efficiently from the input data using graph algorithms included in
Matlab [25]. More precisely, we solve a surrogate maximum flow problem on a
graph representation of the pore space. In Section 3.2.2, we give a thorough
description of this chosen quantity.

Figure 4: Correlation of different physical measures of the underlying
geometry with the computed permeability values $k_{\mathrm{cmp}}$ for
Bentheimer data samples. From left to right: Computed surface area
$A_{\mathrm{cmp}}$, computed porosity $\phi_{\mathrm{cmp}}$, specific surface
area $A_{\mathrm{spec}}$, maximum flow on graph representation
${f_{\mathrm{max}}}$.

Our final network architecture is presented in Table 3 and Figure 5. The
structure is based on the findings of [4], where the hyper-parameters of a
purely convolutional neural network have been optimized using a grid search
algorithm. For our study, we decrease the bottleneck constituted by the final
convolutional layer and improve convergence by adding batch normalization
layers and relaxed cut-offs as nonlinearities (LeakyReLU). Furthermore, the
maximum flow value ${f_{\mathrm{max}}}$ is treated as an additional scalar
physical input quantity which is duplicated to a vector of length $64$ and
subsequently concatenated with the first dense layer of the network. By doing
so, we balance both inputs to the second dense layer which results in a more
uniform weight distribution.

Figure 5: Schematic representation of physics-informed convolutional neural networks (PhyCNNs). The scalar valued physical input circumvents the image convolution layers and is directly fed into the subsequent dense layers. Graphics produced using NN-SVG [46]. block | layers | learnables
---|---|---
input1 | image input 100$\times$100$\times$100 | —
input2 | physics input 1$\times$1$\times$1 | —
conv1 |  conv(32,5) – BN – LeakyReLU(0.1) – MP(5,5) | 4 096
conv2 |  conv(64,5) – BN – LeakyReLU(0.1) – MP(5,5) | 256 192
conv3 |  conv(128,3) – BN – LeakyReLU(0.1) – MP(2,2) | 221 568
dense1 |  dense(64) – LeakyReLU(0.1) | 65 728
dense2 |  dense(32) – LeakyReLU(0.1) | 4 161
output |  regression(1) | —
Table 3: Layer structure of our PhyCNN. The physical scalar input2 is expanded
to dimension $64{\times}1$ and concatenated with the output of ‘dense1’ using
a depth concatenation layer. Nomenclature:
conv($N$,$K$): | convolutional layer with $N$ channels and $K\times K$ kernel size;
---|---
BN: | batch normalization layer;
MP($N$,$P$): | maxPooling Layer, size $N$ stride $P$;
dense($N$): | dense layer with $N$ neurons;
regression($N$): | regression layer with $N$ neurons;
LeakyReLU($\alpha$): | leaky rectified linear unit, slope $\alpha$ on negative inputs.

####  3.2.2 Maximum flow problems on graphs

In the field of optimization, maximum flow problems aim at identifying the
maximum flow rate a network of pipes is able to sustain. Let the network be
described by an undirected graph $G(N,E,\omega)$, where the nodes $N$
correspond to distribution nodes, edges $E$ to connections by pipes, and the
weights $\omega$ encode the maximal capacity of each pipe. For any two nodes
$n_{1}\neq n_{2}\in N$ (source and sink), we can determine the maximum flow
between those nodes allowed by the network. For a thorough introduction to
graph algorithms and maximum flow problems, we refer to [23]. In this sense,
we aim at approximating relevant properties of the physical flow between two
opposite sides of a sample by solving flow problems on a suitable graph
network.

To do so, we consider the pore space of our segmented image data (stemming
from geological specimens) to be a graph by identifying each voxel as a node
and each neighboring relation of voxels via a common face as an edge. This
procedure is illustrated exemplarily in Figure 6 for a 2D structure. For
simplicity, all edge weights $\omega_{i}$ are assumed to be one. Note that
this graph representation still comprises multiple standard characteristics of
the pore space. For example, porosity translates to the number of nodes in the
graph divided by the total number of voxels; the number of surface elements is
approximately determined by $6\operatorname{card}(N)-2\operatorname{card}(E)$.
Moreover, we add a further node $n_{\text{in}}$ connected to each voxel of the
sample’s inflow face, as well as $n_{\text{out}}$ being connected to each
voxel of the sample’s outflow face, cf. Figure 6. As such, we can approximate
the permeability determination problem by a maximum flow problem through $G$
between the nodes $n_{\text{in}},n_{\text{out}}$.

For the samples considered in this paper, the computational effort for
calculating ${f_{\mathrm{max}}}$ using Matlab’s built-in routine maxflow is
reasonably small compared to the solution of stationary Stokes equations as
described in Section 4.3.

Figure 6: Left: Derivation of the graph representation from a pore space
illustrated by a 2D $4{\times}4$ pixel image. Nodes are depicted as circles,
edges as lines. Gray voxels refer to the solid matrix. The maximum flow value
${f_{\mathrm{max}}}$ of the depicted graph is one. Right: More complex
$8{\times}8$ pixel setup exhibiting ${f_{\mathrm{max}}}$ of three w.r.t. the
horizontal axis and zero w.r.t. the vertical axis.

By the min-cut max-flow theorem, cf. [23], the maximum flow problem allows for
another interesting interpretation. In the setup introduced within this
section, ${f_{\mathrm{max}}}$ corresponds to the minimal number of edges that
need to be deleted from $G$ to obtain a disconnected graph such that inflow
and outflow faces of the sample are contained in different connected
components, cf. Figure 6. Hence, we obtain information about the connectivity
of the pore space with respect to the direction of interest. More precisely,
${f_{\mathrm{max}}}$ classifies structures containing thin channels regarding
their restricting effects on fluid flow. This behavior is comprehensible using
the geometries illustrated in Figure 3: Exhibiting ${f_{\mathrm{max}}}$ of
$1230$, the permeability of Subsample 0 is governed by wide channels. On the
other hand, exhibiting ${f_{\mathrm{max}}}$ of $55$, narrow pore throats
hamper the fluid flow of Subsample 9213.

We further note that for moderate tortuosity, ${f_{\mathrm{max}}}$ is
approximately proportional to the volume of the channel needed to transport
the maximum flow through the sample. As such, by exclusively regarding the
subgraph of $G$ actively used by the maximum flow, we cut off parts of $G$
that do not effectively contribute to fluid flow. Hence, the determination of
${f_{\mathrm{max}}}$ can be understood as approximating the porosity with
respect to the pore space participating in fluid transport. That way, we
naturally improve the permeability estimation based on classical porosities.

Plotting ${f_{\mathrm{max}}}$ against the calculated permeabilities
$k_{\mathrm{cmp}}$, we conclude an almost linear relationship between both
quantities in the logarithmic scale. Using mean-square linear regression, we
obtain the approximate functional relation

$\displaystyle k_{\mathrm{cmp}}\sim 50625\cdot(f_{\mathrm{max}})^{1.407}\cdot
10^{-8.183}\quad[$\mathrm{D}$],$ (5)

cf. Figure 4. On the displayed logarithmic scale, this regression yields an
approximation quality of $R^{2}=0.8692$, cf. Section 2.3. Apparently, the
derived quantity is strongly correlated to the target permeability
$k_{\mathrm{cmp}}$. Therefore, the permeability estimation using formula eq. 5
is considered a suitable initial guess, which the CNN is able to improve by
relating to features extracted from the pore geometry.

Figure 7: Computed porosity $\phi_{\mathrm{cmp}}$ and permeability
$k_{\mathrm{cmp}}$ distribution among the total of $10\,000$ Bentheimer
samples. Mean porosity is $26.20\%$ with standard deviation $\sigma=0.067$,
mean permeability is $3930\text{\,}\mathrm{m}\mathrm{D}$ with
$\sigma=$8319\text{\,}\mathrm{m}\mathrm{D}$$.

##  4 PhyCNN training and workflow performance

In this section, information concerning our PhyCNN’s training procedure is
provided. Furthermore, we evaluate the quality of our resulting predictions on
three different types of sandstone as well as a challenging artificially
deformed series of geometries. Finally, we compare the computational
efficiency of our Stokes solver and the PhyCNN.

###  4.1 PhyCNN training

As the data points obtained by the forward computation are quite sparse for
extremely high and low permeability values (Figure 7), we disregarded all data
samples ranging outside the interval of [$50\text{\,}\mathrm{m}\mathrm{D}$,
$50\text{\,}\mathrm{D}$]. As a result, we improve the quality of predictions
within the permeability range where the availability of data is favorable and
reliable. Splitting the remaining data set according to a $90\%/10\%$ key,
$8\,876$ samples are used for training the network, another 987 for
validation. In order to accurately account for the labels covering three
orders of magnitude, a logarithmic transformation was applied before training
as done in [47]. Using a standard mean-squared-error (MSE) loss in the
regression as given in Section 2.3, this results in measuring the error in
relative deviation rather than absolute deviation.

The training was performed over 15 epochs using a stochastic gradient descent
(SGD) optimizer with momentum $0.9$. An almost constant validation loss over
the last epochs indicated convergence of the network. Starting from an initial
learning rate of $0.0020$, the step size has been decreased by $60\%$ every
four epochs.

###  4.2 Prediction quality

The network’s predictive performance is illustrated in Figure 8 via regression
plots. Our PhyCNN achieved coefficient of determination ($R^{2}$) values of
$96.33\%$ on training data and $93.22\%$ on validation data, cf. Table 4. As
such, the network accomplished very high accuracy with only weak tendency to
overfitting. Furthermore, the plots in Figure 8 show a low number of outliers,
underlining the stability of our approach.

Figure 8: PhyCNN accuracy on validation a) and training data b) originating from Bentheimer sandstone. The related $R^{2}$ values are $0.8820$, $0.9289$ in the natural scale and $0.9322$, $0.9633$ in the logarithmic scaling. Images c) and d) show our PhyCNN’s accuracy on the Berea and Castlegate test sets, respectively. | Bentheimer (train.) | Bentheimer (val.) | Berea | Castlegate
---|---|---|---|---
$R^{2}$ | 92.89% | 88.20% | 89.43% | 94.55%
$R^{2}$ log | 96.33% | 93.22% | 94.13% | 94.51%
mean $k_{\mathrm{cmp}}$ [mD] | 3417.7 | 3628.8 | 660.9 | 1661.0
$\sigma(k_{\mathrm{cmp}})$ [mD] | 4948.7 | 5233.2 | 657.8 | 2442.2
Table 4: Overview prediction quality. $R^{2}$ values in natural and
logarithmic scale for the training and validation Bentheimer data as well as
for Berea and Castlegate sandstone. Furthermore, mean permeability and
standard deviation $\sigma(k_{\mathrm{cmp}})$ are listed for the respective
data sets.

In order to further stress the robustness, we additionally validate our PhyCNN
on samples originating from different types of sandstone that were not used
for training. Using the data provided by [26, 27], 200 additional subsamples
were extracted from each Berea and Castlegate µCT scan, see Table 1, using the
same methods as described in Section 2. Achieving $R^{2}$ values of 89.43% and
94.13% (logarithmic) on Berea as well as 94.55% and 94.51% (logarithmic) on
Castlegate, the networks proves excellent generalization properties across
different types of pore geometries, cf. Figure 8, Figure 2, and Table 1. As
such, the network seems in fact to perform slightly better on Berea and
Castlegate data samples than on the validation data set of Bentheimer.
However, Table 4 marks a very high standard deviation for the permeability in
Bentheimer compared to the other sandstone types. As such, the latter data
sets comprise less heterogeneity facilitating the network’s prediction. We
further emphasize that both Berea and Castlegate validation sets have not been
screened to match the trained range [$50\text{\,}\mathrm{m}\mathrm{D}$,
$50\text{\,}\mathrm{D}$]. More precisely, eleven Berea data samples exhibit a
permeability $k_{\mathrm{cmp}}$ below $50\text{\,}\mathrm{m}\mathrm{D}$ as
well as two Castlegate data samples.

Figure 9: Testing generalization ability on artificially distorted data
samples. PhyCNN prediction performance on diluted and reduced pore spaces
exhibiting a large range of porosities $\phi_{\mathrm{cmp}}$. For each
manipulated geometry, Stokes simulation data $k_{\mathrm{cmp}}$ (red) and CNN
prediction $k_{\mathrm{prd}}$ (blue) are compared. Furthermore, we provide
relative prediction errors. Markers referring to the original, not manipulated
data sample are increased in size. The permeability range spanned by our
PhyCNN’s training data set is highlighted in gray.

Finally, we investigate the generalization limits of our PhyCNN by subjecting
it to an artificially distorted data set covering a challengingly wide
permeability and porosity range. To achieve that, we applied a level-set-based
algorithm already used in [48] to expand or shrink the pore space of Subsample
0 displayed in Figure 3. More precisely, the pore space is contracted with a
uniform level-set velocity directed perpendicularly to the pore walls until
all flow channels collapse, i. e. the data sample becomes impermeable. The
same number of deformation steps is also used to expand the pore volume of
Subsample 0. Each of those steps corresponds to a constriction/expansion of
the pore space by a single additional voxel layer on average.

Subsequently, we compute the permeabilities $k_{\mathrm{cmp}}$ using the
Stokes solver and $k_{\mathrm{prd}}$ using the PhyCNN of the series of pore
spaces and compare the results in Figure 9. The data show an almost perfect
match for expanded pore spaces even for porosities up to 70%. Using linear
interpolation, we estimate $k_{\mathrm{cmp}}$ to exceed the trained
permeability range [$50\text{\,}\mathrm{m}\mathrm{D}$,
$50\text{\,}\mathrm{D}$] beyond $\phi_{\mathrm{cmp}}=57.33\%$. Therefore, we
conclude that our PhyCNN is capable of properly characterizing also highly
permeable samples. On the other hand, predictions remain reasonably accurate
for strongly confined geometric structures. Three samples in Figure 9 show
porosities below $9.08\%$, which approximately refers to the lower end of the
trained permeability range for this example. These exhibit
${f_{\mathrm{max}}}$ of 4, 16, and 30, respectively, referring to an almost
disconnected pore space. Since the flow channels narrow in these cases to only
very few voxels in diameter, errors from the µCT scan as well the
discretization of Stokes equations (1) may become non-negligible. As such,
these samples may leave the current operating limits of digital rock physics,
resulting in a significant systematic overestimation of permeability compared
to lab experiments [49].

###  4.3 Computational performance

In this section, we provide computational performance indicators for the
forward simulations as performed using the method described in Section 2.2 as
well as the ones obtained for estimations using the PhyCNN. Subsequently, we
compare the forward simulation run times for the generation of $10\,000$ data
samples including both permeability computation approaches on the same compute
cluster to estimate the actual speed up. All subsequent specifications of
computation times refer to the wall time.

Each of the $10\,000$ forward simulations on Bentheimer sandstone was
performed with classical Taylor–Hood elements on voxels ($\ell{=}1$, cf.
Section 2.2.3) in parallel on 50 compute nodes of the Emmy compute cluster at
RRZE [50], each being equipped with two Intel®Xeon 2660v2 ‘Ivy Bridge’
processors (10 physical cores per chip, i. e., $50{\cdot}2{\cdot}10{=}1000$
cores total, no hyper-threading) and $64\text{\,}\mathrm{G}\mathrm{B}$ of RAM.
The linear solver (MINRES) converged within 969 iterations on average with a
relative tolerance of $1.010-6$ in the norm $\|\cdot\|_{\mathcal{P}^{-1}}$.
Furthermore, meshing, assembly, and solution were performed within $2.35$,
$0.44$, and $33.95$ seconds on average per solver call, respectively. As such,
the labeling procedure of all $10\,000$ data samples was completed within
roughly $100$ hours. We conclude that the forward simulation effort of our
Stokes approach is roughly comparable to the LBM-based simulations applied in
related publications, cf. [4].

Operating on a graphics cluster featuring two Intel®Xeon E5-2620v3 (6 cores
each), two Nvidia® Geforce Titan X GPUs, and 64 GB of RAM, the PhyCNN’s
training process terminated within two hours. Using 10 CPU cores and a single
graphics chip on this machine, the permeabilities of the $10\,000$ elements
data set of Bentheimer sandstone is estimated within $1294$ seconds by our
PhyCNN. More precisely, $867$ seconds on CPU were spent on the graph flow
problems related to ${f_{\mathrm{max}}}$ as described in Section 3.2.2 while
$367$ seconds were required to perform the subsequent network inferencing on
GPU. Estimating the time required to solve the related Stokes problem (1)
extrapolated from a subset of ten data samples indicates a total time of
approximately $1469$ hours. Accordingly, we conclude an acceleration factor of
$4087$ by using our PhyCNN on the stated hardware configuration and data set.

From the data presented above, we infer that the time invested in the
network’s training procedure including the calculation of ${f_{\mathrm{max}}}$
for the complete database is compensated after about $15$ data samples
computed beyond the training and validation data sets. This strongly
underlines the capability of our approach to pose a time-saving yet accurate
alternative to flow simulation-based permeability estimators on large data
sets.

##  5 Summary and conclusions

In this work, we demonstrated the feasibility of direct numerical simulation
(DNS) for flow through porous media to generate a library of computed
permeability labels from 3D images acquired using specimens of natural rocks.
Our distributed-parallel stationary Stokes solver achieved computational
efficiency comparable to classical lattice Boltzmann method (LBM) approaches
while resolving LBM convergence issues that typically appear on complex 3D
geometries.

As a result of computations with the proposed numerically robust forward
simulation algorithm, an unbiased (by computational artifacts) data set was
constructed to train a convolutional neural network for accelerated
permeability predictions. An easily-computable graph-based characteristic
quantity of the pore space, namely the maximum flow value, was introduced to
the neural network as another novel development in our work. As a result, this
methodology enabled our machine learning model to achieve an $R^{2}$ value of
$93.22\%$ on the validation data set. Moreover, similar prediction qualities
were found for types of sandstone rock that are different from the training
samples, as well as for artificially generated voxel sets. These observations
underline the robustness of our artificial intelligence augmented permeability
estimation approach.

In a one-on-one comparison before accounting for the training-investment-
related computational costs, the neural-network-based permeability estimation
approach delivered a speed-up in excess of $4000$ fold. On the other hand,
computational results indicate that the proposed approach exceeds the
performance of a purely DNS-based workflow beyond approximately 15 data
samples when the training-related computational costs are accounted for. This
indicates that the proposed artificial intelligence augmented permeability
estimation workflow is viable for real-life digital rock physics applications.

Further research is required to extend our refined methodology to geological
samples on the REV scale. Recent publications suggest hierarchical multiscale
neural networks [51] to alleviate the memory requirements of the training and
inference procedure. Merging both approaches holds the potential of leveraging
our current results to larger scale geological samples.

## Data availability

The CNN code used in this paper is available as part of the porous media
numerical toolkit RTSPHEM [52] on Github.

## Notation

### Symbols

$A_{\mathrm{cmp}}$ |  Computed interior surface area of µCT scan [$$\mathrm{mm}$^{2}$].
---|---
$A_{\mathrm{spec}}$ |  Specific (w.r.t. material volume) interior surface area [$$\mathrm{mm}$^{-1}$].
$\alpha$ |  Parameter in LeakyReLU nonlinearity.
${f_{\mathrm{max}}}$ |  Maximum flow value.
$k_{\mathrm{ana}}$ |  Analytical permeability (all permeabilities in darcies D, millidarcies mD).
$k_{\mathrm{cmp}}$ |  Computed permeability, cf. (4).
$k_{\mathrm{exp}}$ |  Experimental permeability.
$k_{\mathrm{prd}}$ |  Predicted permeability.
$\ell$ |  Polynomial degree along principle axes ($\mathds{Q}_{\ell}$ is the local space of polynomials of degree at most $\ell$ in each variable).
$\Omega$ |  Pore space, domain $\Omega\subset\mathds{R}^{3}$ of the Stokes equations (1).
$p$ |  Hydrostatic pressure on the pore scale (dimensionless), solution $p:\Omega\rightarrow\mathds{R}$ of the Stokes equations (1).
$\phi_{\mathrm{cmp}}$ |  Computed porosity.
$\phi_{\mathrm{exp}}$ |  Experimental porosity.
$R^{2}$ |  R-squared value, coefficient of determination.
$\mathrm{Re}$ |  Reynolds number (definition is irrelevant here, since $\mathrm{Re}$ cancels out in the computation for $k_{\mathrm{cmp}}$, cf. (4)).
$\sigma$ |  Standard deviation.
${\boldsymbol{u}}$ |  Solenoidal fluid velocity on the pore scale (dimensionless), solution ${\boldsymbol{u}}:\Omega\rightarrow\mathds{R}^{3}$ of the Stokes equations (1).
$n$ |  Node $n\in N$ in an undirected graph $G(N,E,\omega)$, cf. Section 3.2.2.

### Abbreviations

CNN |  Convolutional neural network.
---|---
DNS |  Direct numerical simulation.
DOF |  Degrees of freedom.
FFN |  Feed forward network.
LBM |  Lattice Boltzmann method.
µCT |  Microcomputed tomography.
MSE |  Mean squared error.
PhyCNN |  Physics-informed convolutional neural network.
ReLU |  Rectified linear unit.
REV |  Representative elementary volume.
SGD |  Stochastic gradient descent.

## Acknowledgments

S. Gärttner and A. Meier were supported by the DFG Research Training Group
2339 Interfaces, Complex Structures, and Singular Limits.
N. Ray was supported by the DFG Research Training Group 2339 Interfaces,
Complex Structures, and Singular Limits and the DFG Research Unit 2179
MadSoil.
F. Frank was supported by the Competence Network for Scientific High
Performance Computing in Bavaria (KONWIHR).
We further thank Martin Burger for insightful discussions and Fabian Woller
for assisting with the MFEM implementation.

## References

  * [1] Ying Da Wang, Martin J. Blunt, Ryan T. Armstrong and Peyman Mostaghimi  “Deep learning in pore scale imaging and modeling”  In _Earth-Science Reviews_ 215, 2021, pp. 103555  DOI: 10.1016/j.earscirev.2021.103555
  * [2] Mauricio Araya-Polo, Joseph Jennings, Amir Adler and Taylor Dahlke  “Deep-learning tomography”  In _The Leading Edge_ 37 Society of Exploration Geophysicists, 2018, pp. 58–66  DOI: 10.1190/tle37010058.1
  * [3] Jinlong Wu, Xiaolong Yin and Heng Xiao  “Seeing permeability from images: Fast prediction with convolutional neural networks”  In _Science Bulletin_ 63.18 Elsevier BV, 2018, pp. 1215–1222  DOI: 10.1016/j.scib.2018.08.006
  * [4] Jin Hong and Jie Liu  “Rapid estimation of permeability from digital rock using 3D convolutional neural network”  In _Computational Geosciences_ 24.4 Springer ScienceBusiness Media LLC, 2020, pp. 1523–1539  DOI: 10.1007/s10596-020-09941-w
  * [5] Dorthe Wildenschild and Adrian P. Sheppard  “X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems” 35th Year Anniversary Issue  In _Advances in Water Resources_ 51, 2013, pp. 217–246  DOI: 10.1016/j.advwatres.2012.07.018
  * [6] Tom Bultreys, Wesley De Boever and Veerle Cnudde  “Imaging and image-based fluid transport modeling at the pore scale in geological materials: A practical introduction to the current state-of-the-art”  In _Earth-Science Reviews_ 155, 2016, pp. 93–128  DOI: 10.1016/j.earscirev.2016.02.001
  * [7] Mauricio Araya-Polo, Faruk O. Alpak, Sander Hunter, Ronny Hofmann and Nishank Saxena  “Deep learning-driven permeability estimation from 2D images”  In _Computational Geosciences_ 24.2 Springer ScienceBusiness Media LLC, 2019, pp. 571–580  DOI: 10.1007/s10596-019-09886-9
  * [8] Yu Zhou and Yali Wu  “Analyses on influence of training data set to neural network supervised learning performance”  In _Advances in Computer Science, Intelligent System and Environment_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2011, pp. 19–25  DOI: 10.1007/978-3-642-23753-9_4
  * [9] Arndt Wagner, Elissa Eggenweiler, Felix Weinhardt, Zubin Trivedi, David Krach, Christoph Lohrmann, Kartik Jain, Nikolaos Karadimitriou, Carina Bringedal, Paul Voland, Christian Holm, Holger Class, Holger Steeb and Iryna Rybak  “Permeability estimation of regular porous structures: a benchmark for comparison of methods”  In _Transport in Porous Media_ 138.1 Springer ScienceBusiness Media LLC, 2021, pp. 1–23  DOI: 10.1007/s11242-021-01586-2
  * [10] Xiaofan Yang, Yashar Mehmani, William A. Perkins, Andrea Pasquali, Martin Schönherr, Kyungjoo Kim, Mauro Perego, Michael L. Parks, Nathaniel Trask, Matthew T. Balhoff, Marshall C. Richmond, Martin Geier, Manfred Krafczyk, Li-Shi Luo, Alexandre M. Tartakovsky and Timothy D. Scheibe  “Intercomparison of 3D pore-scale flow and solute transport simulation methods” Pore scale modeling and experiments  In _Advances in Water Resources_ 95 Elsevier BV, 2016, pp. 176–189  DOI: 10.1016/j.advwatres.2015.09.015
  * [11] Quan Liao and Tien-Chien Jen  “Computational Fluid Dynamics”  IntechOpen, 2011  DOI: 10.5772/10585
  * [12] Timm Krüger, Halim Kusumaatmaja, Alexandr Kuzmin, Orest Shardt and Goncalo Silva  “The Lattice Boltzmann Method”  Springer International Publishing, 2016  DOI: 10.1007/978-3-319-44649-3
  * [13] Javier E. Santos, Duo Xu, Honggeun Jo, Christopher J. Landry, Maša Prodanović and Michael J. Pyrcz  “PoreFlow-Net: A 3D convolutional neural network to predict fluid flow through porous media”  In _Advances in Water Resources_ 138 Elsevier BV, 2020, pp. 103539  DOI: 10.1016/j.advwatres.2020.103539
  * [14] P.. Bhatnagar, E.. Gross and M. Krook  “A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems”  In _Physical Review_ 94.3 American Physical Society (APS), 1954, pp. 511–525  DOI: 10.1103/physrev.94.511
  * [15] F.. Alpak, F. Gray, N. Saxena, J. Dietderich, R. Hofmann and S. Berg  “A distributed parallel multiple-relaxation-time lattice Boltzmann method on general-purpose graphics processing units for the rapid and scalable computation of absolute permeability from high-resolution 3D micro-CT images”  In _Computational Geosciences_ 22.3 Springer ScienceBusiness Media LLC, 2018, pp. 815–832  DOI: 10.1007/s10596-018-9727-7
  * [16] Zhaoxia Yang  “Analysis of lattice Boltzmann boundary conditions” http://nbn-resolving.de/urn:nbn:de:bsz:352-opus-35079, 2007
  * [17] Jonas Tölke  “Implementation of a lattice Boltzmann kernel using the compute unified device architecture developed by nVIDIA”  In _Computing and Visualization in Science_ 13.1 Springer ScienceBusiness Media LLC, 2008, pp. 29–39  DOI: 10.1007/s00791-008-0120-2
  * [18] Nishank Saxena, Ronny Hofmann, Faruk O. Alpak, Steffen Berg, Jesse Dietderich, Umang Agarwal, Kunj Tandon, Sander Hunter, Justin Freeman and Ove Bjorn Wilson  “References and benchmarks for pore-scale flow simulated using micro-CT images of porous media and digital rocks”  In _Advances in Water Resources_ 109 Elsevier BV, 2017, pp. 211–235  DOI: 10.1016/j.advwatres.2017.09.007
  * [19] M. Siena, J.. Hyman, M. Riva, A. Guadagnini, C.. Winter, P.. Smolarkiewicz, P. Gouze, S. Sadhukhan, F. Inzoli, G. Guédon and E. Colombo  “Direct numerical simulation of fully saturated flow in natural porous media at the pore scale: a comparison of three computational systems”  In _Computational Geosciences_ 19.2 Springer ScienceBusiness Media LLC, 2015, pp. 423–437  DOI: 10.1007/s10596-015-9486-7
  * [20] R. Anderson, J. Andrej, A. Barker, J. Bramwell, J.-S. Camier, J. V. Dobrev, Y. Dudouit, A. Fisher, Tz. Kolev, W. Pazner, M. Stowell, V. Tomov, I. Akkerman, J. Dahm, D. Medina and S. Zampini  “MFEM: A Modular Finite Element Library”  In _Computers & Mathematics with Applications_ 81, 2021, pp. 42–74  DOI: 10.1016/j.camwa.2020.06.009
  * [21] Moussa Tembely, Ali M. AlSumaiti and Waleed S. Alameri  “Machine and deep learning for estimating the permeability of complex carbonate rock from X-ray micro-computed tomography”  In _Energy Reports_ 7 Elsevier BV, 2021, pp. 1460–1472  DOI: 10.1016/j.egyr.2021.02.065
  * [22] Huafeng Sun, Hasan Al-Marzouqi and Sandra Vega  “EPCI: A new tool for predicting absolute permeability from computed tomography images”  In _Geophysics_ 84.3, 2019, pp. F97–F102  DOI: 10.1190/geo2018-0653.1
  * [23] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest and Clifford Stein  “Introduction to Algorithms”  The MIT Press, 2001
  * [24] H. Viswanathan, Jeffrey Hyman, S. Karra, Daniel O’Malley, Shriram Srinivasan, A. Hagberg and Gowri Srinivasan  “Advancing graph-based algorithms for predicting flow and transport in fractured rock”  In _Water Resources Research_ 54.9 American Geophysical Union (AGU), 2018, pp. 6085–6099  DOI: 10.1029/2017wr022368
  * [25] MATLAB  “version 9.10.0.1602886 (R2021a)”  Natick, Massachusetts: The MathWorks Inc., 2021
  * [26] Rodrigo Neumann, Mariane Andreeta and Everton Lucas-Oliveira  “11 Sandstones: raw, filtered and segmented data”  Digital Rocks Portal, http://www.digitalrocksportal.org/projects/317, 2020  DOI: 10.17612/f4h1-w124
  * [27] Rodrigo F. Neumann, Mariane Barsi-Andreeta, Everton Lucas-Oliveira, Hugo Barbalho, Willian A. Trevizan, Tito J. Bonagamba and Mathias Steiner  “High accuracy capillary network representation in digital rock reveals permeability scaling functions”  In _Scientific Reports_ 11.1 Springer ScienceBusiness Media LLC, 2021  DOI: 10.1038/s41598-021-90090-0
  * [28] Linxian Gong, Lei Nie and Yan Xu  “Geometrical and topological analysis of pore space in sandstones based on X-ray computed tomography”  In _Energies_ 13.15 MDPI AG, 2020, pp. 3774  DOI: 10.3390/en13153774
  * [29] Oleg Sudakov, Evgeny Burnaev and Dmitry Koroteev  “Driving digital rock towards machine learning: Predicting permeability with gradient boosting and deep neural networks”  In _Computers & Geosciences_ 127 Elsevier BV, 2019, pp. 91–98  DOI: 10.1016/j.cageo.2019.02.002
  * [30] John Hopcroft and Robert Tarjan  “Algorithm 447: efficient algorithms for graph manipulation”  In _Communications of the ACM_ 16.6 New York, NY, USA: Association for Computing Machinery (ACM), 1973, pp. 372–378  DOI: 10.1145/362248.362272
  * [31] H.. Elman, D.. Silvester and A.. Wathen  “Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics”, Numerical Mathematics and Scientific Computation  Oxford University Press, 2014
  * [32] A. Ern and J.. Guermond  “Theory and Practice of Finite Elements”, Applied Mathematical Sciences  Springer, 2004
  * [33] Florian Frank, Chen Liu, Faruk Omer Alpak and Béatrice Rivière  “A finite volume / discontinuous Galerkin method for the advective Cahn–Hilliard equation with degenerate mobility on porous domains stemming from micro-CT imaging”  In _Computational Geosciences_ 22.2, 2018, pp. 543–563  DOI: 10.1007/s10596-017-9709-1
  * [34] A.. Wathen  “Preconditioning”  In _Acta Numerica_ 24 Cambridge University Press (CUP), 2015, pp. 329–376  DOI: 10.1017/S0962492915000021
  * [35] Robert D. Falgout, Jim E. Jones and Ulrike Meier Yang  “The design and implementation of hypre, a library of parallel high performance preconditioners”  In _Numerical Solution of Partial Differential Equations on Parallel Computers_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2006, pp. 267–294  DOI: 10.1007/3-540-31619-1_8
  * [36] Michele Benzi, Gene H. Golub and Jörg Liesen  “Numerical solution of saddle point problems”  In _Acta Numerica_ 14 Cambridge University Press (CUP), 2005, pp. 1–137  DOI: 10.1017/s0962492904000212
  * [37] Miroslav Rozložník  “Saddle-Point Problems and Their Iterative Solution”, Neças Center Series  Springer International Publishing, 2018
  * [38] John W. Pearson, Jennifer Pestana and David J. Silvester  “Refined saddle-point preconditioners for discretized Stokes problems”  In _Numerische Mathematik_ 138.2 Springer ScienceBusiness Media LLC, 2017, pp. 331–363  DOI: 10.1007/s00211-017-0908-4
  * [39] Björn Gmeiner, Markus Huber, Lorenz John, Ulrich Rüde and Barbara Wohlmuth  “A quantitative performance study for Stokes solvers at the extreme scale”  In _Journal of Computational Science_ 17 Elsevier BV, 2016, pp. 509–521  DOI: 10.1016/j.jocs.2016.06.006
  * [40] Chen Liu, Florian Frank, Faruk Omer Alpak and Béatrice Rivière  “An interior penalty discontinuous Galerkin approach for 3D incompressible Navier–Stokes equation for permeability estimation of porous media”  In _Journal of Computational Physics_ , 2019, pp. 669–686  DOI: 10.1016/j.jcp.2019.06.052
  * [41] R. Guibert, M. Nazarova, P. Horgue, G. Hamon, P. Creux and G. Debenest  “Computational permeability determination from pore-scale imaging: sample size, mesh and method sensitivities”  In _Transport in Porous Media_ 107.3 Springer ScienceBusiness Media LLC, 2015, pp. 641–656  DOI: 10.1007/s11242-015-0458-0
  * [42] H. Bruus  “Theoretical Microfluidics”, Oxford Master Series in Physics  OUP Oxford, 2008
  * [43] Charu C. Aggarwal  “Neural Networks and Deep Learning”  Springer International Publishing, 2018  DOI: 10.1007/978-3-319-94463-0
  * [44] Lu Lu  “Dying ReLU and initialization: theory and numerical examples”  In _Communications in Computational Physics_ 28.5 Global Science Press, 2020, pp. 1671–1706  DOI: 10.4208/cicp.oa-2020-0165
  * [45] Raphael Schulz, Nadja Ray, Simon Zech, Andreas Rupp and Peter Knabner  “Beyond Kozeny–Carman: Predicting the permeability in porous media”  In _Transport in Porous Media_ 130.2 Springer ScienceBusiness Media LLC, 2019, pp. 487–512  DOI: 10.1007/s11242-019-01321-y
  * [46] Alexander LeNail  “NN-SVG: Publication-Ready Neural Network Architecture Schematics”  In _Journal of Open Source Software_ 4.33 The Open Journal, 2019, pp. 747  DOI: 10.21105/joss.00747
  * [47] Krzysztof M. Graczyk and Maciej Matyka  “Predicting porosity, permeability, and tortuosity of porous media from images by deep learning”  In _Scientific Reports_ 10.1 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1038/s41598-020-78415-x
  * [48] Stephan Gärttner, Peter Frolkovič, Peter Knabner and Nadja Ray  “Efficiency of micro-macro models for surface active two-mineral systems”, 2020  URL: https://www1.am.uni-erlangen.de/research/preprint/pr411.pdf
  * [49] Nishank Saxena, Amie Hows, Ronny Hofmann, Faruk O. Alpak, Jesse Dietderich, Matthias Appel, Justin Freeman and Hilko De Jong  “Rock properties from micro-CT images: Digital rock transforms for resolution, pore volume, and field of view”  In _Advances in Water Resources_ 134, 2019, pp. 103419  DOI: 10.1016/j.advwatres.2019.103419
  * [50] Emmy Compute-Cluster, RRZE, Friedrich-Alexander-Universität Erlangen-Nürnberg, Germany, 2021  URL: https://www.anleitungen.rrze.fau.de/hpc/emmy-cluster/
  * [51] Javier E. Santos, Ying Yin, Honggeun Jo, Wen Pan, Qinjun Kang, Hari S. Viswanathan, Maša Prodanović, Michael J. Pyrcz and Nicholas Lubbers  “Computationally Efficient Multiscale Neural Networks Applied to Fluid Flow in Complex 3D Porous Media”  In _Transport in Porous Media_ Springer ScienceBusiness Media LLC, 2021  DOI: 10.1007/s11242-021-01617-y
  * [52] Stephan Gärttner and Florian Frank  “RTSPHEM – Reactive Transport Solver in Porous Homogenized Evolving Media” Accessed: September 3, 2024, 2021  Department Mathematik, Friedrich-Alexander-Universität Erlangen-Nürnberg  DOI: 10.5281/zenodo.5166669