Abstract
High heterogeneity and nonuniformly distributed multiscale pore systems are two characteristics of the unconventional reservoirs, which lead to very complex transport mechanisms. Limited by inadequate computational capability and imaging field of view, flow simulation cannot be directly performed on complex pore structures. The traditional methods usually coarsen the grid to reduce the computational load but will lead to the missing microstructure information and inaccurate simulation results. To develop a better understanding of flow properties in unconventional reservoirs, this study proposed a new upscaling method integrated gray lattice Boltzmann method (GLBM) and pore network model (PNM), accounting for the fluid flow in heterogeneous porous media. This method can reasonably reduce the computational loads while preserving certain micropore characteristics. Verifications are conducted by comparing the simulation and experimental results on tight sandstones, and good agreements are achieved. The proposed method is proven to be capable of estimating bulk properties in highly heterogenous unconventional reservoirs. This method could contribute to the development of multiscale pore structure characterizations and enhance the understandings of fluid flow mechanisms in unconventional reservoirs.
1. Introduction
With the advanced technology and growing demand for oil and gas, increasing attention has been paid to the development of unconventional oil and gas resources [1]. Unconventional reservoirs have multiple types of pore spaces, complex pore structures, and various mineral compositions and usually feature a high degree of heterogeneity and multiscale characteristics [2–4]. The flow mechanisms in unconventional reservoirs are greatly affected by microscale effects, and it is very challenging to describe the complex pore structures accurately and transmission characteristics of such porous media [5].
At present, there are two main methods to obtain and describe the complex pore structure and morphology of unconventional reservoirs: laboratory methods and high-precision imaging analysis methods. High-precision imaging analysis methods include the methods of scanning electron microscopy (SEM), micro-nano CT X-ray scanning, and dual-beam scanning electron microscopy (FIB-SEM) [6, 7]. However, due to the limitation of the imaging field of view, the entire high-precision pore space of the core cannot be obtained by one imaging; therefore, the multiscale characterization of the flow characteristics should fully consider the interaction between different scale pores in the heterogeneity reservoirs.
Microsimulation methods can accurately describe microscale fluid flow mechanisms. However, it is difficult to consider microscale flow effects in core-scale flow simulation analysis, limited by the method applicability and computational power. Worldwide, scholars have developed simulation methods for highly heterogeneous rocks at different scales, such as the pore network model method, lattice Boltzmann method (LBM), and molecular simulation method [8]. With topology of the pore space, the three-dimensional pore network model can not only characterize the pore structure parameters (pore radius, roaring channel radius, and coordination number) [9, 10] but also provide the basis of research for microscopic flow simulation of porous media. The 3D PNM can reduce the cost of experiments, shorten experimental data obtaining time, and get experimental data that are difficult to measure in the laboratory.
In recent years, the simulation methods of pore network models have also continued to be refined. Martin et al. further divided the pore systems and proposed a generalized PNM method, which eliminated the effects of model parameters on simulation results [11, 12]. Jiang et al. obtained an equivalent multiscale PNM by analyzing the correlation between the statistical characteristics and geometric structure parameter variables of pore network models at different scales [13]. The lattice Boltzmann method is suitable for pore-scale or representative elementary volume (REV) scale studies. However, due to the pore size distribution and connectivity, the LBM involves a sufficient number of meshes and calculations, which requires a large number of computing resources and space storage. It also leads to the limited size of the simulation space. Walsh et al. proposed the gray lattice Boltzmann method based on the partial bounce, which corrected the effect of microstructure on model permeability by adding the net flow particle ratio [14]. Ma et al. improved the convergence of the algorithm by limiting the direction of flow bounce and carried out an ideal model simulation [15]. Chenchen et al. further optimized the GLB method and carried out the 2D heterogeneous model flow simulation. The simulation results have a high consistency with the analytical solution [16, 17].
In fractured and heterogeneous tight cores, the flow patterns of oil and gas in fractures and micro-nano pores are distinct, and single-scale governing equations or flow simulation methods cannot accurately reveal the fluid flow in these porous media [18]. Therefore, scale upgrade methods are currently used for shale structure characterization and flow simulation studies, such as homogenization theory and volume average method [19]. Homogenization theory is mostly used in porous materials with a relatively regular pore structure and periodic distribution [20]. The volume averaging method is generally upgraded on the REV scale using smoothing and spatial averaging formulas [21]. At present, traditional scale upgrading methods are mostly ideal models based on periodic assumptions, which cannot fully consider the differences in real core pore space. They are difficult to apply to the flow simulation of complex and highly heterogeneous shales.
This paper employs multiple imaging techniques to analyze the interaction between pore structure and fluid flow at different scales. At the same time, the pore structure parameters and flow parameters are coupled and iterated using the pore network model and the gray lattice Boltzmann method. The results are of great significance for finely characterizing the complex multiscale pore structure and flow mechanism of unconventional reservoirs.
2. Methodologies
2.1. Multiscale Pore Space Modeling
The pore radius of unconventional reservoir cores is mainly at the nanometer and micrometer scales [22–24], and high heterogeneity exists in pore numbers, volumes, and distributions, which causes anisotropy and greatly increases the complexity of flow simulation [25, 26]. The complexity of pore space mainly lies in the following aspects. First, there are various types of pore structures, which are poorly sortable and difficult to determine the size of REV, including organic pores, inorganic pores, natural microcracks, and artificial cracks. Secondly, the pore size distribution is complex, and nano-micron-millimeter-level pore structures are well developed. The effects of strong compaction, cementation, and clay transformation result in small pore spaces, narrow pore throats, complicated structures, and poor connectivities. The pore sizes in the cores are quite diverse, including large intragranular pores formed by dissolution, and numerous intercrystalline micropores formed by clay minerals, and these differences lead to diverse oil and gas fluid patterns. Due to the existence of multiscale and interconnected pore space, heterogeneity in unconventional reservoirs is very strong (Figure 1). A scanning electron microscope (SEM), field emission electron microscope, atomic force microscope, and nano-CT can visually and intuitively analyze the micro-nano pore morphology, distribution location, connectivity, and degree of pore development in pore space. However, due to the limitation of a single observation field of vision, only combining multiple imaging technologies can effectively observe the complex pore space (Figure 2).
Because the REV scale of digital core image processing is far beyond the computational capability, and the traditional method of directly coarsening the grid will lose numerous microscopic pore structure information, reducing the simulation accuracy, this paper proposed low-precision and high-resolution subscale imaging technology (Figure 3). First, the low-precision imaging data are classified based on the selection of rock physical properties, and a coarse mesh model is established to reduce the calculation requirement. A position in the pore structure is selected for high-precision scanning to obtain the micro-nano pore information of core pore space and for multiscale flow simulation. The CT scanning method uses the attenuation law of X-ray intensity to obtain a three-dimensional structure data of the pore space inside the core [28]. One of the advantages is CT can provide gray data (electron density map) containing massive pore structure information of the sample. The core structure can be divided into pores and matrix by the binary method through the gray map. Moreover, there are much other information that can be obtained from grayscale CT images, such as local porosity and mass concentration of the transported component ([29]):
The local porosity can be calculated using the above formula by analyzing the change of the gray value of two different materials in the pore space, in which is the porosity of () and and represent the saturated and dry gray values of (). Imaging the core using appropriate resolution, the pore distribution pattern can be observed and the gray map can be obtained. On the gray map, the gray value of each point contains certain information of the pore space structure of the region. Since porosity is one of the most important factors affecting the fluid flow, local porosity, which can be calculated by the above equation, is taken as the basis for grouping the core zones. According to different porosity value ranges, core data are divided into groups, each with a different label, and a specific location in each category is selected for high-precision scanning. By means of high-precision scanning, the pore information of the typical pore structure in each type of porosity range can be obtained, consequently establishing and coupling the model at both scales.
2.2. Pore Network Model
Since Fatt (1956) constructed the first pore network model composed of unequal diameter cylindrical capillaries, the refinement of the pore network model has been continuously improved, which can not only describe pore structure characteristics more accurately but also make great progress in flow simulation [30]. In this paper, the medial axis method [31], which is commonly used in the construction of the pore network model, was adopted to obtain the positions of the pores and the pore throats by identifying the intersection of the central axis.
2.2.1. Extraction of the Medial Axis of Pore Space
The medial axis of the pore space is the skeleton curve that connects the center of each pore space, and it could reflect the basic topological characteristics of the real pore space simply and compactly. To get the central axis, the first step is to obtain the Euclidean distance map of pore space by distance transformation. The local maximum value (ridge points) was accounted to represent the topological skeleton of pore space. Secondly, pore points (object points) were defined as 26-adjacency, and matrix points (background points) were defined as 6-adjacency. Based on the Euclidean distance map and calculation of the Euler number, thinning procedure which is iteratively deleting pore points and searching the skeleton points can be used without changing the original image topology. At last, we would obtain the medial axis of pore space which is only a single voxel wide. The voxels in the medial axis can be distinguished as three types according to how many neighbours they have: regular points, endpoints, and junction points. The connectivity would be one of the key considerations during construction of the network model.
2.2.2. Recognization of Nodes (Pore Elements) and Bonds (Throat Elements)
An essential step in constructing the pore network model is to divide the digital core pore space into pores and throats in a general sense. In this paper, the node with a coordination number greater than 2 on the central axis is defined as the pore center, and the pore region is defined by the maximum incised sphere region. Parameters such as pore radius are obtained by the method of spherical isometric expansion. After removing the pore space classified above, the remaining pore space is the pore throat space. Combined with the pore node information, the corresponding throat structural characteristics are statistically analyzed.
2.2.3. Pore Structure Characteristics and Permeability Calculation
As we can assume that relatively homogenous areas exist in the heterogeneous pore space, the pore network model is helpful to better analyze the characteristics of the pore structure in the different areas. General geometric and topological properties of the pores and throats including pore (or throat) size distribution, shape factor, and coordination number can be calculated and analyzed. The pore-scale network model can be used to model flow in pore space and get flow parameters based on mathematical description of conductance and threshold pressure. And the permeability would then be used in the iteration of the GLB model. Meanwhile, each individual lattice in the GLB model is not isolated which means that there were flows between lattices in different directions (Figure 4). Correspondingly, flow capacity of the pore network model should be affected by the flows in directions. Triaxial permeabilities could separately calculate with different settings of the inlet and outlet of the medial axis, which will provide a data basis for subsequent GLB simulation.
2.3. GLBM-Based Multiscale Flow Capacity Characterization
Using the classified grid model as the base model for GLBM simulation, and coupling the microscopic triaxial permeability obtained from PNM, the pore structure upgrading and downgrading can be integrated, so that the computational loads can be reasonably reduced on the premise of preserving certain microscopic pore characteristics. The LBM method is a set of models for simulating the movement of microparticles based on molecular dynamics. It has the characteristics of entirely independent time, space, and interaction. It is a fully parallel algorithm that can carry out numerical simulation of complex fluid problems, widely used in petroleum, chemical, energy, water conservancy, machinery, and other fields. Traditional LBM only distinguishes flow lattice from nonflow lattice (Figure 5), while GLBM thoroughly considered the effects of anisotropy and heterogeneity of the porous media, and the semiflowing gray lattice was added, in which there are collisions between fluid particles and solid matrix, representing the effect of the solid scattering degree in the lattice on the flow. Generally speaking, when the fluids flow through the pore space, part of it is bounced back through the rock particles, which means that the flow field includes the sum of penetrating the pore space and the nonpenetrating bounced back fluid (Figure 5).
The governing equation can be expressed as
Different from the traditional LBM, the simulation of fluid flow can be divided into three processes: (1)Collision: in which is the internal driving force in the lattice [32], , and (2)Practical bounce-back [15]: (3)Flow:
The effective viscosity is
During the partial bounce-back process, represents the bounce-back coefficient when the fluids flow through the porous media, and it is the average scattering density of solid lattice points ranging from 0 to 0.5. It characterizes the solid structure of the porous media and is related to local solid weight in the gray lattice.
It is assumed that Figure 6 is the profile of a heterogeneous theoretical model. There are two solid boundaries in the positive and negative directions of and , and the fluid flows along the -axis. The flow area is divided into upper, middle, and lower parts, in which the upper and lower parts have the bounce-back coefficient , and the bounce-back coefficient of the middle part is . Key parameters are set as follows: the grid size is equal to 1; the timestep is 1; and the values of , , , and are set as 1, 2, 0.5, and , respectively. Since the flow is nonhomogeneous, by using the surface boundary calculation method, the velocity profile along the -axis can be obtained as follows in Figure 2, and it conforms to the basic flow law.
(a)
(b)
Based on the analysis of the flow field of the porous medium, the general N-S equation is deduced in the pore space. The flow field of partial bounce-back coefficients is explained as
By solving the Brinkman equation flow field, is related to the permeability of the local flow field:
During upscaling, the permeability in the typical core pore network model is treated as the permeability of the local flow field in the GLB model grid. By coupling permeability and partial bounce-back coefficient, the relationship between microstructure and complex multiscale porous media flow features can be established, so that the multiscale permeability of the entire core can be characterized and the upgrade is completed.
3. Results and Discussion
Tight sand cores from unconventional reservoirs are selected as the experimental sample for multiscale characterizations. The cylindrical core sample is 25 mm long with a 30 mm diameter. The permeability and porosity is 4.16 mD and 6.98%, respectively. The CT images with different resolutions are shown in Figure 7, the general trend of pore space can be analyzed based on the grayscale distribution, and the few macropores can be seen with low-resolution images. The observational field at this resolution equates with the traditional laboratory test of permeability. However, the connectivity between visible pores is very poor, so it is not suitable of carrying out pore network model modeling and flow simulation. Some large connected pores with the pore radius range of 10 to 70 μm can be extracted at high resolution, and the tiny pores and microfractures which have significant influence on fluid flow in multiscale pore space would be captured when the scanning resolution reaches 0.545 μm.
The raw CT images with a resolution of 24.8 μm were downscaled to build a GLB model, so each lattice contains information of voxels. There were three kinds of lattice: solid, pore, and gray lattice. Based on the grayscale analysis on the heterogeneous pore structure of tight sandstone, the scanning porosity of pore space was roughly classified into 8 ranges (Figure 8). About 40% lattices in the model were defined as solid which means that there were no flows. For the rest of the area, lattices with porosity less than 5% constitute most of the flow area. Therefore, the multiscale GLB model should be built by combining pore structure characteristics at different scales. The gray lattice in the GLB model which represents these areas can add more microscopic properties to help the flow simulation more accurately.
By further high-precision observation on the selected core region, relatively complete and effectively connected pore space can be observed. Each gray lattice is, respectively, represented by subvolume A (size: ) with a resolution of 2.36 μm and subvolume B (size: ) with a resolution of 0.545 μm. Affected by scanning view, the morphology of the large pore can be seen clearer and more complete. So, we used subvolume A to analyze lattices with porosity greater than 30%. These areas usually have good percolation ability and high permeability based on the Poiseuille equation and pore network model. The pore space structure becomes a little more complex when the porosity of lattice is lower than 30%. There are 9 kinds of pore structures extracted from the CT images with high resolution during porosity between 0 and 30%. The pore structure parameters of all these types of pore structures are shown in Table 1. Besides the traditional pore with different sizes, microporous clusters generated by clay mineral dissolution and interconnected microfractures are the most common pore structure. As shown in Table 1, the porosity of porous clusters usually distributes between 5% and 30% with a pore radius of about 1.75 μm; their gray value is relatively focused and slightly lower than pores. There are enough flow channels in the porous cluster so that the permeability could reach hundreds of millidarcies. Meanwhile, permeability is only 0.07 mD in the area where pores poorly develop. Compared with the microfracture network which has higher permeability in all directions, the individual microfractures usually have low permeability on two axes and no flows on the other axis. All of these microscale pore structures impacted the flow simulation in a significant way.
Limited by the porosity measured in the lab, the lattices of the GLB model, respectively, are associated with the corresponding type of the pore network model based on the porosity of the lattice and dispersion degree of the grayscale distribution of the voxels (Figure 9). The upscaling process was carried out through assigning the permeability parameters of each pore network model to the partial bounce-back coefficient of lattice in the GLB model (Figure 2).
(a)
(b)
(c)
The simulation results are shown in Table 2. Besides sample 1# which is used for imaging, we also choose two other samples 2# and 3# from the same block to analyze and model. Their CT images with low resolution have similar distribution range of grayscale map, so the pore network models extracted from sample 1# can match their downscale GLB model directly. Even if there comes an error in either of the simulation and experimental method, the results were basically the same. Using the method described in this paper, multiscale pore structure characteristics can be considered into flow simulation. The 12 types of pore network model were all extracted from connected pore space in the high-resolution images; this may be the reason why the simulation result is higher than the permeability measured in lab.
4. Conclusion
(1)The influence of the core microscopic pore structure and its heterogeneity on core physical property cannot be ignored. With the limited observing field of imaging techniques, the modeling area of pore space only represented a certain range of pore/throat size which is decided by imaging resolution. Based on the combination of the downscale method and upscale method, a new upscaling method is proposed to effectively describe the characterization of multiscale pore space(2)Imaging data with low resolution usually showed disconnected pore space, which means that the areas identified as “skeleton” under low resolution actually have a lot of micropores. To solve the problem, we used gray lattice besides the pore to represent the multiscale pore space and got its microscopic pore structure information by means of high-precision scanning. And based on classification using local porosity, the imaging steps can be able to simplify to a finite number while not affecting the model accuracy(3)By combining the pore network model and the GLB model, the new method can be a practical option to carry out flow simulation in highly heterogeneous unconventional reservoirs. Permeability of gray lattices calculated from the pore network model with high resolution would be assigned to a partial bounce-back coefficient of the gray lattices in the GLB model to solve the problem that connection between different scales of pores can be barely obtained. And we also consider putting the pore network model in the loop iteration of the GLB model to further increase accuracy in future work
Data Availability
The data used to support the findings of this study are available from the first author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
We would like to express appreciation to the following financial support: Shandong Province Natural Science Foundation (ZR2019PEE039 and ZR2018BEE008), National Natural Science Foundation of China (51804328), and Production-study-research Innovation Foundation of Shandong Academy of Sciences (YDLH41903CZ, YDLH41906CZ, and YDLH41901CZ).