Abstract
Water breaks through along fractures is a major concern in tight sandstone reservoirs with a bottom aquifer. Analytical models fail to handle the three-dimensional two-phase flow problem for partially penetrating inclined fractures, so time-consuming numerical simulation are often used for this problem. This paper presents an efficient semianalytical model for this problem considering three-dimensional fractures and two-phase flow. In the model, the hydraulic fracture is handled discretely with a numerical discrete method. The three-dimensional volumetric source function in real space and superposition principle are employed to solve the model analytically for fluid flow in the reservoir. The transient flow equations for flow in three-dimensional inclined fractures are solved by the finite difference method numerically, in which two-phase flow and stress-dependent properties are considered. The eventual solution of the model and transient responses are obtained by coupling the model for flow in the reservoir and discrete fracture dynamically. The validation of the semianalytical model is demonstrated in comparison to the solution of the commercial reservoir simulator Eclipse. Based on the proposed model, the effects of some critical parameters on the characteristics of water and oil flow performances are analyzed. The results show that the fracture conductivity, fracture permeability modulus, inclination angle of fractures, aquifer size, perforation location, and wellbore pressure drop significantly affect production rate and water breakthrough time. Lower fracture conductivity and larger inclination angle can delay the water breakthrough time and enhance the production rate, but the increment tends to decline gradually. Furthermore, water breakthrough will occur earlier if the wellbore pressure drop and aquifer size are larger. Besides, the stress sensitivity and perforation location can delay the water breakthrough time.
1. Introduction
In recent years, tight sandstone reservoirs have been discovered as unconventional oil and gas resources with great potential, and hydraulic fracturing technology has been very critical for improving the well productivity of this kind of reservoir [1–4]. The hydraulic fracture can be inclined if the least principal stress is not perpendicular or parallel to the formation plane [5–8]. Due to the complex fluid distribution characteristics, some tight sandstone reservoirs have bottom aquifer, and water breakthroughs along the inclined fracture are quite common in developing tight sandstone reservoirs [9, 10]. Water invading along the hydraulic fracture will seriously influence well productivity, leading to a rapid decline in bottom hole pressure and production rate. It is difficult to accurately predict the water breakthrough time and analyze two-phase production data for tight sandstone reservoirs with a bottom aquifer. Complex fracture characteristics and two-phase flow behavior add to this challenge. There is no work investigating the transient behavior for this kind of well in tight sandstone reservoirs.
Much research has been studied on the transient behavior of fractured wells in tight sandstone reservoirs. However, most of them are limited to the single-phase flow problem [11–14]. Besides, to simplify the mathematical model, the proposed analytical and semianalytical methods assume that the fluid flow in the fracture is 1D flow. The flux variation along the inclined direction in the fracture is seldom considered in their work. Also, most of the models assume that the fracture both is vertical and fully penetrates the formation. In fact, the microseismic map’s interpretation result illustrates that a sub-horizontal fracture is induced after the fracturing treatment, and the fracture partially penetrates the formation [15–17]. Many theoretical models have recently been proposed to investigate the two-phase flow performance of wells in unconventional reservoirs. Empirical, reservoir simulation and analytical and semianalytical methods have been widely used to predict water breakthrough time and analyze two-phase production data. The empirical technique is too simple and commonly used in fields [18, 19]. In recent studies, numerical simulation methods have been employed to study the two-phase behavior of fractured wells [20–22]. However, to ensure simulation accuracy, the local grid refinement method is generally applied to model the fracture, so the computing efficiency is not satisfying. Analytical models are mainly used to calculate the water influx rate and predict water breakthrough time in bottom water-driven reservoirs [23, 24]. The complex potential theory and potential superposition principle are the common methods to obtain the solution of analytical models. However, these models are generally inappropriate for fractured wells. Semianalytical models capture the two-phase flow behavior in the fracture and have much higher computational efficiency than numerical simulation methods, which have recently attracted increasingly attention [25–29]. However, these methods are mainly employed for simulate early time flowback of fractured wells. Besides, the existing semianalytical models still have several problems in accurately capturing the three-dimensional fracture characteristics and its stress-sensitive. Although much work has been devoted to the two-phase flow models, there are also many challenges in developing semianalytical solutions for fractured wells in tight sandstone reservoirs with a bottom aquifer.
This paper proposed a semianalytical model to analyze production data and predict water breakthrough time in tight sandstone reservoirs with a bottom aquifer. The three-dimensional fracture characteristics and two-phase flow behavior and stress-dependent fracture permeability are all incorporated in the model. The fluid flow in the fracture system is numerically characterized by the finite difference method. In contrast, the fluid flow in the matrix system is analytically characterized based on the three-dimensional source function theory. The solution of the model is solved by coupling the transient flow in the fracture and the matrix. The proposed model is shown to efficiently simulate the behavior of tight sandstone reservoirs with a bottom aquifer. The accuracy of the proposed semianalytical model is validated by using numerical simulation, and the influences of several parameters on production rates and water breakthrough time are also analyzed. The novelty of the new model is in the ability to semianalytically obtain the performance of production wells in tight sandstone reservoirs with a bottom aquifer. Furthermore, it provides an efficient method for modeling bottom water coning, combining three-dimensional fracture characteristics and two-phase flow.
2. Physical Model
This paper mainly focuses on the oil and water production performance of a vertical well intersecting a rectangle-shaped partially penetrating fracture in tight sandstone reservoirs with a bottom aquifer. As shown in Figure 1, the bottom aquifer is fractured and water breaks through along the fracture during depletion development. Because the fracture permeability is much larger than the matrix permeability in tight sandstone reservoirs, it is assumed that water is mobile only in the fracture system to simplify the model. Some simplifying physical model assumptions for the derivation of the governing equation are listed as follows: (1)The reservoir is box-shaped, and the external boundaries of the top and side are assumed to be closed, and the reservoir is assumed to be homogeneous(2)The fracture is located at the center of the reservoir, and the inclination angle of the fracture is (3)This model considers the three-dimensional fracture characteristics, and the fluid flow in the fracture is considered a 2D flow(4)Darcy seepage is applicable in both the fracture and matrix systems, and stress-dependent fracture permeability is considered(5)The vertical well produces constant bottom hole pressure and fluid flow into the wellbore only through the fracture
3. Mathematical Model
As illustrated in the physical model, two distinct flow processes are governed by different mechanisms: partially penetrating inclined fracture flow and matrix flow. In this section, we first describe the theoretical model for fluid flow in the three-dimensional inclined fracture, and the finite difference method is used to obtain the solution of this system. Next, a three-dimensional volumetric source function in real space and superposition principle is applied to solve unsteady state flow in the matrix system. Then, the eventual solution of the model and transient responses are obtained by coupling the two systems with continuity conditions.
3.1. Model for Flow in the Fracture System
With consideration of the stress-sensitivity effect and oil-water two-phase flow behavior, the governing equation of a three-dimensional fracture system can be described as follows:
For the oil phase,
For the water phase,
The relationship between oil and water saturation is given by
The initial condition is
Under constant bottom-hole flowing pressure condition, a well production equation is given by [30] where and are the relative permeability of oil and water phase, respectively; is the absolute permeability of fracture; and are the oil and water viscosity, respectively; represents fracture permeability modulus; is fracture porosity; and are oil and water volume coefficient, respectively; and are oil and water saturation, respectively; is fracture pressure; is initial reservoir pressure; is a source term representing the oil flux entering the fracture from the matrix at a point on the fracture surface per unit volume at standard conditions; and represent the oil and water production rate per unit volume at standard conditions, respectively; is bottom water influx during its coning through the fracture; and represent the direction of the fracture; is the transmissibility conversion factor and for the field unit system; and is wellbore radius and .
We first divide the fracture into several fracture panels, as illustrated in Figure 2 and then apply the finite difference approximation to Equations (1) and (2) at the time level . For each fracture panel (), the implicit finite difference form for the oil and water flow equations can be written as where and ; the transmissibility of oil and water phase between fracture panel and () at time level are defined as
Geometric transmissibility () is calculated using the harmonic average of the properties of the panels, which is expressed as where the geometric transmissibility between the fracture panel () and the interface is given by
The coefficients of , , , and in Equations (7) and (8) depend on the pressure and saturation change of panel () for the oil and water phases at the time level , which are defined as
Bulk volume of the fracture panel () can be expressed as
For a slightly compressible fluid, the coefficients of , , and in Equations (13) and (14) are given by where , , and are the water volume coefficient, oil volume coefficient, and fracture porosity at a reference pressure, respectively; , , and are the oil compressibility, water compressibility, and fracture compressibility, respectively.
Applying the finite difference approximation to the production boundary conditions, Equations (5) and (6) can be rewritten as where is the wellbore pressure and and are the transmissibility of oil and water phase between fracture panel and the wellbore at time level , which are given by
In this paper, the flowing material balance method is used to calculate the water influx rate () in Equation (8). Using the flowing material balance method, the average pressure can be calculated by
Assuming that fracture panels are directly connected to the bottom aquifer, then the cumulative water influx rate is expressed as
Inserting Equation (26) into Equation (25) yields:
The water influx rate of each fracture panel connecting to the bottom aquifer is given by where is the aquifer size, is the initial pressure of the bottom aquifer, represents the water invasion index, and is the water influx rate between the fracture panels and bottom aquifer at each time level.
Substituting the production condition Equation (21) into Equation (7), the implicit finite difference form for the oil flow equation can be rearranged as
Taking Equation (22) and Equation (28) into Equation (8), we can rewrite the water flow equation as
Equations (29) and (30) are written for fracture panels () at time step ; then, we apply these two equations to fracture panels and arrange these fracture flow equations into a matrix format yield: where the matrix and vector in Equation (31) are given by
There are equations to be solved simultaneously for the unknowns (oil saturation and fracture panel pressure) in Equation (31). Here, it should be noted that the system in Equation (31) is nonlinear for the reason that the terms in matrix and vector are the functions of pressure and saturation of the fracture panels.
3.2. Model for Flow in the Matrix System
In our work, the three-dimensional macro fracture is discretized into panels, and the panels are considered plane sources and the withdrawal of fluid into the sinks causes a pressure drop in the matrix system. In this section, the plane source function is obtained by integrating the three-dimensional volumetric source function along the horizontal and inclined direction of the fracture panel, and the point source function is derived by solving the matrix flow model.
For an isotropic permeability reservoir, the matrix flow equation can be expressed as
Initial condition
The outer boundary is assumed to be closed and given by
The top boundary is assumed to be closed
And the bottom boundary is assumed to produce at constant pressure, which is given by
According to Jia et al.’s [31] study, the solution for an isotropic permeability reservoir with closed top and bottom boundaries can be obtained by using the orthogonal transformation method. We derive the solution of Equation (37) with the orthogonal transformation method analytically, and the pressure solution for a three-dimensional volumetric source is expressed as follows: where is the current pressure of the matrix system; is the initial pressure of the matrix system; is the matrix permeability; is the matrix porosity; is the total compressibility of matrix system; is diffusivity; is the unit conversion factor; , , and are dimensions of the reservoir; is the center position of the plane source; is a point source influx from the matrix to the fracture; and is the withdrawal rate at standard condition.
As we can see in Figure 3, the fracture panel is regarded as a plane source. The pressure response at an arbitrary position () at time due to the contribution of a continuous plane source with a withdrawal rate of can be calculated by integrating the three-dimensional volumetric source solution along the horizontal and inclined direction of the fracture panel, which is given by where is the length of the fracture panel (); is the height of the fracture panel (); and , , are the center position of the fracture panel ().
Equation (44) is the pressure response caused by the fracture panel (). Then, for the system in Figure 3, the pressure drop at the point () caused by all fracture panels can be calculated by using the superposition principle, and the equation is given by where where is the oil flow rate from the matrix to the fracture panel per unit fracture area at time step and is the plane source function of fracture panel .
To simplify the form of the equations, let denote
Taking Equations (46) and (47) into Equation (45) yields
Then, the pressure response of fracture segment g caused by all fracture panels at time level can be calculated by
We arrange Equation (49) into a matrix format yields where the matrix and vector are given as
It should be noted that there are equations to be solved simultaneously for the unknowns in Equation (50).
3.3. Coupled Model and Solution
The continuity conditions of pressure and flux along the fracture plane are given by
Based on the continuity condition, we combine the two-phase flow equations in the fracture and single-phase flow equations in the matrix to develop a system of equations to characterize the transient flow behavior of the three-dimensional inclined fracture.
Equation (58) has equations and unknowns, containing the pressure () and saturation () in fracture panels, and the flow rate from the matrix to the fracture panels (). As described above, the terms in matrix and vector are functions of pressure and saturation of fracture panel, and both the saturation and pressure in the fracture panels are unknown parameters. To linearize the equations, we use the upstream weighting method to linearize the nonlinear terms (transmissibility and compressibility ) in space and the simple iteration method is used for linearization in time. Then, the equations can be solved by a computer program.
4. Results and Discussion
4.1. Model Validation
In this section, the proposed semianalytical method is benchmarked against the commercial numerical simulator Eclipse. As shown in Figure 4, a numerical model is built to simulate the transient behavior of fractured wells in tight sandstone reservoirs with a bottom aquifer. The dimensions of the reservoir are , and the grid dimensions used to discrete the reservoir are . To ensure the accuracy of computing, we apply the local grid refinement approach to model the fracture. The number of discretized fracture panels in the semianalytical model is 96. The input parameters for the proposed model and Eclipse are identical, in which the length, width, and height of the fracture are 120 m, 0.01 m, and 10 m, respectively, the fracture conductivity is 5 mD·m, the inclination angle of the fracture is 75°, the fracture permeability modulus is 0.05 MPa-1, the aquifer size is the initial reservoir pressure is 30 MPa, and the oil viscosity is 0.3 mPa·s. The other input parameters are provided in Table 1. The oil-water relative permeability curves in the hydraulic fracture system are shown in Figure 5.
Figure 6 shows the validation results computed by our model and the fully numerical model. As illustrated in Figure 6, the solutions obtained from the proposed method and the Eclipse are almost the same. The matching error is less than 3%. Although there is a slight mismatch in production profiles and water breakthrough time, the match is acceptable within the engineering error. The results suggest that the proposed model can model fluid flow in three-dimensional large-scale fractures and predict the production profiles and water breakthrough time of fractured wells in tight sandstone reservoirs with a bottom aquifer. The calculation speed of the two methods is also compared in our work. In this validation case, the computation time is less than 15 min in our model, while 80 min for the commercial simulator, which indicates the new model is efficient.
(a) Oil production rate
(b) Water production rate
4.2. Sensitivity Analysis
Based on the proposed model, the influences of the following parameters on oil and water production performance and water breakthrough time are analyzed: fracture conductivity, fracture permeability modulus, inclination angle of fracture, aquifer size, perforation location, and wellbore pressure drop. Table 2 shows the values of the parameters used for sensitivity analysis.
4.2.1. The Effect of the Fracture Conductivity
Figure 7 shows the influence of fracture conductivity on oil and water production rates calculated by the semianalytical model. As one can see from this figure, oil and water exhibit higher rates with larger fracture conductivity at an early time. Besides, the oil rate will decline significantly and the water influx rate increases rapidly once the water breaks through along the hydraulic fracture. It is also observed that the fracture conductivity has a relatively small effect on the water influx rate after 200 days and the water breakthrough time is delayed with lower fracture conductivity. It is possible because a lower fracture conductivity leads to a larger pressure drop in the fracture, increasing the flow resistance of both the oil and water phases. Thus, the water influx rate decreases and the water breakthrough time delays with a lower fracture conductivity.
(a) Oil production rate
(b) Water production rate
4.2.2. The Effect of the Stress-Dependent Permeability of the Fracture
It can be seen from Figure 8 that the effect of stress-dependent permeability of fracture on production rates is significant, especially in the early flow period. As the permeability modulus increases, oil and water exhibit lower production rates at the early time, but the oil production rate declines faster and the water influx rate increases rapidly at the late time. Besides, the water breakthrough time is delayed with a larger permeability modulus. This is because the fluid flow will be more difficult and with increasing permeability modulus. Consequently, the production rate declines and the water breakthrough time is delayed with the influence of stress-dependent permeability.
(a) Oil production rate
(b) Water production rate
4.2.3. The Effect of the Inclination Angle of Fracture
Figure 9 illustrates the effect of the inclination angle of the fracture on production rate and water breakthrough time. It can be observed that a larger inclination angle will reduce the decline of the oil production rate. Besides, the water breakthrough time will be delayed as the inclination angle increases. This is because the water breaks through along the fracture will be more difficult when increasing the inclination angle, so the larger the inclination angle, the later the water breakthrough time will be.
(a) Oil production rate
(b) Water production rate
4.2.4. The Effect of the Aquifer Size
Figure 10 depicts the oil and water production dynamic curves with different aquifer sizes. The aquifer size is given 5 m3, 7 m3, and 9 m3 for the case study in this section. As illustrated in Figure 10, water breaks through along the fracture will seriously decrease the oil production rate with a larger aquifer size. It can also be seen from Figure 10 that the water breakthrough time of the three aquifer sizes is 5 days, 50 days, and 120 days, respectively, which indicates that a smaller aquifer size can efficiently delay the water breakthrough time.
(a) Oil production rate
(b) Water production rate
4.2.5. The Effect of the Perforation Location
The effect of the perforation location on production rate and water breakthrough time is analyzed in this part. We denote by the perpendicular distance from the perforation location to the bottom aquifer surface. The modelling results are shown in Figure 11; it should be noticed that the perforation location has a significant effect on production rate and water breakthrough time. With a larger , the oil exhibits a higher production rate in the early period and the water breakthrough time is clearly delayed. This shows that a proper perforation position is excellent for controlling the water breakthrough time and improving the production efficiency during the development of tight sandstone reservoirs with a bottom aquifer.
(a) Oil production rate
(b) Water production rate
4.2.6. The Effect of the Wellbore Pressure Drop
Figure 12 shows the oil and water production rates for different producing pressure drops. As shown from this figure, oil and water exhibit higher rates as the producing pressure drop increases at the early time. It is also observed that the water breakthrough time is delayed with a smaller wellbore pressure drop. This is because the producing pressure drop determines the oil recovery rate and stable production period. A larger wellbore pressure drop will lead to a shorter stable production time and an earlier water breakthrough time. It is believed that rational wellbore pressure drop is an effective way to postpose water breakthrough time and reduce production decline in tight sandstone reservoirs with a bottom aquifer.
(a) Oil production rate
(b) Water production rate
5. Summary and Conclusions
This paper’s main contribution is the provision of a semianalytical method to efficiently simulate two-phase flow behavior in three-dimensional large-scale fractures and accurately predict the production profiles and water breakthrough time of fractured wells in tight sandstone reservoirs with a bottom aquifer. The main conclusions of this work are as follows: (i)The proposed semianalytical model can be used to evaluate the two-phase flow performance of a three-dimensional incline fracture with satisfying precision, and much computational time can be saved compared with fully numerical simulation(ii)The important properties of the fracture system, including fracture conductivity, inclination angle, and stress-dependent permeability, have great effects on the production rates and water breakthrough time. A larger fracture permeability modulus and inclination angle of fracture can delay the water breakthrough time and enhance the oil production rate. A more significant fracture conductivity will lead to an earlier water breakthrough time(iii)Aquifer size has a significant impact on the production rate, and a smaller aquifer size can efficiently delay the water breakthrough time(iv)The perforation location and wellbore pressure drop play an important role in production performance. Reasonable perforation position and producing pressure difference are excellent for delaying the water breakthrough time and improving the production efficiency during the development of tight sandstone reservoirs with a bottom aquifer
Nomenclature
| : | Oil formation volume factor (Sm3/m3) |
| : | Oil flow rate entering the fracture from the matrix per unit volume at standard conditions (1/d) |
| : | Water formation volume factor (Sm3/m3) |
| : | Oil production rate per unit volume at standard conditions (1/d) |
| : | Fracture compressibility (MPa-1) |
| : | Water production rate per unit volume at standard conditions (1/d) |
| : | Oil compressibility (MPa-1) |
| : | Bottom water influx during its coning through the fracture (1/d) |
| : | Total compressibility of matrix system (MPa-1) |
| : | Radial distance of side external boundary (m) |
| : | Water compressibility (MPa-1) |
| : | Wellbore radial (m) |
| : | Geometric transmissibility of fracture panel and panel () |
| : | Oil saturation (dimensionless) |
| : | Formation thickness (m) |
| : | Water saturation (dimensionless) |
| : | Water invasion index |
| : | Time (day) |
| : | Permeability of fracture (mD) |
| : | Transmissibility of oil phase between fracture panel and () |
| : | Initial permeability of fracture (mD) |
| : | Transmissibility of water phase between fracture panel and () |
| : | Permeability of matrix (mD) |
| : | Transmissibility of oil phase between fracture panel and the wellbore |
| : | Relative permeability to oil (dimensionless) |
| : | Transmissibility of water phase between fracture panel and the wellbore |
| : | Relative permeability to water (dimensionless) |
| : | Bulk volume of fracture panel () (m3) |
| : | Number of panels of the discretized fracture |
| : | Aquifer size (m3) |
| : | Number of panels of the discretized fracture directly connected with the bottom aquifer |
| : | Cumulative water influx rate (m3) |
| : | Pressure of the bottom aquifer (MPa) |
| : | , , and coodinates |
| : | Fracture pressure (MPa) |
| : | Center position of the plane source |
| : | Initial formation pressure (MPa) |
| : | Reservoir dimension (m) |
| : | Matrix pressure (MPa) |
| : | Height of fracture panel () |
| : | Wellbore pressure (MPa) |
| : | Length of fracture panel () |
| : | Withdrawal rate at standard condition (m3/d) |
| : | Time step size (d). |
| ϕ: | Porosity (dimensionless) |
| : | Fluid viscosity (mPa·s) |
| : | Inclination angle of fracture (°) |
| : | Differential operator |
| : | Conversion factor |
| : | Difference calculation |
| : | Conversion factor |
| : | Permeability modulus (MPa-1) |
| : | Direction of the fracture |
| : | Direction of the fracture. |
| : | Time level |
| : | Average parameter. |
| e: | External |
| : | Oil |
| : | Fracture |
| : | Relative |
| : | Initial |
| : | Standard condition |
| : | Matrix |
| : | Wellbore. |
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge that this study was partially funded by the CNOOC Research Institute Co., Ltd., Science and Technology Projects (No. CNOOC-YXKJ-QZJC-TJ-2020-01).