# Pressure Transient Behavior of Horizontal Well with Time-Dependent Fracture Conductivity in Tight Oil Reservoirs

Authors: Qihong Feng; Tian Xia; Sen Wang; Harpreet Singh
Journal: Geofluids
 (2017)

---

## Abstract

This work presents a discussion on the pressure transient response of multistage fractured horizontal well in tight oil reservoirs. Based on Green’s function, a semianalytical model is put forward to obtain the behavior. Our proposed model accounts for fluid flow in four contiguous regions of the tight formation by using pressure continuity and mass conservation. The time-dependent conductivity of hydraulic fractures, which is ignored in previous models but highlighted by recent experiments, is also taken into account in our proposed model. We also include the effect of pressure drop along a horizontal wellbore. We substantiate the validity of our model and analyze the different flow regimes, as well as the effects of initial conductivity, fracture distribution, and geometry on the pressure transient behavior. Our results suggest that the decrease of fracture conductivity has a tremendous effect on the well performance. Finally, we compare our model results with the field data from a multistage fractured horizontal well in Jimsar sag, Xinjiang oilfield, and a good agreement is obtained.


---

## Body

## 1. Introduction

Unconventional resources are playing increasingly important roles in the energy industry [1, 2]. We have witnessed a great success in North America [3, 4]. Owing to its extremely low permeability [5], tight oil could not be economically developed via conventional technology [6–8]. Multistage fractured horizontal well (MFHW) is an efficient technique in the development of unconventional reserves [9]. However, analyzing the pressure response of multifractured horizontal wells is challenging because several factors, for example, fracture conductivity, fracture geometry, skew angle between fracture, and horizontal well, are responsible for the pressure transient behavior. Therefore, the variation of pressure as a function of these factors is still ambiguous.Pressure transient analysis is an important tool to estimate the formation characteristics. Significant efforts have been dedicated in proposing models to describe the process. Gringarten and Ramey Jr. analyzed the transient behavior of uniform-flux fracture and infinite-conductivity fracture well through source function method [10]. Cinco-Ley and Samaniego developed a mathematical model for finite-conductivity fracture [11]. Then Cinco-Ley and Meng gave the solution for wells with finite-conductivity fractures in Laplace domain [12]. Ozkan and Raghavan employed point-source solution method to get transient pressure solutions under a variety of conditions [13]. Although these studies only dealt with a single fracture, they have laid a solid foundation for the analysis of pressure behavior in multiple fractures.The tools that are commonly used to make the transient pressure analysis of multiple fractures can be divided into three categories: analytical, semianalytical, and numerical method. Guo et al. developed an analytical method for horizontal well intersecting multiple fractures [14]. However, the interference of the fractures was neglected. Wan and Aziz derived an analytical 3D solution for horizontal well with multiple random fractures by using Fourier analysis to a 2D solution [15]. Ozkan et al. proposed a trilinear flow model, in which the linear flow in outer reservoir, inner reservoir, and hydraulic fractures are included [16]. Brown et al. improved this model to simulate the pressure transient and production behaviors of fractured horizontal wells in unconventional shale reservoirs [17]. Al Rbeawi and Djebbar introduced a new analytical model that can be used to investigate the pressure behavior and flow regimes of a horizontal well with multiple inclined hydraulic fractures and applied it in type curve matching [18].Semianalytical approaches are another important way in analyzing the transient behavior. Horne and Temeng considered the interference among the fractures via the superposition of influence functions [19]. Zerzar and Bettam combined the boundary element method and Laplace transformation to deal with interaction of reservoir flow and fracture flow [20]. Yao et al. presented a method based on Green’s functions and the source/sink method to obtain the transient pressure response for a multifractured horizontal well in a closed box-shaped reservoir [21]. Zhou et al. proposed a semianalytical model to simulate the pressure transient behavior in complex hydraulic fracture networks [22]. Yu combined gas desorption and Zhou et al.’s model to solve the gas production problem in shale gas reservoir [23]. Jia et al. presented a model to solve the transient behavior in complex fracture networks with deep consideration for flow in fractures [24].Numerical approaches overcome many limitations in analytical and semianalytical method in studying unconventional reservoirs. Al-Kobaisi et al. presents a hybrid, numerical-analytical model for the pressure transient response of a finite-conductivity fracture intercepted by a horizontal well [25]. Freeman et al. used numerical sensitivity studies to show the effect of mechanisms and factors on the performance of multifractured horizontal well [26–28]. Olorode et al. employed numerical method to study the effect of fracture angularity and nonplanar fracture configurations on well performance [29]. Yu et al. conducted numerical simulation to investigate the impact of fracture patterns, matrix permeability, cluster spacing, and fracture conductivity [30].In previous studies, the conductivity of hydraulic fractures was often assumed to be uniformly distributed and remained constant with time. However, this assumption contradicted field practice. After hydraulic fracturing, the fractures close rapidly. In order to mitigate the production decrease, proppants are added in the fracturing fluid to prop the fracture and maintain the productivity. However, this operation can only alleviate the decrease rate of the fracture conductivity. The conductivity will eventually decrease [31–34]. Previous studies suggest that the fracture conductivity decreases rapidly during the first couple of days and for the rest of the time the decline degree remains very small (Figure 1).Figure 1
Schematic showing the variation of relative fracture conductivity with time (Montgomery, 1984).Table1 shows some experiments that studied the variation of fracture conductivity with time. McDaniel reported that the sand lost 80% of the conductivity within 15 days, among which the resin coated sand lost 55%, and ceramic proppants lost 25% to 30% [35]. Tests by Cobb and Farrell showed that ceramic proppants lost ~20% of conductivity in 70 days when confined at 10000 psi and the sand lost over 30% when confined at 5000 psi [36]. Handren and Palisch reported decline in the conductivity with sands losing 55% and resin coated sand losing 25% to 30% [37]. Because the production rate is strongly dependent on the fracture conductivity, its variation with time must be taken into account to accurately predict the well performance.Table 1
Background information for the experimental studies on the variation of fracture conductivity with time.
Author
Material
Confining stress
McDaniel (1986)
Sand, resin coated sand, andthree ceramic proppants
8000 psi
Cobb and Farrell (1986)
Ceramic proppants
10000 psi and 5000 psi
Handren and Palisch (2007)
Sand and resin coated sand
6000 psiContrary to the reported observations, many studies have assumed horizontal well to be of infinite conductivity [16, 24, 38]. This assumption is not reliable as it cannot reflect the radial influx, frictional, and acceleration effects. It is necessary to examine the effect of pressure drop within the wellbore on the production performance.The objective of this study is to examine the effect of time-dependent fracture conductivity on the transient behavior of MFHW. As shown in Figure2(a), some hydraulic fractures are not perpendicular to the horizontal wellbore, which contradicts the common assumption in the analytical models. That is, the pressure transient behavior of this complex fracture network is not readily to be analytically accounted for. Therefore, we present a semianalytical model to take into account the time-dependent fracture conductivity and the pressure drop along a horizontal wellbore, as well as the complex fracture networks.Figure 2
(a) Microseismic data showing the complex fractures in Jimsar sag, Xinjiang oilfield. (b) Schematic of a fracture network. The black, red, and blue lines represent the wellbore, hydraulic fractures, and natural fractures, respectively.
(a)
(b)


## 2. Methodology

### 2.1. Mathematical Model

The reservoir after fracturing includes four regions: the matrix, the natural fracture network, hydraulic fractures, and the horizontal wellbore. Therefore, the fluid flow in the reservoir consists of(1) fluid flow from the reservoir to the hydraulic fractures ((i) in Figure 2(b)); (2) fluid flow from the reservoir to the natural fractures ((ii) in Figure 2(b)); (3) fluid flow from natural fracture to hydraulic fractures ((iii) in Figure 2(b)); (4) fluid flow among natural fractures ((iv) in Figure 2(b)); (5) fluid form hydraulic fractures to the horizontal wellbore ((v) in Figure 2(b)); (6) fluid flow in the horizontal wellbore ((vi) in Figure 2(b)). These six types of flow can be categorized into three groups: reservoir flow ((i) and (ii) in Figure 2(b)), fracture flow ((iii), (iv), and (v) in Figure 2(b)), and wellbore flow ((vi) in Figure 2(b)).Following assumptions are made:(1)
Reservoir is isotropic, homogeneous, box-shaped, and of uniform thickness with impermeable boundaries.(2)
Fluid in the reservoir is single-phase and slightly compressible and its compressibility and viscosity are constant.(3)
Fractures are rectangular and vertical. The flux rate is uniform along the fracture.(4)
Horizontal well is parallel to the upper and lower boundary of the reservoir.(5)
Effect of gravity is neglected.First, we define some dimensionless parameters for generality. The dimensionless pressure and time can be defined as(1)pD=2πρkhQμpi-p,tD=ktϕμCtL2.The dimensionless flow rate and influx rate can be defined as(2)qD=qQ,qfD=qfQL.The dimensionless length along x and y directions can be defined as(3)xD=xL,yD=yL.


#### 2.1.1. Reservoir Flow

We use Green’s function that has been frequently applied to solve problems of transient flow since its usage was first explored in well testing by Gringarten and Ramey Jr. [10]. Based on the Newman product method, instantaneous source functions of fracture panels can be obtained; therefore, the pressure response at any point in the reservoir from one fracture panel can be expressed as follows (see Appendix for further details):(4)Δpx,y,z,t=pi-px,y,z,t=1ϕCt∫0tqfjt-τSjx,y,z,τdτ,where Sj(x,y,z,τ) is the instantaneous plane source function of thejth panel and pi is the initial reservoir pressure which is assumed to be uniformly distributed in the reservoir.For fracturej rotated at any horizontal angle to the well, the plane source function is(5)Sjx,y,z=1xeyeze∫xj1xj2∫0ze1+2∑n=1∞exp⁡-n2π2ηxtxe2cos⁡nπx′xecos⁡nπxxe·1+2∑n=1∞exp⁡-n2π2ηytye2cos⁡nπx′tan⁡θjyecos⁡nπyye·1+2∑n=1∞exp⁡-n2π2ηztze2cos⁡nπz′zecos⁡nπzzedx′dz′.By using the superposition principle, pressure at any point in the reservoir at a given time can be given by(6)Δpx,y,z,t=pi-px,y,z,t=1ϕC∑j=1Np∫0tqfjt-τSjx,y,z,τdτ.Then we can get the dimensionless form of (6), (7)pDxD,yD,zD,tD=2π∑j=1Np∫0tDqfjDtD-τDSjDxD,yD,zD,τDdτD.Therefore, pressure map can be obtained at a given time.


#### 2.1.2. Fracture Flow

We assume that each hydraulic fracture is of finite-conductivity and fluid flow inside the fracture is one dimensional [17, 39–41]. We use Darcy’s equation to describe the fluid flow in fractures. For thejth panel, the pressure at any point in the fracture is (see Appendix for further details):(8)pj1-pjm=∫xj1xjmμρkfbfhjqj1+qfjx-xj1dx.Then the dimensionless form is(9)pjmD-pj1D=2πkLkfbf∫xj1DxjmDqj1D+qfjDxD-xj1DdxD.Equations (8) and (9) can be applied to both hydraulic fractures and natural fractures.It is worth noting that fluid flow from the fractures to the horizontal wellbore is radial in the near-well region (Figure3). In order to minimize the error caused by linear flow approximation, a choke flow factor is introduced [42]:(10)sc=h2πxeln⁡h2rw-π2,where rw is the radius of the horizontal well and bf is the hydraulic fracture width.Figure 3
Sketch showing the change of flow pattern from linear to radial.


#### 2.1.3. Time-Dependent Conductivity

Proppants are often pumped into the formation to maintain the fracture conductivity. However, the proppant particles usually break and embed, which causes the fracture conductivity to reduce until equilibrium is established. Montgomery and Steanson suggested that there is a logarithmic relation between hydraulic fracture conductivity and time for 10/20 Sand and 20/40 Sand (Figure4(a)) [43]. Other tests also show logarithmic relation between fracture conductivity and time, including Lanzhou sand from China with diameter of 0.45~0.90 mm [44], and sand from Shanshan oil field in China (Figure 4(c)) [45]. Other kinds of proppants such as ceramic proppants with diameter of 0.45~0.9 mm and Lanzhou sand with diameter of 0.9~1.25 mm show similar correlations, but the slope of the curve, also referred to as the decline coefficient, for each proppant is different. The decline coefficient decreases in the order of Lanzhou sand with diameter of 0.90~1.25 mm, the ceramic proppant, and the Lanzhou sand with smaller diameter.Figure 4
Relation between relative fracture conductivity and time. (a) Montgomery and Steanson (1985) tested two samples with different kind of proppants for 9 months, case 1 denotes 10/20 Sand at 250°F, and case 2 denotes 20/40 Sand at 75°F. (b) Experiment conducted by Yu (1987). (c) Field data from Shanshan oilfield.
(a)
(b)
(c)Based on the trend of observations in Figure4, we use the following model to describe the variation of fracture conductivity as a function of time,(11)Cf=Cf01-βlg⁡tt0,where Cf is the fracture conductivity at time t, Cf0 is the initial conductivity, and β is the decline coefficient. Cf in (11) is substituted by kf·bf in (8) and (9). It is assumed that when the effect of hydraulic fracture disappears totally, the fracture conductivity is equal to the product of formation permeability and fracture width.


#### 2.1.4. Wellbore Flow

Previous studies assumed horizontal wellbore as infinite-conductivity pipe [15, 16], which implicitly assumes no pressure drop along the horizontal well because of uniform pressure distribution along the wellbore. To examine the validity of this approximation, we include the wellbore pressure drop in the present model. The pressure decrease along the wellbore consists of frictional losses and acceleration losses (Figure 5). Based on the Darcy-Weisbach Equation [46], the fictional pressure drop can be expressed as(12)Δpfric=fi2ΔliDiρivi2,where fi is the frictional coefficient and Δli is the length of the wellbore segment. For each segment of the wellbore, fluid flow from the hydraulic fractures to the wellbore will cause an increase of the flow rate. This will result in the change of the momentum of the fluid which leads to the acceleration pressure drop:(13)Δpacce=ρivout2-vin2.Dimensionless form of the pressure drop along the well bore is given by (see Appendix for further details):(14)pD,i-pD,i+1=2πkhQμA2qpD,j2-qpD,j+12-fjΔlj2DqpD,j2.Figure 5
Sketch of microsegment in horizontal well.


#### 2.1.5. Coupling Relationship

Due to the pressure continuity at the center of each fracture, the pressure response obtained from the reservoir flow should be consistent with the fracture flow. Therefore, we have(15)pjcD1=pjcD2.pjcD1 and pjcD2 are given by (7) and (9), respectively. The mass balance is applied to both the intersection nodes between fractures and the intersection nodes of fractures and wellbore (Figure 6). Therefore, for each node, the inflow of fluid must be equal to the outflow of the fluid,(16)qin,D=qout,D.For the intersection nodes of fractures and wellbore, the inflow and outflow are taken into account between the fractures and between the fractures and the wellbore pipe.Figure 6
Sketch of mass balance at intersection nodes. (a) Intersection node between fractures and wellbore. (b) Intersection node between fractures. (Black line denotes wellbore, red lines represent hydraulic fracture, and blue line is natural fracture.)
(a)
(b)


### 2.2. Computational Approach

From (17), we obtain nv equations at the nodes for mass balance, np equations for pressure drop along the fractures, and np equations for pressure continuity in the panel centers, resulting in nv+2np nonlinear equations that need to be solved. Newton’s method has been widely used to solve systems of equations because of its quick convergence. However, it requires the inverse of Hessian matrix at each iteration and the convergence may not be reached if the Hessian matrix is ill-conditioned or nonpositive definite. Therefore, we use Gauss-Newton method, which is an improved version of Newton’s method for finding a minimum of a function [47]. The basic idea is to use the Taylor series expansion to approximate nonlinear regression model and correct the solution through iteration. This algorithm is robust and has a good convergence rate. The iteration equation of Gauss-Newton method can be expressed as(17)xk+1=xk-∇FxkT∇Fxk-1∇FxkTFxk.Pressure at each node, flow rate inside the fracture, and fluxes along the fractures can be obtained through iteration. Thus, dimensionless pressure at any point in the reservoir can be calculated via (7). If the bottom-hole pressure is given, we can get the flow rate in the same way. The flow rate behavior in different conditions can be obtained consequently.


### 2.3. Model Validation

While solving for the numerical solution, fractures are discretized into several panels and to ensure that we obtain correct numerical solution it is necessary to probe the grid independence with respect to the number of panels used. Figure7 shows the grid independence analysis for a horizontal well with three transverse hydraulic fractures, where m is the number of panels that a fracture is discretized into. We can conclude that the result we obtained from our model is independent of segment when m is greater than 1.Figure 7
Effect of discretized number per fracture on the simulation result,m is the number of discrete panels in each fracture.CMG, a commercial reservoir numerical simulator [48], was employed to validate our model. In this article, the IMEX module in CMG, which is a conventional three-phase black-oil simulator, is utilized to make the comparison. Interested readers may find the governing equations of this module from the textbooks on reservoir numerical simulation. The conventional Cartesian grids are employed and the total number of cells is 121 × 121 × 5 and the basic parameters used in the validation are listed in Table 2. The outer boundary of our simulation domain is impermeable. Meanwhile, to be consistent with the assumptions of our proposed model, we maintain the oil production rate constant in the simulation. Figure 8 compares the bottom-hole pressure obtained from CMG simulator and our model and it shows a good match during initial stage. Figure 9 shows the comparison of reservoir pressure from CMG with our model, and except for the pressures near fractures, we see a good match between the two. The reason for this slight difference around the fractures is because in CMG the fluid flows directly from the reservoir to the horizontal well, however, many studies [27, 49, 50] have reported that the correct representation of fluid flow must consider intersection and activation of preexisting natural fractures with hydraulic fractures, which will create a complex fracture network. This representation of complex fracture network can be readily incorporated in our model; however, it is an extremely challenging task to do that in CMG.Table 2
Basic parameters for model validation.
Properties
Value
Reservoir permeabilityk, μm2
1 × 10−4
Reservoir porosity,ϕ
10%
Reservoir length,xe, m
600
Formation thickness,h, m
20
Total compressibility,Ct, MPa−1
4 × 10−3
Initial pressure,pi, MPa
30
Fluid viscosity,μ, mPa·s
20
Fluid density,ρ, kg/m3
900
Fracture conductivity,Cdf, μm2·m
0.5
Fracture half-length,xf, m
75
Fracture spacing,d, m
100
Production rate,Q, t/d
4.32Figure 8
Bottom-hole pressure obtained using our proposed model and CMG.Figure 9
Pressure distribution of the reservoir. (a) CMG simulator. (b) Our model. The horizontal line represents the well.
(a)
(b)


## 2.1. Mathematical Model

The reservoir after fracturing includes four regions: the matrix, the natural fracture network, hydraulic fractures, and the horizontal wellbore. Therefore, the fluid flow in the reservoir consists of(1) fluid flow from the reservoir to the hydraulic fractures ((i) in Figure 2(b)); (2) fluid flow from the reservoir to the natural fractures ((ii) in Figure 2(b)); (3) fluid flow from natural fracture to hydraulic fractures ((iii) in Figure 2(b)); (4) fluid flow among natural fractures ((iv) in Figure 2(b)); (5) fluid form hydraulic fractures to the horizontal wellbore ((v) in Figure 2(b)); (6) fluid flow in the horizontal wellbore ((vi) in Figure 2(b)). These six types of flow can be categorized into three groups: reservoir flow ((i) and (ii) in Figure 2(b)), fracture flow ((iii), (iv), and (v) in Figure 2(b)), and wellbore flow ((vi) in Figure 2(b)).Following assumptions are made:(1)
Reservoir is isotropic, homogeneous, box-shaped, and of uniform thickness with impermeable boundaries.(2)
Fluid in the reservoir is single-phase and slightly compressible and its compressibility and viscosity are constant.(3)
Fractures are rectangular and vertical. The flux rate is uniform along the fracture.(4)
Horizontal well is parallel to the upper and lower boundary of the reservoir.(5)
Effect of gravity is neglected.First, we define some dimensionless parameters for generality. The dimensionless pressure and time can be defined as(1)pD=2πρkhQμpi-p,tD=ktϕμCtL2.The dimensionless flow rate and influx rate can be defined as(2)qD=qQ,qfD=qfQL.The dimensionless length along x and y directions can be defined as(3)xD=xL,yD=yL.


### 2.1.1. Reservoir Flow

We use Green’s function that has been frequently applied to solve problems of transient flow since its usage was first explored in well testing by Gringarten and Ramey Jr. [10]. Based on the Newman product method, instantaneous source functions of fracture panels can be obtained; therefore, the pressure response at any point in the reservoir from one fracture panel can be expressed as follows (see Appendix for further details):(4)Δpx,y,z,t=pi-px,y,z,t=1ϕCt∫0tqfjt-τSjx,y,z,τdτ,where Sj(x,y,z,τ) is the instantaneous plane source function of thejth panel and pi is the initial reservoir pressure which is assumed to be uniformly distributed in the reservoir.For fracturej rotated at any horizontal angle to the well, the plane source function is(5)Sjx,y,z=1xeyeze∫xj1xj2∫0ze1+2∑n=1∞exp⁡-n2π2ηxtxe2cos⁡nπx′xecos⁡nπxxe·1+2∑n=1∞exp⁡-n2π2ηytye2cos⁡nπx′tan⁡θjyecos⁡nπyye·1+2∑n=1∞exp⁡-n2π2ηztze2cos⁡nπz′zecos⁡nπzzedx′dz′.By using the superposition principle, pressure at any point in the reservoir at a given time can be given by(6)Δpx,y,z,t=pi-px,y,z,t=1ϕC∑j=1Np∫0tqfjt-τSjx,y,z,τdτ.Then we can get the dimensionless form of (6), (7)pDxD,yD,zD,tD=2π∑j=1Np∫0tDqfjDtD-τDSjDxD,yD,zD,τDdτD.Therefore, pressure map can be obtained at a given time.


### 2.1.2. Fracture Flow

We assume that each hydraulic fracture is of finite-conductivity and fluid flow inside the fracture is one dimensional [17, 39–41]. We use Darcy’s equation to describe the fluid flow in fractures. For thejth panel, the pressure at any point in the fracture is (see Appendix for further details):(8)pj1-pjm=∫xj1xjmμρkfbfhjqj1+qfjx-xj1dx.Then the dimensionless form is(9)pjmD-pj1D=2πkLkfbf∫xj1DxjmDqj1D+qfjDxD-xj1DdxD.Equations (8) and (9) can be applied to both hydraulic fractures and natural fractures.It is worth noting that fluid flow from the fractures to the horizontal wellbore is radial in the near-well region (Figure3). In order to minimize the error caused by linear flow approximation, a choke flow factor is introduced [42]:(10)sc=h2πxeln⁡h2rw-π2,where rw is the radius of the horizontal well and bf is the hydraulic fracture width.Figure 3
Sketch showing the change of flow pattern from linear to radial.


### 2.1.3. Time-Dependent Conductivity

Proppants are often pumped into the formation to maintain the fracture conductivity. However, the proppant particles usually break and embed, which causes the fracture conductivity to reduce until equilibrium is established. Montgomery and Steanson suggested that there is a logarithmic relation between hydraulic fracture conductivity and time for 10/20 Sand and 20/40 Sand (Figure4(a)) [43]. Other tests also show logarithmic relation between fracture conductivity and time, including Lanzhou sand from China with diameter of 0.45~0.90 mm [44], and sand from Shanshan oil field in China (Figure 4(c)) [45]. Other kinds of proppants such as ceramic proppants with diameter of 0.45~0.9 mm and Lanzhou sand with diameter of 0.9~1.25 mm show similar correlations, but the slope of the curve, also referred to as the decline coefficient, for each proppant is different. The decline coefficient decreases in the order of Lanzhou sand with diameter of 0.90~1.25 mm, the ceramic proppant, and the Lanzhou sand with smaller diameter.Figure 4
Relation between relative fracture conductivity and time. (a) Montgomery and Steanson (1985) tested two samples with different kind of proppants for 9 months, case 1 denotes 10/20 Sand at 250°F, and case 2 denotes 20/40 Sand at 75°F. (b) Experiment conducted by Yu (1987). (c) Field data from Shanshan oilfield.
(a)
(b)
(c)Based on the trend of observations in Figure4, we use the following model to describe the variation of fracture conductivity as a function of time,(11)Cf=Cf01-βlg⁡tt0,where Cf is the fracture conductivity at time t, Cf0 is the initial conductivity, and β is the decline coefficient. Cf in (11) is substituted by kf·bf in (8) and (9). It is assumed that when the effect of hydraulic fracture disappears totally, the fracture conductivity is equal to the product of formation permeability and fracture width.


### 2.1.4. Wellbore Flow

Previous studies assumed horizontal wellbore as infinite-conductivity pipe [15, 16], which implicitly assumes no pressure drop along the horizontal well because of uniform pressure distribution along the wellbore. To examine the validity of this approximation, we include the wellbore pressure drop in the present model. The pressure decrease along the wellbore consists of frictional losses and acceleration losses (Figure 5). Based on the Darcy-Weisbach Equation [46], the fictional pressure drop can be expressed as(12)Δpfric=fi2ΔliDiρivi2,where fi is the frictional coefficient and Δli is the length of the wellbore segment. For each segment of the wellbore, fluid flow from the hydraulic fractures to the wellbore will cause an increase of the flow rate. This will result in the change of the momentum of the fluid which leads to the acceleration pressure drop:(13)Δpacce=ρivout2-vin2.Dimensionless form of the pressure drop along the well bore is given by (see Appendix for further details):(14)pD,i-pD,i+1=2πkhQμA2qpD,j2-qpD,j+12-fjΔlj2DqpD,j2.Figure 5
Sketch of microsegment in horizontal well.


### 2.1.5. Coupling Relationship

Due to the pressure continuity at the center of each fracture, the pressure response obtained from the reservoir flow should be consistent with the fracture flow. Therefore, we have(15)pjcD1=pjcD2.pjcD1 and pjcD2 are given by (7) and (9), respectively. The mass balance is applied to both the intersection nodes between fractures and the intersection nodes of fractures and wellbore (Figure 6). Therefore, for each node, the inflow of fluid must be equal to the outflow of the fluid,(16)qin,D=qout,D.For the intersection nodes of fractures and wellbore, the inflow and outflow are taken into account between the fractures and between the fractures and the wellbore pipe.Figure 6
Sketch of mass balance at intersection nodes. (a) Intersection node between fractures and wellbore. (b) Intersection node between fractures. (Black line denotes wellbore, red lines represent hydraulic fracture, and blue line is natural fracture.)
(a)
(b)


## 2.1.1. Reservoir Flow

We use Green’s function that has been frequently applied to solve problems of transient flow since its usage was first explored in well testing by Gringarten and Ramey Jr. [10]. Based on the Newman product method, instantaneous source functions of fracture panels can be obtained; therefore, the pressure response at any point in the reservoir from one fracture panel can be expressed as follows (see Appendix for further details):(4)Δpx,y,z,t=pi-px,y,z,t=1ϕCt∫0tqfjt-τSjx,y,z,τdτ,where Sj(x,y,z,τ) is the instantaneous plane source function of thejth panel and pi is the initial reservoir pressure which is assumed to be uniformly distributed in the reservoir.For fracturej rotated at any horizontal angle to the well, the plane source function is(5)Sjx,y,z=1xeyeze∫xj1xj2∫0ze1+2∑n=1∞exp⁡-n2π2ηxtxe2cos⁡nπx′xecos⁡nπxxe·1+2∑n=1∞exp⁡-n2π2ηytye2cos⁡nπx′tan⁡θjyecos⁡nπyye·1+2∑n=1∞exp⁡-n2π2ηztze2cos⁡nπz′zecos⁡nπzzedx′dz′.By using the superposition principle, pressure at any point in the reservoir at a given time can be given by(6)Δpx,y,z,t=pi-px,y,z,t=1ϕC∑j=1Np∫0tqfjt-τSjx,y,z,τdτ.Then we can get the dimensionless form of (6), (7)pDxD,yD,zD,tD=2π∑j=1Np∫0tDqfjDtD-τDSjDxD,yD,zD,τDdτD.Therefore, pressure map can be obtained at a given time.


## 2.1.2. Fracture Flow

We assume that each hydraulic fracture is of finite-conductivity and fluid flow inside the fracture is one dimensional [17, 39–41]. We use Darcy’s equation to describe the fluid flow in fractures. For thejth panel, the pressure at any point in the fracture is (see Appendix for further details):(8)pj1-pjm=∫xj1xjmμρkfbfhjqj1+qfjx-xj1dx.Then the dimensionless form is(9)pjmD-pj1D=2πkLkfbf∫xj1DxjmDqj1D+qfjDxD-xj1DdxD.Equations (8) and (9) can be applied to both hydraulic fractures and natural fractures.It is worth noting that fluid flow from the fractures to the horizontal wellbore is radial in the near-well region (Figure3). In order to minimize the error caused by linear flow approximation, a choke flow factor is introduced [42]:(10)sc=h2πxeln⁡h2rw-π2,where rw is the radius of the horizontal well and bf is the hydraulic fracture width.Figure 3
Sketch showing the change of flow pattern from linear to radial.


## 2.1.3. Time-Dependent Conductivity

Proppants are often pumped into the formation to maintain the fracture conductivity. However, the proppant particles usually break and embed, which causes the fracture conductivity to reduce until equilibrium is established. Montgomery and Steanson suggested that there is a logarithmic relation between hydraulic fracture conductivity and time for 10/20 Sand and 20/40 Sand (Figure4(a)) [43]. Other tests also show logarithmic relation between fracture conductivity and time, including Lanzhou sand from China with diameter of 0.45~0.90 mm [44], and sand from Shanshan oil field in China (Figure 4(c)) [45]. Other kinds of proppants such as ceramic proppants with diameter of 0.45~0.9 mm and Lanzhou sand with diameter of 0.9~1.25 mm show similar correlations, but the slope of the curve, also referred to as the decline coefficient, for each proppant is different. The decline coefficient decreases in the order of Lanzhou sand with diameter of 0.90~1.25 mm, the ceramic proppant, and the Lanzhou sand with smaller diameter.Figure 4
Relation between relative fracture conductivity and time. (a) Montgomery and Steanson (1985) tested two samples with different kind of proppants for 9 months, case 1 denotes 10/20 Sand at 250°F, and case 2 denotes 20/40 Sand at 75°F. (b) Experiment conducted by Yu (1987). (c) Field data from Shanshan oilfield.
(a)
(b)
(c)Based on the trend of observations in Figure4, we use the following model to describe the variation of fracture conductivity as a function of time,(11)Cf=Cf01-βlg⁡tt0,where Cf is the fracture conductivity at time t, Cf0 is the initial conductivity, and β is the decline coefficient. Cf in (11) is substituted by kf·bf in (8) and (9). It is assumed that when the effect of hydraulic fracture disappears totally, the fracture conductivity is equal to the product of formation permeability and fracture width.


## 2.1.4. Wellbore Flow

Previous studies assumed horizontal wellbore as infinite-conductivity pipe [15, 16], which implicitly assumes no pressure drop along the horizontal well because of uniform pressure distribution along the wellbore. To examine the validity of this approximation, we include the wellbore pressure drop in the present model. The pressure decrease along the wellbore consists of frictional losses and acceleration losses (Figure 5). Based on the Darcy-Weisbach Equation [46], the fictional pressure drop can be expressed as(12)Δpfric=fi2ΔliDiρivi2,where fi is the frictional coefficient and Δli is the length of the wellbore segment. For each segment of the wellbore, fluid flow from the hydraulic fractures to the wellbore will cause an increase of the flow rate. This will result in the change of the momentum of the fluid which leads to the acceleration pressure drop:(13)Δpacce=ρivout2-vin2.Dimensionless form of the pressure drop along the well bore is given by (see Appendix for further details):(14)pD,i-pD,i+1=2πkhQμA2qpD,j2-qpD,j+12-fjΔlj2DqpD,j2.Figure 5
Sketch of microsegment in horizontal well.


## 2.1.5. Coupling Relationship

Due to the pressure continuity at the center of each fracture, the pressure response obtained from the reservoir flow should be consistent with the fracture flow. Therefore, we have(15)pjcD1=pjcD2.pjcD1 and pjcD2 are given by (7) and (9), respectively. The mass balance is applied to both the intersection nodes between fractures and the intersection nodes of fractures and wellbore (Figure 6). Therefore, for each node, the inflow of fluid must be equal to the outflow of the fluid,(16)qin,D=qout,D.For the intersection nodes of fractures and wellbore, the inflow and outflow are taken into account between the fractures and between the fractures and the wellbore pipe.Figure 6
Sketch of mass balance at intersection nodes. (a) Intersection node between fractures and wellbore. (b) Intersection node between fractures. (Black line denotes wellbore, red lines represent hydraulic fracture, and blue line is natural fracture.)
(a)
(b)


## 2.2. Computational Approach

From (17), we obtain nv equations at the nodes for mass balance, np equations for pressure drop along the fractures, and np equations for pressure continuity in the panel centers, resulting in nv+2np nonlinear equations that need to be solved. Newton’s method has been widely used to solve systems of equations because of its quick convergence. However, it requires the inverse of Hessian matrix at each iteration and the convergence may not be reached if the Hessian matrix is ill-conditioned or nonpositive definite. Therefore, we use Gauss-Newton method, which is an improved version of Newton’s method for finding a minimum of a function [47]. The basic idea is to use the Taylor series expansion to approximate nonlinear regression model and correct the solution through iteration. This algorithm is robust and has a good convergence rate. The iteration equation of Gauss-Newton method can be expressed as(17)xk+1=xk-∇FxkT∇Fxk-1∇FxkTFxk.Pressure at each node, flow rate inside the fracture, and fluxes along the fractures can be obtained through iteration. Thus, dimensionless pressure at any point in the reservoir can be calculated via (7). If the bottom-hole pressure is given, we can get the flow rate in the same way. The flow rate behavior in different conditions can be obtained consequently.


## 2.3. Model Validation

While solving for the numerical solution, fractures are discretized into several panels and to ensure that we obtain correct numerical solution it is necessary to probe the grid independence with respect to the number of panels used. Figure7 shows the grid independence analysis for a horizontal well with three transverse hydraulic fractures, where m is the number of panels that a fracture is discretized into. We can conclude that the result we obtained from our model is independent of segment when m is greater than 1.Figure 7
Effect of discretized number per fracture on the simulation result,m is the number of discrete panels in each fracture.CMG, a commercial reservoir numerical simulator [48], was employed to validate our model. In this article, the IMEX module in CMG, which is a conventional three-phase black-oil simulator, is utilized to make the comparison. Interested readers may find the governing equations of this module from the textbooks on reservoir numerical simulation. The conventional Cartesian grids are employed and the total number of cells is 121 × 121 × 5 and the basic parameters used in the validation are listed in Table 2. The outer boundary of our simulation domain is impermeable. Meanwhile, to be consistent with the assumptions of our proposed model, we maintain the oil production rate constant in the simulation. Figure 8 compares the bottom-hole pressure obtained from CMG simulator and our model and it shows a good match during initial stage. Figure 9 shows the comparison of reservoir pressure from CMG with our model, and except for the pressures near fractures, we see a good match between the two. The reason for this slight difference around the fractures is because in CMG the fluid flows directly from the reservoir to the horizontal well, however, many studies [27, 49, 50] have reported that the correct representation of fluid flow must consider intersection and activation of preexisting natural fractures with hydraulic fractures, which will create a complex fracture network. This representation of complex fracture network can be readily incorporated in our model; however, it is an extremely challenging task to do that in CMG.Table 2
Basic parameters for model validation.
Properties
Value
Reservoir permeabilityk, μm2
1 × 10−4
Reservoir porosity,ϕ
10%
Reservoir length,xe, m
600
Formation thickness,h, m
20
Total compressibility,Ct, MPa−1
4 × 10−3
Initial pressure,pi, MPa
30
Fluid viscosity,μ, mPa·s
20
Fluid density,ρ, kg/m3
900
Fracture conductivity,Cdf, μm2·m
0.5
Fracture half-length,xf, m
75
Fracture spacing,d, m
100
Production rate,Q, t/d
4.32Figure 8
Bottom-hole pressure obtained using our proposed model and CMG.Figure 9
Pressure distribution of the reservoir. (a) CMG simulator. (b) Our model. The horizontal line represents the well.
(a)
(b)


## 3. Results and Discussions

### 3.1. Flow Regime Analysis

The transient behavior can be shown by type curves, which is employed to obtain the characteristic of formation and the reservoir fluids and to figure out different flow regimes. Figure10 shows the dimensionless well pressure and its derivative for a horizontal well with three fractures. The transient behavior can be divided into several flow periods.Figure 10
Pressure response of multistage fractured horizontal well (fracture numbern=3, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=8, and dimensionless reservoir length xeD=240).(
1) Linear Flow (Figure11(a)). Fluid flows linearly from the reservoir to the fracture and each fracture behaves independently. Both the slope of dimensionless pressure and derivative is 1/2 in this stage.Figure 11
Flow regimes for a multistaged fractured horizontal well: (a) linear flow; (b) early radial flow, (c) biradial flow, and (d) pseudo radial flow.
(a)
(b)
(c)
(d)(
2) Early Radial Flow (Figure11(b)). An early radial flow occurs around each fracture after the linear flow. This period mainly depends on the fracture length and fracture spacing. Besides, in this period, fractures still behave independently. The characteristic of this stage is a horizontal line of 1/(2N) in pressure derivative curve (N is the fracture stage). We can figure out that the value of the horizontal line plateau is 1/6.(
3) Biradial Flow (Figure11(c)). Fractures interact with each other and flow becomes elliptical to the wellbore. The slope of pressure derivative is 0.36.(
4) Pseudo Radial Flow (Figure11(d)). Fluid flows to the fracture-well system appear to be radial and flow across the outermost elements plays important part. The pressure derivative curve demonstrates a horizontal line of 0.5.(
5) Boundary Dominated Flow. In the closed system, the flow will reach pseudo-steady state. The pressure curve and derivate tend to merge and the slope equals 1.


### 3.2. Effect of Initial Fracture Conductivity

Figure12 illustrates the effect of initial conductivity on the pressure transient response and production rate. This figure shows that increasing hydraulic fracture conductivity results in an increase in well productivity; however, the incremental benefit decreases as the fracture conductivity increases. Figure 12 also shows that variations in the dimensionless pressure and the production rate tend to disappear as the flow regime approaches pseudo radial flow. The reason why the variations in these two parameters disappear is because the flow at late times is dominated by the outer reservoir, and therefore, the effect of fracture conductivity will be extremely small. Therefore, to utilize the potential of the hydraulic fractures, we must delay the occurrence of pseudo radial flow and that can be done by optimal well placement.Figure 12
Effect of initial fracture conductivity on pressure and rate behavior (fracture numbern=2, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=0.67, and dimensionless reservoir length xeD=120). (a) Pressure and derivative curve. (b) Dimensionless production rate.
(a)
(b)Although it is well documented that the fracture conductivity decreases rapidly at first and tends to converge asymptotically [31–34, 37], previous studies always assumed that hydraulic fractures do not vary with time. Below, we analyze the difference between time-independent and time-dependent fractures. As discussed earlier, the conductivity, and hence the productivity of the horizontal well, decreases with time. Therefore, the pressure response of the well is larger than the time-independent situation, as suggested by Figure 13(a). The flow rate of a horizontal well with time-independent conductivity fractures exceeds that with time-dependent conductivity fractures in the initial stage (Figure 13(b)), and that could be attributed to the dominant role of fractures. During late times of production, the fluid flow in the outer reservoir becomes dominant and negligible difference can be observed between time-dependent and time-independent cases. Therefore, for the optimal exploitation of an unconventional reservoir it is necessary to maintain the fracture conductivity, especially in the initial stage of the production.Figure 13
Effect of time-dependent fracture conductivity on pressure and rate behavior (fracture numbern=2, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=0.67, dimensionless reservoir length xeD=120, and decline coefficient β=0.11). (a) Dimensionless pressure and derivative curve. (b) Dimensionless production rate.
(a)
(b)Figure14 shows the effect of decline coefficient on the production performance of the horizontal well, where a larger decline coefficient results in smaller production rate and if the decline coefficient is too large the fracture permeability tends to quickly approach the formation permeability, which results in sharp decline in the production rate.Figure 14
Effect of decline coefficient on the production performance of the horizontal well (fracture numbern=2, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=0.67, and dimensionless reservoir length xeD=120).


### 3.3. Effect of Fracture Distribution

The pressure difference created during hydraulic fracturing allows the proppants to stay in the fractures away from the horizontal well heel [51], and as a result, the fracture conductivities are different that may further have an impact on the pressure behavior. Two scenarios of even and uneven proppant distribution can be seen as schematics in Figures 15(a) and 15(b), respectively, which are used to study the effect of proppant distribution.Figure 15
Sketch of proppant distribution in different fracture: (a) uneven distribution; (b) even distribution.
(a)
(b)For unevenly distributed fractures, the fracture nearest to the well toe has the largest conductivity, whereas the conductivity of the intermediate fracture equals that of the evenly distributed fractures.Figure16(a) shows the dimensionless pressure drop of a horizontal well (at constant production rate multistage fractured well) with evenly distributed and unevenly distributed fractures. The pressure response with unevenly distributed proppants is slightly larger than that with evenly distributed proppants at initial period but the gap disappears after that. Result of Figure 16(b) suggests that horizontal well with evenly distributed proppants produces with higher rate initially than with unevenly distributed proppants. Therefore, these two results indicate that well with evenly distributed proppants performs better than that with unevenly distributed proppants.Figure 16
Effect of unevenly distributed fracture conductivity on the pressure and production behavior (fracture numbern=3, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=0.67, dimensionless reservoir length xeD=120. The conductivity of evenly distributed fracture is 0.375 μm2·m; for unevenly distributed fracture, the conductivities are 0.25 μm2·m, 0.375 μm2·m, and 0.5 μm2·m.). (a) Dimensionless pressure curve. (b) Dimensionless production rate curve.
(a)
(b)Figure17 shows the dimensionless flux rate of three fractures that have evenly and unevenly distributed proppants. When the fracture conductivities (or proppant distribution) are evenly distributed, the dimensionless flux rate of Fractures 1 and 3 increases with time, but the flux rate of Fracture 2 decreases sharply with time. At the initial stage, each fracture behaves independently; however, with time the fractures begin to interfere with each other and the flux rate of Fractures 1 and 3 tends to increase, whereas owing to the symmetry, the flux rate of intermediate fracture, Fracture 2, is hindered. If the fracture conductivities are unevenly distributed, the flux rate of the third fracture will first decrease and then increase with the time. However, production rate of the first fracture increases steadily with time and it crosses the conductivity of the second fracture at some point.Figure 17
Dimensionless flux rate of each fracture. (a) Fracture conductivity evenly distributed. (b) Fracture conductivity unevenly distributed.
(a)
(b)


### 3.4. Effect of Fracture Geometry

In hydraulic fracturing, the existent propped fractures result in the redistribution of local earth stresses. Moreover, microseismic measurements have proved that there is mechanical-stress interference between multiple transverse fractures. The stress-shadow effect can restrict the growth of the fracture in the middle section while promoting the growth of the fractures at the heel or the toe [49, 52]; therefore, different fracture geometry may exist in the reservoir. Three geometries shown in Figure 18, equilong type, spindle-shaped type, and dumbbell-shaped type, are studied.Figure 18
Sketch of different fracture geometries. (a) Equilong type. (b) Spindle-shaped type. (c) Dumbbell-shaped type.
(a)
(b)
(c)Figure19 shows the dimensionless flux rate of three fractures from three fracture geometries of Figure 18, respectively. For spindle-shaped geometry, the flux rate of the first and third fracture will first decrease and then slightly increase before levelling off at long time (Figure 19(b)). However, the flux of the second fracture always remains smaller than the other two fractures, although it has a larger fracture length. Figure 19(c) shows that, for dumbbell-shaped geometry, the flux rate of fractures near the heel and toe will slightly increase with time, whereas the flux of the intermediate fracture will sharply decrease with time. From these results it is evident that no matter what the fracture geometry is, the flux of the intermediate fracture always decreases with time. Figure 20 suggests that if the bottom-hole pressure is constant, the production rate of equilong fracture type is higher than the others because it has the longest effective interference fracture length. Apparently, this is the reason why equilong fracture type is usually preferred over other fracture geometries.Figure 19
Dimensionless flux rate of each fracture. (a) Equilong type,xf1=xf2=xf3=100 m. (b) Spindle-shaped type, xf1=75 m, xf2=150 m, and xf3=75 m. (c) Dumbbell-shaped type, xf1=120 m, xf2=60 m, and xf3=120 m.
(a)
(b)
(c)Figure 20
Dimensionless production rate for different fracture configurations.


### 3.5. Effect of Horizontal Wellbore Pressure Drop

Many previous studies considered the horizontal wellbore as an infinite wellbore; however, it has been shown that wellbore pressure drop exists in the production [53, 54]. We investigated the effect of horizontal wellbore pressure drop on pressure behavior as shown in Figure 21, which shows that the effect of considering wellbore pressure drop is negligible. This observation can be explained by the low permeability and low flow rate in ultratight reservoir. Therefore, we conclude that the horizontal wellbore pressure drop can be neglected in ultratight reservoir.Figure 21
Effect of horizontal wellbore pressure drop on pressure behavior. (a) Pressure curve of multistage fractured horizontal well. (b) An enlarged view of pressure curve within green circle.
(a)
(b)


### 3.6. Complex Fracture Network

It is inappropriate to simulate the production with biwing fracture when the characteristic of complex fracture network is evident [55–57]. Figure 22 shows the schematic of a complex fracture network, the red lines represent the hydraulic fractures, the blue lines are natural fractures, and black line denotes horizontal well. Parameters of the reservoir, fluid, and the fractures are shown in Table 3. Figure 23 shows the pressure distribution of the complex fracture network at different times. With increasing time, the area that contributes to production as well as the elastic production in a specific volume increases.Table 3
Basic parameters for complex fracture network case.
Properties
Value
Reservoir permeability,k, μm2
1 × 10−4
Reservoir porosity,ϕ
10%
Reservoir length,xe, m
600
Formation thickness,h, m
20
Total compressibility,Ct, MPa−1
2.76 × 10−3
Initial pressure,pi, MPa
30
Fluid viscosity,μ, mPa·s
20
Fluid density,ρ, kg/m3
900
Hydraulic fracture conductivity,Cdf, μm2·m
0.2
Natural fracture conductivity,Cdn, μm2·m
0.04
Production rate,Q, t/d
1.92Figure 22
Schematic of a complex fracture network. The black, red, and blue lines represent the wellbore, hydraulic fractures, and natural fractures, respectively.Figure 23
Pressure distribution of complex fracture network at different production time. (a)t=100 d, (b) t=500 d, and (c) t=1000 d.
(a)
(b)
(c)


## 3.1. Flow Regime Analysis

The transient behavior can be shown by type curves, which is employed to obtain the characteristic of formation and the reservoir fluids and to figure out different flow regimes. Figure10 shows the dimensionless well pressure and its derivative for a horizontal well with three fractures. The transient behavior can be divided into several flow periods.Figure 10
Pressure response of multistage fractured horizontal well (fracture numbern=3, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=8, and dimensionless reservoir length xeD=240).(
1) Linear Flow (Figure11(a)). Fluid flows linearly from the reservoir to the fracture and each fracture behaves independently. Both the slope of dimensionless pressure and derivative is 1/2 in this stage.Figure 11
Flow regimes for a multistaged fractured horizontal well: (a) linear flow; (b) early radial flow, (c) biradial flow, and (d) pseudo radial flow.
(a)
(b)
(c)
(d)(
2) Early Radial Flow (Figure11(b)). An early radial flow occurs around each fracture after the linear flow. This period mainly depends on the fracture length and fracture spacing. Besides, in this period, fractures still behave independently. The characteristic of this stage is a horizontal line of 1/(2N) in pressure derivative curve (N is the fracture stage). We can figure out that the value of the horizontal line plateau is 1/6.(
3) Biradial Flow (Figure11(c)). Fractures interact with each other and flow becomes elliptical to the wellbore. The slope of pressure derivative is 0.36.(
4) Pseudo Radial Flow (Figure11(d)). Fluid flows to the fracture-well system appear to be radial and flow across the outermost elements plays important part. The pressure derivative curve demonstrates a horizontal line of 0.5.(
5) Boundary Dominated Flow. In the closed system, the flow will reach pseudo-steady state. The pressure curve and derivate tend to merge and the slope equals 1.


## 3.2. Effect of Initial Fracture Conductivity

Figure12 illustrates the effect of initial conductivity on the pressure transient response and production rate. This figure shows that increasing hydraulic fracture conductivity results in an increase in well productivity; however, the incremental benefit decreases as the fracture conductivity increases. Figure 12 also shows that variations in the dimensionless pressure and the production rate tend to disappear as the flow regime approaches pseudo radial flow. The reason why the variations in these two parameters disappear is because the flow at late times is dominated by the outer reservoir, and therefore, the effect of fracture conductivity will be extremely small. Therefore, to utilize the potential of the hydraulic fractures, we must delay the occurrence of pseudo radial flow and that can be done by optimal well placement.Figure 12
Effect of initial fracture conductivity on pressure and rate behavior (fracture numbern=2, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=0.67, and dimensionless reservoir length xeD=120). (a) Pressure and derivative curve. (b) Dimensionless production rate.
(a)
(b)Although it is well documented that the fracture conductivity decreases rapidly at first and tends to converge asymptotically [31–34, 37], previous studies always assumed that hydraulic fractures do not vary with time. Below, we analyze the difference between time-independent and time-dependent fractures. As discussed earlier, the conductivity, and hence the productivity of the horizontal well, decreases with time. Therefore, the pressure response of the well is larger than the time-independent situation, as suggested by Figure 13(a). The flow rate of a horizontal well with time-independent conductivity fractures exceeds that with time-dependent conductivity fractures in the initial stage (Figure 13(b)), and that could be attributed to the dominant role of fractures. During late times of production, the fluid flow in the outer reservoir becomes dominant and negligible difference can be observed between time-dependent and time-independent cases. Therefore, for the optimal exploitation of an unconventional reservoir it is necessary to maintain the fracture conductivity, especially in the initial stage of the production.Figure 13
Effect of time-dependent fracture conductivity on pressure and rate behavior (fracture numbern=2, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=0.67, dimensionless reservoir length xeD=120, and decline coefficient β=0.11). (a) Dimensionless pressure and derivative curve. (b) Dimensionless production rate.
(a)
(b)Figure14 shows the effect of decline coefficient on the production performance of the horizontal well, where a larger decline coefficient results in smaller production rate and if the decline coefficient is too large the fracture permeability tends to quickly approach the formation permeability, which results in sharp decline in the production rate.Figure 14
Effect of decline coefficient on the production performance of the horizontal well (fracture numbern=2, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=0.67, and dimensionless reservoir length xeD=120).


## 3.3. Effect of Fracture Distribution

The pressure difference created during hydraulic fracturing allows the proppants to stay in the fractures away from the horizontal well heel [51], and as a result, the fracture conductivities are different that may further have an impact on the pressure behavior. Two scenarios of even and uneven proppant distribution can be seen as schematics in Figures 15(a) and 15(b), respectively, which are used to study the effect of proppant distribution.Figure 15
Sketch of proppant distribution in different fracture: (a) uneven distribution; (b) even distribution.
(a)
(b)For unevenly distributed fractures, the fracture nearest to the well toe has the largest conductivity, whereas the conductivity of the intermediate fracture equals that of the evenly distributed fractures.Figure16(a) shows the dimensionless pressure drop of a horizontal well (at constant production rate multistage fractured well) with evenly distributed and unevenly distributed fractures. The pressure response with unevenly distributed proppants is slightly larger than that with evenly distributed proppants at initial period but the gap disappears after that. Result of Figure 16(b) suggests that horizontal well with evenly distributed proppants produces with higher rate initially than with unevenly distributed proppants. Therefore, these two results indicate that well with evenly distributed proppants performs better than that with unevenly distributed proppants.Figure 16
Effect of unevenly distributed fracture conductivity on the pressure and production behavior (fracture numbern=3, dimensionless fracture half-length xfD=1, dimensionless fracture spacing dD=0.67, dimensionless reservoir length xeD=120. The conductivity of evenly distributed fracture is 0.375 μm2·m; for unevenly distributed fracture, the conductivities are 0.25 μm2·m, 0.375 μm2·m, and 0.5 μm2·m.). (a) Dimensionless pressure curve. (b) Dimensionless production rate curve.
(a)
(b)Figure17 shows the dimensionless flux rate of three fractures that have evenly and unevenly distributed proppants. When the fracture conductivities (or proppant distribution) are evenly distributed, the dimensionless flux rate of Fractures 1 and 3 increases with time, but the flux rate of Fracture 2 decreases sharply with time. At the initial stage, each fracture behaves independently; however, with time the fractures begin to interfere with each other and the flux rate of Fractures 1 and 3 tends to increase, whereas owing to the symmetry, the flux rate of intermediate fracture, Fracture 2, is hindered. If the fracture conductivities are unevenly distributed, the flux rate of the third fracture will first decrease and then increase with the time. However, production rate of the first fracture increases steadily with time and it crosses the conductivity of the second fracture at some point.Figure 17
Dimensionless flux rate of each fracture. (a) Fracture conductivity evenly distributed. (b) Fracture conductivity unevenly distributed.
(a)
(b)


## 3.4. Effect of Fracture Geometry

In hydraulic fracturing, the existent propped fractures result in the redistribution of local earth stresses. Moreover, microseismic measurements have proved that there is mechanical-stress interference between multiple transverse fractures. The stress-shadow effect can restrict the growth of the fracture in the middle section while promoting the growth of the fractures at the heel or the toe [49, 52]; therefore, different fracture geometry may exist in the reservoir. Three geometries shown in Figure 18, equilong type, spindle-shaped type, and dumbbell-shaped type, are studied.Figure 18
Sketch of different fracture geometries. (a) Equilong type. (b) Spindle-shaped type. (c) Dumbbell-shaped type.
(a)
(b)
(c)Figure19 shows the dimensionless flux rate of three fractures from three fracture geometries of Figure 18, respectively. For spindle-shaped geometry, the flux rate of the first and third fracture will first decrease and then slightly increase before levelling off at long time (Figure 19(b)). However, the flux of the second fracture always remains smaller than the other two fractures, although it has a larger fracture length. Figure 19(c) shows that, for dumbbell-shaped geometry, the flux rate of fractures near the heel and toe will slightly increase with time, whereas the flux of the intermediate fracture will sharply decrease with time. From these results it is evident that no matter what the fracture geometry is, the flux of the intermediate fracture always decreases with time. Figure 20 suggests that if the bottom-hole pressure is constant, the production rate of equilong fracture type is higher than the others because it has the longest effective interference fracture length. Apparently, this is the reason why equilong fracture type is usually preferred over other fracture geometries.Figure 19
Dimensionless flux rate of each fracture. (a) Equilong type,xf1=xf2=xf3=100 m. (b) Spindle-shaped type, xf1=75 m, xf2=150 m, and xf3=75 m. (c) Dumbbell-shaped type, xf1=120 m, xf2=60 m, and xf3=120 m.
(a)
(b)
(c)Figure 20
Dimensionless production rate for different fracture configurations.


## 3.5. Effect of Horizontal Wellbore Pressure Drop

Many previous studies considered the horizontal wellbore as an infinite wellbore; however, it has been shown that wellbore pressure drop exists in the production [53, 54]. We investigated the effect of horizontal wellbore pressure drop on pressure behavior as shown in Figure 21, which shows that the effect of considering wellbore pressure drop is negligible. This observation can be explained by the low permeability and low flow rate in ultratight reservoir. Therefore, we conclude that the horizontal wellbore pressure drop can be neglected in ultratight reservoir.Figure 21
Effect of horizontal wellbore pressure drop on pressure behavior. (a) Pressure curve of multistage fractured horizontal well. (b) An enlarged view of pressure curve within green circle.
(a)
(b)


## 3.6. Complex Fracture Network

It is inappropriate to simulate the production with biwing fracture when the characteristic of complex fracture network is evident [55–57]. Figure 22 shows the schematic of a complex fracture network, the red lines represent the hydraulic fractures, the blue lines are natural fractures, and black line denotes horizontal well. Parameters of the reservoir, fluid, and the fractures are shown in Table 3. Figure 23 shows the pressure distribution of the complex fracture network at different times. With increasing time, the area that contributes to production as well as the elastic production in a specific volume increases.Table 3
Basic parameters for complex fracture network case.
Properties
Value
Reservoir permeability,k, μm2
1 × 10−4
Reservoir porosity,ϕ
10%
Reservoir length,xe, m
600
Formation thickness,h, m
20
Total compressibility,Ct, MPa−1
2.76 × 10−3
Initial pressure,pi, MPa
30
Fluid viscosity,μ, mPa·s
20
Fluid density,ρ, kg/m3
900
Hydraulic fracture conductivity,Cdf, μm2·m
0.2
Natural fracture conductivity,Cdn, μm2·m
0.04
Production rate,Q, t/d
1.92Figure 22
Schematic of a complex fracture network. The black, red, and blue lines represent the wellbore, hydraulic fractures, and natural fractures, respectively.Figure 23
Pressure distribution of complex fracture network at different production time. (a)t=100 d, (b) t=500 d, and (c) t=1000 d.
(a)
(b)
(c)


## 4. Field Example

Well J172-H is a multistage fractured horizontal well in Jimsar sag, which is a typical tight oil region located in the southeast of the Junggar Basin of China. Tight oil in the sag is mainly in Permian Lucaogou Formation, which is divided into the first member and the second member from bottom to top (Figure24). Sweet spots are developed in the lower parts of both members, that is, the upper sweet spot and the lower sweet spot. J172-H is placed in the upper sweet spot, and the reservoir is characterized by alternating layers of fine dolarenite, dolomitic siltstone, and fine dolomitic mudstone in the vertical direction [58]. The interpretation of nuclear magnetic log shows that the porosity of the upper sweet spot lies between 0.061 and 0.258, with an average of 0.101, whereas the permeability of the upper sweet spot is in the range of 0.001 × 10−3μm2 to 0.284 × 10−3μm2, with an average of 0.012 × 10−3μm2. The permeability of over 90 percent of the samples is less than 0.1 × 10−3μm2.Figure 24
Tectonic location and comprehensive logging evaluation of Lucaogou Formation in Jimsar sag.
(a)
(b)J172-H went into production in September of 2012, and during the initial stage the production rate of J172-H was approximately 15 times that of the adjacent vertical well. A build-up well test was conducted in May of 2013 and Figure25 shows the pressure and pressure derivative data of that build-up test, which was conducted for a short period. Thus, it can be concluded that the data shown in Figure 25 comes from very beginning of the build-up test.Figure 25
Well build-up test data of J172-H.Due to the short test and ultralow reservoir permeability, some characteristic responses in the middle and late time period, corresponding to the biradial flow and boundary dominated flow regimes featured by the slope of 0.36 and 1 in the pressure derivative curve, cannot be found from the well test data. In the regimes of wellbore storage and early linear flow, the slopes of pressure derivative curves are equal to 1 and 0.5, respectively. We matched the curves with straight lines of slope 1 and 0.5 and identified the flow regime from the data as wellbore storage period and linear flow period. Note that the effect of wellbore storage is not taken into account in this study, because our model is proposed on the basis of Green’s function in real domain and it is difficult to obtain the derivative of pressure with respect to time in the inner boundary condition. However, the wellbore storage effect can be accounted for by transforming the equations from real domain to the Laplace domain by using the following equation:(18)p¯wD,storage,skin=sp¯wD+Ss+CDs2sp¯wD+S,where p¯wD,storage,skin stands for the dimensionless pressure including wellbore storage and skin effects (in Laplace space), p¯wD is the dimensionless pressure without these effects (in Laplace space), s is the Laplace variable, and CD and S are the wellbore storage coefficient and skin factor, respectively. This is the subject of future studies. Therefore, data from the linear flow period is employed in our analysis. A physical experiment was conducted to estimate the decline coefficient of hydraulic fractures, and the resulting value is 0.106. Figure 26 shows the match for linear flow regime, from which we estimate the fracture and reservoir parameters. In Table 4, we compare the fitting parameters from our proposed model and the field data from microseismic monitoring. Our estimated average half-length of these fractures approximates to 148 m, which lies in the range of 110–230 m. The fracture spacing obtained from our match, 81 m, approximates to the average value from microseismic map. Because the reservoir heterogeneity is not taken into account in our proposed model, we only get the mean permeability of the formation, which is fairly comparable to the well logging results.Table 4
Comparison between fitting parameter and the field data.
Parameter
Our model
Field data
Fracture half-length, m
148
110~230
Fracture spacing, m
81
78.5 (in average)
Formation permeability, 10−3 μm2
0.0107
0.001~0.284Figure 26
Curve match for linear flow.


## 5. Conclusions

A semianalytical model was proposed to analyze the pressure behavior of multistage fractured horizontal well in tight oil reservoirs. Factors that influence the pressure behavior of multistage fractured horizontal well were analyzed to provide a deep understanding of the pressure transient behavior. Following conclusions were reached.(1) Higher hydraulic fracture conductivity will result in an increase in well productivity. Therefore, hydraulic fracture conductivities should be optimized according to the well performance and investment.(2) The stimulated region has more significant importance on the transient pressure and rate behavior than the outer reservoir region. In order to make use of the hydraulic fractures, the existence of pseudo radial flow must be delayed through optimal well placement.(3) The comparison among different factors suggests that fracture conductivity is the most important factor when exploiting the full potential of a horizontal well. The decline of fracture conductivity has a tremendous influence on the well performance. Therefore, the proppant selection should be given the first priority to maintain the fracture conductivity.(4) The interference among fractures of equilong type is the strongest. This geometry is favorable for both fracture flux rate and total production rate.(5) Due to the low permeability and flow rate of ultratight reservoir, the effect of wellbore pressure drop is negligible; therefore, the horizontal wellbore pressure drop can be ignored in ultratight reservoir.



---

*Source: 5279792-2017-09-28.xml*