ORIGINAL PAPER
Development of Kriging-approximation simulated annealing
optimization algorithm for parameters calibration of porous media
flow model
Ming-Che Hu1• Chia-Hui Shen1•Shao-Yiu Hsu1 •Hwa-Lung Yu1•Krzysztof Lamorski2• Cezary Sławin´ski2
Published online: 18 January 2019
Ó Springer-Verlag GmbH Germany, part of Springer Nature 2019
Abstract
This research proposes an innovative Kriging-approximation simulated annealing (KASA) optimization algorithm to increase optimization efficiency and reduce the computation time for the parameter calibration of simulation model. The newly developed KASA optimization algorithm utilizes the simulated annealing algorithm to search the global optimum; meanwhile, Kriging approximation, a statistical estimation method, is incorporated with simulated annealing to interpolate unknown objective values in solution spaces. Furthermore, this research establishes a network-based porous media flow simulation (NET-PFS) model and then KASA is applied in the calibration of the NET-PFS representative pore network. NET-PFS is a pore network based model constructing a representative pore network to approximate soil characteristics and pore geometry with limited access to pore-scale imaging processes. NET-PFS is applied to estimate the permeability of a sand-packing porous media. NET-PFS establishes a framework for simplifying the pore network but remaining the same hydraulic conductivity and the flow status of pore networks with limited information about the pore structure. In the case study, a quartz sand-packing porous media is scanned by X-ray micro computed tomography. The NET-PFS model is applied to estimate the hydraulic conductivity and flow velocity distribution from the original pore network. The results demonstrate the proposed KASA algorithm effectively calibrated the NET-PFS model; in addition, a representative pore network and the determined flow status in the pore network is presented by NET-PFS.
Keywords Kriging approximation Simulated annealing Network-based porous media flow Representative pore network List of symbols Indices i Unobserved points, i¼ 1; 2; 3; . . .; m j1; j2; k Observed points, j1; j2; k¼ 1; 2; 3; . . .; n r1; r2 Pores, r1; r2¼ 1; 2; 3; . . .; q1 r3 Inner pores, r3¼ 1; 2; 3; . . .; q3 r4 Boundary pores, r4¼ 1; 2; 3; . . .; q4
Parameters and variables
alr1;r2 Length of the arc between pores j1 and j2 arr1;r2 Diameter of the arc between pores j1 and j2 dv Dynamic viscosity
dj1;j2 Distance between points j1 and j2 objnew New objective value
objcurrent Current objective value
p Number of simulated annealing iterations prr4 Given node pressure at node r4 of the NET-PFS
model
qa Number of arcs in the pore network q1 Number of pores in the pore network q3 Number of inner pores in the pore network q4 Number of boundary pores in the pore network temp Temperature of the simulated annealing
algorithm
tk Computational time of the Kriging approximation for the NET-PFS model
& Shao-Yiu Hsu syhsu@ntu.edu.tw
1 Department of Bioenvironmental Systems Engineering,
National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan
2 Institute of Agrophysics, Polish Academy of Sciences,
Dos´wiadczalna 4, 20-290 Lublin, Poland https://doi.org/10.1007/s00477-018-01646-y(0123456789().,-volV)(0123456789().,-volV)
tn Computational time of the NET-PFS model vEi Kriging estimation of a random variable at
unobserved point i vT
i Realization of a random variable at unobserved
point i
vj1 Realization of a random variable at observed point j1
wj1 Kriging weight at point j1
xfr1;r2 Arc flow from nodes r1 to r2 of the NET-PFS model
xpr1 Node pressure at node r1 of the NET-PFS model xpin Pore pressure at the inlet boundary
xpout Pore pressure at the outlet boundary
z Random variable with a uniform distribution in the interval of [0,1]
cðdj1;j2Þ Semivariogram function of distance dj1;j2
k Lagrange multiplier
l Expected values of the random variables r2 Variance of the random variables
Functions
COV Covariance function
ERRORi Estimation error at unobserved point i EXP Expected value
EXPO Exponential function
LAGRi Lagrange function at unobserved point i Qm Total flow rate
VAR Variance function
1 Introduction
This study proposes an innovative Kriging-approximation simulated annealing (KASA) algorithm. The KASA algo-rithm incorporates simulated annealing with Kriging approximation to save computational time of iterative objective function evaluation. The proposed KASA algo-rithm applies Kriging approximation to interpolate the objective value in the solution space rather than to estimate random variables in physical space. The simulated annealing algorithm is often used for water resources management, model calibration, decision making, etc. (Bechler et al.2013; Cao and Ye 2013; Cieniawski et al. 1995; Dougherty and Marryott 1991; Jiang et al. 2018; Marryott et al. 1993). It is a heuristic optimization algo-rithm that seeks the global optimum by simulating the cooling process of heated metals. The algorithm evaluates a new solution if the new solution has a better objective value than the current solution. To jump out of the local optimum, the algorithm still accepts the new solution on the basis of the acceptance probability, if the new solution has a worse objective value than the current solution. The
acceptance probability is negatively proportional to the metal temperature. Since the temperature decreases during the cooling process, the acceptance probability function decreases and it leads to the convergence of the global optimum. Most importantly, the newly proposed KASA algorithm utilizes Kriging approximation, a statistical interpolation method, to reduce the evaluation time of the objective function for simulated annealing search. The Kriging method is often used for statistical interpolation in hydrology, atmospheric processes, geophysics, etc. (de Lavenne et al.2016; Verdin et al.2015; Laaha et al.2014; Oliver and Webster 1990; Yeh 1986). It is a spatial sta-tistical estimation method that can estimate an unknown variable at a given point by calculating the weighted average of known values at given points in the neighbor-hood. The method is also a linear interpolation method that provides an optimal and unbiased estimation of unknown random variables.
In this research, a case study of the KASA optimization algorithm was conducted for parameters calibration of the representative pore network. Then, a porous media flow model was formulated to analyze the flow velocity and pore pressure distribution in a packed quartz sand medium. Traditionally, owing to the complexity of the pore geom-etry, the flows in porous media were treated macroscopi-cally and are described by Darcy’s law. The uncertainty and the heterogeneity of the subsurface geological structure is one of the key issues of flow and mass transportation simulation of the groundwater contamination. In the last three decades, pore-scale modeling has developed rapidly from a pure scientific-oriented method for studying fluid displacement to a predictive tool used in the oil industry (Blunt et al.2013). The rapid development is based on the more sophisticated construction of models for the pore geometry and on the predictive power for the dynamics of multiphase flow (Blunt et al. 2002; Keehm et al. 2004; Oren et al.1998; Piri and Blunt2005; Raeini et al.2015; Steefel et al.2015). Continuum and discrete approaches are commonly used to compute the fluid properties in the pore space. The first is to discretize the voids using a grid derived from a binarized three-dimensional (3D) pore-scale image and then compute the flow and transport on this grid. Nevertheless, most of the existing techniques are compu-tationally demanding. Network modeling is an alternative method, which first extracts a representative network on the basis of the binarized pore-scale image (Blunt et al.2013). Then, the flow and transport are computed, usually semi-analytically, through this network. The computational demand of the pore-network model is much lower than that of the first method. The abovementioned pore-scale simu-lation can provide not only the details of the flow behaviors but also the pore structure properties or the soil charac-teristics, i.e., the parameters of the water retention curves,
which, traditionally, can only be obtained through time-consuming laboratory experiments. An overview of these imaging techniques applied to the earth sciences is pro-vided by Ketcham and Carlson (2001). For example, if a 3D image is available (experimental or synthetically gen-erated), pore networks can be constructed directly from this image, assuming what is considered to be a pore and what to be a throat between pores. (Miao et al.2017) used neural networks to predict the dimensionless hydraulic conduc-tance based on the circularity, convexity and elongation of each pore element. However, binarizing the pore imaging and extracting the pore network are also computationally demanding (Xiong et al.2016). Abovementioned methods are still limited to samples with centimeter scales, which might have contained more than 10,000 pores and pore throats (Dong and Blunt 2009). The need for up-scaling from the typically small imaged volume to larger domains leads to the need for construction of statistically repre-sentative networks.
Statistical methods have been used to reconstruct the equivalent pore-network structure (Xiong et al. 2016). Yeong and Torquato (1998) developed a stochastic method based on simulated annealing to minimize the objective function and they obtained the correct porosity by moving pore space voxels around. Okabe and Blunt (2004) devel-oped a multi-point statistical method aiming to construct the 3D representative pore network from the measured 2D images of the porous media with similar pore structure statistics. Bryant and co-workers pioneered the use of geologically realistic networks based on random close packing of equally-sized spheres (Bryant and Blunt1992; Bryant et al. 1993a, b). (Bakke and Oren 1997; Lerdahl et al.2000; Øren and Bakke 2002; Oren et al.1998) sim-ulated the packing of spheres of different size based on the grain size distribution from analysis of thin sections of the rock of interest. The works of Pillotti, and Coehlo and co-workers have demonstrated a method to simulate the deposition of grains of non-spherical shapes, in particular, enabling the reconstruction technique to be applied more generally (Coelho et al.1997; Pilotti2000). Also, in order to reduce the computational time, reconstructing a repre-sentative pore network based on statistically information extracted from the pore-scale image to keep the pore-scale information is of interest to the up-scaling work.
Adding pore structure information based on pore-scale imaging into the statistical method might enhance the reconstruction process. Therefore, we proposed a hybrid method combining a heuristic/statistical method and lim-ited information provided by direct mapping method to reconstruct an equivalent pore-network incorporating information of real pore structure. This study applied a network-based porous media flow simulation (NET-PFS) model to explore pore geometry and analyze flows in
porous media. The model is capable of simulating a rep-resentative pore network and analyzing the flow status with limited information about the pore structure.
However, constructing a representative pore network by parameter calibration is computationally demanding. An effective search algorithm is therefore needed. Hence, the newly proposed KASA algorithm was developed to cali-brate the pore network of the target sample for the NET-PFS model.
This study makes three major contributions. First, an innovative KASA algorithm based on simulated annealing and Kriging approximation was developed to increase opti-mization searching efficiency. Second, this research built a representative pore network with limited information about the soil samples. A framework combining pore network simulation and KASA algorithm was developed to generate the equivalent pore network based on the pore-structure statistical information. Third, this study successfully con-ducted a case study of the KASA algorithm for calibrating the parameters of the equivalent pore network of a sand packing porous medium for NET-PFS model. The rest of the paper is organized as follows. The methodology of the KASA algorithm is developed in Sect.2. An application of the KASA algorithm for calibrating the NET-PFS model with a packed quartz sand medium is demonstrated in Sect.3. Results and discussion are presented in Sect.4. Finally, the conclusion is presented in Sect.5.
2 The KASA algorithm
The newly developed KASA algorithm is an efficient opti-mization algorithm especially designed to optimize problems with a complicated objective function. The KASA algorithm combines simulated annealing with Kriging interpolating approximation. More specifically, the KASA search proce-dures mainly follow the simulated annealing algorithm; in addition, the Kriging approximation method is utilized to approximate the evaluation of the complicated objective function during the running of the KASA algorithm. The detailed KASA algorithm is derived in the following, and a flowchart of the algorithm is proposed in Fig.1.
The search procedure of the KASA algorithm is similar to that of simulated annealing. The simulated annealing algorithm is a heuristic optimization algorithm that simu-lates the cooling process of heated metals. In the cooling process, while the metal temperature decreases, the metal structure gradually converges to a stable condition. Simi-larly, this optimization algorithm approaches the global optimum step by step with a probabilistic technique. At each iteration of the algorithm, the KASA algorithm takes the steps from the current solution to the new solution when the new solution has a better objective value than the
current solution. Suppose the problem seeks a minimal objective value. The new solution is better if it has a smaller objective value, i.e., objnew objcurrent; it would then be accepted without any exception. In contrast, if the new solution is worse, i.e., objnew[ objcurrent, the algorithm accepts the new solution on the basis of a probability dis-tribution function that is negatively proportional to the metal temperature. The algorithm generates a random number z with a uniform distribution in the interval of [0, 1]. The new solution is accepted only if z is less than or
equal to an exponential function
EXPO obj new objcurrent=temp in Eq. (1); otherwise, the new solution is rejected when Eq. (2) holds.
objnewis accepted
if z EXPO obj new objcurrent=temp ð1Þ objnewis rejected
if z [ EXPO objnew objcurrent
=temp
ð2Þ
Notice that the evaluation of the complicated objective function could be time consuming for a large-scale and complex model. Therefore, the newly developed KASA algorithm estimates the objective value by using Kriging Simulated Annealing Kriging approximaon method
Generate the inial soluon as the current soluon
Calculate the objecve value of the current soluon by using the
Kriging method
Generate a new soluon
If the new soluon is beer than the current soluon
If the new soluon is worse than the current soluon Calculate the objecve value of
the new soluon by using the Kriging method
Accept the new soluon as the current soluon
Accept the new soluon using the acceptance
probability funcon
Decline the new soluon using the acceptance
probability funcon
First stopping criteria
Approximaon to the global opmum Second stopping criteria
Calculate objecve value of new soluon by running NET-PFS model No
Yes
No
Yes
Kriging-approximaon Simulated Annealing (KASA) algorithm
Calculate the semivariogram to esmate the spaal dependence
Calculate the Kriging weights
Perform the Kriging approximaon for the objecve
funcon
Fig. 1 Flowchart of the Kriging-approximation simulated annealing (KASA) algorithm
approximation. The more approximation is used, the greater is the computational efficiency provided by the KASA algo-rithm. Kriging is a statistical method that is used for spatial interpolation; it is used to predict the best linear unbiased estimation (BLUE) of unknown random variables. It estimates an unknown variable at a given point i by calculating the weighted average of known values at given points j1 in the neighborhood where i¼ 1; 2; 3; . . .; m and j1 ¼ 1; 2; 3; . . .; n. vTi and vj1are the realization of random variables at point i and point j1, respectively. Suppose vT
i and vj1 have the same
expected values l and the same variance r2. The estimation vEi of random variable vTi at point i can be approximated by the weighted average of other known values, i.e. vE
i ¼
P
j1 wj1 vj1
h i
. Then expected estimation of vT i is
unbiased as long asPj1hwj1i¼ 1. Furthermore, the expected square error EXP vE
i vTi
2
h i
is computed in Eq. (3). To find the minimal expected square error, we need to minimize the Lagrange objective function LAGRiin Eq. (4). Hence, the partial derivatives of LAGRiwith respect to decision variables wkand k should be equal to zero. The optimal weight wj1and the optimal Lagrange multiplier k are determined to minimize the expected square error of the Kriging method. The Kriging method uses a semivariogram function to estimate the spatial dependence among points. The semivariogram cðdj1;j2Þ is
defined as half of the expected square difference between two field values at points j1 and j2. Notice that the semivariogram is a function that depends only on the distance between sam-ples. Since covariance COV½vj1; vj2 is equal to l2 cðdj1;j2Þ,
the optimality conditions becomes Eq. (5). Similarly, the Kriging weights and the Lagrange multiplier can be computed in matrix form of Eq. (6).
EXP ERRORi2¼ EXP½ vE i v T i 2 ¼ VAR ½ vE i v T i h i ¼X j1;j2 wj1 wj2 COV vj1; vj2 h i h i þ r2 2 X j1 wj1 COV vj1; v T i h i h i 8i ð3Þ MIN LAGRi¼X j1;j2 wj1 wj2 COV vh j1; vj2i h i þ r2 2 X j1 wj1 COV vj1; v T i h i h i þ 2 k X j1 wj1 h i 1 ! 8i ð4Þ X j1 wj1 c d k;j1 h i k ¼ c d k;i 8i; k ð5Þ w1 w2 . . . wn k 2 6 6 6 4 3 7 7 7 5¼ c d1;1 c d2;1 . . . c dn;1 1 c d1;2 c d2;2 . . . c dn;2 1 . . . . . . . . . . . . . . . c d1;n c d2;n . . . c dn;n 1 1 1 . . . 1 0 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 1 c d1;i c d2;i . . . c dn;i 1 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 8i ð6Þ This research incorporates a Kriging approximation method into the simulated annealing optimization algo-rithm. The KASA algorithm is constructed on the basis of the simulated annealing search mechanism, and, addition-ally, Kriging approximation is used to enhance the evalu-ation of the complicated objective function of the KASA algorithm. Since the calibration of the NET-PFS simulation model needs to be simulated for every optimization step, the KASA algorithm successfully increases the optimiza-tion efficiency by using the Kriging method to estimate the simulation results instead of running the NET-PFS model.
3 Case study
3.1 The NET-PFS model
The permeability of a porous medium can be estimated by simulate a steady-laminar single-phase flow in the porous medium. This research applied the NET-PFS model, which can be used to simulate a pore network and analyze porous media flows. Given the total pore volume and distribution of the pore size, the number, size, and location of pores are simulated by using Monte Carlo simulation. Then, a pore network is constructed by simulating a pore-throat con-nection. Assume the pore network has qa arcs (pore throats) and q1 nodes (pore spaces). Denote xpr1as a node pressure at node r1, where index r1¼ 1; 2; 3; . . .; q1. Fur-thermore, xfr1;r2 represents the arc flow delivered from r1 to r2; index r2 is also a node index that is the same as r1 (r2¼ 1; 2; 3; . . .; q1). In addition, nodes can be categorized into q3 inner nodes and q4 boundary nodes, i.e., q3þ q4 ¼ q1. In addition, the inner and boundary nodes are denoted by r3 and r4, respectively, where r3¼ 1; 2; 3; . . .; q3 and r4¼ 1; 2; 3; . . .; q4.
The governing equations of the NET-PFS model contain mass balance equations and Poiseuille’s equations. The details of the physical meanings and mathematical derivations can be found in Dullien (1992). The mass balance equation in Eq. (7) defines that, for each pore, the total inflow into the pore equals to total outflow from the pore. X r1 xfr1;r3X r1 xfr3;r1¼ 0 8r3 ¼ 1; 2; 3; . . .; q3 ð7Þ Poiseuille’s equation in Eq. (10) is a physical law that describes an incompressible fluid flowing in a pipe with a constant cross section under a pressure drop. In Eq. (8), the pair-wise pressure equations that describe the arc flow rate xfr1;r2 from nodes r1 to r2 are proportional to the node pressure drop, xpr1 xpr2, of those two nodes (see Fig.2).
xpr1 xpr2¼ 8lalr1;r2xfr1;r2 ½p ar= r1;r24
8r1 ¼ 1; 2; 3; . . .; q1; r2 ¼ 1; 2; 3; . . .; q1
ð8Þ
where l is the viscosity of the fluid, al is the length of the arc, ar is the radius of the arc.
A fixed pressure boundary condition is imposed on the boundary nodes as shown in Eq. (9).
xpr4¼ prr4 8r4 ¼ 1; 2; 3; . . .; q4 ð9Þ
The NET-PFS model solves the system of equations in Eqs. (7)–(9). The system of equations has (qaþ q3 þ q4) equations including q3 mass balance equations, qa Poi-seuille’s equations, and q4 boundary conditions. There are qaþ q3 þ q4 unknown variables including qa arc flow and q1 (¼ q3 þ q4) node pressure variables. By solving the system of equations of the NET-PFS model, the flow rate and pore pressure of each arc and node of the porous media are determined.
The single pair pore network structure (see Fig.2) can be extended to two or three-dimensional pore networks. As the dimension of the pore network increases, the number of the pores and throats exponentially increase as well as the size of the system of equations of the NET-PFS model. The NET-PFS model with m pores and n throats forms a system of equations with nþ m equations based on Eqs. (7)–(9). By solving the system of equations, the distributions of the
pressure in the pores and flow rate in the arc are determined as shown in the Fig.3.
The hydraulic conductivity K of the pore network is found from Darcy’s law, K¼lQtL
ADU, where l is the viscosity
of the fluid. The total flow rate through the network, Qt, is found by imposing a potential drop DU across its length L, with A being the cross-sectional area of the model. Potential is defined as U¼ P qgh, where P is the pres-sure, q is the fluid density, g is the gravitational constant and h is the height above datum.
3.2 Pore-scale imaging and model parameters
of a packed quartz sand medium
Since the number of control variables and equations of the NET-PFS model is large, we found that calibrating the model is computationally demanding. Hence, the case study used the NET-PFS model to simulate the flow rate, pore pressure, and permeability of a porous medium in a pore-network based on the pore structure of a porous medium with a sand-packing. The KASA algorithm was applied to seek the optimal indeterminate parameters of the NET-PFS model. The properties of the sand-packing soil were scanned and investigated by Computational X-ray microtomography (lCT), and, then, the NET-PFS model simulated the porous media flow in the pore network.
lCT analyses of sand samples were carried out using the GE Nanotom S device at the microtomography facility of the Institute of Agrophysics of the Polish Academy of Sciences (IAPAS) microtomography facility (Lamorski 2016). The voxel size achieved in the scan was 2 lm. The parameters of the X-ray source were: accelerating voltage of 90 kV and cathode current of 120 lA. A tungsten exit window was used for the X-ray generation. For the purpose of noise reduction, each of the two-dimensional (2D) images recorded during the scan was an average of 15
Fig. 2 Schematic of arc flow between two pores
Fig. 3 Arc flows and pressure distribution in a pore network. The colors in the nodes indicate the fluid pressure from high (red) to low (blue), and the arrows show the flow directions and magnitudes
images. These images were used for the reconstruction of the 3D representation of the sample.
During the scan, the temperature in the lCT chamber was stabilized, but the sample itself was heated by X-rays, which could cause thermal changes in the sample dimen-sions, which, in turn, could impact the reconstruction process in case of high-resolution scans. To overcome this problem, we performed a 30-min pre-scan prior to the main scan. During this pre-scan, the sample temperature was stabilized. After the pre-scan, the proper scan was started instantly.
After the scan, volume reconstruction and image pro-cessing were performed. 3D volume reconstruction was done using the DatosX 2.0 software. For the 3D image processing and visualization, the VG Studio Max software was used. The image processing procedure consisted of region-of-interest selection and 3D median filtering, with the kernel diameter equal to 3 pixels.
Figure4shows a 2D cross-sectional image of a packed quartz sand medium and a 2D pore network extracted from the selected section of the pore structure (red square). The reconstructed 3D image provides the basic information about the pore structure of the quartz sand medium; dis-tributions of the radii of the pores (354 pores) and throats (443 throats) based on the maximum ball method (Al-Kharusi and Blunt 2007) follow gamma distribution as shown in Fig.5. The dashed lines in Fig.5are the gamma distributions with the fitted parameters a and b.
The distributions of the radius of the pores and throats show that the pore structure of the packed quartz sand medium consists of pore spaces of a uniform radius with a few large pores/throats. However, the quality of the pore network of the quartz sand medium was highly dependent on the limitation of the resolution of the X-ray CT scanner and on the quality of the image. No pore throat/body with a radius smaller than 1 lm was identified. The parameters of the model based on the limited information about the pore
structure were also extracted from the micro X-ray CT. The hydraulic conductivity was directly measured by Darcy’s method.
Based on the probability distributions of the pore radius in Fig.5, pore soil ratio, arc length, arc radius, and con-nectivity ratio among the pores of a real sand-packing porous medium, the pore network structures of the NET-PFS model were constructed by Monte Carlo simulation. A simulation of the size and location of the pore body is shown in Fig. 6. The connectivity of the nodes was simu-lated by the k nearest neighbor method. Furthermore, a pore network was simulated and is plotted in Fig.7.
4 Results and discussion
In the case study, the NET-PFS model simulated the porous media flow in a 216 mm3 (¼ 6 mm 6 mm 6 mm) packed quartz sand medium. Figures5 and 6 show an example of a randomly generated pore network for an arc radius of r¼ 0:17 (mm) and a connecting distance of d ¼ 1:2 (mm). The generated pores were classified into three groups: isolated pores, boundary pores, and inner pores. In this example, there were 15 pores at the inlet boundary, 31 pores at the outlet boundary, and 303 inner pores. The volumes of the pores were generated on the basis of the pore volume distribution measured from the micro X-ray experiment. In the model, the constant pressure conditions at the inlet and outlet boundaries were xpin¼ 104ðPaÞ and
xpout ¼ 0 Pað Þ, respectively; dynamic viscosity was dv¼ 1 Pa sð Þ. With the randomly generated pore network, we simulated the distributions of the water pressures at the pore bodies and the flow rates at the throats. In this case, the NET-PFS model generated a pore network model with 354 pores (303 inner pores and 46 boundary pores) and 443 pore throats. Accordingly, the model had 746 variables
Fig. 4 A two-dimensional image of the micro X-ray images of the packed quartz sand medium, and a network structure extracted from the selected red square area
(303 pressure variables, xpr1, for the inner pores and 443 flow variables, xfr1;r2, for the pore throats) and 746 equa-tions [303 mass balance equaequa-tions in Eq. (9) and 443
Poiseuille’s equations in Eq. (10)]. By solving the system of equations, we can determine the pore pressure and throat flow rate. See Fig.8 for the distribution of the pore
Fig. 5 Distributions of a pore radius and b throat radius for the packed quartz sand. The dashed lines are the gamma
distributions with the fitted parameters a and b 0 1 2 3 4 5 6 0 2 4 6 0 1 2 3 4 5 6 Y (mm) X (mm) Z( m m ) 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
Fig. 6 Simulation of the pore size and location of the network by using the Monte Carlo method
0 1 2 3 4 5 6 0 2 4 6 0 1 2 3 4 5 6 Y (mm) X (mm) Z( m m ) 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
Fig. 7 Simulation of a pore network by using the Monte Carlo method
Fig. 8 Water pressures (pa) of the pore bodies
Fig. 9 Relationship between the simulation error, the throat radius, and the connecting distance
pressure. We also calculated the total flow rate, Qm, which
is the summation of flow rates in all pore throats of the pore structure.
By adjusting the value of the arc radius and the con-necting distance, we were able to plot the relationship between the simulation error, the throat radius, and the connecting distance. Figure9 shows the deviation (simu-lation error) between the modeled flow rates of the pore networks (based on different values of r and d) from the measured total flow rates. The pore throat radius and connecting distance were significant factors for the porous media flow simulation; a sensitivity analysis of the throat radius and connecting distance was conducted. The results showed that a larger radius and a shorter connecting dis-tance increased the total flow rate in the model. In Fig.10, the total flow rate increased as we increased the throat radius under the same pressure gradient since the viscous pressure drop decreased as the arc radius increased on the basis of Poiseuille’s law. When the connecting distance was increased, each pore connected to more pores, and the connectivity of the pore network also increased. As a consequence, the hydraulic conductivity increased. There-fore, on the basis of Darcy’s law, the total flow rate increased as we increased the connecting distance, as shown in Fig.11.
Furthermore, the computational efficiency of the newly developed KASA algorithm and that of the original simu-lated annealing algorithm were compared. An NET-PFS error minimization problem is formulated as follows. Equation (10) computes objective function of simulation error between estimated flow rate, Qm, and true flow rate,
Qr. The optimization model seeks minimal error by
adjusting arc radius of r and connecting distance of d; the objective function is subject to constraints of NET-PFS model. MINr;d Qm xpr1ðr; dÞ; xfr1;r2ðr; dÞjr; d Qr ð10Þ
Figure12 shows that both algorithms converged to optimal solutions. However, KASA can save a significant calculation time by using Kriging approximation instead of running the NET-PFS model. Then, the computation time and cost were reduced because of Kriging approximation. The KASA algorithm was able to provide a higher com-putational efficiency for the optimization procedures. Suppose the original running time of the model is tn, and Kriging approximation takes time tk. For every p iteration of the KASA algorithm, pð 1Þ Kriging approximation is performed and the NET-PFS model is run only once. Equation (11) shows that the reduction of the computa-tional time is approximately pð 1Þ tn tkð Þ.
p ð Þ tnð Þ
½ p 1½ð Þ tkð Þ þ tn ¼ p 1ð Þ tn tkð Þ ð11Þ For the NET-PFS simulation of the packed quartz sand medium, Fig.13shows that the KASA algorithm reduced the calculation time from 20,000 (s) to 7000 (s). The computational efficiency can be improved further if the simulation model is more complicated. In addition, a sen-sitivity analysis of the throat radius and connecting dis-tance in terms of computational time was also conducted. Since the throat radius does not change the numbers of control variables and equations of the NET-PFS model, Fig.14 shows that the sensitivity analysis of the throat radius had less impact on the computational time. How-ever, a larger connecting distance led to greater network complexity, which caused an increase in the computational time for the NET-PFS model, as shown in Fig.15.
5 Conclusions
In this research, we proposed a framework combining pore network simulation and KASA algorithm to construct the equivalent pore network based on the pore-structure R² = 0.9893 0 100 200 300 400 500 600 700 800 0 0.05 0.1 0.15 0.2 T o tal flow rate (m m 3/s) Throat radius (mm) data
trend line (Order 5)
Fig. 10 Sensitivity analysis of the pore throat radius in relation to total flow rate
y = -259.14x3+ 2092.1x2- 2911.7x + 1137.8 R² = 0.9765 0 200 400 600 800 1000 1200 1400 1600 1800 0 0.5 1 1.5 2 T o tal flow rate (m m 3/s) Connecting distance (mm) data
trend line (polynomial)
Fig. 11 Sensitivity analysis of the connecting distance in relation to total flow rate
statistical information from the X-Ray CT image. The case study showed that the proposed method can effectively find a representative pore network for a packed quartz sand
medium. The hydraulic conductivity of the representative pore network was the same as that of the complicated pore structure of the packed quartz sand medium.
Fig. 12 Optimization steps of the KASA algorithm and the original simulated annealing algorithm
KASA prepara
1 2 3 4 5 6 (hr)
0
Fig. 13 Calculation time of the KASA algorithm and the original simulated annealing algorithm y = 2.2643x + 89.352 R² = 0.0002 0 20 40 60 80 100 120 140 0 0.05 0.1 0.15 0.2 Com p utational tim e (sec) Throat radius (mm) data
trend line (linear)
Fig. 14 Relationship between throat radius and computational time (d = 1.2 mm) y = 0.0668e5.684x R² = 0.9718 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 0.5 1 1.5 2 Com p utational tim e (sec) Connecting distance (mm) data
trend line (exponential)
Fig. 15 Relationship between connecting distance and computational time (r = 0.15 mm)
In addition, this research developed a new KASA opti-mization algorithm. The KASA algorithm uses Kriging approximation to avoid the evaluation of the complicated objective function. The proposed KASA algorithm applies Kriging approximation to interpolate the objective value in the solution space rather than to estimate variables in physical space. The advantage of KASA is to reduce the evaluation time of the objective function. Similar to sim-ulated annealing, the KASA algorithm searches for the optimal solution by simulating the cooling process of heated metals. The KASA algorithm is also a heuristic optimization method and can reach the global optimum with a probabilistic technique. Notice that the proposed KASA algorithm interpolates the objective value in the solution space and the uncertainty estimation of Kriging approximation can be determined. The current KASA algorithm does not use uncertainty information in optimum searching procedure. Uncertainty estimation could be incorporated with solutions searching strategy of simulated annealing for future study. In this research, KASA was applied to calibrate the computationally demanding NET-PFS model. The results showed the computational effi-ciency of the KASA optimization algorithm.
Compared with those of previous studies, the significant contribution of this research includes the development of both the NET-PFS simulation model/framework and an innovative KASA optimization algorithm. Nevertheless, in this study, the control variables were limited to pore throat radius and connecting distance. Future studies should fur-ther consider the variations in pore body locations, pore throats, and arc radius. Moreover, in the study, the distri-butions of the pore volume and arc radius were given. In reality, micro X-ray images of the pore structure might actually be difficult to obtain. The inputs of the optimiza-tion framework should be extended to the informaoptimiza-tion that can be obtained with conventional soil analysis experi-ments without a micro X-ray image, i.e., pore size distri-bution. Then, the applicability of the optimization framework can be significantly extended. Determined representative pore networks can further be used to esti-mate the relative permeability, residual phase saturation, and other soil characteristics for single- or multi-phase flows in porous geo-materials.
Acknowledgements The authors thank the editors and anonymous referees for thoughtful comments and suggestions. The authors are responsible for the opinions and comments. This research was funded by the Taiwanese Ministry of Science and Technology
(MOST-105-2627-M-002-037; 105T612C502; MOST
106-2628-M-002-009-MY3) and National Taiwan University
(NTU-CCP-106R891007; NTU-CC-107L892607).
References
Al-Kharusi AS, Blunt MJ (2007) Network extraction from sandstone and carbonate pore space images. J Pet Sci Eng 56(4):219–231 Bakke S, Øren PE (1997) 3-D pore-scale modelling of sandstones and
flow simulations in the pore networks. Spe J 2(02):136–149 Bechler A, Romary T, Jeannee N, Desnoyers Y (2013) Geostatistical
sampling optimization of contaminated facilities. Stoch Environ Res Risk Assess 27:1967
Blunt MJ, Jackson MD, Piri M, Valvatne PH (2002) Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv Water Resour 25(8):1069–1089
Blunt MJ, Bijeljic B, Dong H, Gharbi O, Iglauer S, Mostaghimi P, Paluszny A, Pentlanda C (2013) Pore-scale imaging and modelling. Adv Water Resour 51:197–216
Bryant S, Blunt M (1992) Prediction of relative permeability in simple porous media. Phy Rev A 46(4):2004
Bryant SL, King PR, Mellor DW (1993a) Network model evaluation of permeability and spatial correlation in a real random sphere packing. Transp Porous Media 11(1):53–70
Bryant SL, Mellor DW, Cade CA (1993b) Physically representative network models of transport in porous media. AIChE J 39(3):387–396
Cao K, Ye X (2013) Coarse-grained parallel genetic algorithm applied to a vector based land use allocation optimization problem: the case study of Tongzhou Newtown, Beijing, China. Stoch Environ Res Risk Assess 27:1133
Cieniawski SE, Eheart JW, Ranjithan S (1995) Using genetic algorithms to solve a multiobjective groundwater monitoring problem. Water Resour Res 31(2):399–409
Coelho D, Thovert JF, Adler PM (1997) Geometrical and transport properties of random packings of spheres and aspherical particles. Phy Rev E 55(2):1959
de Lavenne A, Skøien JO, Cudennec C, Curie F, Moatar F (2016) Transferring measured discharge time series: large-scale com-parison of top-Kriging to geomorphology-based inverse model-ing. Water Resour Res 52(7):5555–5576
Dong H, Blunt MJ (2009) Pore-network extraction from micro-computerized-tomography images. Phy Rev E 80(3):036307 Dougherty DE, Marryott RA (1991) Optimal groundwater
manage-ment: 1. Simulated annealing. Water Resour Res
27(10):2493–2508
Dullien FA (2012) Porous media: fluid transport and pore structure. Academic Press
Jiang X, Lu W, Na J, Hou Z, Wang Y, Chi B (2018) A stochastic optimization model based on adaptive feedback correction process and surrogate model uncertainty for DNAPL-contami-nated groundwater remediation design. Stoch Environ Res Risk Assess 32:3195.https://doi.org/10.1007/s00477-018-1559-4
Keehm Y, Mukerji T, Nur A (2004) Permeability prediction from thin sections: 3D reconstruction and Lattice-Boltzmann flow simu-lation. Geophys Res Lett 31:L04606. https://doi.org/10.1029/ 2003GL018761
Ketcham RA, Carlson WD (2001) Acquisition, optimization and interpretation of X-ray computed tomographic imagery: appli-cations to the geosciences. Comput Geosci 27(4):381–400 Laaha G, Skøien JO, Blo¨schl G (2014) Spatial prediction on river
networks: comparison of top-Kriging with regional regression. Hydrol Process 28(2):315–324
Lamorski K (2016) X-ray computational tomography facility. Insti-tute of Agrophysics PAS. http://tomography.ipan.lublin.pl. Accessed 18 Sept 2016
Lerdahl TR, Oren PE, Bakke S (2000) A predictive network model for three-phase flow in porous media. In: SPE/DOE improved oil recovery symposium. Society of Petroleum Engineers
Marryott RA, Dougherty DE, Stollar RL (1993) Optimal groundwater management: 2. Application of simulated annealing to a field-scale contamination site. Water Resour Res 29(4):847–860 Miao X, Gerke KM, Sizonenko TO (2017) A new way to
parame-terize hydraulic conductances of pore elements: A step towards creating pore-networks without pore shape simplifications. Adv Water Resour 105:162–172
Okabe H, Blunt MJ (2004) Prediction of permeability for porous media reconstructed using multiple-point statistics. Phy Rev E 70(6):066135
Oliver MA, Webster R (1990) Kriging: a method of interpolation for geographical information systems. Int J Geogr Inf Sci 4(3):313–332
Øren PE, Bakke S (2002) Process based reconstruction of sandstones and prediction of transport properties. Transp Porous Media 46(2–3):311–343
Oren P-E, Bakke S, Arntzen OJ (1998) Extending predictive capabilities to network models. SPE J 3(04):324–336
Pilotti M (2000) Reconstruction of clastic porous media. Transp Porous Media 41(3):359–364
Piri M, Blunt MJ (2005) Three-dimensional mixed-wet random pore-scale network modeling of two- and three-phase flow in porous media. I. Model description. Phys Rev E 71(2):026301
Raeini AQ, Bijeljic B, Blunt MJ (2015) Modelling capillary trapping using finite-volume simulation of two-phase flow directly on micro-CT images. Adv Water Resour 83:102–110
Steefel CI, Beckingham LE, Landrot G (2015) Micro-continuum approaches for modeling pore-scale geochemical processes. Rev Mineral Geochem 80:217–246
Verdin A, Rajagopalan B, Kleiber W, Funk C (2015) A Bayesian Kriging approach for blending satellite and ground precipitation observations. Water Resour Res 51(2):908–921
Xiong Q, Baychev TG, Jivkov AP (2016) Review of pore network modelling of porous media: experimental characterisations, network constructions and applications to reactive transport. J Contam Hydrol 192:101–117
Yeh WWG (1986) Review of parameter identification procedures in groundwater hydrology: the inverse problem. Water Resour Res 22(2):95–108
Yeong C, Torquato S (1998) Reconstructing random media. Phys Rev E 57(1):495
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
