Abstract
The objective for this work is to develop a datadriven proxy to highfidelity numerical flow simulations using digital images. The proposed model can capture the flow field and permeability in a large verity of digital porous media based on solid grain geometry and pore size distribution by detailed analyses of the local pore geometry and the local flow fields. To develop the model, the detailed pore space geometry and simulation runs data from 3500 twodimensional highfidelity Lattice Boltzmann simulation runs are used to train and to predict the solutions with a high accuracy in much less computational time. The proposed methodology harness the enormous amount of generated data from highfidelity flow simulations to decode the often underutilized patterns in simulations and to accurately predict solutions to new cases. The developed model can truly capture the physics of the problem and enhance prediction capabilities of the simulations at a much lower cost. These predictive models, in essence, do not spatiotemporally reduce the order of the problem. They, however, possess the same numerical resolutions as their Lattice Boltzmann simulations equivalents do with the great advantage that their solutions can be achieved by significant reduction in computational costs (speed and memory).
Velocity dataset: https://doi.org/10.6084/m9.figshare.7889759
Geometry dataset: https://doi.org/10.6084/m9.figshare.7889741
1 Introduction
1.1 Imagedbased flow simulations
Darcy’s principles (Darcy, 1856) describe the fluid flow of singlephase fluids in porous media at low Reynolds numbers, which is of significant importance in earth sciences, hydrology, and petroleum engineering. According to Darcy equation, pressure gradients are linearly proportional to the fluid rate; the proportionality constant is permeability, which is merely a function of pore space topology of porous media irrespective of the fluid type. In numerical flow simulators for porous media, permeability values are obtained based on the data collected from the field and experiments. An accurate quantification of permeability is difficult due to the variations in pore space morphology characteristics. Permeability have been obtained from experiments and also from analytical and empirical expressions that relate permeability to some attributes of the porous media, such as porosity and pore size distribution. The analytical expressions are, however, only approximations for ideal cases while the empirical expressions have utility only in media similar to scenarios for which they were obtained and thus, are inaccurate when applied to a wide range of other media. Experimental approaches are generally preferred if it is not possible to account for all relevant physics by an equation or model; however, they tend to be time consuming and expensive. Furthermore, they do not capture the effect of pore space morphology characteristics on the flow field and thus, on the permeability. For certain properties, such as permeability, tortuosity, and inertial factors of the porous media, highfidelity numerical simulations using digital images have become a credible alternative, enabled by improvements of 3D imaging techniques, numerical methods, and computing power (TakbiriBorujeni et al., 2013; Sanematsu et al., 2015). Appealing aspects of this approach include the ability to probe porescale physics at a level not possible with traditional experiments and the ability to perform an endless set of numerical tests without degrading or altering the sample. There are considerations that can limit this digital approach including whether the imaging technique can resolve all relevant characteristic scales in the pore space and whether numerical algorithms can accurately model the physical processes. Higher resolution, however, mandates higher computation power. In highfidelity numerical simulation models, the expensive computational costs, the intensive memory requirements, and the poor scaling performances have traditionally prevented their applications beyond toy or smallscale problems, even using the modern high performance computing systems.
In imagebased porescale modeling, the domain is discretized into nodes, voxels, or volume elements, and the resulting grid is used to numerically approximate the relevant partial differential equations for flow, namely computational fluid dynamics (CFD). There is a group of numerical modeling techniques that can utilize the voxel data from Xray tomography or similar methods as the numerical grid. This gridding approach has become widely used in porous media studies in conjunction with the Lattice Boltzmann (LB) simulations, and has been proved to be highly effective for simulating fluid flow through porous media
(TakbiriBorujeni et al., 2013).LB simulations have been applied to flow simulations of realistic porous media (Succi et al., 1989; Ferreol and Rothman, 1995; Jin et al., 2004) with the advantage of being flexible in the specification of variables on the complex boundaries in terms of simple particle bounce back and reflection. This flexibility has opened up the potential in its use for modeling and simulating flow in complex media, such as porous rocks. Challenges for applying LB to real problems include finitesize effects and relaxation time dependence of noflow boundaries. In imagebased simulations, the accuracy of the calculated macroscopic properties depends on the spatial resolution of the rock image (Ginzbourg and Adler, 1994; Manwart et al., 2002). However, there is always a tradeoff between image resolution and computational power. Furthermore, in all digital samples, there is a resolution threshold, below which certain flow characteristics, such as recirculation, are not resolved (Maier and Bernard, 2010).
The main advantage of porescale flow simulations is that explicit influence of each impacting factor can be studied by isolating the effect of other parameters. Attempts of this tabulation of all these impacts have not been manageable yet since such a multidimensional parametric study requires comprehensive efforts and time. In this respect, this work aspires to change the status quo and make a transformative leap by combining porescale modeling with physicsbased ML to develop proxy models, which can be used to determine the flow fields at very little additional cost. It is also important to note that, using trained and validated physicsbased datadriven proxy (PBDDP) model will give us a luxury of performing porescale flow simulations, in which computational expenses are not restrictive. Recently, there have been numerous studies on application of machine learning machine learning (ML) in CFD, most of which are limited to building interpretable reduced order models (ROMs). In ROMs, where the number of variables is reduced to simplify the governing equations and the relationships between inputs and outputs, some details are inevitably overlooked. On the other hand, the widespread success of MLbased predictive modeling in other disciplines, such as autonomous cars, suggests a great opportunity to advances in the stateoftheart by combining conventional CFD simulation techniques with predictive capabilities of PBDDP models to truly capture the physics of the problem and enhance prediction capabilities of the simulations at a much lower cost. These predictive models, in essence, do not spatiotemporally reduce the order of the problem. They, however, possess the same numerical resolutions as their CFD equivalents do with the great advantage that their solutions can be achieved by significant reduction in computational costs (speed and memory). Essentially, the predictive models learn the nature of communications among grid cells and decode the spatial correlations between them (auto and crosscorrelations) in the entire computational domain and can accurately predict solutions to completely new sets of simulation runs, from beginning to end.
Recently, deep convolutional networks (CNNs) with hierarchical feature learning capability have outperformed the state of the art in many computer vision tasks, including image classification
(Simonyan and Zisserman, 2014), segmentation (Long et al., 2015), and synthesis (Goodfellow et al., 2014). Despite in classification tasks, where the network predict a single class label for an input image, in many visual tasks, the desired output could be a class label, or a continuous value, assigned to each pixel of the input image.Ciresan et al. (2012) predicted the class label of each pixel by training a network in a slidingwindow fashion which takes a patch around each pixel. This network, then, is able to localize and also is more robust to overfitting the training data, i.e., generated patches, is much larger than the number of training images. However, this framework is quite slow due to separate processing of each patch, which results in a lot of redundancy on overlapping patches. Moreover, such network should deal with the tradeoff between the localization and context. Large patches need many pooling layers that can reduce the localization performance, while small patches only incorporate little context information in the final decision. More recent studies (Long et al., 2015; Seyedhosseini et al., 2013) proposed to fuse the fine to coarse features from multiple layers in different depth. This enables the network to achieve an accurate localization while having a large receptive field (context) at the same time. In the work performed by Ronneberger et al. (2015)
, the authors introduced UNet which employed contracting path in its AutoEncoder architecture to capture context and enable precise localization. Furthermore, training a very deep neural network is quite a challenging task. More specifically, it is hard for a deep network to find an optimal solution compared to shallower counterparts. One of the main issue in training a deep network is the vanishing gradient problem, making it difficult to tune the parameters of the early layers in the network
(Glorot and Bengio, 2010). In the past couple of years, multiple training strategies have been proposed to train a deep neural network effectively, including deep supervision in hidden layers (Lee et al., 2015), initialization scheme (Glorot and Bengio, 2010), and batch normalization
(Ioffe and Szegedy, 2015). He et al. (2016)introduced residual connections in which they employ additive merging of signals to improve the training speed, and gradient flow through the networks.
1.2 Application of datadriven modeling in engineering problems
Applications of ML have gained lots of popularity in the past few years throughout various industries. This rise in popularity is due to new technologies, such as sensors and highperformance computing services, e.g., Apache Hadoop and NoSQL that enable bigdata acquisition and storage in different fields of studies, such as social networks, stock market, natural language processing, oil and gas industry, and automotive industry. The application of ML in CFD has gained considerable interest recently, mostly to build reduced order models (ROMs), where the number of variables is reduced to define simpler relationships between inputs and outputs. However, in such applications of ML in CFD, it is inevitable to overlook some details. On the other hand, predictive ML techniques suggest a greater opportunity, when the conventional CFD simulation techniques are combined with predictive capabilities of PBDDP models. Such approaches can truly capture the physics of the problem and enhance prediction capabilities of the simulations at a much lower cost.
Unlike automotive industry, the application of Artificial Intelligence (AI) in CFD has been limited to interpretable models from data
(Shelton et al., 1993; Tracey et al., 2013; Kutz, 2017), and predictive models are yet to be employed. The widespread success of predictive modeling in complex problems suggests a great opportunity to advances in the stateoftheart by combining conventional CFD simulation techniques with ML predictive modeling to truly capture the physics of the problem and enhance prediction capabilities of the simulations at a much lower cost. This can be achieved by developing physically interpretable spatiotemporal simulations of complex CFD problems and introducing significant reduction in computational cost (speed and memory).2 Samples
In this work, LB simulations are performed over computer generated structures, which provide a number of advantages for testing porescale modeling algorithms. The most intuitive advantage is the ability to fully control the pore structure. Another advantage related to imagebased modeling is that the geometricbased data, e.g., locations and sizes of solid grains in a random packing can be converted to voxel data at any desired image resolution without segmentation error. Computergenerated packings have been widely used to simulate granular materials. In some cases, unconsolidated sphere packs have been modified using procedures that mimic diagenetic processes, thus producing consolidated materials (Bosl et al., 1998; Zhan et al., 2010; Kameda, 2004). We generated twodimensional circle pack images of size 256x256 pixels, each consisting of 8 grains (circles) with 50 pixels diameter with random positions. A total of 3550 images are generated for LB simulation runs to determine the permeability.
A body force approach, which is an alternative to specifying pressure values at the inlet and outlet of the domain, is used. Periodic boundary condition is applied on all the external faces. The bounceback boundary scheme is used to implement the noflow boundary conditions at the voidsolid interfaces. Having a reasonably large pore and grainsizes of the porous media, calculation of the permeability are done without substantial numerical errors (finitesize errors and relaxationtime dependence of the noflow boundaries) (TakbiriBorujeni et al., 2013; TakbiriBorujeni and Tyagi, 2016).
3 Lattice Boltzmann Method
The LB simulation method is based on kinetic theory and can be used to simulate many hydrodynamic systems (Sukop, 2006). In this method, positions of particles are limited to nodes of a lattice with equal spacing. The LB equation with streaming and single velocity relaxation operator (LBGK) collision is
(1) 
Directions i in which fluid particles can move Discrete distribution functions in i direction in velocity space Discrete equilibrium distribution functions in i direction [G]Relaxation time in which are directions in which fluid particles can move, is the relaxation time, are the discrete distribution functions in velocity space, and and are the equilibrium distributions,
(2) 
Weight factors for i direction in Equation 2 Sound speed in the fluid where are weight factors specific to different directions, is the speed of sound in the fluid, is the fluid velocity, and is the fluid density. In this work, the
model (two dimensions and nine directions of fluid movement) is used. Velocity vectors for this model are described below,
(3) 
In LB simulations, parameterized values of the lattice constants and fluid in lattice units are used in simulation while correspondence between the real physical system being simulated and the parameterized simulation is achieved through Reynolds number (principle of dimensional similarity) (Chukwudozie and Tyagi, 2013).
LBM simulations in this study are performed using the Parallel Lattice Boltzmann Solver (PALABOS, 2012). The PALABOS (Parallel Lattice Boltzmann Solver) code is used for solving the flow problems in this study (Latt, 2009).
3.1 Permeability.
Permeability is calculated from
(4) 
Permeability tensor of the porous medium
Dynamic pressure gradient in the fluid in which is the permeability tensor of the porous medium, is average velocity of the fluid in the domain, is viscosity of the fluid, and is dynamic pressure gradient of the fluid. Velocity values in each node are computed in all directions using the LB simulations to determine the permeability tensor.4 Deep convolutional neural network
In order to achieve an accurate and efficient model, we employ a deep convolutional neural network (DCNN) based on contracting paths and residual blocks. Since the network consists of only convolutional layers, it can take any arbitrarysized image as input and generate an output of similar size. For downsampling, we use convolutional layers with increased stride instead of pooling layers. After a series of successive strided convolution, the spatial size of feature maps become much smaller than that of the input image. To increase the computational capacity of the network, the generated feature maps by the last strided convolution is followed by multiple residual blocks before upscaling to the same size as the input image. The residual connections improve the gradients flow at the training time. Finally, to rescale the feature maps to the size of input image, we exploit Nearest Neighbor (NN) upsampling followed by a convolutional layer, instead of deconvolutional layers
(Zeiler et al., 2010) to prevent checkerboard artifacts. Generally, as we go deeper in a DCNN, the size of receptive field increases, which means the learned feature maps represent more abstract and global contextual features. However, the information about the exact local structure of the image may be lost. On the other hand, the feature maps in early layers, which have smaller receptive fields, preserve the local structure information. This information is critical for effective velocity field predictions. Consequently, to preserve the local structure information, high resolution features from the contracting path (downsampling) are combined with the output of the NN upsampling layer. Then, the subsequent convolution learns to produce a more precise output based on this information. Exploiting the learned discriminative features by the proposed DCNN, we can produce an accurate prediction of velocity fields.4.1 Architecture setup
Figure 1
shows the architecture of the proposed network. It consists of 6 strided convolutions which reduce the size of input by a factor of 64, followed by four residual blocks. Each residual block consists of two 3x3 convolutional layers. Finally, In order to obtain the final prediction map, we add six subsequent upsampling blocks on top of the residual blocks. The input to each upsampling block is the feature maps of the previous layer concatenated in depth with those of contracting downsampling path. As mentioned earlier each upsampling block comprises successive NNupsampling and 3x3 convolution with unit stride. Note that all the convolutions are followed by a Batch Normalization and Relu activation function. Since our input (geometry input image) is the simulated velocity fields we employ reflection 1x1 padding for all the convolutions.
To train the network, we first normalize the velocity maps. To reduce overfitting, we employed data augmentation to enlarge the training data samples. To this end, we randomly flip the maps vertically with a probability of
. The input density and their corresponding velocity maps are used to update the parameters of the network minimizing the norm error:(5) 
Calculated xdirection velocity by LB simulations Estimated xdirection velocity by the network input map (image) norm error where is the input map, is the estimated velocity map by the network, is the ground truth velocity map, and is number of samples in the training dataset. Adam optimization technique is used with a learning rate of , and an weight decay of
. The network was implemented in Pytorch running on a NVIDIA TITAN X GPU. The network is trained for 100 epochs and the model with the minimum error on validation is selected.
5 Results
To develop the model, the velocity values for the entire images are normalized between zero and one (the minimum value the velocity values is transformed into zero, the maximum value is transformed into one, and every other value is transformed into a decimal between 0 and 1). All the simulation cases are divided into two sections. The first section with 85% of the data (3,000 image pairs out of 3,550 total pairs) is used to train the model while the remaining 15% of the data are used as test data. The training and validation loss curve for the training process is depicted in Figure 2.
The test portion of the data, which is not used in the training process, is only used to examine the predictive capabilities and the robustness of the model. All the simulations are tested to verify that they have reached steady state conditions, where kinetic energy of the system becomes constant. The binary images are used for the LB simulations (Figure 3a). The regions away from the solidpore interfaces exhibit higher velocity values compared the ones adjacent to the interfaces. Contour plots of xdirection velocity for the developed model (Figure 3b) and those taken from the LB simulations (Figure 3c) show similar behavior. The contour plots in Figure 3c are adjusted based on the minimum and maximum values obtained in contours in Figure 3b. The behavior in Figure 3b and Figure 3c are fairly similar, except for the high (hotcolored) values in Figure 3c. The positions of the circles (zerovelocity valued pixels in Figure 3b and Figure 3c) are accurately predicted. The error plots (Figure 3d) exhibit errors in the domain. The error for most of the cases are bounded within 20% errors (Figure 3e), confirming the plausibility of the approach to replicate costly numerical simulations.
In this work, a pointbypoint comparison of the predicted normalized velocity fields using DDP model and the LB simulations is performed. Velocity profiles in a vertical and horizontal crosssection of the simulation domain for one of the test cases are depicted in Figure 4. By inspection of these plots, one can see that DDP model mimics the LB simulation results with negligible errors at points in the solid grain regions (circles). Please note that the velocity values at these regions is not zero due to the fact that the minimum value the velocity value in the dataset in transferred to zero, resulting the zerovalued velocities to be transformed into a decimal value (0.07) in the normalized form. The DDP predicted velocity values along the vertical crosssection exhibit deviations from the LB simulation results (Figure 4b). However, at regions close to the interfaces (noflow boundary), the predictions are close to the LB simulation results. In horizontal crosssection in Figure 4c is selected in a high velocity region. The error in predicted velocity values at inside the grains is negligible. However, the deviations increase as velocity increases. In this region, DDP predicts lower values that those calculated using LB simulations.
The predicted mean of xdirection velocity values for all test cases are shown in Figure 5. All the points are projected along the unitslope line, which shows that predicted values are fairly close to the LB simulation results.
Based on the results, the predicted permeability results for the twodimensional domains (images) are predicted using the PBDDP model with high accuracy.
6 Conclusions
A datadriven proxy to highfidelity numerical flow simulations is presented by employing a deep convolutional neural network based on contracting paths and residual blocks. The network consists of only convolutional layers and can take any arbitrarysized image as input and generate an output of similar size. The developed model captures the flow field and permeability for samples at the grid level that had not been used in the development of the model. Based on the results, the predicted permeability results for the twodimensional domains (images) are predicted using the PBDDP model with high accuracy. This work aspires to make a transformative leap by combining fluid flow modeling with physicsbased machine learning to develop proxy models, which can be used to determine the flow fields at very little additional cost.
References
 Bosl et al. (1998) Bosl, W. J., J. Dvorkin, and A. Nur (1998). A study of porosity and permeability using a lattice boltzmann simulation. Geophysical Research Letters 25(9), 1475–1478.
 Chukwudozie and Tyagi (2013) Chukwudozie, C. and M. Tyagi (2013). Pore scale inertial flow simulations in 3d smooth and rough sphere packs using lattice boltzmann method. AIChE Journal 59(12), 4858–4870.

Ciresan et al. (2012)
Ciresan, D., A. Giusti, L. M. Gambardella, and J. Schmidhuber (2012).
Deep neural networks segment neuronal membranes in electron microscopy images.
In Advances in neural information processing systems, pp. 2843–2851.  Darcy (1856) Darcy, H. (1856). Les fontaines publiques de la ville de dijon, dalmont. Paris: Dalmont.
 Ferreol and Rothman (1995) Ferreol, B. and D. H. Rothman (1995). Latticeboltzmann simulations of flow through fontainebleau sandstone. In Multiphase flow in porous media, pp. 3–20. Springer.
 Ginzbourg and Adler (1994) Ginzbourg, I. and P. Adler (1994). Boundary flow condition analysis for the threedimensional lattice boltzmann model. Journal de Physique II 4(2), 191–214.
 Glorot and Bengio (2010) Glorot, X. and Y. Bengio (2010). Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 249–256.
 Goodfellow et al. (2014) Goodfellow, I., J. PougetAbadie, M. Mirza, B. Xu, D. WardeFarley, S. Ozair, A. Courville, and Y. Bengio (2014). Generative adversarial nets. In Advances in neural information processing systems, pp. 2672–2680.

He
et al. (2016)
He, K., X. Zhang, S. Ren, and J. Sun (2016).
Deep residual learning for image recognition.
In
Proceedings of the IEEE conference on computer vision and pattern recognition
, pp. 770–778.  Ioffe and Szegedy (2015) Ioffe, S. and C. Szegedy (2015). Batch normalization: Accelerating deep network training by reducing internal covariate shift. arXiv preprint arXiv:1502.03167.
 Jin et al. (2004) Jin, G., T. W. Patzek, D. B. Silin, et al. (2004). Direct prediction of the absolute permeability of unconsolidated and consolidated reservoir rock. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
 Kameda (2004) Kameda, A. (2004). Permeability evolution in sandstone: Digital rock approach. Ph. D. thesis, Stanford University.
 Kutz (2017) Kutz, J. N. (2017). Deep learning in fluid dynamics. Journal of Fluid Mechanics 814, 1–4.
 Latt (2009) Latt, J. (2009). Palabos, parallel lattice boltzmann solver. FlowKit, Lausanne, Switzerland.
 Lee et al. (2015) Lee, C.Y., S. Xie, P. Gallagher, Z. Zhang, and Z. Tu (2015). Deeplysupervised nets. In Artificial Intelligence and Statistics, pp. 562–570.
 Long et al. (2015) Long, J., E. Shelhamer, and T. Darrell (2015). Fully convolutional networks for semantic segmentation. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 3431–3440.
 Maier and Bernard (2010) Maier, R. and R. Bernard (2010). Latticeboltzmann accuracy in porescale flow simulation. Journal of Computational Physics 229(2), 233–255.
 Manwart et al. (2002) Manwart, C., U. Aaltosalmi, A. Koponen, R. Hilfer, and J. Timonen (2002). Latticeboltzmann and finitedifference simulations for the permeability for threedimensional porous media. Physical Review E 66(1), 016702.
 Ronneberger et al. (2015) Ronneberger, O., P. Fischer, and T. Brox (2015). Unet: Convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computerassisted intervention, pp. 234–241. Springer.
 Sanematsu et al. (2015) Sanematsu, P., Y. Shen, K. Thompson, T. Yu, Y. Wang, D.L. Chang, B. Alramahi, A. TakbiriBorujeni, M. Tyagi, and C. Willson (2015). Imagebased stokes flow modeling in bulk proppant packs and propped fractures under high loading stresses. Journal of Petroleum Science and Engineering 135, 391–402.
 Seyedhosseini et al. (2013) Seyedhosseini, M., M. Sajjadi, and T. Tasdizen (2013). Image segmentation with cascaded hierarchical models and logistic disjunctive normal networks. In Proceedings of the IEEE International Conference on Computer Vision, pp. 2168–2175.
 Shelton et al. (1993) Shelton, M., B. Gregory, S. Lamson, H. Moses, R. Doughty, and T. Kiss (1993). Optimization of a transonic turbine airfoil using artificial intelligence, cfd and cascade testing. In ASME 1993 International Gas Turbine and Aeroengine Congress and Exposition, pp. V03AT15A012–V03AT15A012. American Society of Mechanical Engineers.
 Simonyan and Zisserman (2014) Simonyan, K. and A. Zisserman (2014). Very deep convolutional networks for largescale image recognition. arXiv preprint arXiv:1409.1556.
 Succi et al. (1989) Succi, S., E. Foti, and F. Higuera (1989). Threedimensional flows in complex geometries with the lattice boltzmann method. EPL (Europhysics Letters) 10(5), 433.
 Sukop (2006) Sukop, M. (2006). DT Thorne, Jr. Lattice Boltzmann Modeling Lattice Boltzmann Modeling. Springer.
 TakbiriBorujeni et al. (2013) TakbiriBorujeni, A., N. Lane, K. Thompson, and M. Tyagi (2013). Effects of image resolution and numerical resolution on computed permeability of consolidated packing using lb and fem porescale simulations. Computers & Fluids 88, 753–763.
 TakbiriBorujeni and Tyagi (2016) TakbiriBorujeni, A. and M. Tyagi (2016). Multiscale modeling of permeability and nondarcy factor in propped fractures. Hydraulic Fracturing Journal 3(2), 64–74.

Tracey
et al. (2013)
Tracey, B., K. Duraisamy, and J. Alonso (2013).
Application of supervised learning to quantify uncertainties in turbulence and combustion modeling.
In 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, pp. 259.  Zeiler et al. (2010) Zeiler, M. D., D. Krishnan, G. W. Taylor, and R. Fergus (2010). Deconvolutional networks.
 Zhan et al. (2010) Zhan, X., L. M. Schwartz, M. N. Toksöz, W. C. Smith, and F. D. Morgan (2010). Porescale modeling of electrical and fluid transport in berea sandstone. Geophysics 75(5), F135–F142.