*Article* **Numerical Simulation of Multiarea Seepage in Deep Condensate Gas Reservoirs with Natural Fractures**

**Lijun Zhang 1, Wengang Bu 2,\*, Nan Li 1, Xianhong Tan <sup>1</sup> and Yuwei Liu <sup>2</sup>**


**\*** Correspondence: binbeiu@163.com

**Abstract:** Research into condensate gas reservoirs in the oil and gas industry has been paid much attention and has great research value. There are also many deep condensate gas reservoirs, which is of great significance for exploitation. In this paper, the seepage performance of deep condensate gas reservoirs with natural fractures was studied. Considering that the composition of condensate gas changes during the production process, the component model was used to describe the condensate gas seepage in the fractured reservoir, modeled using the discrete fracture method, and the finite element method was used to conduct numerical simulation to analyze the seepage dynamic. The results show that the advancing speed of the moving pressure boundary can be reduced by 55% due to the existence of threshold pressure gradient. Due to the high-speed flow effect in the near wellbore area, as well as the high mobility of oil, the condensate oil saturation near the wellbore can be reduced by 42.8%. The existence of discrete natural fractures is conducive to improving the degree of formation utilization and producing condensate oil.

**Keywords:** condensate gas reservoir; condensate oil; discrete fracture

#### **1. Introduction**

The development of petroleum industrial technologies and the quick growth in demand for oil and gas resources have meant that unconventional oil and gas resources have gradually become a hot spot in the industry, and the development of deep oil and gas reservoirs is gradually increasing. Of the newly discovered recoverable reserves of global condensate oil in 2020, deep and ultra-deep condensate oil accounts for 95% [1], and 100 billion cubic meters of large, deep condensate gas fields exist in Bohai [2]. However, problems such as their great depth and the complex hydrocarbon fluid phase state bring difficulties to their development.

Deep, fractured condensate gas reservoirs have obvious retrograde condensate characteristics in the development process, forming a multiscale, multiphase flow system with matrix fractures. Syed et.al [3] reported a comprehensive review of the complex nature of unconventional reservoirs. Apart from tight reservoirs, deep natural gas reservoirs and geo-pressurized zones, which usually exist under high pressure and high temperature, are also referred to as unconventional reservoirs. The low, dual porosity distribution and ultra-low permeability of those reservoirs are evident. According to an evaluation of core samples [4], the pore size of deep reservoirs can be tiny and heterogeneous. The matrix of deep formation is tight and has a low permeability, and the permeability difference between it and its fractures is large, which causes different flow characteristics to form [5,6]. For low-permeability reservoirs, many research results have shown the impact of threshold pressure gradient on seepage [7–10]. Many studies focus on oil, water, or oil–water phases [11,12], while practice has proved that there is a threshold pressure gradient in the development of low-permeability gas reservoirs [13–15]. Tian et al. [16] pointed out that the threshold pressure gradient of a tight sandstone condensate gas reservoir is positively related to the liquid saturation, and the threshold pressure gradient is significantly

**Citation:** Zhang, L.; Bu, W.; Li, N.; Tan, X.; Liu, Y. Numerical Simulation of Multiarea Seepage in Deep Condensate Gas Reservoirs with Natural Fractures. *Energies* **2023**, *16*, 10. https://doi.org/10.3390/ en16010010

Academic Editors: Xingguang Xu, Kun Xie and Yang Yang

Received: 16 November 2022 Revised: 12 December 2022 Accepted: 14 December 2022 Published: 20 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

increased when fractures have not developed. Martyushev et.al [17] reported that the properties of fractured reservoirs cannot fully recover after damage caused by the stress state increasing. Therefore, there is a risk of a production decline due to reservoir pressure drops and fracture closures [18]. For the study of fractured reservoirs, the dual medium model, equivalent model, and discrete fracture model are mainly used. Due to its advantages with regard to simulation and accuracy, the discrete fracture model is widely used [19–21]. The discrete fracture method is popular for dealing with the matrix-fracture flow of shale gas. Wang et al. [22] established a discrete fracture model to analyze the impact of hydraulic fractures on shale gas development, taking into account non-Darcy flow, and this model can be used for ultra-low permeability reservoirs. Zhao et al. [23] established a discrete fracture network including micro fractures and hydraulic fractures, coupled with the gas exchange between the matrix and fracture, and analyzed the impact of fracture geometry and distribution on production.

The purpose of this paper is to study the seepage performance of deep condensate gas reservoirs with natural fractures and to establish a mathematical model for this subject. The model can be applied to calculate the formation pressure distribution and the expansion dynamic of the moving pressure boundary at different development stages, and to analyze the changes in condensate oil in different areas of the formation and the impact of natural discrete fractures on the seepage field.

#### **2. Model and Methods**

#### *2.1. Physical Model*

The seepage problem present in deep fractured condensate gas reservoirs was considered. The porosity and permeability of the formation were assumed to be homogeneous, and the temperature was constant. The seepage in the formation was unsteady; that is, the formation pressure was a function of time and position. The physical model of formation with natural fractures was built using the discrete fracture method, as shown in Figure 1. The model was divided into matrix area (gray entity) and fracture area (black line).

**Figure 1.** Formation model diagram.

#### *2.2. Mathematical Description*

The seepage field structure of condensate gas reservoirs can be considered as fitting the four-zone flow model [24], as shown in Figure 2. From the boundary to the bottom of the well, the formation is divided into a single-phase gas area, a transition area, an outer two-phase flow area, and an inner two-phase flow area. The formation pressure in the single-phase gas area is higher than the dew pressure, and no condensate is produced. The formation pressure in the transition zone decreases gradually, and there is no flowing

condensate oil, only gas-phase flows. Gas and oil two-phase flow exists in both the outer and inner two-phase flow areas, which are divided by Reynolds number.

**Figure 2.** Schematic diagram of four-zone flow model.

#### 2.2.1. Continuity Equation

The flow of condensate gas belongs to a multicomponent, multiphase flow, and the continuity equation of component *i* is:

$$\nabla \cdot \left( \mathbf{x}\_{i} \rho\_{o} v\_{o} + y\_{i} \rho\_{\mathcal{S}} v\_{\mathcal{S}} \right) + \frac{\partial}{\partial t} \left[ \phi\_{\hat{\jmath}} \left( \mathbf{x}\_{i} \rho\_{o} \mathbf{S}\_{o} + y\_{i} \rho\_{\mathcal{S}} \mathbf{S}\_{\mathcal{S}} \right) \right] - q\_{i} = 0,\tag{1}$$

where *xi* and *yi* are the mole fractions of component *i* in the oil phase and gas phase respectively; *<sup>t</sup>* is the time; *<sup>φ</sup>* is the porosity; *<sup>ρ</sup>* is the fluid density; *<sup>S</sup>* is the saturation; *v* is the velocity vector; *qi* is the source and sink item; subscript *j* = *m* and *f*, representing matrix and fracture, respectively; and subscripts *o* and *g* denote oil and gas phases, respectively.

The dimension of the discrete fracture model can be reduced [25], and the continuity equation of the fracture can be written as:

$$\nabla \cdot \left[ d\_f \left( \mathbf{x}\_i \rho\_o \mathbf{v}\_o + \mathbf{y}\_i \rho\_\mathcal{g} \mathbf{v}\_\mathcal{g} \right) \right] + d\_f \frac{\partial}{\partial t} \left[ \phi\_\mathcal{j} \left( \mathbf{x}\_i \rho\_o \mathbf{S}\_o + \mathbf{y}\_i \rho\_\mathcal{g} \mathbf{S}\_\mathcal{g} \right) \right] - d\_f q\_i = 0,\tag{2}$$

where *df* is the fracture width.

#### 2.2.2. Motion Equation

In the single-phase gas area, the transition area, and the outer two-phase flow area, fluid flow can be described by Darcy's law modified by threshold pressure gradient:

$$v\_n = -k\_m \frac{k\_j}{\mu\_n} (\nabla p\_n - G)\_\prime \tag{3}$$

where *v* is the velocity; *kr* is the relative permeability; *μ* is the viscosity; *p* is the pressure; *G* is the threshold pressure gradient; and subscript *n* = *o* and *g*.

In the inner two-phase flow area, the influence of the near-well high-speed non-Darcy effect shall be considered [26]. The non-Darcy term was introduced to describe the flow law and is written as:

$$-\left(\nabla p\_n - \mathcal{G}\right) = \frac{\mu\_n}{k\_{rn}k\_j}v\_n + \beta \rho\_n v\_{n\prime}^2\tag{4}$$

where *β* is the inertial coefficient. The inertial coefficient reflects the deviation of seepage from Darcy's law. When oil and gas two-phases coexist in the near-wellbore area, the inertial coefficient can be calculated using [27]:

$$\beta = \frac{3.3 \times 10^{-9}}{\phi^{0.75} k^{1.25} s\_{\text{\textdegree}}^{4.5}},\tag{5}$$

The inner and outer two-phase flow areas are divided by Reynolds number, which can be calculated by Kajiahoff Equation (6), and the critical value is 0.2~0.3. The area where Reynolds number is greater than the critical value is the inner two-phase flow area, and the opposite is the outer two-phase flow area.

$$Re = \frac{v\rho\sqrt{k}}{17.5\mu\phi^{1.5}},\tag{6}$$

2.2.3. Constraint Equation

The constraint condition of the component is:

$$
\sum \mathbf{x}\_i = \sum \mathbf{y}\_i = \mathbf{1},
\tag{7}
$$

The constraint condition of saturation is:

$$S\_o + S\_{\mathcal{K}} + S\_w = 1,\tag{8}$$

where *Sw* is the connate water saturation.

The total molar amount of component *i* is:

$$z\_i = L\mathbf{x}\_i + V\mathbf{y}\_{i\prime} \tag{9}$$

where *L* and *V* are the mole fractions of the oil phase and gas phase in the system, respectively, as shown in Equation (10).

$$L = \frac{\rho\_o S\_o}{\rho\_o S\_o + \rho\_\mathcal{g} S\_\mathcal{g}}, V = \frac{\rho\_\mathcal{g} S\_\mathcal{g}}{\rho\_o S\_o + \rho\_\mathcal{g} S\_\mathcal{g}} \, \tag{10}$$

2.2.4. Phase Equilibrium Equation

$$f\_i^L = f\_i^V{}\_{\;\prime} \tag{11}$$

where *f <sup>L</sup> <sup>i</sup>* and *<sup>f</sup> <sup>V</sup> <sup>i</sup>* represent the liquid fugacity and gas fugacity of component *i*, respectively.

#### 2.2.5. State Equation

The state equation of rock is:

$$
\phi = \phi\_{int} [1 + \mathbb{C}\_s(p - p\_{int})],
\tag{12}
$$

where *Cs* is the rock compression coefficient; *φint* is the original porosity; and *pint* is the original pressure.

The fluid phase density is calculated using:

$$\rho\_n = \frac{pM\_n}{Z\_nRT} \tag{13}$$

where *M* is the average molecular weight; *Z* is the deviation factor; *R* is the ideal gas constant; and *T* is the temperature.

The PVT properties are described by Peng-Robinson equation [28]:

$$p = \frac{RT}{V\_{\rm r} - b} - \frac{a(T)}{V\_{\rm r}(V\_{\rm r} + b) + b(V\_{\rm r} - b)},\tag{14}$$

$$\alpha(T\_c) = 0.45724 \frac{\left(RT\_c\right)^2}{p\_c} \,\tag{15}$$

$$b = 0.0778 \frac{RT}{p\_c} \,\text{s}\tag{16}$$

where *Vr* is the molar volume; *pc* is the contrast pressure; and *Tc* is the contrast temperature.

2.2.6. Boundary Conditions and Initial Conditions

The closed outer boundary condition is:

$$-\mathfrak{n} \cdot (\rho \mathfrak{v}) = 0,\tag{17}$$

where *n* is the boundary normal vector.

The specified rate condition for the inner boundary is:

$$\left.v\right|\_{\Omega\_{\rm in}} = \frac{q\_0}{2\pi r\_w h'}\tag{18}$$

where Ω*in* represents the inner boundary, and *q*<sup>0</sup> is the production rate under formation condition.

The initial conditions can be written as:

$$p|\_{t=0} = p\_{int\prime} \tag{19}$$

$$S\_o|\_{t=0} = 0.\tag{20}$$

#### **3. Model Validation**

*3.1. Static Parameters*

The modeled typical deep fractured condensate gas reservoir had a depth of 4550.6 m, and a fracture width between 0.1 mm and 0.4 mm, with east-west strike fractures,. The static parameters of the formation are shown in Table 1.

**Table 1.** Formation statics parameters.


#### *3.2. Fluid Properties*

The content of C1, C2~C6+CO2, and C7+ in the condensate gas of the gas reservoir accounted for 71.97%, 18.76%, and 9.27%, respectively. The dew pressure was 47.56 MPa at the formation temperature of 178 ◦C. The P-T phase diagram of formation fluid is shown in Figure 3, in which 0%, 5%, 10%, and 20% represent the phase envelope lines at different condensate oil volumes. The components of condensate gas are complex, and calculations using the complete component will greatly reduce the calculation efficiency. Therefore, C2~C5 components and C6~C10 components were each combined into pseudo components. The PVT fitting was performed and dew pressure and relative volume fitting results are shown in Figures 4 and 5. The results show that the fitting effect of the pseudo components is good and can be used for calculation.

#### *3.3. Validation*

The equations were solved using the finite element numerical method. The measured production performance data of a production well within 480 days are matched with the above model, and the result is shown in Figure 6. In practice, the production rate is always changing. Therefore, in the early period of production, the measured production pressure data are unstable, and there are discrepancies between measured data and calculated data. After 390 days, the errors almost disappear, meaning that the model fits well.

**Figure 3.** P-T phase diagram of formation fluid.

**Figure 4.** Dew pressure fitting diagram.

**Figure 5.** Relative volume fitting diagram.

**Figure 6.** Comparison of measured and calculated pressure data.

#### **4. Simulation Results and Discussion**

In this section, a production simulation was conducted using the specified rate condition with a value of 14.27 × <sup>10</sup><sup>4</sup> <sup>m</sup>3/d from the 480th day to the 1500th day.

#### *4.1. Formation Pressure Distribution*

Figure 7 shows the radial distribution nephogram of formation pressure at different times. With the development of production, the pressure propagation range increases, meaning that the producing range of the formation gradually expands. As can be seen from the pressure curves shown in Figure 8, the advancing speed of the pressure moving boundary gradually decreases. From 500 days to 750 days, the moving pressure boundary advances 20 m, while from 1250 days to 1500 days, the moving pressure boundary advances only 9 m, showing a relative decrease of 55% compared with the former. The formation energy decreases when the production is conducted, and at the same time, the pressure propagation also needs to overcome the threshold pressure gradient, so the moving boundary moves slowly.

#### *4.2. Condensate Oil Distribution*

Figure 9 shows the distribution of condensate oil from the wellbore to the far side of formation at different times. Generally speaking, the closer to the wellbore, the lower the bottom hole pressure is, the higher the condensate oil saturation is. However, high saturation always brings high liquidity. With the development of production, the two-phase flow area expands, and the condensate oil is separated and gradually forms continuous phases and flows. In addition, the high-speed gas flow will carry part of condensate oil out of the reservoir, so condensate oil saturation near the wellbore decreases. When the production reaches 1500 days, the condensate oil saturation near the wellbore decreases from 0.07 to 0.04, forming a decrease of 42.8%. The pressure drop spreads to the far side of the formation, thus the transition area also gradually expands. From the 500th day to 1500th day, the transition area boundary expands from 94.4 m to 128.3 m.

#### *4.3. Influence of Discrete Natural Fractures*

Figure 10 shows the nephogram of formation pressure distribution under different fracture permeability. The discrete natural fractures and matrix are spatially continuous, so the pressure field is continuous. Compared with the initial fracture permeability, the advancing distance of the moving pressure boundary increases from 120 m to 140 m, with a relative increase of 16.7% after the fracture permeability doubled. Figure 11 shows the variation of condensate oil saturation in the presence of fractures with different permeabilities. With an increase in fracture permeability, the condensate oil near the wellbore is further

exploited, and the condensate oil saturation decreases. As the moving pressure boundary advances farther, the oil–gas two-phase area is further expanded, and more condensate oil is produced in the middle of the formation (101–135 m).

**Figure 7.** Formation pressure distribution at different times. (**a**) 500 d, (**b**) 1000 d, (**c**) 1250 d, (**d**) 1500 d.

**Figure 8.** Radial pressure distribution at different times.

**Figure 9.** Radial oil saturation distribution at different times.

**Figure 10.** Formation pressure distribution under different fracture permeabilities. (**a**) *kf* = 1.79 mD, (**b**) *kf* = 3.58 mD.

**Figure 11.** Radial condensate oil saturation under different fracture permeability.

#### **5. Conclusions**

In this paper, the component equation and discrete fracture method was applied to establish a seepage model of deep condensate gas reservoirs with natural fractures by considering the seepage laws of different flow areas, and the numerical simulation study was carried out to analyze the seepage dynamic. The following conclusions are drawn:


**Author Contributions:** Conceptualization, L.Z. and W.B.; Methodology, L.Z., W.B. and Y.L.; Validation, N.L. and X.T.; Writing—original draft, W.B., N.L. and Y.L.; Writing—review and editing, L.Z. and X.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
