Abstract
Characterizing internal microscopic structures of porous media is of vital importance to simulate fluid and electric current flow. Compared to traditional rock mechanics and geophysical experiments, digital core and pore network modeling is attracting more interests as it can provide more details on rock microstructure with much less time needed. The axis extraction algorithm, which has been widely applied for pore network modeling, mainly consists of a reduction and burning algorithm. However, the commonly used methods in an axis extraction algorithm have the disadvantages of complex judgment conditions and relatively low operating efficiency, thus losing the practicality in application to large-scale pore structure simulation. In this paper, the updated algorithm proposed by Palágyi and Kuba was used to perform digital core and pore network modeling. Firstly, digital core was reconstructed by using the Markov Chain Monte Carlo (MCMC) method based on the binary images of a rock cutting plane taken from heavy oil reservoir sandstone. The digital core accuracy was verified by comparing porosity and autocorrelation function. Then, we extracted the central axis of the digital pore space and characterize structural parameters through geometric transformation technology and maximal sphere method. The obtained geometric parameters were further assigned to the corresponding nodes of pore and throat on the central axis of the constructed model. Moreover, the accuracy of the new developed pore network model was measured by comparing pore/throat parameters, curves of mercury injection, and oil-water relative permeability. The modeling results showed that the new developed method is generally effective for digital core and pore network simulation. Meanwhile, the more homogeneity of the rock, which means the stronger “representative” of binary map the rock cutting plane, the more accurate simulated results can be obtained.
1. Introduction
Porous media, such as metal, wood, soil, and rock, is one of the most common substances in daily life [1–3]. The common character of all porous materials is that the internal pore structure usually governs the flow of fluid and electric current inside the media [4–6]. Even for the two materials with equal porosity, the fluid flow behaviors may be still different because of the different spatial distribution of internal pore structures, the connection types, and the shape or size of the pore and throat [7–9]. Therefore, in order to comprehensively study porous media, characterizing pore network structures is of significant importance. This is typically true for petroleum/reservoir engineering, as pore network structures determine the capacity of hydrocarbon accumulation and also the flow behavior.
As we are approaching the end of conventional hydrocarbon resources, a giant potential of unconventional resources of energy has been manifested. However, there are difficulties with unconventional resources in exploitation, such as the very low permeability of tight oil/gas reservoir and the high viscosity of heavy oil in heavy oil reservoir. So much attention and research have been paid to heavy oil reservoirs, as it accounts for a large proportion in oil and gas resources. Since reservoir macroscopic properties (i.e., permeability and capillary pressure) are controlled by its microstructure, especially for heavy oil reservoirs where the oil viscosity is so high, it is necessary to examine the characteristics of pore networks at microscale and eventually improve oil and gas recovery. Thus, in order to provide efficient production from heavy oil reservoirs, sufficiently reliable information about the pore structure of the hosting rock is required.
In fact, due to specially designed experimental setups and complicated operation procedures, traditional rock mechanics and geophysical experiments are usually difficult and time-consuming methods to extract pore network [10]. Meanwhile, the rapid development of computer technology and the explosive growth of computing power make it possible to perform large-scale simulations using rock digital geometric data [11, 12]. The digital core modeling can provide detailed information on rock internal pore network, and the simulation results can be also compared with laboratory observations to comprehensively understand rock microstructures.
The current computational study method focuses on rock microstructures mainly and consists of (i) digital core model and (ii) pore network model. For the digital core model, it can accurately capture the pore network properties based on the numerical reconstruction of data obtained from rock geophysics tests. Specifically, the reconstructed digital core using data from X-ray CT scanning could present almost the same statistical characteristics as the real core. Therefore, the modeling results via this type of digital core can truly reflect the rock microscopic properties. Given the natural features for capturing all detailed microstructures, the reconstructed pore networks by X-ray CT scanning are very complex. This raises some difficulties for core-scale numerical simulation (i.e., seepage simulation and acoustic and electrical characteristic simulation), as it needs very high computation to reflect every detail of microstructures while has very limited contribution for accuracy improvement [13–15]. Therefore, urgent needs are calling to develop a new method which can retain the topological geometry properties of real core pore space and with far less computing power to improve modeling efficiency.
The pore network model is one of such promising models that is efficient and accurate [16, 17]. By reasonably simplifying the microscopic details of pore space, the pore network model can achieve similar topological and geometric structures as digital core while requiring much less computing power [18–20]. Since the first two-dimensional pore network model introduced by Fatt [21–23] in 1956, a wide range of pore network models have been developed, which are mainly divided into two categories: (1) regular topological pore network model and (2) real topological pore network model. The regular topological pore network model consists of pores and throats that are regularly arranged in plane or space. By assigning different values of pore/throat shape and size, this model could present a certain degree of heterogeneity. However, the uniform distributions of pores and throats in the regular topological network cannot represent the high complexity of pore structures within the real rock, which constrains the application in numerical simulation. On the other hand, the real topological pore network model is developed on the basis of digital core. Therefore, it has the equivalent topological structure of digital/real pore space and is more accurate for pore network modeling. Previously, a variety of methods were proposed to construct this type of model. For example, Zhao et al. [24] proposed a multidirectional slice scanning method to characterize the values of internal pore/throat. They defined the located in situ minimum value as throat while pores were hardly detected. Lindquist et al. [25], Sheppard et al. [26], and Prodanovic et al. [27] used the pore central axis method to characterize the topological structure of pore space, as the central axis can accurately represent the characteristics of topological networks. The internal pores were defined as the nodes on the central axis, and the nodes with minimum area were defined as the throat. Later, Delerue and Perrier [28] developed a method that can build Voronoi polyhedron in any type of pore space and subsequently form pore networks based on the established polyhedron. This method was applied by Okabe and Blunt [29] to establish pore networks for Berea sandstones, but they found that the generated topological networks were poorly structured and concluded that the Voronoi polyhedron method is not a suitable tool for the pore network modeling of digital core. Silin et al. [30] established a pore network model on Fontainebleau sandstones by using the maximum sphere method, and they calculated the mercury injection capillary pressure (MICP). They reported that although the developed model with the maximum sphere method seems to be reasonable, the connectivity number (the number of throats connected with a pore) of pore/throat was still relatively higher than the actual value.
It is worth mentioning that the combination of the central axis method and maximum sphere method has been widely applied to characterize pore space topological structure and identify pores and throats. The central axis can be obtained through thinning algorithm [31, 32] and burning algorithm [33]. However, the complicated judgment conditions of these algorithms usually present relatively low running efficiency. Meanwhile, the central axis extraction algorithm proposed by Palágyi and Kuba [34] only requires to judge the neighborhood relationship of pore voxel, so that it dramatically reduces the system calculation amount. In the discrete system, a voxel has no more than 6, 18, and 26 neighboring voxels, when the distance between neighboring voxels is defined as no more than 1, , and unit lengths. Therefore, when extracting the central axis with the algorithm proposed by Palágyi and Kuba, only 27 neighboring voxels at most rather than all points or many points of local area need to be analyzed in one iteration. Furthermore, six deletion templates in the algorithm simplify the analysis of neighboring voxels, thus making this central axis method computationally effective. While this method has been widely used in the study of heart and cerebral vessels [35–37] and pulmonary vessels [38, 39] of human, early embryonic mouse heart [40], and plants [41] and shows satisfactory results, to the best of our knowledge, it has not ever been applied petroleum engineering for core-scale pore network simulation. Therefore, in this paper, we try to use this central axis extraction algorithm to extract the pore topological structure of reservoir rock and compare the modeling results with the curves of MICP and oil-water relative permeability. Firstly, four images of binary cutting plane from four individual sandstones were used for the reconstruction of digital core with the MCMC (Markov Chain Monte Carlo) method [42]. Then, pore axis networks were extracted by applying the refinement algorithm proposed by Palágyi and Kuba. Finally, pore networks were obtained by adding the identified pore throat information on the pore central axis network. In order to construct the digital core and pore network model, Visual Basic 6.0 was used.
2. Methodology
2.1. Reconstruction of Digital Core
Markov chain describes a state sequence that has lack-of-memory or memoryless properties, where the value of each state depends on the previous finite states and is independent of other states. The probability of this state is called the transfer probability. The calculating transfer probability is the core in the application of Markov chain, but the calculation procedures are very complex. To address this issue, Wu et al. [43] introduced the neighborhood concept, namely, that the state of a certain point only depends on the states of several nearby points. In this work, we use a traversal scanning algorithm (based on the MCMC method) to generate sufficient image size to ensure the formed digital core has the same statistical characteristics with the original image. In general, the MCMC method requires three different two-dimensional images as input files to build digital cores on the corresponding three coordinate axis planes. But when the core presents high homogeneity or lacks images of good quality, the identical image can be also used to construct a 3D digital core.
The procedure’s steps of 3D digital core construction through the MCMC method are listed as follows (Figures 1 and 2); the -neighborhood, which consists of one voxel being traversed, traversed voxels, and one untraversed voxel adjacent (component of partial -neighborhoods) to the voxel being traversed, denotes that the number of voxels involved in one iteration of simulation is . (1)Use calculated porosity from core scanning images as the conditional probability of the first voxel state(2)Generate voxels on the first row of the first layer along axial direction. The second voxel state in the first row is simulated by the 2-neighborhood, and the next voxel is simulated by 3-neighborhood. Their condition probabilities are derived from the 6-neighborhood of core scanning image in plane(3)Repeat step (2) to form voxels in axial direction. The edge voxels are simulated by 3-neighborhood and 4-neighborhood, whereas the inner voxels are simulated by 5-neighborhood and 6-neighborhood. Their condition probabilities are also derived from the 6-neighborhood of core scanning image in plane(4)Repeat step (3) to form voxels in axial direction. From the second row of the second layer, the edge voxels are simulated by 9-neighborhood and 10-neighborhood, and the inner voxels are simulated by 14-neighborhood and 15-neighborhood. Their condition probabilities are obtained from the combination of 6-neighborhood systems of three planar core scanning images
2.2. Correction of Digital Core
Considering the intrinsic properties of the image segmentation process and reconstruction algorithm, isolated skeleton and pores often exist in the reconstructed digital core [44]. The isolated skeleton needs to be removed as it would form redundant pores and faulty connecting channels, which actually do not exist. For the isolated pores, they may exist in real cases, but the contribution of these pores to fluid flow is negligible. Besides, these isolated pores could affect the simulated results of core pore topological structure, the calculation of pore structure parameters, and the simulation of microscopic seepage. Therefore, these pores also need to be corrected.
In this paper, we used the Seed Filling Method to delete isolated skeleton and unreasonable isolated pores. The adjacent skeleton voxels or pore voxels are firstly put into a set . Then, we count the number of elements in set and compare the number with threshold to determine whether the set belongs to an isolated skeleton voxel set or an isolated pore voxel set. The detailed procedures of removal of isolated pores (similar steps can be used for isolated skeleton correction) are explained as the following: (1)Set the threshold (the value of varies with target cores)(2)Travers all pore voxels of digital core until finding the pore voxel A without a label(3)Take A as seed and assign it with a label. Then, push all pore voxels of the seed but without a label, back to the stack(4)Pop up the top element B, assign it with the same label, and store it in the set . Then, put all pore voxels in the 26-neighborhood of element B but without the label in the set (5)If the number of elements in set is greater than the threshold , go back to step (2). Otherwise, repeat step (4) until the stack is empty. Then, a connected pore space containing A in the digital core can be found, and the pore voxels in this area are all marked with the label(6)If the number of elements in set is less than the threshold , delete all the pore voxels (converted into skeleton voxels) belonging to the set . Otherwise, return to step (1) until all pore voxels are marked with the label
2.3. Extraction of Pore Axis
The theory and criterion proposed by Bakken and Eliassen [45], Kong and Rosenfeld [46], and Ma [47] lay the foundation for the pore axis extraction. Later, Palágyi and Kuba [34] developed a refinement algorithm based on delete template. Besides, Wildenschild and Sheppard [32] developed a fully parallel refinement algorithm which is on the basis of theory of Ma [47], and Lee et al. [48] proposed another refinement algorithm (the LKC algorithm) with consideration of the octree algorithm. Among these extraction algorithms, the method developed by Palágyi and Kuba only requires to consider the neighborhood relationship information, which is more efficient than other algorithms.
The core concept of this algorithm is to transform the black points (pore voxels) that match deletion conditions into white points (skeleton voxels), while keeping other original white points the same. Whether the black points match the deletion conditions or not is mainly determined by the definition of connectivity in a discrete medium. In a 3D discrete system, voxels of different phases (skeleton and pore in this work) should have maximum and minimum number of connected neighboring voxels in the same phase correspondingly to avoid situations where voxels of both phases are disconnected or connected. When voxels of one phase have the maximum number of connected voxels, i.e., 26, they may connect with each other in the form of a point or a line, which do not occur between solid-phase voxels. Considering rock skeleton is solid and the simple point described by Bakken and Eliassen [45], in this work, we define that the pore voxel is connected with 26 nearby pore voxels at most and the rock skeleton voxel is connected with 6 nearby skeleton voxels at most. Therefore, the condition to delete black point depends on the relationship in its 26-neighborhood voxels. Palágyi and Kuba [34] introduced six deletion templates to judge whether the black point should be deleted or not, and these templates actually describe six types of 26-neighborhood relations. These templates vary with deletion direction (, , , , , and ), whereas the deletion direction is determined by the 6-neighborhood relationship. For example, if the black point lacks of 6-neighborhood neighbors in the direction, then, the deletion direction is . Figure 3 shows the deletion templates M1 to M6 in the direction. M1 to M6 combined with the rotation around the direction axis (angles of rotation are 90°, 180°, and 270°) form the final deletion templates along the direction. The meaning of each symbol in the figure is as follows: at least one “★” is a black point, “1” is a black point, “0” is a white point, and “♦” is a black or white point.
The detailed procedures for black points processing based on the aforementioned deletion principle are shown as follows: (1)Import binary image of the digital core into the voxel set (2)Select a black point from set to judge whether it matches deletion condition along the direction. If satisfies one or more neighborhood relations along the direction, then delete ; otherwise, keep . Repeat this process until all black points are traversed(3)Change the deletion direction and select a black point to judge whether it matches the deletion condition in that direction. If satisfies one or more neighborhood relations in this direction, then delete ; otherwise, keep . Repeat this process until all black points are traversed(4)Repeat step (2) until all black points are traversed(5)Repeat step (3) until the deletion process traversed along all directions(6)If the black points in set do not change anymore, the deletion process ends, and all black points in set are the central axis voxels; otherwise, repeat step (2) until black points in set do not change (no more black points need to be deleted)
2.4. Optimization of Pore Axis
The aforementioned refinement algorithm would extract the central axis of pores very well, but there are some unreasonable outcomes, caused by surface noise, limit size of digital core, and complex pore structure, that could still exist. These unwanted outcomes are namely unreasonable short branches, redundant branches on the boundary, and multiple central axis nodes in a single pore. The existing of unreasonable microstructures cannot accurately reflect the topological structure of pore space. Therefore, we need to further modify these unreasonable microstructures. In this paper, the processing method proposed by Jiang et al. [49] was applied for pore axis correction, which is explained as the following: (1)Delete short branch on the central axis: given the influence of surface noise, some unreasonable short branches may appear in the upper part of the pore (as shown in Figure 4(a), the black point represents the pore point, and the red point represents the central axis point). It is necessary to delete these short-axis branches where the number of black points is less than a certain threshold or the pore node radius(2)Correct the boundary central axis: the size limitation of digital or real core could lead to extra generation of unreasonable boundary pores and boundary central axis, as shown in Figure 4(b). To eliminate the effect of undesired central axis, we first add several layers that have the same properties of the outermost layer of digital cores, then extract the central axis, and remove the added parts(3)Correct central axis node: Figure 4(c) shows the complex pore structure could trigger the multiexisting of several nodes in one pore. To ensure the accurate and reasonable corresponding relationship between pore node and core pore, it is necessary to appropriately merge multiple nodes which exist in the same pore
(a)
(b)
(c)
2.5. Construction of the Pore Network Model
To properly develop the pore network model, the key step is to assign the geometric parameters of the pore and throat nodes on the central axis at the corresponding position; thus, the central axis network of the pore space is transformed into the pore network model. During this construction process, the pore space geometric parameters need to be correctly identified.
The maximum sphere method is one of the most commonly used methods to analyze pore and throat parameters. However, the length of the segmented pore is usually way too large and affects the accuracy of parameters of other pores and throats. Therefore, in this paper, we applied a geometric transformation method [50] combined with the maximum sphere method [51] to segment pore space. Meanwhile, the OTSU method [52], a method to obtain the threshold value leading to the maximum between-cluster variance, was used to determine the pore/throat length and throat shape factor.
To effectively determine the pore length, geometric transformation has been widely used to cut pore space in a certain angle at the pore point (central axis node) to form a series of cutting planes, as shown in Figure 5. Then, on all cut planes, rays are emitted from the pore point at a certain angle. As Figure 6(a) shows, the emitted ray would keep extending until reaching the skeleton voxel and then record the length of all line segments. These line segments can form a single set, and the pore length can be determined by the OTSU method using optimal segmentation.
(a)
(b)
(c)
(d)
The pore radius can be determined by the maximum sphere method. Taking pore point as the sphere center, the radius of the sphere is kept increasing until reaching the skeleton voxel and forms an inscribed sphere. Therefore, the pore radius is equivalent to the radius of the largest sphere, as shown in Figure 6(b).
Figure 6(c) shows that throat length is given by the difference of the distance between adjacent pore points and their respective pore radius and , i.e., . The throat radius can be also determined through the maximum sphere method which is similar to pore radius estimation. The difference is that the throat radius is the minimum value of the maximum sphere radius within the throat length.
The last parameter is the throat shape factor. We first get several cutting planes of the throat along different directions using image transformation, as shown in Figure 6(d). The perimeter of the cutting planes can be obtained by counting boundary voxels, and the corresponding area can be obtained by counting all voxels on the plane. The shape factor of a single cutting plane can be then calculated via the shape factor calculation formula. Consequently, the shape factor of the throat can be calculated by counting obtained shape factors of all cutting planes using the OTSU method.
Finally, we construct the pore network model by assigning the aforementioned pore structure parameters to the corresponding position on the central axis of the pore space.
3. Validation
As aforementioned, although the digital core reconstructed by CT scanning has almost the same pore structure characteristics as the real core, the expensive cost and long consuming time still limit its application. On the other hand, the numerical reconstruction based on series of pore structure parameters or cutting plane pictures can effectively avoid these issues while presenting satisfactory running results.
In this paper, four individual core samples of heavy oil reservoir sandstones (named as sandstone a, sandstone b, sandstone c, and sandstone d) were used in X-ray CT scanning experiment. The lab equipment used for CT scanning is the Xradia MicroXCT-400 CT scanner. Basic parameters of the five cores of the CT scanning are shown in Table 1, and the binarized scanning results of four core samples are shown in Figure 7.
(a)
(b)
(c)
(d)
The initial input files for the digital core modeling method used in this work are four binary images (as shown in Figures 8(a)–8(d)) from the CT scanning result of sandstone a, sandstone b, sandstone c, and sandstone d. And all the parameters of the control group (sandstones a, b, c, and d) of digital cores reconstructed later are all processed from the CT scanning results.
Figures 9(a), 10(a), 11(a), and 12(a) show the reconstructed digital cores with the size of voxel3 using the MCMC method, i.e., digital core a, digital core b, digital core c, and digital core d, respectively. Given the intrinsic properties of the image segmentation process and reconstruction algorithm, isolated skeleton and pores may exist in the reconstructed digital core. We then use the Seed Filling Method to remove these unreasonable isolated skeletons isolated pores (the modified sections are indicated by blue circles in Figures 9(b), 10(b), 11(b), and 12(b)).
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
In fact, the removal of redundant isolated skeleton pore can make digital core more concise and realistic. The corrected skeleton voxels and pore voxels can transform into each other. Besides, the proportion of the modified skeleton/pore voxel points in the whole digital core is very low so that the porosity before and after correction should roughly remain the same value. For example, the porosity of sample A of the real core, digital core, and corrected digital core is 26.90%, 27.70%, and 28.10%, respectively, which only gives 0.4% difference before and after correction (Table 2). If only referring to porosity data, the reconstructed and corrected digital core are in line with the original core. However, as said, the internal microstructures of two porous media with exactly the same porosity could still behave differently. To further examine the conformity between constructed core and real core, we used autocorrelation function to characterize pore space properties. The autocorrelation function [53] is widely used to evaluate image structural properties. It represents the probability that any two points in the image are in the same phase (i.e., skeleton or pore):
where is the value of autocorrelation function, is the coordinate of point , is the observation value at the coordinate , is the observation value at coordinate , and is the distance between two observation coordinates in μm.
Figure 13 shows the calculated autocorrelation function varies with . It can be clearly seen that the variation trends of autocorrelation functions of reconstructed cores and original cores for all sandstones are very similar, indicating a high similarity of pore structure characteristics between reconstructed cores and real cores. With consideration of similar porosity of real and digital cores, it can be concluded that the digital core reconstructed in this paper has high accuracy compared to raw core, and the reconstruction algorithm and correction process are trustworthy.
(a)
(b)
(c)
(d)
The next step is to construct the final pore network model. We first extract the central axis of pore space based on the reconstructed digital core and identify corresponding geometric parameters using geometric transformation technology, maximum sphere method, and OTSU method. Then, the desired pore network model can be obtained once assigning geometric parameters to the corresponding position of central axis node. Figures 14(a) and 14(b) show the central axis network of digital core a before and after correction, respectively.
(a)
(b)
Figures 15(a)–15(d) show the final pore network models for all candidate cores. It can be qualitatively seen that core a has relatively larger pore size compared to the rest of the cores. For core b, a small portion of pores is large, but in general, they seem to be smaller and more densely distributed than core a. Core c has the worst uniformity of pore distribution, where the pore size is generally small but large pores could be locally present. The average pore size of core d is the smallest among all samples, and it is in fine distribution. Considering the high match of aforementioned digital core characteristics with the properties of the raw cutting plane binary map, we believe that the extracted pore network model can accurately reflect the microstructure of the real core.
(a)
(b)
(c)
(d)
To further validate the accuracy of the extracted pore network, we compare microstructure parameters of the pore network model with the same parameters of the original core, including throat radius, throat length, pore radius, and connectivity number, as shown in Figures 16–19. It is worth noting that the biggest difference among all microstructure properties between digital cores and real core is the connectivity number. The main reason is that compared to the rest of parameters, the connectivity number is spatially more sensitive, so that the three-dimensional space position changes would significantly affect the value of connectivity number. Since the initial input file for reconstruction modeling is a two-dimensional binary cutting plane image, the information contained in this image will largely determine the microstructure parameters of reconstructed digital core and pore networks. Given one two-dimensional image cannot completely and accurately reflect the pore throat coordination number which is significantly affected by the spatial position, the obtained coordination numbers of digital cores are therefore slightly lower than the original cores’ values.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
For core a and d, the throat radius, throat length, and pore radius are almost identical; for core b, the throat length of digital core is slightly concentrated at relatively smaller sizes than sandstone, while throat radius and pore radius are almost identical; for core c, there is the biggest difference in the length of throat among these 4 cores, while throat radius and pore radius are almost identical. In general, except for coordination number, the other three parameters of digital core are almost in line with the real core. Among them, pore and throat radius have higher consistency than throat length. Particularly the throat length, core c shows the largest difference between digital core and real core, which is followed by core b. The rest two cores have very well consistency of throat lengths between digital core and real core.
From the above results, we can observe that the parameters of reconstructed digital core, such as throat radius, throat length, pore radius, and connectivity number, are generally in line with the real core. Therefore, we conclude that the modeling method used can capture the same microstructure properties from real core to reconstructed digital core or/and pore network model. On the other hand, the main purpose of digital core and pore network model construction is to replace the physical experiment of the real core with numerical simulation. However, the accuracy of the established digital core and pore network model is yet to be totally assured of this replacement. To achieve this, the digital core and pore network model should have the same micropore structures and internal flow behaviors as real core. Therefore, to further validate the proposed modeling method, we performed simulations on mercury injection and oil-water two-phase flow based on the extracted pore network model.
The basis of flow (both single-phase and multiphase) simulations in the pore network model is the Hagen–Poiseuille model for hydraulic conductance between two pores: where is the flux between pores connected via throat, is the fluid conductance of the pore-throat-pore, is the length of path between pore body centers, and is the pressure difference between pores. The computation of the hydraulic conductance for any 2D cross-section (or fluid menisci for multiphase flow) is based on the dimensionless conductance model of Patzek and Silin [54]: where is the area of a cross-section, is the dimensionless hydraulic conductance (rigorously established for circle, square, and triangle within a circle-triangle-square model), and is the fluid viscosity.
As shown in Figures 20–23, it can be observed that the fluid flow behaviors in digital cores and pore network models are basically consistent with the corresponding real cores by comparing the mercury injection curves and oil-water relative permeability of digital cores and real cores. Mercury injection curves show better consistency than the relative permeability curves of oil-water two-phase flow, which could be explained by the more complex force in multiphase flow than that in single-phase flow. To be more specific, the sequence of consistency of relative permeability curves between digital cores and real cores is , which is in line with the results from pore and throat microstructure parameters. In other words, when the other parameters show high consistency (except connectivity number), the throat length of digital core b is slightly less consistent with real core b and worse for core c compared to core a and core d. Nevertheless, we observed high accuracy between digital cores and real cores. Therefore, we can conclude that the method used for digital core/pore network reconstructions in this paper can fully reflect the microstructure characteristics of real core, which gives convenience and paves the road for large-scale fluid flow simulation in oil/gas reservoirs.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
4. Conclusions
In this work, we developed a new method for digital core and core network modeling based on the algorithm proposed by Palágyi and Kuba and also the Markov Chain Monte Carlo (MCMC) method. (1)The modeling results show that the digital core of heavy oil reservoir sandstone reconstructed with this digital core modeling method can accurately reflect the microstructure characteristics of the real core. Moreover, the accuracy could be further improved when three representative images on the three vertical coordinate planes are available in the digital core modeling(2)Considering the low initial input requirement, which is one binary image of the core cutting plane for the least and the image can be easily obtained from X-ray CT scanning, casting thin section imaging or scanning electron microscope, etc., and high accuracy of simulation results, the new proposed method for digital core/pore network modeling in rock engineering should have wide applied future(3)As the microstructure properties of core cutting plane binary image, the only input in the modeling directly governs the corresponding pore network characteristics of the reconstructed digital core/pore network model. Therefore, the quality of input image highly determines the accuracy of reconstruction results. In general, the more homogeneous the rock behaves, the more precise the produced binary map section will be and the modeling results will be more accurate. As in this paper, the sequence of homogeneity of sandstones is , so the reconstructed digital core a and d presents higher accuracy than digital core b and c(4)In this work, we extracted the pore network model based on the obtained digital core. The results of pore/throat microstructure parameters and flow behavior simulation between extracted pore network and real core show a high consistency. Therefore, the modeling method of pore network in rock engineering presented in this paper is accurate and can be used to extract the pore network model of any other digital core
Data Availability
The data used to support the finding of this study are included within the article.
Conflicts of Interest
The author declares no conflict of interest regarding the publication of this manuscript.
Acknowledgments
This research was funded by the National Natural Science Foundation of China (No. 51874249), the National Key Research and Development Program of China (No. 2019YFC0312304-4), and the National Science and Technology Major Project of China (2016ZX05015-002) Open Fund (PLN201914) of State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation (Southwest Petroleum University).