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

In recent years, convolutional neural networks (CNNs) have experienced an increasing interest in their ability to perform a 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 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. The robustness of this approach is reflected by very high prediction accuracy, which is observed for a variety of sandstone samples from archetypal rock formations.


Introduction
Artificial neural networks can accelerate or replace classical methods for estimating various hydrological and petrophysical properties of artificial and natural rock [50].In [6], deep neural networks were successfully used to approximate tomography operators to reconstruct wave velocity models from seismic data.Moreover, convolutional neural networks (CNNs) were used in [30] to determine effective electrical parameters in a homogenized model for electric conductivity.Likewise, CNNs have proven their ability to provide fast predictions of scalar permeability values directly from images of the pore space of geological specimens [23,35,53].Finally, CNNs were successfully used to replace standard solution schemes to predict effective permeabilities in 2D multiscale flow simulations [18].This paper aims at contributing to the deployment of machine learning techniques in effective permeability prediction 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.We use the term Phy-CNN throughout this paper to underline the incorporation of an additional physically motivated input parameter to the network.This approach is to be carefully distinguished from implementing physical laws into the loss-function as performed in [11] to estimate flow velocity fields.
Specimens of natural rock are typically obtained by microcomputed tomography (µCT) scanning [10,52] or indirectly from 2D colored images [5].As an important characteristic quantity of porous media, permeability measures the ability 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, i. e., the accuracy of the forward simulations in our case, as shown in [56].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 [49].
In the literature, various methods are available to compute the stationary flow on complex geometries as constituted by the pore space of porous media, a thorough comparison of which is found in [54].Most commonly, lattice Boltzmann methods (LBM) are used to solve the flow problem on a discrete modeling basis, cf.[27].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 [25].Yet, since a global equilibrium has to be reached, complex geometries containing thin channels may cause LBM to converge slowly or even diverge [23,37], especially in the case of very constricted and therefore typically underresolved pore morphologies.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 [8].However, the novel multiple-relaxation-time (MRT) scheme addresses this problem to a great extent (e. g., [3]) at the cost of a moderate additional computational overhead.Having stated that, residual discretization errors stemming from the explicit nature of time integration in LBM schemes still remain in the numerical solution.Especially for under-resolved pore morphologies, these errors may be amplified causing divergence of the numerical flow solution.We further note that classical bounceback rules used for implementing boundary conditions tend to develop oscillations that might locally dominate fluid behavior in scenarios with thin channels, cf.[55].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 [37].On the other hand, industrial-grade state-of-theart LBM solvers are typically proprietary without public access to the code base for academic research purposes (e. g., [3,47]).Moreover, a comprehensive permeability computation benchmarking study [39] 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 [23].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 [42].More precisely, our forward simulation is based on a distributed-parallel Stokes solver utilizing the finite element library MFEM [4].As studied in [49] for simple cylindrical obstacles, such DNS approaches (in this case FEM) deliver permeability values comparable to the ones from LBM simulation.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 [46,53].As shown in [44], carefully chosen specific quantities derived from the pore space such as connectivity indices can deliver reasonable approximation quality for permeability estimation.As illustrated in [6], 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 [12].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

Data labeling
Perform forward simulation to determine permeability labels kcmp.

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.additionally provide our CNN with information that reflects possible thin, channel-like structures.As discussed in [48], 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.

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 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 [31].

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 φ exp = 22.64% and permeability k exp = 386 mD, provided by [32,33].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 [20].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 derived from a 1000×1000×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  In the first step, we extract subsamples of 100×100×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 [43].More precisely, we sweep a 100×100×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 [23].Furthermore, by rotating the original 100×100×100 voxel subsamples by 90 • 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 erosion and dilation 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•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 using natural rock sample for training, cf.[23,37].However, using artificially generated pore-space geometries, even larger data sets including 90 000 data point are available in the literature, cf.[35].More precisely, these data are obtained from stochastic computer models, evading the necessity of expensive imaging and segmentation of real rock.
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 (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 [24].
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 kB per subsample.Accordingly, our whole library of data samples allocates only 1.25 GB of disk space .Therefore, it is manageable on standard personal computers and file read/write access is reasonably cheap.
We note that a subsample of size 225 µm is too small to be a representative elementary volume (REV) for most sandstone types, cf.[20,33].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.Moreover, hierarchical neural networks have proven to be a powerful tool to leverage the performance of CNNs predicting permeability to samples of REV-size (in the order of 500 3 voxels), cf.[38].Therefore, we believe that our proposed methodology also benefits existing REV-scale models.

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 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 Ω for the fluid velocity u and pressure p.The discretization uses arbitrary order (stable) Taylor-Hood or reduced-order stabilized Taylor-Hood mixed finite elements.In this paper, we consider voxel meshes only, yet the usage of unstructured grids by using obvious respective discrete spaces is straight forward.
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 cmp to train and validate our PhyCNN.
We note that-despite the finite resolution of the µCT images-the voxel scans are henceforth considered to represent the actual (ground-truth) pore geometry.Hence, errors arising from imaging are neglected as they pose an independent part of the total methodology that is beyond the scope of this paper.

Computation of the flow field
We consider the stationary Stokes equation for a Newtonian fluid in the nondimensionalized form, where Ω ⊂ (0, 1) 3 is a domain that consists of the union of fluid voxels of a considered 100×100×100 voxel subsample (i.e. the subsample is inscribed into the unit cube).In ( 1 ] is invariant with respect to 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 where Γ N := (x, y, z) ∈ ∂Ω | x ∈ {0, 1} and e x denoting the unit vector in x direction.On the remaining boundary Γ D := ∂Ω \ Γ N , no-slip boundary conditions are prescribed, The weak formulation of ( 1) is discretized by generalized Taylor-Hood pairs of spaces, Q 3 +1 /Q , where Q denotes the local space of polynomials of degree at most in each variable x, y, z, cf.[13,14].For = 0, the pressure space Q 0 = P 0 is discontinuous and consists of elementwise constants.The respective "reduced Taylor-Hood pair" Q 3  1 /P 0 requires stabilization (see below) due to a lack of discrete inf-sup-stability [13].Let velocity magnitude i = 1, . . ., m denote the basis functions for the global discrete spaces for velocity and pressure, respectively.The discrete velocity u h = u h (x, y, z) and discrete pressure p h = p h (x, y, z) then have the representation with degree-of-freedom vectors x u ∈ R n , x p ∈ R m , which are unique solutions of the linear system The sparse blocks in A are the vector-Laplacian matrix A ∈ R n,n and the divergence matrix and C is a stabilization matrix that is required only for lowest order = 0 to guarantee the full rank of A.
We choose C as in [13], (3.84), in which case, C can be interpreted as a pressure-Laplacian discretized by cell-centered finite volumes [16].For > 0, C is set to zero.In either case, 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 [51].The chosen precondition operator for A in (3) is the symmetric and positively definite block-diagonal matrix with W being the pressure-mass matrix [W k,l ] := Ω ψ k ψ l .The action of the inverse 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 [15]).Since W is spectrally equivalent to the (negative) Schur complement BA −1 B T of A (for > 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 [7,34,36] (however, it highly depends on the geometry of the domain).MINRES with preconditioning as described above belongs to the stateof-the art Stokes solvers in the high-performance computing context [19].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.

Permeability estimation
We deduce the permeability value k cmp of interest from the previously calculated Stokes velocity 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 and thus disjoint subsets of Γ N .From the solution (u, p) of (1), we compute the approximated permeability k cmp [m 2 ] by the classical Darcy law, which reads in nondimensionalized form, with (dimensionless) volume flow rate Q, and (dimensionless) area-averaged inflow and outflow pressures, P in , P out , given by [28] Q From (4a), k cmp is derived from the Darcy number where L c is the characteristic length, in our case, the edge length of a 100 3 -voxel subsample, i. e., L c = 225 µm.We choose Re equal to one, since Re does not influence the permeability value k cmp (a rescaling of u by Re −1 implies a rescaling of Q by Re −1 and therefore, 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 [22], 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.

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 cmp obtained from (4a) to the analytical value k ana .In order to have analytical results available, we restrict our considerations to viscous flow through rectangular channels Ω a,b .This is supposed to constitute a sufficient benchmark, since the approximation quality of the computed permeability k 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 of width a ∈ 0, 1 2 and height b ∈ 0, 1 2 .An analytical expression for its permeability k ana is [9] , where denoting the limit of K j for j → ∞ by K ∞ .By application of the triangular inequality, 0 < min(a, b) ≤ max(a, b), and | tanh(x)| < 1 ∀x ∈ R, we obtain the following approximation error bound: Piecewise application of Jensen's inequality to the convex function (2x − 1) −5 finally yields the estimate: As such, we obtain at least six significant digits using K 10 for the calculation of the analytical reference permeability k ana .) and mesh refinement levels.As in (2), m and n are the numbers of DOF for velocity u h and pressure p h , respectively.The relative tolerance of MINRES is set to 1.0E−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 cmp obtained from (4a) on Ω 0.06,0.03using different polynomial orders 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 Q 3  2 /Q 1 elements ( = 1) on the original grid.In particular, this choice poses a significant improvement in accuracy over the use of lowest-order spaces Q 3  1 /Q 0 ( = 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 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 ( = 1) for permeability labeling in our forward simulations, cf.Section 2.2.1.We emphasize that our scheme does not require mesh refinement, which is typically required in LBM simulations of natural rock.
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.

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 σ by for a number N of real-valued data t i and arithmetic mean value t, measuring the expected deviation from t.In order to quantify the correlation between different measures of a data set, we further introduce the coefficient of determination 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

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 [2].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.

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 ∈ {0, . . ., 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 ≤ i ≤ L on input data x i−1 ∈ R Ni−1 is given as an affine-linear transformation, followed by application of a nonlinearity (activation function) σ i : We refer to W i ∈ R Ni,Ni−1 as the weight matrices and to b i ∈ R Ni 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. [2].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 α ∈ (0, 1), the LeakyReLU(α) function is defined by LeakyReLU(α; x) := max(0, x) − α 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 [29], which arise formally by setting α = 0. Furthermore, backpropagation during the learning phase is facilitated.Our implementation is conducted using the Deep Learning Toolbox in Matlab R2021a [31].

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.

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.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 alone may even provide sufficient information to the network to accurately perform predictions, cf.[6], underlining the potential of physics-informed strategies.
This approach has been successfully applied to the prediction of permeabilities from pore-space geometries: In [37], the Euclidean distance map, maximum inscribed sphere radii, and time of flight maps are inserted as inputs to determine the overall fluid field.In [53], porosity and surface area have been used to increase predictions performance for 2D image data.Similar observations have been made on 3D data where the additional consideration of porosity and tortuosity significantly reduced the number of outliers, especially when training on small data sets, cf.[45].The specific choice of these parameters is related to the Kozeny-Carman formula estimating permeability from porosity and surface area / tortuosity.As indicated in [41], the proposed relation deteriorates in quality for increasingly complex geometries.However, [46] 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 [46].
As indicated by our data set statistics in Figure 4, porosity and surface area are only weakly correlated with the computed permeability k cmp for our complex 3D setting.For a given surface area A cmp , we observe data samples covering four orders of magnitude in their respective permeability values, about half the range is observed for fixed porosities φ 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 [31].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.Our final network architecture is presented in Table 3 and Figure 5.The structure is based on the findings of [23], 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 cutoffs as nonlinearities (LeakyReLU).Furthermore, the maximum flow value f 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.

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, ω), where the nodes N correspond to distribution nodes, edges E to connections by pipes, and the weights ω encode the maximal capacity of each pipe.For any two nodes n 1 = n 2 ∈ 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 [12].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 ω 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 card(N ) − 2 card(E).Moreover, we add a further node n in connected to each voxel of the sample's inflow face, as well as n 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 in , n out .
For the samples considered in this paper, the computational effort for calculating f max using Matlab's built-in routine maxflow with about one second per sample on a single CPU core, also cf.Section 4, is reasonably small compared to the solution of stationary Stokes equations as described in Section 4.3.By the min-cut max-flow theorem, cf.[12], the maximum flow problem allows for another interesting n in n out interpretation.In the setup introduced within this section, f 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 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 max of 1230, the permeability of Subsample 0 is governed by wide channels.On the other hand, exhibiting f max of 55, narrow pore throats hamper the fluid flow of Subsample 9213.We further note that for moderate tortuosity, f 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 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 max against the calculated permeabilities k 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 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 cmp .Therefore, the permeability estimation using formula ( 5) is considered a suitable initial guess, which the CNN is able to improve by relating to features extracted from the pore geometry.

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.

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 mD, 50 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 [21].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.

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.In order to further stress the robustness, additionally validate our PhyCNN on samples originating from different types of sandstone that were not used for training.Using the data provided by [32,33], 200 additional subsamples were extracted from each Berea and Castlegate µCT scan, see   4 lists the prediction quality indices of other related CNNs for comparison.To clearly and consistently separate the different approaches investigated, we refer to the PhyCNN described in Section 3.1 more precisely as PhyCNN(f max ) within this comparison.In case of the PhyCNN(φ-σ), the general structure depicted in Figure 5 is maintained while replacing the maxflow input variable by the porosity φ and surface area σ, and retraining the network, cf. Figure 4.As the data of Table 4 show, the PhyCNN(φ-σ) cannot reach the accuracy of our PhyCNN(f max ).Moreover, robustness appears significantly lower since accuracy on other sandstone types than the training data decreases rapidly.Especially for Berea sandstone, the results are poor, which might be a direct consequence of the significantly different pore-space characteristics in comparison to Bentheimer, cf.Table 1.Furthermore, we also compare to the analogous network without any additional input quantities (plain CNN), resulting in slightly worse predictions than the PhyCNN(φ-σ).Finally, we note that both of these CNNs provide less accurate predictions on validation samples than f max alone (plain f max ) being fitted to the Bentheimer training data set using a power law approach and measured with respect to R 2 -log.As such, we conclude that f max already contains highly relevant information regarding the samples' permeability delivering more accurate predictions in the logarithmic measure than standard CNN approaches.However, our PhyCNN(f max ) improves the results of pure f max by extracting further information directly from the pore-space, cf. Figure 5.
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 [18] to erode or dilate 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 cmp using the Stokes solver and k 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 cmp to exceed the trained permeability range [50 mD, 50 D] beyond φ 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 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 [40].

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 ( =1, cf.Section 2.2.3) in parallel on 50 compute nodes of the Emmy compute cluster at RRZE [1], each being equipped with two Intel®Xeon 2660v2 'Ivy Bridge' processors (10 physical cores per chip, i. e., 50•2•10=1000 cores total, no hyper-threading) and 64 GB of RAM.The linear solver (MINRES) converged within 969 iterations on average with a relative tolerance of 1.0E−6 in the norm • 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.[23].
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 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 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.

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) implementations while successfully addressing convergence challenges that typically appear on complex 3D geometries including poorly resolved subdomains.
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.
Finally, we note that the refined methodology presented in this paper may generalize to larger sample sizes approaching the REV-scale.Recent publications suggest hierarchical multiscale neural networks [38] 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.

Figure 1 :
Figure 1: Overall workflow in flow chart representation.
), u = u(x, y, z) denotes the (dimensionless) fluid velocity, p = p(x, y, z) the (dimensionless) pressure and Re := ρ U c L c /µ the Reynolds number of the system with characteristic length L c [m], characteristic velocity U c [m s −1 ], fluid density ρ [kg m −3 ], and fluid viscosity µ [Pa s].Note that the permeability k cmp [m 2

Figure 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.

Figure 4 :
Figure 4: Correlation of different physical measures of the underlying geometry with the computed permeability values kcmp for Bentheimer data samples.From left to right: Computed surface area Acmp (R 2 = 0.055), computed porosity φcmp (R 2 = 0.547), specific surface area Aspec (R 2 = 0.455), maximum flow on graph representation fmax (R 2 = 0.869).The given R 2 values refer to a linear regression model.

Figure 5 :
Figure5: Schematic representation of physics-informed convolutional neural networks (PhyCNNs).Numbers at the layers denote their dimension (number of neurons).Note that the cubic convolutional kernels and their dimension are illustrated in 2D for clarity.The scalar valued physical input circumvents the image convolution layers and is directly fed into the subsequent dense layers.Graphics produced using NN-SVG[26].

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

Figure 8 :
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.

Figure 9 :
Figure 9: Testing generalization ability on artificially distorted data samples.PhyCNN prediction performance on eroded and dilated pore spaces exhibiting a large range of porosities φcmp.For each manipulated geometry, Stokes simulation data kcmp (red) and CNN prediction k 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.
µm yielding an overall cube side length of 2.25 mm.

Table 2 :
Computed permeabilities kcmp of a channel with 6×3 voxel rectangular cross section for different polynomial orders of the pressure space (cf.Section 2.2.1

Table 3 :
Layer structure of our PhyCNN.The physical scalar input2 is expanded to dimension 64×1 and concatenated with the output of 'dense1' using a depth concatenation layer.Nomenclature: LeakyReLU(α): leaky rectified linear unit, slope α on negative inputs.

Table 1 ,
using the same methods as described in Section 2. Achieving R 2 values of 89.43% and 94.13% (logarithmic)

Table 4 :
Overview prediction quality.In the first block, mean permeability and standard deviation σ(kcmp) are listed for the respective data sets.Subsequently, R 2 values in natural and logarithmic scale for the training and validation Bentheimer data as well as for Berea and Castlegate sandstone are presented as achieved by our PhyCNN, i. e., the maxflow-informed CNN (PhyCNN(fmax)), the porosity/surface-area-informed CNN (PhyCNN(φ-σ)), and a plain CNN without additional inputs.Finally, the additional input quantity fmax is considered in a power-law (linear fitting applied to logarithmic variables) regression model (plain fmax).onBereaaswellas94.55%and94.51%(logarithmic)onCastlegate, the networks proves excellent generalization properties across different types of pore geometries, cf.Figure8, Figure2, and Table1.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, Table4marks 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 mD, 50 D].More precisely, eleven Berea data samples exhibit a permeability k cmp below 50 mD as well as two Castlegate data samples.Moreover, Table