Comprehensive Investigation of the Petrophysical and Two-Phase Flow Properties of the Tight Sandstone in Yanchang Formation, Ordos Basin, China: Insights from Computed Tomography Imaging and Pore Scale Modelling

To meet for the increasing global demands of energy, the development of unconventional reservoirs, such as the tight sandstone reservoirs, shale oil/gas, natural gas hydrate, and coalbed methane, has attracted more and more attentions in petroleum industry [1–7]. The proved reserves of the unconventional reservoir are abundant in China, and the geological resources of tight oil/gas in Ordos, Songliao, and Junggar Basins have become an important replacement area for China’s oil production capacity construction [8]. Due to the sedimentation and diagenesis process, the tight sandstone reservoirs are generally characterized by complex pore-throat structure with low porosity and permeability, which makes it significantly different from the conventional reservoirs [9–12]. Routine core analysis relies on core drilling of natural samples; however, the strong heterogeneity and complex pore-throat structure of tight sandstone reservoirs will lead to a large discrepancy of test results from sample to sample, even under a consistent test environment. Besides, the conventional laboratory experiments are mostly destructive (e.g., compressive test, mercury intrusion porosimetry, and adsorption experiments) which will cause damage and pollute the core samples over time and restrict the reuse of the core samples for the parallel analysis [13, 14]. As an alternative to laboratory tests, the digital rock petrophysics (DRP) based on CT imaging and model reconstruction has been regarded as an efficiency approach to compute and predict petrophysical parameters over the last decades. The DRP approaches can predict the petrophysical, acoustical, electrical, and thermal properties, by simulating the corresponding physical process to the laboratory tests on the reconstructed pore or matrix model of rocks. The DRP technique owes the advantages on model repeatable, realizing the coupling of multiple physical responses, as well as the quantitatively controlling of microeffecting factors and model parameters [15–18]. There has been a continuous increasing attention on the DRP research and application in tight sandstone reservoirs [19–22].

The reconstruction of digital models from the pore structure images and corresponding datasets (e.g., CT, SEM, and FIB-SEM) is a key step of the DRP analysis and is becoming relevant. According to the different reconstruction algorithms, the digital model can be divided into two types, i.e., physical-based model and statistical-based model [23]. The physical-based model is reconstructed from 3D volume datasets such as CT and FIB-SEM images, which is a physical reappearance of rock structure. The statistical-based model (e.g., multipoint geological statistics method (MPS) and process-based method) is reconstructed from 2D images (e.g., thin section and SEM image) or statistical parameters derived from routine core analysis (e.g., MIP) [24–28]. The main challenge that is associated with the physical-based model is how to balance the optimal scanning resolution with the sample size [13, 29]. Compared to the physical-based model, the challenges of the statistical-based model are how to find a more suitable evaluation function to improve the accuracy, the multiscale, and multicomponent issues [18]. The inherent limitations of the algorithms make it difficult to model those samples with complex structure and strong heterogeneity, e.g., tight sandstone and carbonate [30–33].

Using the reconstructed pore scale models, the absolute and relative permeability of the core sample can be predicted. There are mainly four types of numerical solvers for permeability prediction based on CT images, including PNM [29], lattice Boltzmann method (LBM) [34], voxel-based solver (VBS) [35], and CFD [36]. The PNM solver is operated on the extracted connected pore-throat network with topologically representative [37], which calculates the conductance between adjacent pores and throats by the Hagen-Poiseuille equation and the permeability by Darcy’s law [38]. In the last three decades, PNM has been widely used in porous media-related industries (e.g., petroleum, food, material, and chemistry) to study the complex porous flow phenomenon, including phase change problem, reactive flow, non-Newtonian fluid, and particle transport [39–47]. However, during the pore network extraction, approximations are made concerning the pore-throat geometry of rock. The LBM solver solves the Boltzmann equation on the segmented CT image and treats the flow process as a continuous motion and collision process of fluid particles in mesoscale [48]. The VBS solver directly solves the Stokes and Navier-Stokes (NS) equation using the finite volume method (FVM) on the segmented CT image. During computation, the adaptive grid is used to reduce the cell number rather the regular grid [49, 50]). Both the LBM and VBS treat the pixels of the segmented CT image as the lattice or element for numerical calculation, which are suitable to solve porous flow case with complex geometries and internal structure. However, the LBM solver is easy to code, but time-consuming even when using a massively parallel implementation, especially for the tight sandstones [51]. As for the VBS solver, the convergence speed depends on the complexity of the model, and it is also a challenge to simulate the slip boundary conditions [52]. In recent years, the combination of CFD-based solver and phase-interface tracking algorithms (e.g., VOF, level set (LS), coupled of VOF and LS, and phase field (PH)) has been regarded as a practical approach to simulate single- and multiphase flow in porous media with a wide range of viscosity and density ratio due to its mass conservation principle and superior ability of the numerical computing [53–55]. The CFD solver computes the N-S equation on the grids of pore region. The tetrahedron meshes adopted in CFD computation can be generated by most commercial code (e.g., Mimics, Simpleware, and Avizo), which are prior to preserve the natural geometry of samples and describe the transient phase distribution in the flow channel (pores) while considering capillary and viscous forces. However, the challenges in CFD modelling lie in how to balance the number and quality of the meshes, to promote both the convergence and accuracy. In recent years, the DRP approaches have been widely adopted to estimate the petrophysical properties and simulate the pore scale physical processes of unconventional reservoirs (e.g., tight sandstone and carbonate reservoirs) [32, 33, 56–59]. A comprehensive analysis combination with the micro-CT images and multi-DRP approaches is conducted to characterize the petrophysical and fluid transporting characteristics of the tight sandstone reservoirs. DRP has made abundant progress, especially in microscopic flow mechanism. However, limited by the inherent shortages of the DRP model, which cannot be directly scaled up the physical model, one can only scale up the single- or multiphase flow properties from pore scale to core scale and finally reservoir scale by developing various analytical and numerical methods like the homogenization method ([60]). Compared to flow upscaling, the upscaling research is far advanced in geomechanics and rock mechanics ([61–63]). Aiming at the scaling-up studies on flow properties, several works both in numerical methods and analytical algorithm were presented in literatures ([64–68]). Directly scaling up the digital model from pore scale to core scale and then linking the estimation properties to reservoir scale could be the most promising way of the DRP.

This paper analyzes the petrophysical and two-phase flow properties of the tight sandstones drilled from the Chang-7 member in Yanchang Formation in China, coupling with the micro-CT technology, image analysis, PNM, and CFD simulations. Compared to previous studies, the effect of the pore structure heterogeneity on the fluid flow properties and the residual fluid distribution was comprehensively investigated, especially the evolution of the fluid interface during the simulation of the water flooding process using the VOF model by CFD in this work. The main contents include three aspects: the micro-CT imaging is conducted to capture the pore-throat system of the tight sandstones. Then, the pore-throat structure is quantitatively characterized and analyzed by using PNM including pore geometry, pore-throat size distribution, connectivity, and tortuosity. Finally, the pore scale single- and two-phase flow simulations are performed on reconstructed pore structure to investigate the permeability anisotropy and residual fluid trapping mechanism of the tight sandstone.

The tight sandstones used in this study are drilled from Chang-7 member in Yanchang Formation, in the Ordos Basin, which has become the most important petroleum basin in China (Figure 1(a)) [62, 63, 69]. The basin can be divided into 6 tectonic units, as shown in Figure 1(b), i.e., Yimeng Uplift, Yishan Slope, Weibei Uplift, western thrust belts, Tianhua Depression, and Jinxi fold belt.

In this study, two fresh core samples collected from two different wells (Z93 and Z194) are selected and drilled from uniform core profiles (Figure 1(c)). The basic parameters and well information of the samples are listed in Table 1.

The minicore plugs with the size of 2-5 mm in diameter and 5-8 mm in length are drilled from standard core plug (⁠$25∗50mm$⁠) and prepared for CT scanning [70]. The sample is scanned by X-ray from different directions (with rotation of a sample table step-by-step with a specific angle), and the detector is used to acquire the radiation quantity of the rock components (e.g., pore and minerals) with different density. The projection data of each section is then collected by detector to reconstruct the 3D structure of core sample by image reconstruction algorithm. The schematic of the standard scanning process is illustrated in Figure 2(a). The Phoenix nanotom M type scanner (Figure 2(b)) is adopted in this study for CT scanning with the highest theoretical spatial resolution of 0.5 μm.

A series of image processing operations are performed on original CT image to improve the quality including denoising, contrast enhancement, and filtering (Figure 3). Then, for the purpose to reduce the effect of exterior region on image processing, as well as for computation efficiency, the representative elementary volume (REV) is analyzed to determine the region of interest for subsequent computation and analysis (Figures 3(b) and 3(e)). The workflow of image processing is illustrated in Figure 3. After contrast enhancement operation, the pore phase can be distinguished from matrix, and the gray value distribution (Figures 3(c) and 3(f)) can be used as the benchmark for segmentation. The statistical parameters of the original CT images are listed in Table 2. The dimensions of the original CT images are $1257∗1257∗1385$ pixels and $1350∗1350∗814$ pixels for samples Z93 and Z124, respectively. The pixel resolutions are 0.6 and 1 μm; hence, the physical sizes of the samples are $0.7542∗0.7542∗0.831mm$ and $1.35∗1.35∗0.814mm$⁠, respectively.

The binarization can be operated on the processed CT image to extract the pore structure and corresponding PNM. In this study, the watershed segmentation algorithm is used for multicomponents’ identification of pore, matrix, and high-density minerals. The watershed algorithm simulates the flooding from a set of labeled regions in a 2D or 3D image. It expands the regions according to a priority map, until the regions reach at the watershed lines. The process can be regarded as a progressive immersion in a landscape. More details can be referred to Pan et al. [71].

Once the CT image is segmented, the maximum ball (MB) algorithm is applied on the extracted pore space phase to generate the corresponding PNM. The MB algorithm was first proposed by Silin and Jin [72] and improved by Blunt et al. [52] and has been widely used in porous media characterization and modelling. The PNM can be regarded as a series of maximal inscribed spheres of pore bodies to represent the pore space, and the searching process of the maximum inscribed sphere of the pore body is illustrated in Figure 4. And the throat between two adjacent pores is replaced by a cylinder. The details for PNM extraction can be referred to Dong [37]. Pore structure can be quantitatively characterized using PNM analysis including pore-throat size distribution, pore geometry (shape factor), coordination number, and tortuosity.

Only the connected pores contribute to the fluid flow; thus, connectivity analysis is performed on the extracted pore space prior to permeability computation. The connectivity analysis is conducted in Avizo software using the built-in pixel-based algorithm, which determines the interconnection of two pixels if any vertex of pore pixels is connected. The extracted connected pore space of the tight sandstones is shown in Figure 5.

The immiscible two-phase flow is simulated using the VOF model by CFD to investigate the waterflooding development characteristics of tight sandstone reservoir. The effect of pore-throat heterogeneity on waterflooding interfacial dynamics and oil recovery efficiency is investigated and analyzed. The mesh model of the connected pore space is reconstructed using the methods been proposed and validated in our previous study [14, 75]. To reduce the mesh number and promote the computation efficiency, a subdomain of 350³ voxels is extracted for mesh generation, and the reconstructed finite element mesh model of the tight sandstone is shown in Figure 6.

The PNM of the tight sandstone samples is extracted from the segmented CT image by using MB algorithm. The extracted PNM of the tight sandstones is shown in Figure 7. The upper row shows the extracted pore space and corresponding PNM of the sample Z93. The lower row is the extracted pore space and corresponding PNM of the sample Z124, in which the size of the spheres is more homogeneous than the sample Z93. The parameters of the cubic subvolume extracted for PNM are listed in Table 3. Based on the PNM, the statistical parameters of the pore-throat structure are computed and analyzed.

Based on the extracted PNM, the pore-throat size distributions of the tight sandstones are calculated and listed in Figure 8. The pore diameter of the sample Z93 is mainly distributed between 5 and 10 μm (the mean value is 4.17 μm), and the sample Z124 is between 10 and 15 μm (the mean value is 13.41 μm). The throat diameter of the samples Z93 and Z124 is mainly distributed between 2~4 μm and 3~5 μm, and the corresponding mean values are 3.24 μm and 4.02 μm, respectively. Other parameters including the pore volume, throat area, and throat length are also obtained and listed in Table 4. The results indicate that the pore structure of the tight sandstones shows strong heterogeneity on the pore-throat size distribution, as well as the pore volume.

The coordination number is introduced to describe the connectivity which represents the number of throats connected to the pore ([82]). In general, the larger the coordination number is, the higher the connectivity is. Based on the extracted PNM, the coordination number of the tight sandstones is calculated and shown in Figure 10. The coordination numbers of the sample Z93 are mainly distributed between 2 and 5 with the maximal and mean values of 22 and 3.97, respectively. The coordination numbers of the sample Z124 are mainly distributed between 3 and 8 with the maximal and mean values of 38 and 4.09, respectively.

The tortuosity is used to describe the degree of circuity of the flow path in porous media (e.g., rock and soil), which is also a key parameter to describe the complexity of the pore structure. It is firstly proposed by Epstein [83] and defined as the ratio of the actual flow path length to the apparent length of the flow channel (i.e., the physical size of the medium) through the porous media, as shown in Figure 11. The tortuosity is difficult to determine by physical experiments. It is usually estimated by empirical relations combined with other lab-testable petrophysical parameters such as porosity, pore-throat size, surface area, and capillary pressure [49, 50]). Estimation of the tortuosity by empirical relations is a complex and low-accuracy way since there are plenty of parameters. In this study, an open-source code (TauFactor) developed by Cooper et al. [84] is used to calculate the tortuosity value in three directions based on the CT images of the tight sandstones. More details can be referred to our previous study [85]. Both the results calculated by simulations and estimated by the empirical relation are listed in Table 5. The calculation results indicate that the tortuosity shows great deviation in the different directions, which is believed to be the main cause of the heterogenous pore structure. Besides, the estimation value is much lower than the calculation value, since the empirical relation used in this study adopts the total porosity value to estimate the tortuosity, which may weaken the influence of the complexity and heterogeneity of the pore structure.

In general, the pore structure of the tight sandstone is extremely complex and develops many discrete and small pores with diverse morphology. Besides, there is a transition region in the CT image with gray value between matrix and pore phase which is probably some kinds of argillaceous filling during diagenetic evolution. It makes the pore structure more complicated, and the segmentation of this part is also challenging. To further quantify the pore structure, the obtained petrophysical parameters are statistically analyzed, as shown in Figure 12 and Table 6.

As shown in Figure 12, the porosity shows a strong positive correlation with the pore-throat size. While due to the limit of sample number, the correlation between pore-throat size and tortuosity (coordination number) is not obvious. Anyway, it provides a hint that using basic parameters like porosity and pore-throat size to estimate tortuosity is not rigorous. And the permeability estimated by such parameters will also cause larger deviation which agrees well with our previous findings [13]. Those empirical relations apply only to some certain cases.

The results indicate that the pore structure of the tight sandstone exhibits strong heterogeneity in spatial connectivity which may have a strong effect on the accurate description of the effective flow channels and/or reserve estimation at the reservoir scale. Besides, the strong heterogeneity of the pore structure could lead to an uneven distribution of the hydrocarbon in the reservoir, which is also an indirectly evidence to confirm the necessity of large-scale hydraulic fracturing operation to enhance the connectivity of the reservoirs during the exploitation of the tight oil and gas.

Based on the extracted pore structure, the single-phase flow simulation is conducted in three directions on the segmented CT images by using VBS solver to investigate the anisotropic permeability. The periodic pressure boundary is adopted in the simulation with an inlet pressure of $1.3e+05Pa$ and outlet pressure of $1e+05Pa$⁠. The other boundaries and pore-grain interfaces are set as nonslip boundary (i.e., the $u=0$⁠), and the simulated fluid viscosity is 0.001 Pa·s.

The normalized velocity-based streamline patterns and the corresponding pressure fields of the sample Z93 are shown in Figure 13. As shown in the streamline patterns, the velocity distribution varies with the directions. For the $x$ and $y$ directions, there is a narrow throat in the middle of the sample, which restricts the flow capacity, and the corresponding pressure nephogram also shows obvious mutation in this area. The pressure distribution in the $z$ direction parallel to the lamellar area shows a relatively uniform trend, and the distribution density of streamlines is also more uniform. The calculated permeabilities are 0.0157 mD, 0.0253 mD, and 0.0402 mD in the $x$⁠, $y$⁠, and $z$ directions, respectively. The permeability value in the $z$ direction is larger than the other two directions and is also closer to the measured value.

The corresponding pressure fields of the sample Z124 are shown in Figure 14. Compared to the sample Z93, the distributions of streamlines and pressure of the sample Z124 are more homogeneous, even so it still shows great differences in the different directions. The calculated permeabilities are 0.0837 mD, 0.0305 mD, and 0.0595 mD in the $x$⁠, $y$⁠, and $z$ directions, respectively. The permeability value in the $x$ direction is larger than the other two directions and closer to the measured value. The simulated results indicate that the permeability of tight sandstones is anisotropic. The comparison of permeability between the computed and measured values is listed in Table 7.

The simulation results indicate that the permeability of the tight sandstone is a strong anisotropy in spatial distribution, which provides a hint for drawing up the reason considering the preferential flow channels (high-permeable direction or region). Besides, it is also necessary to make dynamic adjustment of the injection-production plan (e.g., water shut-off and profile control) in the high water-cut period to avoid water channeling and other problems.

The pore scale two-phase flow characteristics of the tight sandstone are investigated by simulating the immiscible waterflooding process in the reconstructed mesh model using the VOF model by CFD. During simulation, a two-pore volume (2PV) water is injected to simulate the waterflooding process. In the simulation, the fluid flow direction is set along $z$ direction, and the other four faces are set as wall boundary. The simulation is assumed to be laminar flow, and the SIMPLEC method is used to operate the velocity-pressure coupling [86]. The default relaxation factor is adopted, and the absolute convergence condition is set to $1e−5$⁠. The linear upwind scheme is adopted to discrete the convective term. The uniform wettability is assumed in the simulation, and the contact angle is set to 150° for waterflooding. The related parameters of the fluids used in the simulation are listed in Table 8. To visualize the phase distribution and the displacement front evolution process, taking sample Z93 for example, the distribution of the water-oil saturation at different time steps is shown in Figure 15. Once the simulation is completed, the relative permeability can be obtained by calculating the water saturation (⁠$Sw$⁠) and corresponding flux of different phases at the outlet in different time steps, and the relative permeability curve is shown in Figure 16.

In the simulation, the model is set as oil wet. Thus, the injected water is preferentially advanced along the middle region of the larger pores with the smaller resistance, which is in good agreement with the single-phase flow simulation. As shown in Figure 15, due to the strong heterogeneity of the pore-throat size, the whole displacement process exhibits obvious fingering phenomenon, which causes hysteresis in the displacement process. The injected water is preferentially reaching to the outlet along several main pore-throat chains, which reduces the sweeping efficiency of injected water. When the displacing water breaks through along the channel with smaller resistance, the residual oil phase between the preferential flow channels is trapped. This phenomenon mainly occurs in the small pore-throat group surrounded by the large pores, which also provides the theoretical basis for water plugging and profile control in oilfield operation. In addition, due to the low sweeping efficiency, after water breakthrough, there are still many residual oil blocks trapped in the corner of the model and the dead-end pores, which forms another occurrence state of residual oil. The pore scale visualization of water-oil distribution reveals the mechanism that oil in the narrow throats and microfractures is impossible to be displaced during waterflooding due to the tremendous capillary resistance, which also provides a hint and theoretical support for polymer flooding to adjust the fluid interfacial force state in the field practice.

As shown in Figure 16, due to the complexity of pore-throat system and poor connectivity, the ultimate oil recovery and flooding efficiency of two sets of models are quiet low, and the residual oil saturation is nearly 40%. Due to the strong heterogeneity and flow resistance (large capillary force), the flow capacity of the injected water increases slowly at the early stage. With the water injection, the relative permeability of the water increases significantly when the main flow region is formed along several large pore channels in the model, but the increase is still limited and reduces the sweeping efficiency of the injected water.

Compared to the sample Z93, the porosity, pore-throat size, and coordination number of the Z124 are larger than Z93, which indicates that the seepage conditions of the pore-throat of the sample Z124 are better to some extent. The equal point of the relative permeability of the sample Z124 is reached at a lower water saturation. Besides, the co-seepage area of water and oil of the sample Z124 is also larger.

Although the DRP is a powerful approach for petrophysical parameter estimation and flow characteristics modelling, since the physical dimensions of digital model reconstructed from image data are usually small (~mm scale), the researchers and industry are still confused whether the results derived from DRP can accurately describe the reservoirs. It is true that limited by the computing capacity, as well as to balance the contradiction between the data size and image resolution, scanning and modelling of the whole core at a ~μm or ~nm scale is not a feasible or cost-effective way for petrophysical study. Therefore, currently, the DRP is mostly used for mechanism characterization and study, including heterogeneity analysis, sedimentation environment, and mineral calibration, and the results of multiphase flow simulation can be used for model calibration and further prediction.

Image-based digital rock analysis and pore scale simulation are of great significance and beneficial to reduce the subsurface uncertainties in petrophysical properties prediction compared to routine core-scale laboratory tests, especially for those samples with strong heterogeneity and anisotropy. In this paper, the pore structure and fluid flow properties of the tight sandstone reservoirs are comprehensively analyzed by using the micro-CT imaging, PNM analysis, and pore scale flow simulation. The following conclusions can be achieved:

• (1)

  The micro-CT image-based PNM analysis indicates that the pore structure of the tight sandstone is mainly composed of discontinuous tiny pores and thin-lamellar pores with poor connectivity and diverse pore morphology. In addition, the tortuosity of the tight sandstones is calculated and compared to the estimated value. It is reasonable to assume that the variety of tortuosity in different directions is believed to be one of the main causes of the anisotropy of rocks, and it can also be regarded as an indicator to evaluate the anisotropic permeability

• (2)

  The effects of the pore structure on permeability anisotropy are investigated and analyzed by using CT image-based DNS. Due to the complexity of pore structure, the tight sandstone exhibits obvious anisotropy

• (3)

  The immiscible two-phase flow simulation is conducted on the reconstructed mesh model by CFD modelling. The strong heterogeneity of pore-throat size can lead obvious fingering during waterflooding, which can reduce the sweeping degree of the injected water and affect the ultimate recovery

Compared to routine core analysis, the work presented in this paper approves that the image-based DRP approach is a time-saving and cost-effective alternative way to evaluate petrophysical and hydraulic properties for tight sandstone reservoirs. However, the balance between the minimum resolution and the computation accuracy to capture the representative features of rocks becomes relevant and should be paid more attention in future studies.

All the data and materials that support the findings used in this paper are available from the corresponding author upon reasonable request.

The authors declare that there are no competing financial interests with any other people or groups regarding the publication of this manuscript.

The Changqing Oilfield was greatly acknowledged for providing the core samples used in this study. This work was financially supported by the Natural Science Foundation of SWUST (Grant No. 20zx7129) and the National Natural Science Foundation of China (Grant No. 51909225).