Seepage Mechanism of Tight Sandstone Reservoir Based on Digital Core Simulation Method

: Recently, tight sandstone oil has played an increasingly important role in the energy strate-gies of countries around the world. However, the understanding of a microscopic mechanism is still not clear enough, which has been affecting the improvement of the recovery of tight sandstone oil. In this article, a digital core model was established to simulate the pore network of a physical core with CT scan and difference equations were veriﬁed by Fourier counting. Then, a combination of orthogonal tests and cubic digital cores was used to experimentally investigate various parameters including pressure, length, permeability, viscosity, and time. By combining the physical experiments with the digital core methods, it can be observed that the state of the micro-crack affects the conductivity of the core, which may be the decisive reason for changing the pressure gradient. The orthogonal test showed that the sensitivity of the parameters was pressure, length, permeability, time, and viscosity in order. The results of the numerical simulations showed that this method can reveal the seepage mechanism of a tight sandstone reservoir, greatly shortening the experimental time and improving ﬂexibility.


Introduction
Reservoirs with a permeability of less than 0.1 mD are generally referred to as tight sandstone reservoirs [1]. From 2000 to 2017, production of tight sandstone oil has increased. In 2019, accounting for 53% of a total crude oil product in the United States [2]. The multistage and multi-cluster fracturing technologies in horizontal wells have become the main means while developing tight sandstone oil [3]. When fractured, enormous micro-cracks form in the reservoir around the bottom of the hole, causing the formation of small pieces and improving the conductivity of the reservoir.
At present, there are many macroscopic studies on the seepage problem of a reservoir matrix-fracture dual system, but research studies on the law of matrix failure at a microscopic level are rare [4,5]. Moreover, the study of the seepage law of the matrix now is of great importance for developing tight sandstone oil reservoirs effectively. Liu et al. suggested that distribution of tight sandstones could be observed by Octamethylcyclotetrasiloxane [6]. Several researchers noted that a microscopic mechanism in seepage could be observed by the combination of an indoor experiment and a lithography physical model [2,7]. Taken together, many researchers have done a lot of work and achieved remarkable achievement [3,4,8,9]. However, due to strict experimental conditions and a preparation cycle, many of the above experiments cannot be repeated, or have a high cost.
In this paper, a more convenient method combining the physical experiment with a digital core is designed to reveal seepage mechanism at a microscopic level and greatly shorten the experimental time. First, the digital core for the simulation calculation is established for the subsequent simulated annealing algorithm to build the matrix block simulation. Then, the system simulation results of the core are compared with the indoor experimental results to verify the accuracy of the simulation model. The simulation shows that the opening and closing state of the micro-crack affects the conductivity of the core, which may be the reason for the starting pressure gradient [10]. The pros and cons of methods such as Gaussian field method [11], simulated the annealing method [12], and the Markov chain-Monte Carlo method [13] are discussed. Finally, the pore network model established by a simulated annealing algorithm is used to construct three-dimensional matrix cubes with different parameters. The orthogonal test method [14] was used to combine different levels of various parameters. The factors affecting the degree of exhaustion were analyzed and finally sorted to guide production practices.
The rest of this paper is listed as follows. Section 1 introduces the digital core modeling methods. Section 2 proposes the different equations for digital cores. Section 3 presents and analyzes the stability of different equations. Finally, the results and discussion are drawn in Section 4.

Modeling of the Digital Core
Indoor experimental analysis methods can determine the physical properties of cores by using the fixed instruments and methods [2], which we use to analyze and optimize the elliptical flow around the well bottom shown in Figure 1. The fixed instruments mainly include pumps and core holders. When using this method individually, the cycle and cost will be tremendous. Additionally, it may cost several days or even longer to retrieve complete data of recovery versus time. However, this method needs a lot of time and cost. Scanning electron microscopy imaging by CT is done to excite a beam of electrons to the surface of the sample, under the region where the electrons in the inner layer are bombarded by foreign electrons [15,16]. The characteristic X-ray spectrum with the mineral component information is received by the optical signal receiver. Thus, a point-by-point scan of the core surface provides an image of the mineral distribution on the surface of the core sample [13]. This type of simulation technology can estimate core parameters using minimal cost in a very short time [17]. Since it can model the pore structure inside the core, the microscopic seepage mechanism inside the core can be visually revealed in the simulation process.
In this paper, a more convenient method combining the physical experiment with a digital core is designed to reveal seepage mechanism at a microscopic level and greatly shorten the experimental time. First, the digital core for the simulation calculation is established for the subsequent simulated annealing algorithm to build the matrix block simulation. Then, the system simulation results of the core are compared with the indoor experimental results to verify the accuracy of the simulation model. The simulation shows that the opening and closing state of the micro-crack affects the conductivity of the core, which may be the reason for the starting pressure gradient [10]. The pros and cons of methods such as Gaussian field method [11], simulated the annealing method [12], and the Markov chain-Monte Carlo method [13] are discussed. Finally, the pore network model established by a simulated annealing algorithm is used to construct three-dimensional matrix cubes with different parameters. The orthogonal test method [14] was used to combine different levels of various parameters. The factors affecting the degree of exhaustion were analyzed and finally sorted to guide production practices.
The rest of this paper is listed as follows. Section 1 introduces the digital core modeling methods. Section 2 proposes the different equations for digital cores. Section 3 presents and analyzes the stability of different equations. Finally, the results and discussion are drawn in section 4.

Modeling of the Digital Core
Indoor experimental analysis methods can determine the physical properties of cores by using the fixed instruments and methods [2], which we use to analyze and optimize the elliptical flow around the well bottom shown in Figure 1. The fixed instruments mainly include pumps and core holders. When using this method individually, the cycle and cost will be tremendous. Additionally, it may cost several days or even longer to retrieve complete data of recovery versus time. However, this method needs a lot of time and cost. Scanning electron microscopy imaging by CT is done to excite a beam of electrons to the surface of the sample, under the region where the electrons in the inner layer are bombarded by foreign electrons [15,16]. The characteristic X-ray spectrum with the mineral component information is received by the optical signal receiver. Thus, a point-by-point scan of the core surface provides an image of the mineral distribution on the surface of the core sample [13]. This type of simulation technology can estimate core parameters using minimal cost in a very short time [17]. Since it can model the pore structure inside the core, the microscopic seepage mechanism inside the core can be visually revealed in the simulation process.  The method used in rock mass scanning is to accurately control the parameters of the CT testing machine through a computer shown in Figure 2. The X-ray generator (Zeiss) then emits X-rays that penetrate the sample fixed and to be rotated on the sample holder. The other end is an X-ray detector and the data acquisition system of the CT machine collects the data acquired by the detector [18]. Finally, the holder was rotated at a certain speed and was re-recorded. The cumulative rotation is 360 • . The CT scanning and modeling process are shown in Figure 3.
The plunger-shaped core is placed into the CT scanner to obtain a CT scan sectional view. By modeling the scanned images of different sections and angles, a 3D pore network model of the plug core is obtained as shown in Figure 4.
The method used in rock mass scanning is to accurately control the parameters of the CT testing machine through a computer shown in Figure 2. The X-ray generator (Zeiss) then emits X-rays that penetrate the sample fixed and to be rotated on the sample holder. The other end is an X-ray detector and the data acquisition system of the CT machine collects the data acquired by the detector [18]. Finally, the holder was rotated at a certain speed and was re-recorded. The cumulative rotation is 360°. The CT scanning and modeling process are shown in Figure 3.  The plunger-shaped core is placed into the CT scanner to obtain a CT scan sectional view. By modeling the scanned images of different sections and angles, a 3D pore network model of the plug core is obtained as shown in Figure 4.  The method used in rock mass scanning is to accurately control the parameters of the CT testing machine through a computer shown in Figure 2. The X-ray generator (Zeiss) then emits X-rays that penetrate the sample fixed and to be rotated on the sample holder. The other end is an X-ray detector and the data acquisition system of the CT machine collects the data acquired by the detector [18]. Finally, the holder was rotated at a certain speed and was re-recorded. The cumulative rotation is 360°. The CT scanning and modeling process are shown in Figure 3.  The plunger-shaped core is placed into the CT scanner to obtain a CT scan sectional view. By modeling the scanned images of different sections and angles, a 3D pore network model of the plug core is obtained as shown in Figure 4.   The method used in rock mass scanning is to accurately control the parameters of the CT testing machine through a computer shown in Figure 2. The X-ray generator (Zeiss) then emits X-rays that penetrate the sample fixed and to be rotated on the sample holder. The other end is an X-ray detector and the data acquisition system of the CT machine collects the data acquired by the detector [18]. Finally, the holder was rotated at a certain speed and was re-recorded. The cumulative rotation is 360°. The CT scanning and modeling process are shown in Figure 3. The plunger-shaped core is placed into the CT scanner to obtain a CT scan sectional view. By modeling the scanned images of different sections and angles, a 3D pore network model of the plug core is obtained as shown in Figure 4.

Difference Equations for Digital Cores
For the core numerical model, the difference equations for pressure propagation need to be specified [19]. It is assumed that the n th time step physical quantity P of the problem is known and the n+1 time step P value is to be sought. Choose any time and space point (i, j, k, n + 1). P(x, y, z, t) is substituted for the second-order difference quotient of (x, y, z) and the backward difference quotient of t at any point (x, y, z, t) to obtain difference equations at any point (i, j, k, n + 1).
Among them: For the symbols appearing above, the physical meaning of each symbol is shown in Table 1.

Stability Analysis of Difference Equations
Using the Fourier analysis [20,21] method, let the error be ε and get the error equation as: According to the principle of Fourier series [19][20][21], the error term is written as a simple harmonic (complex) form.
The growth factor A was obtained as follows.
It is known that |A| ≤ 1. Therefore, the difference equation is unconditionally stable and the numerical model can be directly used to solve the equations.

Parameter Selection
Dividing the indoor experiments into four groups, the downstream pressures were 0 MPa, 10 MPa, 20 MPa, and 30 MPa, respectively, and the initial value of the upstream pressure remained unchanged at 43 MPa. By comparing the output cumulative yield and pressure distribution at different outlet pressures at the same time, the difference between the digital core model and the physical experiment can be verified, thereby, verifying the validity of the numerical model. The parameters of the example mentioned next are shown in Table 2.

D Matrix Block Modeling
The pore network model in Figure 4c is modeled by a scanned image. In order to study the 3-dimensional space rather than the one-dimensional seepage law, it is necessary to re-establish the dense matrix model with a fast pore network in 3D space. Extending the dimension makes it satisfy the statistical structure of the pore structure of the original numerical model. Current mainstream statistical modeling methods include Gaussian field method simulated annealing and the Markov chain-Monte Carlo method. Compared with the Gaussian field method, the simulated annealing method can take more rock information into account, making the model closer to the real core [22,23]. Compared with the Markov-Monte Carlo method, the complexity of simulated annealing is lower.
We only regard the core model as a two-phase system of rock skeletons and pores. Since the space where the fluid is the pore, the calculation of each characteristic function value focuses on the pore space. The modeling objective function is composed of the autocorrelation function of the reconstruction system and the reference system. The objective function of the system at a certain time can be defined as follows.
where α i , β i are the weight value of the different constraint functions and E is the total energy of the system. S(r) represents the probability that two points of any distance r in the system are simultaneously pore phases [15]. It has the property S j 2 (0) = ϕ j , and j represents the j th phase. L(r) is the connectivity function of the pores and r is the distance between the two points. For a one-way system with a porosity of ϕ, L(0) = ϕ.
The digital core modeling steps are as follows: (1) Binarize the actual core CT image or electron microscope scan image, collect model statistics ϕ, S(r), L(r). (2) Randomly generate a three-dimensional digital core with a porosity of 0.3 and calculate the initial model objective function E to enter the optimization process. (3) Take a hole point and a skeleton point arbitrarily and interchange their positions to generate a new system, and then calculate the objective function E. (4) The new system is determined by the Metropolis criterion if it is good enough to replace the old system in order (5) to decide whether to terminate the cycle according to the inner and outer loop termination condition. (6) Finally, output the digital core simulation results and end the process [16,24]. The modeling results are shown in Figure 5

Orthogonal Test Analysis
After verifying the validity of the numerical model, the scope of theoretical research can be extended to a three-dimensional space. We study the sensitivities of different factors regarding the efficiency of a dense matrix failure. Figure 6 is an example of a pressure field simulation of a depletion process.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 6 of 10 olis criterion if it is good enough to replace the old system in order (5) to decide whether to terminate the cycle according to the inner and outer loop termination condition. (6) Finally, output the digital core simulation results and end the process [16,24]. The modeling results are shown in Figure 5.

Orthogonal Test Analysis
After verifying the validity of the numerical model, the scope of theoretical research can be extended to a three-dimensional space. We study the sensitivities of different factors regarding the efficiency of a dense matrix failure. Figure 6 is an example of a pressure field simulation of a depletion process. Moreover, theoretical studies were carried out using orthogonal experiments. A total of 16 orthogonal experiments were performed and a total of 5 factors and 4 levels (0 MPa, 10 MPa, 20 MPa, and 30 MPa) were selected with results of the selection being shown in Table 3 below. In addition, the boundary pressure is set to 0, while the remaining settings are the same with Table 2. Finally, output the digital core simulation results and end the process [16,24]. The modeling results are shown in Figure 5.

Orthogonal Test Analysis
After verifying the validity of the numerical model, the scope of theoretical research can be extended to a three-dimensional space. We study the sensitivities of different factors regarding the efficiency of a dense matrix failure. Figure 6 is an example of a pressure field simulation of a depletion process. Moreover, theoretical studies were carried out using orthogonal experiments. A total of 16 orthogonal experiments were performed and a total of 5 factors and 4 levels (0 MPa, 10 MPa, 20 MPa, and 30 MPa) were selected with results of the selection being shown in Table 3 below. In addition, the boundary pressure is set to 0, while the remaining settings are the same with Table 2. Moreover, theoretical studies were carried out using orthogonal experiments. A total of 16 orthogonal experiments were performed and a total of 5 factors and 4 levels (0 MPa, 10 MPa, 20 MPa, and 30 MPa) were selected with results of the selection being shown in Table 3 below. In addition, the boundary pressure is set to 0, while the remaining settings are the same with Table 2.

Results and Discussion
By comparing the simulated value with the physical one, mutual verification of digital and physical models can be performed.

Results and Discussion
By comparing the simulated value with the physical one, mutual verification of digital and physical models can be performed. In the comparison, the upstream pressures under four gradients are taken. They are 0 MPa, 10 MPa, 20 MPa, and 30 MPa, respectively. The time to get the last data point is 30,000 s, 18,000 s, 25,200 s, and 25,200 s, so the numerical simulation is the same as the measured time in the laboratory. A comparison of the recovery ratio of the theoretical and measured curves under upstream pressure was made and shown in Figure 7. From the comparison of Figure 7, it can be observed that the established numerical model is very close to the actual core physical properties. The final exhaustion yields are basically coincidental, which indicates the effectiveness of digital core technology.
(1) The theoretical and measured recovery factors have a large error in the early stage of exhaustion mining. However, as the mining time increases, the theoretical model and the measured recovery data begin to approach. At the end of exhaustion, the theoretical and measured results are basically coincident. The reason for this is that the digital core ignores interfacial tension between gas and oil, thus, resulting in a different flow state at the end of a digital one compared to a physical one.
(2) The relative error between the numerical model and physical experiment is not large, but the four sets of contrast curves show the high pressure of the previous theoretical model. The reason is that, in the early stage of exhaustion mining, the pressure drop of the theoretical model has not yet propagated upstream. Therefore, the upstream pressure remains unchanged during the next period of mining start.
Additionally, during the numerical simulation, an interesting phenomenon in which some cracks have both open and closed states is observed. This phenomenon gives a microscopic interpretation clearly of the starting pressure gradient of a tight reservoir. As From the comparison of Figure 7, it can be observed that the established numerical model is very close to the actual core physical properties. The final exhaustion yields are basically coincidental, which indicates the effectiveness of digital core technology.
(1) The theoretical and measured recovery factors have a large error in the early stage of exhaustion mining. However, as the mining time increases, the theoretical model and the measured recovery data begin to approach. At the end of exhaustion, the theoretical and measured results are basically coincident. The reason for this is that the digital core ignores interfacial tension between gas and oil, thus, resulting in a different flow state at the end of a digital one compared to a physical one.
(2) The relative error between the numerical model and physical experiment is not large, but the four sets of contrast curves show the high pressure of the previous theoretical model. The reason is that, in the early stage of exhaustion mining, the pressure drop of the theoretical model has not yet propagated upstream. Therefore, the upstream pressure remains unchanged during the next period of mining start.
Additionally, during the numerical simulation, an interesting phenomenon in which some cracks have both open and closed states is observed. This phenomenon gives a microscopic interpretation clearly of the starting pressure gradient of a tight reservoir. As shown in Figure 8a, a plan view of the model at a pressure differential of 10 MPa (b) corresponding to a plan view of the model at a pressure of 40 MPa. It can be observed that the micro-cracks within the marked range are in a closed state when the pressure difference of failure is 10 MPa. However, it is turned on under the depletion pressure difference of 40 MPa. It explains, at a microscopic level, why there is a starting pressure gradient in a tight sandstone reservoir in the case of a low pressure difference.
Furthermore, 16 calculations of the orthogonal test were performed according to the selected five factors and four levels, and the simulation results are shown in Table 4.
The calculation results on the orthogonal test were analyzed by considering the absolute deviation of each parameter, and the results are listed in Table 5.
It can be seen from the results in Table 5 that the change of pressure has the largest influence on tight oil recovery and the change of viscosity has the smallest one. Furthermore, the parameters in order of sensibility are pressure, length, permeability, time, and viscosity. shown in Figure 8a, a plan view of the model at a pressure differential of 10 MPa (b) corresponding to a plan view of the model at a pressure of 40 MPa. It can be observed that the micro-cracks within the marked range are in a closed state when the pressure difference of failure is 10 MPa. However, it is turned on under the depletion pressure difference of 40 MPa. It explains, at a microscopic level, why there is a starting pressure gradient in a tight sandstone reservoir in the case of a low pressure difference. Furthermore, 16 calculations of the orthogonal test were performed according to the selected five factors and four levels, and the simulation results are shown in Table 4. The calculation results on the orthogonal test were analyzed by considering the absolute deviation of each parameter, and the results are listed in Table 5.

Conclusions
The research method, which combines a physical experiment with digital core methods, can greatly shorten the experimental period. In addition, after verifying the accuracy of the 3D digital core, the scale shape and production system of the research object can be arbitrarily changed. By this way, the approach greatly increases the flexibility of the research and all of them complete the experiment in a short period. This promising research method is applied to the matrix failure process of tight reservoirs, which can reveal the seepage mechanism of a tight sandstone reservoir as follows.