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

Authors: Suran Wang; Yuhu Bai; Bingxiang Xu; Ling Chen; Wenlan Li; Yanzun Li
Journal: Geofluids
 (2021)

---

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


---

## Body

## 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 Figure1, 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 fractureFigure 1
Schematic diagram of the physical model used in this work.


## 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,(1)∂∂ξβkFie−γFpi−pFkroμoBo∂pF∂ξ+∂∂ηβkFie−γFpi−pFkroμoBo∂pF∂η+q~~SCFo+q~~SCWo=∂∂tϕF1−SwBo.For the water phase,(2)∂∂ξβkFie−γFpi−pFkrwμwBw∂pF∂ξ+∂∂ηβkFie−γFpi−pFkrwμwBw∂pF∂η+q~~SCWw+q~~SCwI=∂∂tϕFSwBw.The relationship between oil and water saturation is given by(3)Sw+So=1.The initial condition is(4)pFξ,η,t=0=pi.Under constant bottom-hole flowing pressure condition, a well production equation is given by [30]
(5)q~~SCWo=2πβkFikrowFe−γFpi−pFμoBolnreq/rWpF−pW,(6)q~~SCWw=2πβkFikrwwFe−γFpi−pFμwBwlnreq/rWpF−pW,where kro and krw are the relative permeability of oil and water phase, respectively; kFi 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; Bo and Bw are oil and water volume coefficient, respectively; So and Sw are oil and water saturation, respectively; pF is fracture pressure; pi 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 rW is wellbore radius and req=0.14ΔlF2+ΔhF2.We first divide the fracture into several fracture panels, as illustrated in Figure2 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
(7)∑l∈ψi,jTol,i,jn+1pFln+1−pFi,jn+1+qSCFoi,jn+1+qSCWoi,jn+1=Copi,jn+1pFi,jn+1−pFi,jn+Cowi,jn+1Swi,jn+1−Swi,jn,(8)∑l∈ψi,jTwl,i,jn+1pFln+1−pFi,jn+1+qSCWwi,jn+1+qSCwIi,jn+1=Cwpi,jn+1pFi,jn+1−pFi,jn+Cwwi,jn+1Swi,jn+1−Swi,jn,where i=1,⋯Nx and j=1,⋯Ny; the transmissibility of oil and water phase between fracture panel l and (i,j) at time level n+1 are defined as
(9)Tol,i,jn+1=Gl,i,jkrol,i,jn+11μoBol,i,jn+1,(10)Twl,i,jn+1=Gl,i,jkrwl,i,jn+11μwBwl,i,jn+1.Figure 2
Discretization of the three-dimensional fracture along both the horizontal axis and the inclined axis.Geometric transmissibility (Gl,i,j) is calculated using the harmonic average of the properties of the panels, which is expressed as
(11)Gl,i,j=γlγi,jγl+γi,j,where the geometric transmissibility between the fracture panel (i,j) and the interface is given by
(12)γi,j=βkFwFΔhFΔlFi,j.The coefficients ofCopi,jn+1, Cwpi,jn+1, Cowi,jn+1, and Cwwi,jn+1 in Equations (7) and (8) depend on the pressure and saturation change of panel (i,j) for the oil and water phases at the time level n+1, which are defined as
(13)Copi,jn+1=Vbi,jΔt1−Swi,jnϕFi,jn1Boi,j′+1Boi,jnϕFi,j′,(14)Cwpi,jn+1=Vbi,jΔtSwi,jnϕFi,jn1Bwi,j′+1Bwi,jnϕFi,j′,(15)Cowi,jn+1=Vbi,jΔtϕFBoi,jn+1,(16)Cwwi,jn+1=Vbi,jΔtϕFBwi,jn+1.Bulk volume of the fracture panel (i,j) can be expressed as
(17)Vbi,j=ΔlFΔhFwFi,j.For a slightly compressible fluid, the coefficients of1/Boi,j′, 1/Bwi,j′, and ϕFi,j′ in Equations (13) and (14) are given by
(18)1Boi,j′=coBo°,(19)1Bwi,j′=cwBw°,(20)ϕFi,j′=ϕF°cF,where Bw°, Bo°, and ϕF° are the water volume coefficient, oil volume coefficient, and fracture porosity at a reference pressure, respectively; co, cw, and cF 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
(21)qSCWoln+1=−Tol,Wn+1pFln+1−pWn,l∈ψW,(22)qSCWwln+1=−Twl,Wn+1pFln+1−pWn,l∈ψW,where pWn is the wellbore pressure and Tol,Wn+1 and Twl,Wn+1 are the transmissibility of oil and water phase between fracture panel l and the wellbore at time level n+1, which are given by
(23)Tol,Wn+1=2πβkFiwFe−γFpi−pFlnreq/rWkroln+11μoBoln+1,(24)Twl,Wn+1=2πβkFiwFe−γFpi−pFlnreq/rWkrwln+11μwBwln+1.In this paper, the flowing material balance method is used to calculate the water influx rate (qSCwIi,jn+1) in Equation (8). Using the flowing material balance method, the average pressure can be calculated by
(25)p^=pe−WpVpcw.Assuming thatNw fracture panels are directly connected to the bottom aquifer, then the cumulative water influx rate is expressed as
(26)Wp=∑σ=1Nt∑g=1NwqSCwIgΔtσ.Inserting Equation (26) into Equation (25) yields:
(27)p^=pe−∑σ=1Nt∑g=1NwqSCwIgΔtσVpcw.The water influx rate of each fracture panel connecting to the bottom aquifer is given by(28)qSCwIln+1=Jwpe−∑σ=1Nt∑g=1NwqSCwIgΔtσVpcw−pFi,jn+1,where Vp is the aquifer size, pe is the initial pressure of the bottom aquifer, Jw represents the water invasion index, and qSCwIg 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
(29)∑l∈ψi,jToli,jn+1pFln+1−∑l∈ψi,jToli,jn+1+Copi,jn+1+Toi,j,Wn+1pFi,jn+1+Cowi,jn+1Swi,jn+1=−qSCFoi,jn+1−Toi,j,Wn+1pWn−Copi,jn+1pFi,jn−Cowi,jn+1Swi,jn.Taking Equation (22) and Equation (28) into Equation (8), we can rewrite the water flow equation as
(30)∑l∈ψi,jTwl,i,jn+1pFln+1+∑l∈ψi,jTwl,i,jn+1+Cwpi,jn+1+Twi,j,Wn+1+JwpFi,jn+1+Cwwi,jn+1Swi,jn+1=−Twi,j,Wn+1pWn−Jwpen+JwWpVpcw−Cwpi,jn+1pFi,jn−Cwwi,jn+1Swi,jn.Equations (29) and (30) are written for fracture panels (i,j) at time step n+1; then, we apply these two equations to NF fracture panels and arrange these fracture flow equations into a matrix format yield:
(31)B⋅X→=C+q→.where the matrix and vector in Equation (31) are given by
(32)X→=pFi,jn+1Swi,jn+1,(33)q→=−qSCFoi,jn+10,(34)Βli,jl≠i,j=Toli,jn+10Twli,jn+10,l∈ψi,j,0000,else,(35)Bi,j,i,j=−∑l∈ψi,jTol,i,jn+1−Copi,jn+1−Toi,j,Wn+1−Cowi,jn+1−∑l∈ψi,jTwl,i,jn+1−Cwpi,jn+1−Twi,j,Wn+1−Jw−Cwwi,jn+1,(36)Ci,j=−Toi,j,Wn+1pWn−Copi,jn+1pFi,jn−Cowi,jn+1Swi,jn−Twi,j,WnpWn−Jwpen+JwWpVpcw−Cwpi,jn+1pFi,jn−Cwwi,jn+1Swi,jn.There are2NF equations to be solved simultaneously for the 2NF 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.


### 3.2. Model for Flow in the Matrix System

In our work, the three-dimensional macro fracture is discretized intoNF 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(37)∂2pm∂x2+∂2pm∂y2+∂2pm∂z2+q~δx−xw,y−yw,z−zw=ϕmμctmkm∂pm∂t.Initial condition(38)pmx,y,z,0=pi.The outer boundary is assumed to be closed and given by(39)∂pm0,y,z,t∂x=0,∂pmxe,y,z,t∂x=0,(40)∂pmx,0,z,t∂y=0,∂pmx,ye,z,t∂y=0.The top boundary is assumed to be closed(41)∂pmx,y,ze,t∂z=0.And the bottom boundary is assumed to produce at constant pressure, which is given by(42)pmx,y,0,t=pi.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:
(43)p⌢m−pmx,y,z,t=2αqSCϕmctmxeyeze∫0t1+2∑n=1∞exp−n2π2χt−τxe2cosnπxwxecosnπxxe⋅1+2∑n=1∞exp−n2π2χt−τye2cosnπywyecosnπyye⋅1+2∑n=1∞exp−2n−12π2χt−τ4ze2cos2n−1πzw2zecos2n−1πz2zedτ,where pm is the current pressure of the matrix system; p⌢m is the initial pressure of the matrix system; km is the matrix permeability; ϕm is the matrix porosity; ctm is the total compressibility of matrix system; χ is diffusivity; α is the unit conversion factor; xe, ye, and ze are dimensions of the reservoir; xw,yw,zw is the center position of the plane source; q~ is a point source influx from the matrix to the fracture; and qSC is the withdrawal rate at standard condition.As we can see in Figure3, the fracture panel is regarded as a plane source. The pressure response at an arbitrary position (x,y,z) at time t due to the contribution of a continuous plane source with a withdrawal rate of qsc 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
(44)p⌢m−pmx,y,z,t=2αqSCϕmctmxeyeze∫τ1τ2∫−ΔlFi,j/2ΔlFi,j/2∫−ΔhFi,j/2ΔhFi,j/21+2∑n=1∞exp−n2π2χt−τxe2cosnπxi,jxecosnπxxe⋅1+2∑n=1∞exp−n2π2χt−τye2cosnπyi,jyecosnπyye⋅1+2∑n=1∞exp−2n−12π2χt−τ4ze2cos2n−1πzi,j2zecos2n−1πz2zedτdξdη,where ΔlFi,j is the length of the fracture panel (i,j); ΔhFi,j is the height of the fracture panel (i,j); and xi,j, yi,j, zi,j are the center position of the fracture panel (i,j).Figure 3
Basic parameters of 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 (xo,yo,zo) caused by all fracture panels can be calculated by using the superposition principle, and the equation is given by
(45)Δpmxo,yo,zo,tn+1=p⌢m−pmxo,yo,zo,tn+1=αϕmctm∑k=1n+1∑g=1NFq~SCmgk∫tk−1tkSxo,yo,zo,tn+1−τ,xg,yg,zgdτ,where
(46)Sxo,yo,zo,tn+1−τ,xg,yg,zg=∫−ΔlFi,j/2ΔlFi,j/2∫−ΔhFi,j/2ΔhFi,j/21+2∑n=1∞exp−n2π2χt−τxe2cosnπxi,jxecosnπxxe⋅1+2∑n=1∞exp−n2π2χt−τye2cosnπyi,jyecosnπyye⋅1+2∑n=1∞exp−2n−12π2χt−τ4ze2cos2n−1πzi,j2zecos2n−1πz2zedτdξdη,where q~SCmgk is the oil flow rate from the matrix to the fracture panel g per unit fracture area at time step k and Sxo,yo,zo,tn+1−τ,xg,yg,zg is the plane source function of fracture panel g.To simplify the form of the equations, letUo,gn+1,k denote
(47)Uo,gn+1,k=αϕmctm∫tk−1tkSxo,yo,zo,tn+1−τ,xg,yg,zgdτ.Taking Equations (46) and (47) into Equation (45) yields
(48)Δpmxo,yo,zo,tn+1=p⌢m−pmxo,yo,zo,tn+1=∑k=1n+1∑g=1NFq~SCmgk⋅Uo,gn+1,k.Then, the pressure response of fracture segment g caused by all fracture panels at time leveln+1 can be calculated by
(49)pmi,jn+1=p⌢m−∑g=1NFq~SCmgk⋅Ui,j,gn+1,n+1−∑k=1n∑g=1NFq~SCmgk⋅Ui,j,gn+1,k.We arrange Equation (49) into a matrix format yields
(50)U⋅q~→=−p→+p⌢→−u→,where the matrix and vector are given as
(51)U=U1,1,1,1n+1,n+1⋯U1,1,gn+1,n+1⋯U1,1,NF,NFn+1,n+1⋮⋮⋮Ui,j,1,1n+1,n+1⋯Ui,j,gn+1,n+1⋯Ui,j,NF,NFn+1,n+1⋮⋮⋮UNF,NF,1,1n+1,n+1⋯UNF,NF,gn+1,n+1⋯UNF,NF,NF,NFn+1,n+1,(52)u→=∑k=1n∑g=1NFq~SCmgk⋅H1,1,gn+1,k⋮∑k=1n∑g=1NFq~SCmgk⋅Hi,j,gn+1,k⋮∑k=1n∑g=1NFq~SCmgk⋅HNF,NF,gn+1,k,(53)p⌢→=p⌢mp⌢m⋯p⌢mT,(54)p→=pm1,1n+1⋯pmi,jn+1⋯pmNF,NFn+1T,(55)q~→=q~SCm1,1n+1⋯q~SCmi,jn+1⋯q~SCmNF,NFn+1T.It should be noted that there areNF equations to be solved simultaneously for the NF unknowns in Equation (50).


### 3.3. Coupled Model and Solution

The continuity conditions of pressure and flux along the fracture plane are given by(56)pFi,jn+1=pmi,jn+1,(57)qSCFoi,jn+1=ΔlFi,jΔhFi,jwFi,jq~SCmi,jn+1.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.(58)B⋅X→=C+q→,U⋅q~→=−p→+p⌢→−u→.Equation (58) has 3NF equations and 3NF unknowns, containing the pressure (pFi,jn+1) and saturation (Swi,jn+1) in fracture panels, and the flow rate from the matrix to the fracture panels (qSCFoi,jn+1). As described above, the terms in matrix B and vector C 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 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.


## 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,(1)∂∂ξβkFie−γFpi−pFkroμoBo∂pF∂ξ+∂∂ηβkFie−γFpi−pFkroμoBo∂pF∂η+q~~SCFo+q~~SCWo=∂∂tϕF1−SwBo.For the water phase,(2)∂∂ξβkFie−γFpi−pFkrwμwBw∂pF∂ξ+∂∂ηβkFie−γFpi−pFkrwμwBw∂pF∂η+q~~SCWw+q~~SCwI=∂∂tϕFSwBw.The relationship between oil and water saturation is given by(3)Sw+So=1.The initial condition is(4)pFξ,η,t=0=pi.Under constant bottom-hole flowing pressure condition, a well production equation is given by [30]
(5)q~~SCWo=2πβkFikrowFe−γFpi−pFμoBolnreq/rWpF−pW,(6)q~~SCWw=2πβkFikrwwFe−γFpi−pFμwBwlnreq/rWpF−pW,where kro and krw are the relative permeability of oil and water phase, respectively; kFi 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; Bo and Bw are oil and water volume coefficient, respectively; So and Sw are oil and water saturation, respectively; pF is fracture pressure; pi 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 rW is wellbore radius and req=0.14ΔlF2+ΔhF2.We first divide the fracture into several fracture panels, as illustrated in Figure2 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
(7)∑l∈ψi,jTol,i,jn+1pFln+1−pFi,jn+1+qSCFoi,jn+1+qSCWoi,jn+1=Copi,jn+1pFi,jn+1−pFi,jn+Cowi,jn+1Swi,jn+1−Swi,jn,(8)∑l∈ψi,jTwl,i,jn+1pFln+1−pFi,jn+1+qSCWwi,jn+1+qSCwIi,jn+1=Cwpi,jn+1pFi,jn+1−pFi,jn+Cwwi,jn+1Swi,jn+1−Swi,jn,where i=1,⋯Nx and j=1,⋯Ny; the transmissibility of oil and water phase between fracture panel l and (i,j) at time level n+1 are defined as
(9)Tol,i,jn+1=Gl,i,jkrol,i,jn+11μoBol,i,jn+1,(10)Twl,i,jn+1=Gl,i,jkrwl,i,jn+11μwBwl,i,jn+1.Figure 2
Discretization of the three-dimensional fracture along both the horizontal axis and the inclined axis.Geometric transmissibility (Gl,i,j) is calculated using the harmonic average of the properties of the panels, which is expressed as
(11)Gl,i,j=γlγi,jγl+γi,j,where the geometric transmissibility between the fracture panel (i,j) and the interface is given by
(12)γi,j=βkFwFΔhFΔlFi,j.The coefficients ofCopi,jn+1, Cwpi,jn+1, Cowi,jn+1, and Cwwi,jn+1 in Equations (7) and (8) depend on the pressure and saturation change of panel (i,j) for the oil and water phases at the time level n+1, which are defined as
(13)Copi,jn+1=Vbi,jΔt1−Swi,jnϕFi,jn1Boi,j′+1Boi,jnϕFi,j′,(14)Cwpi,jn+1=Vbi,jΔtSwi,jnϕFi,jn1Bwi,j′+1Bwi,jnϕFi,j′,(15)Cowi,jn+1=Vbi,jΔtϕFBoi,jn+1,(16)Cwwi,jn+1=Vbi,jΔtϕFBwi,jn+1.Bulk volume of the fracture panel (i,j) can be expressed as
(17)Vbi,j=ΔlFΔhFwFi,j.For a slightly compressible fluid, the coefficients of1/Boi,j′, 1/Bwi,j′, and ϕFi,j′ in Equations (13) and (14) are given by
(18)1Boi,j′=coBo°,(19)1Bwi,j′=cwBw°,(20)ϕFi,j′=ϕF°cF,where Bw°, Bo°, and ϕF° are the water volume coefficient, oil volume coefficient, and fracture porosity at a reference pressure, respectively; co, cw, and cF 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
(21)qSCWoln+1=−Tol,Wn+1pFln+1−pWn,l∈ψW,(22)qSCWwln+1=−Twl,Wn+1pFln+1−pWn,l∈ψW,where pWn is the wellbore pressure and Tol,Wn+1 and Twl,Wn+1 are the transmissibility of oil and water phase between fracture panel l and the wellbore at time level n+1, which are given by
(23)Tol,Wn+1=2πβkFiwFe−γFpi−pFlnreq/rWkroln+11μoBoln+1,(24)Twl,Wn+1=2πβkFiwFe−γFpi−pFlnreq/rWkrwln+11μwBwln+1.In this paper, the flowing material balance method is used to calculate the water influx rate (qSCwIi,jn+1) in Equation (8). Using the flowing material balance method, the average pressure can be calculated by
(25)p^=pe−WpVpcw.Assuming thatNw fracture panels are directly connected to the bottom aquifer, then the cumulative water influx rate is expressed as
(26)Wp=∑σ=1Nt∑g=1NwqSCwIgΔtσ.Inserting Equation (26) into Equation (25) yields:
(27)p^=pe−∑σ=1Nt∑g=1NwqSCwIgΔtσVpcw.The water influx rate of each fracture panel connecting to the bottom aquifer is given by(28)qSCwIln+1=Jwpe−∑σ=1Nt∑g=1NwqSCwIgΔtσVpcw−pFi,jn+1,where Vp is the aquifer size, pe is the initial pressure of the bottom aquifer, Jw represents the water invasion index, and qSCwIg 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
(29)∑l∈ψi,jToli,jn+1pFln+1−∑l∈ψi,jToli,jn+1+Copi,jn+1+Toi,j,Wn+1pFi,jn+1+Cowi,jn+1Swi,jn+1=−qSCFoi,jn+1−Toi,j,Wn+1pWn−Copi,jn+1pFi,jn−Cowi,jn+1Swi,jn.Taking Equation (22) and Equation (28) into Equation (8), we can rewrite the water flow equation as
(30)∑l∈ψi,jTwl,i,jn+1pFln+1+∑l∈ψi,jTwl,i,jn+1+Cwpi,jn+1+Twi,j,Wn+1+JwpFi,jn+1+Cwwi,jn+1Swi,jn+1=−Twi,j,Wn+1pWn−Jwpen+JwWpVpcw−Cwpi,jn+1pFi,jn−Cwwi,jn+1Swi,jn.Equations (29) and (30) are written for fracture panels (i,j) at time step n+1; then, we apply these two equations to NF fracture panels and arrange these fracture flow equations into a matrix format yield:
(31)B⋅X→=C+q→.where the matrix and vector in Equation (31) are given by
(32)X→=pFi,jn+1Swi,jn+1,(33)q→=−qSCFoi,jn+10,(34)Βli,jl≠i,j=Toli,jn+10Twli,jn+10,l∈ψi,j,0000,else,(35)Bi,j,i,j=−∑l∈ψi,jTol,i,jn+1−Copi,jn+1−Toi,j,Wn+1−Cowi,jn+1−∑l∈ψi,jTwl,i,jn+1−Cwpi,jn+1−Twi,j,Wn+1−Jw−Cwwi,jn+1,(36)Ci,j=−Toi,j,Wn+1pWn−Copi,jn+1pFi,jn−Cowi,jn+1Swi,jn−Twi,j,WnpWn−Jwpen+JwWpVpcw−Cwpi,jn+1pFi,jn−Cwwi,jn+1Swi,jn.There are2NF equations to be solved simultaneously for the 2NF 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.


## 3.2. Model for Flow in the Matrix System

In our work, the three-dimensional macro fracture is discretized intoNF 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(37)∂2pm∂x2+∂2pm∂y2+∂2pm∂z2+q~δx−xw,y−yw,z−zw=ϕmμctmkm∂pm∂t.Initial condition(38)pmx,y,z,0=pi.The outer boundary is assumed to be closed and given by(39)∂pm0,y,z,t∂x=0,∂pmxe,y,z,t∂x=0,(40)∂pmx,0,z,t∂y=0,∂pmx,ye,z,t∂y=0.The top boundary is assumed to be closed(41)∂pmx,y,ze,t∂z=0.And the bottom boundary is assumed to produce at constant pressure, which is given by(42)pmx,y,0,t=pi.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:
(43)p⌢m−pmx,y,z,t=2αqSCϕmctmxeyeze∫0t1+2∑n=1∞exp−n2π2χt−τxe2cosnπxwxecosnπxxe⋅1+2∑n=1∞exp−n2π2χt−τye2cosnπywyecosnπyye⋅1+2∑n=1∞exp−2n−12π2χt−τ4ze2cos2n−1πzw2zecos2n−1πz2zedτ,where pm is the current pressure of the matrix system; p⌢m is the initial pressure of the matrix system; km is the matrix permeability; ϕm is the matrix porosity; ctm is the total compressibility of matrix system; χ is diffusivity; α is the unit conversion factor; xe, ye, and ze are dimensions of the reservoir; xw,yw,zw is the center position of the plane source; q~ is a point source influx from the matrix to the fracture; and qSC is the withdrawal rate at standard condition.As we can see in Figure3, the fracture panel is regarded as a plane source. The pressure response at an arbitrary position (x,y,z) at time t due to the contribution of a continuous plane source with a withdrawal rate of qsc 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
(44)p⌢m−pmx,y,z,t=2αqSCϕmctmxeyeze∫τ1τ2∫−ΔlFi,j/2ΔlFi,j/2∫−ΔhFi,j/2ΔhFi,j/21+2∑n=1∞exp−n2π2χt−τxe2cosnπxi,jxecosnπxxe⋅1+2∑n=1∞exp−n2π2χt−τye2cosnπyi,jyecosnπyye⋅1+2∑n=1∞exp−2n−12π2χt−τ4ze2cos2n−1πzi,j2zecos2n−1πz2zedτdξdη,where ΔlFi,j is the length of the fracture panel (i,j); ΔhFi,j is the height of the fracture panel (i,j); and xi,j, yi,j, zi,j are the center position of the fracture panel (i,j).Figure 3
Basic parameters of 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 (xo,yo,zo) caused by all fracture panels can be calculated by using the superposition principle, and the equation is given by
(45)Δpmxo,yo,zo,tn+1=p⌢m−pmxo,yo,zo,tn+1=αϕmctm∑k=1n+1∑g=1NFq~SCmgk∫tk−1tkSxo,yo,zo,tn+1−τ,xg,yg,zgdτ,where
(46)Sxo,yo,zo,tn+1−τ,xg,yg,zg=∫−ΔlFi,j/2ΔlFi,j/2∫−ΔhFi,j/2ΔhFi,j/21+2∑n=1∞exp−n2π2χt−τxe2cosnπxi,jxecosnπxxe⋅1+2∑n=1∞exp−n2π2χt−τye2cosnπyi,jyecosnπyye⋅1+2∑n=1∞exp−2n−12π2χt−τ4ze2cos2n−1πzi,j2zecos2n−1πz2zedτdξdη,where q~SCmgk is the oil flow rate from the matrix to the fracture panel g per unit fracture area at time step k and Sxo,yo,zo,tn+1−τ,xg,yg,zg is the plane source function of fracture panel g.To simplify the form of the equations, letUo,gn+1,k denote
(47)Uo,gn+1,k=αϕmctm∫tk−1tkSxo,yo,zo,tn+1−τ,xg,yg,zgdτ.Taking Equations (46) and (47) into Equation (45) yields
(48)Δpmxo,yo,zo,tn+1=p⌢m−pmxo,yo,zo,tn+1=∑k=1n+1∑g=1NFq~SCmgk⋅Uo,gn+1,k.Then, the pressure response of fracture segment g caused by all fracture panels at time leveln+1 can be calculated by
(49)pmi,jn+1=p⌢m−∑g=1NFq~SCmgk⋅Ui,j,gn+1,n+1−∑k=1n∑g=1NFq~SCmgk⋅Ui,j,gn+1,k.We arrange Equation (49) into a matrix format yields
(50)U⋅q~→=−p→+p⌢→−u→,where the matrix and vector are given as
(51)U=U1,1,1,1n+1,n+1⋯U1,1,gn+1,n+1⋯U1,1,NF,NFn+1,n+1⋮⋮⋮Ui,j,1,1n+1,n+1⋯Ui,j,gn+1,n+1⋯Ui,j,NF,NFn+1,n+1⋮⋮⋮UNF,NF,1,1n+1,n+1⋯UNF,NF,gn+1,n+1⋯UNF,NF,NF,NFn+1,n+1,(52)u→=∑k=1n∑g=1NFq~SCmgk⋅H1,1,gn+1,k⋮∑k=1n∑g=1NFq~SCmgk⋅Hi,j,gn+1,k⋮∑k=1n∑g=1NFq~SCmgk⋅HNF,NF,gn+1,k,(53)p⌢→=p⌢mp⌢m⋯p⌢mT,(54)p→=pm1,1n+1⋯pmi,jn+1⋯pmNF,NFn+1T,(55)q~→=q~SCm1,1n+1⋯q~SCmi,jn+1⋯q~SCmNF,NFn+1T.It should be noted that there areNF equations to be solved simultaneously for the NF unknowns in Equation (50).


## 3.3. Coupled Model and Solution

The continuity conditions of pressure and flux along the fracture plane are given by(56)pFi,jn+1=pmi,jn+1,(57)qSCFoi,jn+1=ΔlFi,jΔhFi,jwFi,jq~SCmi,jn+1.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.(58)B⋅X→=C+q→,U⋅q~→=−p→+p⌢→−u→.Equation (58) has 3NF equations and 3NF unknowns, containing the pressure (pFi,jn+1) and saturation (Swi,jn+1) in fracture panels, and the flow rate from the matrix to the fracture panels (qSCFoi,jn+1). As described above, the terms in matrix B and vector C 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 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.


## 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 Figure4, 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×100m, and the grid dimensions used to discrete the reservoir are 20×20×2m. 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×107m3, 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 4
Schematic diagram of fracture distribution and discretizing for numerical simulation.Table 1
Parameters used for model validation.
ParametersSymbolUnitValueFracture permeabilitykFmD500Fracture porosityϕFDimensionless0.4Fracture half-lengthΔlFm70Fracture half-heightΔhFm50Fracture widthwFm0.01Inclination angle of fractureθ°75Fracture compressibilitycFMPa-11×10−3Permeability modulesγFMPa-10.05Matrix permeabilitykmmD0.1Matrix porosityϕmDimensionless0.1Formation thicknesshm100Aquifer sizeVpm31 × 107Oil volume coefficientBom3/m31.6Oil viscosityμomPa·s0.3Water volume coefficientBwm3/m31.0Water viscosityμwmPa·s0.8Figure 5
Oil-water relative permeability of the hydraulic fracture system.Figure6 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.Figure 6
Comparison of the results of the numerical simulator and the semianalytical method.
(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. Table2 shows the values of the parameters used for sensitivity analysis.Table 2
Parameters used for sensitivity analysis.
ParametersSymbolUnitValueFracture conductivityCFmD·m5, 10, 50Permeability modulesγFMPa-10.05, 0.1, 0.2Inclination angle of fractureθ°45, 60, 75Aquifer sizeVpm37×105, 7×107, 7×109Perforation locationΔhm20, 60Wellbore pressure dropΔpMPa5, 10, 15


#### 4.2.1. The Effect of the Fracture Conductivity

Figure7 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.Figure 7
Effect of fracture conductivity on oil and water production rates.
(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 Figure8 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.Figure 8
Effect of fracture permeability modulus on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


#### 4.2.3. The Effect of the Inclination Angle of Fracture

Figure9 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.Figure 9
Effect of the inclination angle on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


#### 4.2.4. The Effect of the Aquifer Size

Figure10 depicts the oil and water production dynamic curves with different aquifer sizes. The aquifer size is given 7×105 m3, 7×107 m3, and 7×109 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.Figure 10
Effect of the aquifer size on oil and water production rates.
(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Δ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.Figure 11
Effect of the perforation position on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


#### 4.2.6. The Effect of the Wellbore Pressure Drop

Figure12 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.Figure 12
Effect of the wellbore pressure drop on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


## 4.1. Model Validation

In this section, the proposed semianalytical method is benchmarked against the commercial numerical simulator Eclipse. As shown in Figure4, 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×100m, and the grid dimensions used to discrete the reservoir are 20×20×2m. 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×107m3, 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 4
Schematic diagram of fracture distribution and discretizing for numerical simulation.Table 1
Parameters used for model validation.
ParametersSymbolUnitValueFracture permeabilitykFmD500Fracture porosityϕFDimensionless0.4Fracture half-lengthΔlFm70Fracture half-heightΔhFm50Fracture widthwFm0.01Inclination angle of fractureθ°75Fracture compressibilitycFMPa-11×10−3Permeability modulesγFMPa-10.05Matrix permeabilitykmmD0.1Matrix porosityϕmDimensionless0.1Formation thicknesshm100Aquifer sizeVpm31 × 107Oil volume coefficientBom3/m31.6Oil viscosityμomPa·s0.3Water volume coefficientBwm3/m31.0Water viscosityμwmPa·s0.8Figure 5
Oil-water relative permeability of the hydraulic fracture system.Figure6 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.Figure 6
Comparison of the results of the numerical simulator and the semianalytical method.
(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. Table2 shows the values of the parameters used for sensitivity analysis.Table 2
Parameters used for sensitivity analysis.
ParametersSymbolUnitValueFracture conductivityCFmD·m5, 10, 50Permeability modulesγFMPa-10.05, 0.1, 0.2Inclination angle of fractureθ°45, 60, 75Aquifer sizeVpm37×105, 7×107, 7×109Perforation locationΔhm20, 60Wellbore pressure dropΔpMPa5, 10, 15


### 4.2.1. The Effect of the Fracture Conductivity

Figure7 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.Figure 7
Effect of fracture conductivity on oil and water production rates.
(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 Figure8 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.Figure 8
Effect of fracture permeability modulus on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


### 4.2.3. The Effect of the Inclination Angle of Fracture

Figure9 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.Figure 9
Effect of the inclination angle on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


### 4.2.4. The Effect of the Aquifer Size

Figure10 depicts the oil and water production dynamic curves with different aquifer sizes. The aquifer size is given 7×105 m3, 7×107 m3, and 7×109 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.Figure 10
Effect of the aquifer size on oil and water production rates.
(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Δ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.Figure 11
Effect of the perforation position on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


### 4.2.6. The Effect of the Wellbore Pressure Drop

Figure12 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.Figure 12
Effect of the wellbore pressure drop on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


## 4.2.1. The Effect of the Fracture Conductivity

Figure7 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.Figure 7
Effect of fracture conductivity on oil and water production rates.
(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 Figure8 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.Figure 8
Effect of fracture permeability modulus on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


## 4.2.3. The Effect of the Inclination Angle of Fracture

Figure9 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.Figure 9
Effect of the inclination angle on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


## 4.2.4. The Effect of the Aquifer Size

Figure10 depicts the oil and water production dynamic curves with different aquifer sizes. The aquifer size is given 7×105 m3, 7×107 m3, and 7×109 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.Figure 10
Effect of the aquifer size on oil and water production rates.
(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Δ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.Figure 11
Effect of the perforation position on oil and water production rates.
(a)
Oil production rate(b)
Water production rate


## 4.2.6. The Effect of the Wellbore Pressure Drop

Figure12 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.Figure 12
Effect of the wellbore pressure drop on oil and water production rates.
(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



---

*Source: 5549411-2021-04-27.xml*