A Semianalytical Model for Simulating Fluid Flow in Tight Sandstone Reservoirs with a Bottom Aquifer

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 threedimensional 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.


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][2][3][4]. The hydraulic fracture can be inclined if the least principal stress is not perpendicular or parallel to the formation plane [5][6][7][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][12][13][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][16][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 twophase 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][21][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][26][27][28][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 threedimensional 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 threedimensional 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 threedimensional fracture characteristics and two-phase flow.

Physical Model
This paper mainly focuses on the oil and water production performance of a vertical well intersecting a rectangleshaped partially penetrating fracture in tight sandstone reser-voirs 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

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.

Geofluids
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 k ro and k rw are the relative permeability of oil and water phase, respectively; k Fi is the absolute permeability of fracture; μ o and μ w are the oil and water viscosity, respectively; γ F represents fracture permeability modulus; ϕ F is fracture porosity; B o and B w are oil and water volume coefficient, respectively; S o and S w are oil and water saturation, respectively; p F is fracture pressure; p i is initial reservoir pressure; ẽ q SCFo 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; ẽ q SCWo and ẽ q SCWw represent the oil and water production rate per unit volume at standard conditions, respectively; ẽ q SCwI is bottom water influx during its coning through the fracture; ξ and η represent the direction of the fracture; β is the transmissibility conversion factor and β = 0:0864 for the field unit system; and r W is wellbore radius and r eq = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0:14ðΔl 2 F + Δh 2 F Þ q . 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 n + 1. For each fracture panel (i, j), the implicit finite difference form for the oil and water flow equations can be written as where i = 1, ⋯N x and j = 1, ⋯N y ; the transmissibility of oil and water phase between fracture panel l and (i, j) at time level n + 1 are defined as T n+1 Geometric transmissibility (G l,ði,jÞ ) is calculated using the harmonic average of the properties of the panels, which is expressed as where the geometric transmissibility between the fracture panel (i, j) and the interface is given by The coefficients of C n+1 (7) and (8)  3 Geofluids change of panel (i, j) for the oil and water phases at the time level n + 1, which are defined as Bulk volume of the fracture panel (i, j) can be expressed as For a slightly compressible fluid, the coefficients of (13) and (14) are given by 1 where B°w, B°o, and ϕ°F are the water volume coefficient, oil volume coefficient, and fracture porosity at a reference pressure, respectively; c o , c w , and c F 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 p n W is the wellbore pressure and T n+1 o l,W and T n+1 w l,W are the transmissibility of oil and water phase between fracture panel l and the wellbore at time level n + 1, which are given by T n+1 In this paper, the flowing material balance method is used to calculate the water influx rate (q n+1 SCwI i,j ) in Equation (8). Using the flowing material balance method, the average pressure can be calculated bŷ Assuming that N w 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 V p is the aquifer size, p e is the initial pressure of the bottom aquifer, J w represents the water invasion index, and q SCwI g 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 (i, j) at time step n + 1; then, we apply these two equations 4 Geofluids to N F 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 2N F equations to be solved simultaneously for the 2N F 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 B and vector C are the functions of pressure and saturation of the fracture panels.

Model for Flow in the Matrix
System. In our work, the three-dimensional macro fracture is discretized into N F 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 p m is the current pressure of the matrix system; p _ m is the initial pressure of the matrix system; k m is the matrix permeability; ϕ m is the matrix porosity; c tm is the total compressibility of matrix system; χ is diffusivity; α is the unit conversion factor; x e , y e , and z e are dimensions of the reservoir; ðx w , y w , z w Þ is the center position of the plane source;q is a point source influx from the matrix to the fracture; and q SC 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 5 Geofluids (x, y, z) at time t due to the contribution of a continuous plane source with a withdrawal rate of q sc 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 cos nπy i,j y e cos nπy y e ( ) where Δl F i,j is the length of the fracture panel (i, j); Δh F i,j is the height of the fracture panel (i, j); and x i,j , y i,j , z i,j are the center position of the fracture panel (i, j).
Equation (44) is the pressure response caused by the fracture panel (i, j). Then, for the system in Figure 3, the pressure drop at the point (x o , y o, z o ) caused by all fracture panels can be calculated by using the superposition principle, and the equation is given by where whereq k SCm g is the oil flow rate from the matrix to the fracture panel g per unit fracture area at time step k and Sðx o , y o , z o , t n+1 − τ, x g , y g , z g Þ is the plane source function of fracture panel g.
To simplify the form of the equations, let U n+1,k o,g 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 n + 1 can be calculated by

Geofluids
It should be noted that there are N F equations to be solved simultaneously for the N F unknowns in Equation (50).

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 twophase 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 threedimensional inclined fracture.
Equation (58) has 3N F equations and 3N F unknowns, containing the pressure (p n+1 F i,j ) and saturation (S n+1 w i,j ) in fracture panels, and the flow rate from the matrix to the fracture panels (q n+1 SCFo i,j ). As described above, the terms in matrix B and vector C are functions of pressure and saturation of frac-ture 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 T and compressibility C) in space and the simple iteration method is used for linearization in time. Then, the equations can be solved by a computer program.

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 1000 × 1000 × 100 m, and the grid dimensions used to discrete the reservoir are 20 × 20 × 2 m. 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 1 × 10 7 m 3 , the initial reservoir pressure is 30 MPa, and the oil viscosity is    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 threedimensional 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.

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.  Figure 7 shows the influence of fracture conductivity on oil and water production rates calculated by the semianalytical model. As    9 Geofluids 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.

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.

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.

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 7 × 10 5 m 3 , 7 × 10 7 m 3 , and 7 × 10 9 m 3 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.

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 Δh 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 Δh, 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.
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  11 Geofluids 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.

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