A Novel Model Incorporating Geomechanics for a Horizontal Well in a Naturally Fractured Reservoir

Fracture aperture of a fractured reservoir can be affected by both matrix elasticity and fracture compressibility when the reservoir pressure decreases, namely stress sensitivity. An elasticity parameter coupling Young’s modulus and Poisson’s ratio was introduced to reflect this geomechanical behavior, and a new model incorporating geomechanics was developed to analyze the flow behavior of a horizontal well in a naturally fractured reservoir. Pressure solutions for two cases—uniform-flux and infinite-conductivity—were derived, respectively. For the uniform-flux case, the effect of dimensionless elasticity parameter on the pressure-drop profile becomes stronger with continuing production, and the profile may be like a bow. Nine flow regimes can be observed on the transient response of the infinite-conductivity case. Stress sensitivity mainly affects the late-flow period and a larger dimensionless elasticity parameter causes a greater pressure drop. Due to stress sensitivity, the pressure derivative curve exhibits an upward tendency in the pseudo-radial flow regime, and the slope is greater than “1” in the pseudo-steady flow regime. For KT-I formation in the North Truva field, its elasticity parameter decreases with the increase of Young’s modulus or Poisson’s ratio and ranges from 8× 10−8 Pa−1 to 1.1× 10−7 Pa−1. Meanwhile, the transient response of H519 has a slight negative correlation with Young’s modulus and Poisson’s ratio in the pseudo-steady flow regime.


Introduction
Horizontal wells have been an important technique to improve oil recovery, and naturally fractured reservoirs are one of the potential targets to use this kind of well type [1][2][3][4].Reservoir pressure depletion in the naturally fractured reservoirs changes the stress balance of reservoir rock and then changes fracture permeability, namely stress sensitivity, which is a typical geomechanical behavior of naturally fractured reservoirs [5][6][7][8][9].Therefore, it is of great significance to investigate the flow behavior of a horizontal well in the naturally fractured reservoir, especially after a long period of production.
The dual-porosity model, including a natural fracture system and a matrix system, is always used to simulate the fluid flow in a naturally fractured medium [10,11].The geometric characteristics of natural fractures may vary over a range of scales.They can increase reservoir capacity and conductivity, but also can act like a barrier to fluid flow [12,13].Based on the characteristics of the dual-porosity model, three classical models are developed to deal with the fluid-flow problems in naturally fractured reservoirs, including the sugar-cube model presented by Warren and Root [14], the slab-shaped matrix model presented by Kazemi [15] and the sphere-shaped matrix model presented by de Swaan [16].Thus, we can easily further investigate the transient flow of various well types in naturally fractured Energies 2018, 11, 2584 2 of 28 reservoirs.However, in these models, stress sensitivity is neglected and petrophysical parameters of natural fracture system are considered as static properties during the production, which obviously does not match the characteristics of naturally fractured reservoirs.
As mentioned before, fracture porosity and fracture permeability are variable since they are related to reservoir pressure.Several approaches have been proposed to incorporate geomechanical effects into the simulation and characterization of elastic media in the literature [17][18][19][20][21][22].Nur and Yilmaz [17] defined a stress sensitivity coefficient to account for permeability variation, and Kikani and Pedrosa [18] further rewrote it as a form of exponential function.This exponential relationship becomes one of the most common expressions and is widely used to describe stress-sensitive reservoirs.On the basis of this relationship, Ceils et al. [23] and Samaniego and Villalobos [24] investigated the effect of stress sensitivity on the transient flow behavior of naturally fractured reservoirs.Similarly, Yao et al. [25] presented a new semi-analytical model for hydraulic fractures with stress-sensitive conductivity, and Zhang et al. [26] also developed an improved model for a finite-conductivity fracture with stress-sensitive fracture permeability.These proposed models highly depend on the exponential relationship between fracture permeability and reservoir pressure.However, stress sensitivity coefficient is mainly obtained by fitting the experimental data measured from the laboratory and it may have a certain error due to the difference between the experimental environment and the actual reservoir conditions [8,17,18].There are still other approaches to deal with stress sensitivity.Raghavan and Chin [20] established a numerical model coupling geomechanical and fluid-flow aspects to examine the fracture permeability reduction from the partial closure of natural fractures.Rosalind [27] proposed an approach that avoids using pseudo-pressure function to handle stress sensitivity, but this approach is limited to a linear variation of permeability and porosity.
Generally, in terms of stress sensitivity, most proposed models depends on either an exponential relationship or a numerical analysis.Depending on these models, many researchers investigate the transient pressure of vertical wells or fractured wells in stress-sensitive reservoirs [20,[23][24][25][26].However, they do not take the reservoir geomechanical parameters into account.For example, Young's modulus and Poisson's ratio, two key parameters reflecting the reservoirs' geomechanical properties, are not incorporated into these models [28].By using the theory of geomechanics, Jabbari et al. [22] proposed a fracture permeability model to figure out the dependency of fracture properties on reservoir pressure change, which couples fracture porosity, fracture compressibility, Young's modulus and Poisson's ratio.This model provides a completely new method to deal with stress sensitivity and consider geomechanical properties.To be our best knowledge, no research has been done upon horizontal well incorporating nonlinear geomechanics with the model proposed by Jabbari et al. [22], and the corresponding analytical solution for this fluid-flow problem is not available.
The objective of this work was to develop a model incorporating geomechanics for a horizontal well in a naturally fractured reservoir.Firstly, by using a Pedrosa transform, perturbation technique and spatial integral method, analytical pressure solutions for two cases: uniform flux and infinite conductivity, were derived, respectively.Secondly, a special case was employed to verify this uniform-flux pressure solution and the effect of dimensionless elasticity parameter on this horizontal wellbore's pressure drop profiles was also discussed.Thirdly, based on the infinite-conductivity pressure solution, transient pressure, flow regimes and parameters sensitivity of this model were investigated in detail.Finally, by taking KT-I formation of the North Truva field as an example, the effects of Young's modulus and Poisson's ratio of this formation on the elasticity parameter and the transient response of H519 were analyzed.

Physical Model
Figure 1 shows a sketch of a horizontal well in a naturally fractured reservoir with stress sensitivity.In order to simplify the mathematical model, basic assumptions are made as follows: Energies 2018, 11, 2584 3 of 28 (1) This formation consists of two systems, including natural fracture system and matrix system, and can be represented by Warren-Root model.Two kinds of porous media are homogeneous, and the closure of natural fractures caused by reservoir pressure depletion is taken into account.(2) This reservoir is circular and closed at the external boundary of side and bounded by upper and lower impermeable formation, and has a uniform thickness.(3) Flow in this reservoir is considered to be a slightly compressible and single-phase fluid with constant viscosity and obeys Darcy's law.Pseudo-steady flow from matrix system to natural fracture system is assumed.(4) This horizontal well is produced at a constant rate and no fluid is assumed to flow at the tip of the horizontal wellbore.(5) Gravity and capillary effects are neglected.
Energies 2018, 11, x FOR PEER REVIEW 3 of 29 (1) This formation consists of two systems, including natural fracture system and matrix system, and can be represented by Warren-Root model.Two kinds of porous media are homogeneous, and the closure of natural fractures caused by reservoir pressure depletion is taken into account.
(2) This reservoir is circular and closed at the external boundary of side and bounded by upper and lower impermeable formation, and has a uniform thickness.(3) Flow in this reservoir is considered to be a slightly compressible and single-phase fluid with constant viscosity and obeys Darcy's law.Pseudo-steady flow from matrix system to natural fracture system is assumed.(4) This horizontal well is produced at a constant rate and no fluid is assumed to flow at the tip of the horizontal wellbore.(5) Gravity and capillary effects are neglected.

Fracture Permeability vs. Mechanical Elasticity
Warren and Root [14] first proposed a sugar-cube model to describe the naturally fractured medium (Figure 2).Based on this conceptual model, we can obtain the fracture porosity for a naturally fractured reservoir:

Fracture Permeability vs. Mechanical Elasticity
Warren and Root [14] first proposed a sugar-cube model to describe the naturally fractured medium (Figure 2).Based on this conceptual model, we can obtain the fracture porosity for a naturally fractured reservoir: (2) This reservoir is circular and closed at the external boundary of side and bounded by upper and lower impermeable formation, and has a uniform thickness.(3) Flow in this reservoir is considered to be a slightly compressible and single-phase fluid with constant viscosity and obeys Darcy's law.Pseudo-steady flow from matrix system to natural fracture system is assumed.(4) This horizontal well is produced at a constant rate and no fluid is assumed to flow at the tip of the horizontal wellbore.(5) Gravity and capillary effects are neglected.

Fracture Permeability vs. Mechanical Elasticity
Warren and Root [14] first proposed a sugar-cube model to describe the naturally fractured medium (Figure 2).Based on this conceptual model, we can obtain the fracture porosity for a naturally fractured reservoir:  Considering that the fracture aperture is much smaller than the block size, Equation (1) can be approximated as: In the sugar-cube model, assuming that the ratio of fracture width to block height is very small and fluid flows only in y-direction in both the vertical and horizontal fractures, the flow paths in Figure 2 can be approximated by flow through rectangular ducts.Janna [29] proposed the following equation to describe the average velocity of laminar flow through a rectangular duct: For a duct, the average velocity also can be described as a function of flow rate and fracture dimensions: Equating Equations ( 3) and ( 4), we obtain: Thus, the total flow through multiple fractures becomes: Further, Darcy's Law for linear horizontal flow through a porous medium can be applied: Combining Equation (6) with Equation ( 7), fracture permeability with sugar cube geometry is yielded [22]: Obviously, regarding Equation ( 8), fracture permeability is a stronger function of the fracture width, b, than the matrix block size, a. Jabbari et al. [22] take the derivative of Equation ( 8) with respect to reservoir pressure and let fracture width be a function of reservoir pressure, and then obtain the following equation: During the closing of fracture aperture, a positive change in the length of a matrix block corresponds to a negative change in the fracture width, thus, we get a relationship: ∆a = −∆b.Substituting this formula into Equation (9), it yields: Jabbari et al. [22] came up with the idea that the total change in fracture width, ∆b, caused by the change of reservoir pressure, is the sum of the change caused by fracture compressibility and matrix elasticity.On the basis of the general relationship between effective stress, total stress, and fracture pore pressure [30,31], Jabbari et al. [22] presented a formula of the change in fracture width caused by fracture compressibility: where ∆p f = p − p i is the change in fracture pore pressure.Similarly, using the theory of geomechanics, Jabbari et al. [22] obtained the change in fracture width that comes from matrix elasticity: where E is Young's modulus, and v is Poisson's ratio.Note that as a feasible assumption, the change in overburden pressure is neglected.Hence, the total change in fracture width is: where p ob is the overburden pressure.
The differential form of Equation ( 13) can be written as: Further, substituting Equation ( 14) into Equation (10) and neglecting any change in overburden pressure, Jabbari et al. [22] developed an expression of fracture permeability incorporating geomechanics: If overburden pressure is held constant, Equation ( 15) is reduced to the following form: Equation ( 16) is a representation of fracture permeability dependent on reservoir pressure, which considers two significant geomechanical parameters, Young's modulus and Poisson's ratio.

Governing Equations
According to the law of conservation of mass, we can have the continuity equations for natural fracture system and matrix system: respectively.
Energies 2018, 11, 2584 6 of 28 Pseudo-steady flow has been assumed between natural fracture system and matrix system, and the inter-porosity flow equation can be given as: where α is a parameter related to the geometry of natural fractures.With Equation (19), Equations ( 17) and ( 18) can be transformed into: The porosity change of natural fracture and matrix during reservoir production can be described by the following equations, respectively: where φ f 0 and φ m0 are the initial porosity of natural fracture system and matrix system; c f and c m are the compressibility of natural fracture system and matrix system.If the connate-water saturation in the matrix medium is assumed to be negligible, we can obtain a reasonable estimate expression of fracture compressibility: where c o is the oil compressibility.Using Equations ( 22) and ( 23), we can obtain: where φ f 0 c f and φ m0 c m are the storage of natural fracture system and matrix system, respectively.Substituting Equations ( 25) and (26) into Equations ( 20) and ( 21), we have the governing equations for the point sink model in the naturally fractured reservoir: Considering that fracture permeability depends on reservoir pressure, we can transform the term on the left-hand side of Equation ( 27) into: Energies 2018, 11, 2584 7 of 28 With Equation ( 16), the second term on the right-hand side of Equation ( 29) can be rewritten as: Similarly, by introducing Equation ( 24), the last term on the right-hand side of Equation ( 29) can be further rewritten as: Substituting Equations ( 30) and (31) into Equation ( 29), we have: To simplify the model derivation, we define a new elasticity parameter: Obviously, this parameter depends on oil compressibility, fracture porosity, fracture compressibility, Young's modulus and Poisson's ratio, thus, Equation (32) can be simplified as: In the cylindrical coordinate system, we have: Combining Equations ( 34)- (36) with Equation ( 27), the continuity equation for natural fracture system becomes:

Initial Conditions and Boundary Conditions
According to the physical model above, we can obtain the corresponding constraint conditions for the point sink model as follows: Initial conditions: Inner boundary condition: where q point (t) is the production rate of point sink.
Outer boundary condition: Upper and lower boundary conditions: In the cylindrical system, Equations ( 28) and ( 37)-( 41) compose a point sink model incorporating geomechanical property in the naturally fractured reservoir.

Point Sink Model Solution
In order to unify and solve the model easily, the dimensionless variables defined in Table 1 are used to deal with the above equations.

Parameters Definitions
Dimensionless production rate of point sink q pointD (t D ) = q point (t) q Dimensionless wellbore storage coefficient Based on these dimensionless variables, Equations ( 28) and ( 37)-( 41) can be transformed into the following dimensionless forms: Energies 2018, 11, 2584 9 of 28 Governing equations: Initial conditions: Inner boundary condition: Outer boundary condition: Upper and lower boundary conditions:

.1. Pedrosa Transform and Perturbation Transform
Due to the consideration of geomechanics in the governing equation, Equation (42), we obtain a strong nonlinear differential system in the cylindrical system, which causes difficulty in solving the above point sink model.Here we take Pedrosa transform to weaken the strong nonlinearity in Equation ( 42) and further apply the perturbation theory in this case to obtain an approximate solution.Pedrosa [32] provided a useful transform to deal with the nonlinear equation: Using Equation (48), we can obtain a series of differential relationships in the cylindrical system as follows: (49) Energies 2018, 11, 2584 10 of 28 Substituting Equations ( 48)-(53) into Equations ( 42)-(47), we have: The strong nonlinearity in Equation ( 42) is substantially weakened in Equation ( 54), and by taking perturbation transform [33,34], we can have the approximation of point sink model's full solution: For small ε D , these higher-order terms in the series become successively smaller and the zero-order perturbation solution can satisfy the accuracy requirement.Then, Equations ( 54)-( 59) can be transformed into:

Laplace Transform and Fourier Transform
By employing a Laplace transform of Equations ( 61)-(66) with respect to t D , we can obtain the following equations in the Laplace domain: Combining Equation (67) with Equation (68), we can obtain: Applying the Taylor series expansion to Equation (48), it becomes: Similarly, for small ε D the zero-order perturbation solution can satisfy the accuracy requirement, thus, we have in Laplace domain: Substituting Equation (74) into Equation ( 72), the governing equation can be simplified as: where Equations ( 69)-( 71) and (75) compose the zero-order mathematical model of a point sink in the naturally fractured reservoir in Laplace domain.Essentially, this zero-order model is the same as the conventional point-source model in the closed cylindrical system presented by Ozkan [2].
In order to obtain the zero-order point sink solution, we can take the following finite Fourier cosine transform with respect to z D for this model: The inverse of this finite Fourier cosine transform with respect to z D is: where Applying the finite Fourier cosine transform to Equation (75) and combining with Equation (71), we obtain: Similarly, taking a Fourier transform for the boundary conditions, Equations ( 69) and (70) become:

Model Solution
The governing equation, Equation (78), is a modified Bessel equation.Combining with Equations ( 78)-( 80), the general solution for the zero-order point sink model can be obtained: where With the inverse of finite Fourier cosine transform given by Equation (77), the general solution of a point sink in Laplace domain can be expressed as: Equation ( 82) is the point sink model solution in the closed cylindrical system.As shown in Figure 1, the horizontal wellbore is deployed along the x-axis (y wD = 0), so we can get the following relationship: r D = (x D − x wD ) 2 + y 2 D .Further, substituting this relationship into Equation (82), we can rewrite it as:

Horizontal Well Model Incorporating Geomechanics
When fluid flow along the horizontal wellbore is uniform, by integrating with respect to x wD from −1 to 1 in Equation ( 83) and letting y D = 0 and z D = z wD + r wD L D [2,35,36], we can obtain the pressure distribution along the horizontal wellbore face in the naturally fractured reservoir: Due to the assumption of uniform flux, the pressure along the horizontal wellbore is not uniform and changes with the position and time.To represent the pressure behavior of a horizontal wellbore, no pressure drop in its interior is always assumed, which means that pressure is uniform along the wellbore, and this kind of well is said to have infinite conductivity.This inner-boundary condition poses a very difficult boundary-value problem.Generally, the infinite-conductivity solution can be obtained by selecting a fixed point along the wellbore in the corresponding uniform-flux solution.Clonts and Ramey [37] proposed that the uniform-flux solution with a dimensionless distance of 0.732 was a good approximation of the infinite-conductivity case, and this fixed point was taken by Ozkan [2] to evaluate the pressure for the infinite-conductivity horizontal well.However, the uniform-flux solution at a dimensionless distance of 0.7 or 0.68, presented by Daviau et al. [38] and Rosa and de Carvalho [4], is also used to simulate the flow behavior of infinite-conductivity horizontal well.Therefore, there is no definite result on the selection of the fixed point at present.In this study, the pressure averaging technique based on the uniform-flux solution, proposed by Kuchuk et al. [39], was used to approximate the infinite-conductivity solution.Firstly, integrate along the horizontal wellbore (x-axis) with respect to x D from −1 to 1 in Equation (84).Then, divide the integral result by the dimensionless wellbore length.Finally, we obtain an average wellbore pressure for the infinite-conductivity horizontal well: Further, the effects of wellbore storage and skin can be easily incorporated into the uniform-flux solution (Equation ( 84)) and the infinite-conductivity solution (Equation ( 85)), and we can get the following expressions [40,41]: and: respectively.By applying the Stehfest numerical algorithm [42] to Equations ( 86) and (87), the pressure solution in real domain can be obtained.Then, with Pedrosa transform, Equation (48), we can get the pressure response, including the uniform-flux and infinite-conductivity case, for a horizontal well incorporating geomechanics in the naturally fractured reservoir: and: According to the definition of dimensionless elasticity parameter ε D , when ε D = 0 it means no consideration of stress sensitivity in the model.We can obtain the pressure solution under this condition by calculating the limit in Equations ( 88) and (89) as ε D approaches zero: and:

Verification of the Solution
de Carvalho and Rosa [3] presented an infinite-conductivity solution for a horizontal well in the naturally fractured reservoir, but did not consider the effect of the geomechanics, which just corresponds to ε D = 0 in our model.Note that their pressure was calculated at a fixed point along the wellbore, x D = 0.68, by using the uniform-flux solution.They presented a series of type curves for a horizontal well in a dual-porosity reservoir.Thus, we can let x D = 0.68 in our uniform-flux solution, Equation (90), to verify our uniform-flux solution.
Due to the difference of the reference length in dimensionless variables between two models, a reasonable transformation of the parameters value should be conducted before the verification.Base on the data presented in the work of de Carvalho and Rosa [3], the parameter values used for the verification are listed in Table 2.Meanwhile, the reservoir boundary in their model is assumed infinite, which is different from our outer boundary.Therefore, we can only select the result in the early time (t D ≤ 10 3 ) to verify our solution.As shown in Figure 3, there is an excellent agreement on pressure response at x D = 0.68 between our result and the work of de Carvalho and Rosa [3] for different parameters value.It is also verified that the method in this work can calculate accurately.

Parameters in this Work
Values Parameters in Figure 5

Effect of Dimensionless Elasticity Parameter on Horizontal Wellbore Pressure-Drop Profiles
For the uniform-flux case, pressure drop is not uniform along the wellbore face.By using the data in Table 3, Figure 4 discusses the effect of dimensionless elasticity parameter on the wellbore

Effect of Dimensionless Elasticity Parameter on Horizontal Wellbore Pressure-Drop Profiles
For the uniform-flux case, pressure drop is not uniform along the wellbore face.By using the data in Table 3, Figure 4 discusses the effect of dimensionless elasticity parameter on the wellbore pressure-drop profiles at different time.and 0.05 should also be like a bow, but their shapes are not very obvious due to a large scale of wD p .
Compared with Figure 4c, as D ε increases, the pressure drop at the same wellbore position D x increases more obviously at this time.Overall, we can conclude that the effect of dimensionless elasticity parameter on the pressure-drop profiles becomes stronger with the stress sensitivity increasing and the production going on, for the reason of a larger fracture permeability reduction from the closure of natural fractures.
Table 3. Basic data used in Figure 4.  Figure 4a,b present the pressure drop profiles with respect to horizontal wellbore position x D from −1 to 1 for different ε D at t D = 10 3 , 10 5 , respectively.We can find that the curves for different ε D almost overlap together in each figure, which means the effect of stress sensitivity is very weak and can be neglected at small time.Meanwhile, the pressure drop profiles in Figure 4a,b are close to a horizontal line and only the pressure drop at both ends of the wellbore is slightly lower than that at the rest of the wellbore.It demonstrates that the ends of the wellbore cause a lower pressure drop under the uniform-flux condition.Figure 4c shows the horizontal wellbore's pressure-drop profiles at t D = 10 7 .It can be seen from Figure 4c that when the dimensionless production time increases to 10 7 , the pressure drop at the same wellbore position x D slightly increases with the increase of ε D , so a stronger stress sensitivity will lead to a larger pressure drop at this time.All the p wD − x D curves in Figure 4c are symmetric around the vertical line, x D = 0, and form a bow shape, and a farther position away from the center of the wellbore corresponds to a lower pressure drop.Figure 4d presents the pressure drop profiles for different ε D at t D = 10 9 .Like the profiles in Figure 4c, the p wD − x D curve for ε D = 0.1 in Figure 4d also looks like a bow.The p wD − x D curves for ε D = 0, 0.01 and 0.05 should also be like a bow, but their shapes are not very obvious due to a large scale of p wD .Compared with Figure 4c, as ε D increases, the pressure drop at the same wellbore position x D increases more obviously at this time.Overall, we can conclude that the effect of dimensionless elasticity parameter on the pressure-drop profiles becomes stronger with the stress sensitivity increasing and the production going on, for the reason of a larger fracture permeability reduction from the closure of natural fractures.

Comparison with a Model Not Considering Geomechanics
Since we have established a model of a horizontal well incorporating geomechanics, it is necessary to make a comparison with a model not considering geomechanics.In this study, let ε D equal to 0.1 and 0, which represent the model incorporating and not incorporating geomechanics, respectively.Based on the infinite-conductivity solutions, Equation ( 89) and (91), Figure 5 presents the comparison result of the wellbore pressure and pressure derivative curves of these two models.

Comparison with a Model Not Considering Geomechanics
Since we have established a model of a horizontal well incorporating geomechanics, it is necessary to make a comparison with a model not considering geomechanics.In this study, let D ε equal to 0.1 and 0, which represent the model incorporating and not incorporating geomechanics, respectively.Based on the infinite-conductivity solutions, Equation ( 89) and (91), Figure 5 presents the comparison result of the wellbore pressure and pressure derivative curves of these two models.As shown in Figure 5, in the early-and middle-flow periods (

Flow Regime Recognition on Transient Wellbore Pressure Curves
The identification of flow regimes is one of the most important procedures of transient pressure analysis.Similarly, we use the infinite-conductivity solution to analyze this case's flow regimes.Figure 6 shows the transient pressure and pressure derivative responses of a horizontal well incorporating geomechanics in a naturally fractured reservoir ( D ε = 0.05) and nine typical flow regimes can be observed in Figure 6: (1) Pure wellbore storage regime.At small dimensionless time, flow is mainly governed by the wellbore storage effect, and the transient pressure and pressure derivative curves align in an upward straight line with a unit slope.(2) First transition flow.This regime is marked by a hump on the pressure derivative curve.
(3) Radial flow in the natural fracture system.A short radial flow may occur after a transition flow if the reserve ratio in the natural fracture system is large enough.As presented in Figure 6, the As shown in Figure 5, in the early-and middle-flow periods (t D ≤ 10 −1 ), the transient pressure almost overlap together, and the pressure derivative curves almost overlap together, too.As the production time gets longer, especially in the late-flow period, the pressure drop of the model incorporating geomechanics becomes larger than that not incorporating geomechanics, so does the pressure derivative.It can be interpreted by that for the model incorporating geomechanics, the closure of natural fractures gets more and more serious, and more reservoir energy will be consumed to produce the same amount of oil.

Flow Regime Recognition on Transient Wellbore Pressure Curves
The identification of flow regimes is one of the most important procedures of transient pressure analysis.Similarly, we use the infinite-conductivity solution to analyze this case's flow regimes.Figure 6 shows the transient pressure and pressure derivative responses of a horizontal well incorporating geomechanics in a naturally fractured reservoir (ε D = 0.05) and nine typical flow regimes can be observed in Figure 6: (1) Pure wellbore storage regime.At small dimensionless time, flow is mainly governed by the wellbore storage effect, and the transient pressure and pressure derivative curves align in an upward straight line with a unit slope.(2) First transition flow.This regime is marked by a hump on the pressure derivative curve.
(3) Radial flow in the natural fracture system.A short radial flow may occur after a transition flow if the reserve ratio in the natural fracture system is large enough.As presented in Figure 6, the pressure derivative curve in this regime is a horizontal line with a value of "1/(4L D )" and the corresponding streamline distribution in the reservoir is shown in Figure 7a.(4) Inter-porosity flow.This regime represents the cross flow from matrix system to natural fracture system.The pressure derivative curve in this regime shows a concave-shaped segment, which is the main characteristic of naturally fractured reservoirs.(5) Radial flow in the compound system.When the inter-porosity flow finishes, fluid flow in two systems reaches a dynamic balance, and a second radial flow begins to occur in the compound system.A horizontal line with a value of "1/(4L D )" can be also observed on the pressure derivative curve in this regime and its streamline distribution is the same with that in the third flow regime (as shown in Figure 7a).( 6) Linear flow in the compound system.In this regime, the dominating flow is linear flow in the horizontal plane, and the streamlines are parallel to the upper and lower boundary (as shown in Figure 7b).Its pressure derivative curve shows an upward straight line with a half slope.( 7) Second transition flow.As the drainage area expands, a second transition flow may occur in the reservoir before it reaches a pseudo-radial flow in the horizontal plane.( 8) Pseudo-radial flow in the compound system.When reservoir radius is large enough, a pseudo-radial flow in the horizontal plane may be formed after the second transition flow (as shown in Figure 7c).Unlike the conventional reservoir, due to the effect of stress sensitivity, the pressure derivative curve is not a horizontal line with a value of "0.5" but exhibits an upward tendency.(9) Pseudo-steady flow in the compound system.When the transient wave reaches the outer boundary, both the transient pressure curve and pressure derivative curves go up rapidly, and their slopes are greater than "1" due to stress sensitivity (as shown in Figure 6).
Energies 2018, 11, x FOR PEER REVIEW 18 of 29 pressure derivative curve in this regime is a horizontal line with a value of " ( ) L " and the corresponding streamline distribution in the reservoir is shown in Figure 7a.(4) Inter-porosity flow.This regime represents the cross flow from matrix system to natural fracture system.The pressure derivative curve in this regime shows a concave-shaped segment, which is the main characteristic of naturally fractured reservoirs.( 5) Radial flow in the compound system.When the inter-porosity flow finishes, fluid flow in two systems reaches a dynamic balance, and a second radial flow begins to occur in the compound system.A horizontal line with a value of " ( ) L " can be also observed on the pressure derivative curve in this regime and its streamline distribution is the same with that in the third flow regime (as shown in Figure 7a).( 6) Linear flow in the compound system.In this regime, the dominating flow is linear flow in the horizontal plane, and the streamlines are parallel to the upper and lower boundary (as shown in Figure 7b).Its pressure derivative curve shows an upward straight line with a half slope.( 7) Second transition flow.As the drainage area expands, a second transition flow may occur in the reservoir before it reaches a pseudo-radial flow in the horizontal plane.(8) Pseudo-radial flow in the compound system.When reservoir radius is large enough, a pseudo-radial flow in the horizontal plane may be formed after the second transition flow (as shown in Figure 7c).Unlike the conventional reservoir, due to the effect of stress sensitivity, the pressure derivative curve is not a horizontal line with a value of "0.5" but exhibits an upward tendency.(9) Pseudo-steady flow in the compound system.When the transient wave reaches the outer boundary, both the transient pressure curve and pressure derivative curves go up rapidly, and their slopes are greater than "1" due to stress sensitivity (as shown in Figure 6).

Analysis of Parameters' Influence on Transient Pressure Behavior
In this section, we investigated several important parameters' influence on the transient response of the proposed model in this work.Note that these investigations are based on the infinite-conductivity solution proposed in this work.The influence of wellbore storage coefficient, skin factor, interporosity flow coefficient and fracture storage capacity, has been discussed in the literature, and will not be further discussed here.

Dimensionless Elasticity Parameter
Figure 8 shows the influence of dimensionless elasticity parameter D ε on the transient response.Its influence is mainly concentrated in the late-flow period.For the case not incorporating geomechanics ( 0 D ε = ), the pressure derivative curve is a horizontal line with a value of 0.5 in the pseudo-radial regime, and is an upward straight line with a unit slope in the pseudo-steady flow regime.However, when the geomechanics is incorporated, the pressure derivative curve bends upward in the late-flow period.It can be concluded that with the increase of D ε , the slope of the derivative curve in the late-flow period becomes larger, which means a greater pressure drop.

Analysis of Parameters' Influence on Transient Pressure Behavior
In this section, we investigated several important parameters' influence on the transient response of the proposed model in this work.Note that these investigations are based on the infinite-conductivity solution proposed in this work.The influence of wellbore storage coefficient, skin factor, interporosity flow coefficient and fracture storage capacity, has been discussed in the literature, and will not be further discussed here.

Dimensionless Elasticity Parameter
Figure 8 shows the influence of dimensionless elasticity parameter ε D on the transient response.Its influence is mainly concentrated in the late-flow period.For the case not incorporating geomechanics (ε D = 0), the pressure derivative curve is a horizontal line with a value of 0.5 in the pseudo-radial regime, and is an upward straight line with a unit slope in the pseudo-steady flow regime.However, when the geomechanics is incorporated, the pressure derivative curve bends upward in the late-flow period.It can be concluded that with the increase of ε D , the slope of the derivative curve in the late-flow period becomes larger, which means a greater pressure drop.

Analysis of Parameters' Influence on Transient Pressure Behavior
In this section, we investigated several important parameters' influence on the transient response of the proposed model in this work.Note that these investigations are based on the infinite-conductivity solution proposed in this work.The influence of wellbore storage coefficient, skin factor, interporosity flow coefficient and fracture storage capacity, has been discussed in the literature, and will not be further discussed here.

Dimensionless Elasticity Parameter
Figure 8 shows the influence of dimensionless elasticity parameter D ε on the transient response.Its influence is mainly concentrated in the late-flow period.For the case not incorporating geomechanics ( 0 D ε = ), the pressure derivative curve is a horizontal line with a value of 0.5 in the pseudo-radial regime, and is an upward straight line with a unit slope in the pseudo-steady flow regime.However, when the geomechanics is incorporated, the pressure derivative curve bends upward in the late-flow period.It can be concluded that with the increase of D ε , the slope of the derivative curve in the late-flow period becomes larger, which means a greater pressure drop.9 presents the effect of dimensionless reservoir radius r eD on the transient pressure behavior.Due to the closed boundary, this parameter also mainly influences the oil flow in the late-flow period, and a smaller dimensionless reservoir radius will cause the final pseudo-steady flow to occur earlier.We also can see in Figure 9 that in the pseudo-steady flow regime, at the same production time, in terms of the pressure and pressure derivative, their differences between the models incorporating and not incorporating geomechanics, are increasing significantly with the decrease of dimensionless reservoir radius.It can be interpreted by that, with the production going on, the continuous fracture closure is more serious for the reservoir with a smaller dimensionless reservoir radius, and the effect of stress sensitivity on the transient response becomes stronger, too.

Dimensionless Reservoir Radius
Figure 9 presents the effect of dimensionless reservoir radius eD r on the transient pressure behavior.Due to the closed boundary, this parameter also mainly influences the oil flow in the late-flow period, and a smaller dimensionless reservoir radius will cause the final pseudo-steady flow to occur earlier.We also can see in Figure 9 that in the pseudo-steady flow regime, at the same production time, in terms of the pressure and pressure derivative, their differences between the models incorporating and not incorporating geomechanics, are increasing significantly with the decrease of dimensionless reservoir radius.It can be interpreted by that, with the production going on, the continuous fracture closure is more serious for the reservoir with a smaller dimensionless reservoir radius, and the effect of stress sensitivity on the transient response becomes stronger, too.Obviously, it has a great influence on the early-flow period, including the following regimes: first transition flow, radial flow in the natural fracture system, inter-porosity flow and radial flow in the compound system, but has no influence on the middle-and late-flow periods.There is an implicit relationship between dimensionless wellbore radius, wellbore length and reservoir thickness, In this study, let D h equal to 200, 250 and 400, and the corresponding D L is 300, 400 and 500, respectively.As shown in Figure 10, the concave-shaped segment gradually uplifts with the increase of dimensionless reservoir thickness and the value of the pressure derivative in two radial flow regimes is 1/80, 1/64 and 1/40, respectively.However, the change of dimensionless reservoir thickness does not enhance or weaken the effect of stress sensitivity like dimensionless radius reservoir and its effect will not extend to the late-flow period.

Dimensionless Reservoir Thickness
Figure 10 demonstrates the effect of dimensionless reservoir thickness h D on transient pressure dynamics of a horizontal well in the naturally fractured reservoir.Obviously, it has a great influence on the early-flow period, including the following regimes: first transition flow, radial flow in the natural fracture system, inter-porosity flow and radial flow in the compound system, but has no influence on the middle-and late-flow periods.There is an implicit relationship between dimensionless wellbore radius, wellbore length and reservoir thickness, r wD L D h D = 1.In this study, let h D equal to 200, 250 and 400, and the corresponding L D is 300, 400 and 500, respectively.As shown in Figure 10, the concave-shaped segment gradually uplifts with the increase of dimensionless reservoir thickness and the value of the pressure derivative in two radial flow regimes is 1/80, 1/64 and 1/40, respectively.However, the change of dimensionless reservoir thickness does not enhance or weaken the effect of stress sensitivity like dimensionless radius reservoir and its effect will not extend to the late-flow period.Figure 11 shows the effect of dimensionless horizontal wellbore length L D on transient pressure behavior.As presented in Figure 11, the influence of dimensionless wellbore length is similar to that of dimensionless reservoir thickness on the transient pressure behavior and this parameter also mainly affects the early-flow period.Their differences lie in that the concave-shaped segment gradually descents with the increase of dimensionless wellbore length.Figure 12 shows the effect of dimensionless wellbore vertical position z wD on transient response.This parameter only influences the duration and occurrence of radial flow in the compound system.Similarly, the change of wellbore vertical position will not enhance or weaken the effect of stress sensitivity and its impact does not extend to the late-flow period.
Energies 2018, 11, x FOR PEER REVIEW 22 of 29 6.5.5.Dimensionless Wellbore Vertical Position Figure 12 shows the effect of dimensionless wellbore vertical position wD z on transient response.This parameter only influences the duration and occurrence of radial flow in the compound system.Similarly, the change of wellbore vertical position will not enhance or weaken the effect of stress sensitivity and its impact does not extend to the late-flow period.

Analysis of the Influence of Geomechanical Parameters in North Truva Field
Here we take KT-I formation of North Truva field in Kazakhstan, a typical fractured formation, as an example to understand the influence of geomechanical parameters on the elasticity parameter and transient response.Note that the elasticity parameter ε is dimensional and its unit is Pa −1 .

Influence of Geomechanical Parameters on the Elasticity Parameter
Considering the data of KT-I formation presented in Table 4, we focused on two important geomechanical parameters, Young's modulus E and Poisson's ratio ν , to investigate their influence on the elasticity parameter ε .According to the definition of elasticity parameter ε in Equation ( 28), we calculate the valve of elasticity parameter for different Young's modulus at f φ = 0.015, 0.035 and 0.055 respectively, and then plot the curves of E ε − (as shown in Figure 13a).It can be observed that for a given fracture

Analysis of the Influence of Geomechanical Parameters in Truva Field
Here we take KT-I formation of North Truva field in Kazakhstan, a typical fractured formation, as an example to understand the influence of geomechanical parameters on the elasticity parameter and transient response.Note that the elasticity parameter ε is dimensional and its unit is Pa −1 .

Influence of Geomechanical Parameters on the Elasticity Parameter
Considering the data of KT-I formation presented in Table 4, we focused on two important geomechanical parameters, Young's modulus E and Poisson's ratio ν, to investigate their influence on the elasticity parameter ε.According to the definition of elasticity parameter ε in Equation ( 28), we calculate the valve of elasticity parameter for different Young's modulus at φ f = 0.015, 0.035 and 0.055 respectively, and then plot the curves of ε − E (as shown in Figure 13a).It can be observed that for a given fracture porosity, as the Young's modulus of reservoir rock increases, the elasticity parameter first decreases sharply, and then decreases slowly, and finally decreases very slowly and even this slight decrease can be neglected.Similarly, we plot the curves of ε − ν at φ f = 0.015, 0.035 and 0.055 (as shown in Figure 13b).We can see that the elasticity parameter linearly decreases with the increase of Poisson's ratio at a given fracture porosity.Combining Figure 13a with Figure 13b, it can be concluded that under the same values of Young's modulus and Poisson's ratio, the elasticity parameter will decrease with the increase of fracture porosity.Using the properties of KT-I formation, its elasticity parameter ε ranges from 8 × 10 −8 Pa −1 to 1.1 × 10 −7 Pa −1 .
Energies 2018, 11, x FOR PEER REVIEW 23 of 29 porosity, as the Young's modulus of reservoir rock increases, the elasticity parameter first decreases sharply, and then decreases slowly, and finally decreases very slowly and even this slight decrease can be neglected.Similarly, we plot the curves of ε ν − at f φ = 0.015, 0.035 and 0.055 (as shown in Figure 13b).We can see that the elasticity parameter linearly decreases with the increase of Poisson's ratio at a given fracture porosity.Combining Figure 13a with Figure 13b, it can be concluded that under the same values of Young's modulus and Poisson's ratio, the elasticity parameter will decrease with the increase of fracture porosity.Using the properties of KT-I formation, its elasticity parameter ε ranges from    As demonstrated in Figure 14, for H519 in North Truva field, when Young's modulus increases from

Conclusions
This paper presents a new model of a horizontal well incorporating geomechanics in the naturally fractured reservoir.The analytical pressure solutions for the uniform-flux and infinite-conductivity horizontal wells including Young's modulus and Poisson's ratio are derived by using Pedrosa transform, perturbation technique and spatial integral method.Based on this study, several important conclusions are obtained: (1) For a uniform-flux horizontal well, the effect of dimensionless elasticity parameter on the pressure drop profiles becomes stronger as the production time gets longer, due to a larger fracture permeability reduction caused by the closure of natural fractures.When the effect of dimensionless elasticity parameter is strong enough, the profile looks like a bow, and a farther position away from the center of the horizontal wellbore corresponds to a lower pressure drop.(2) For an infinite-conductivity horizontal well, while the production going on, especially in the late-flow period, the pressure drop of the model incorporating geomechanics becomes larger than that not incorporating geomechanics, and more reservoir energy will be consumed to produce the same amount of oil.(3) Nine typical flow regimes can be observed on the transient response of the infinite-conductivity horizontal well incorporating geomechanics.The differences between the models incorporating and not incorporating this property lie in the final flow regimes.Due to the stress sensitivity, the pressure derivative curve is not a horizontal line with a value of "0.5" but exhibits an upward tendency in the pseudo-radial flow in the compound system, and in the pseudo-steady flow regime, both the transient pressure curve and the pressure derivative curve rise quickly and their slope is greater than "1".(4) Analysis of parameters' influence on the transient pressure behavior shows that dimensionless elasticity parameter and dimensionless reservoir radius mainly affect the late-flow period, but dimensionless reservoir thickness, dimensionless wellbore length and dimensionless wellbore vertical position mainly influence the early-flow period.There is a larger pressure drop with the increase of dimensionless elasticity parameter or the decrease of dimensionless reservoir radius in the late-flow period.

Figure 1 .
Figure 1.Horizontal well model in the closed cylindrical system.

Figure 2 .
Figure 2. Fracture porosity definition in a naturally fractured reservoir.

Figure 1 .
Figure 1.Horizontal well model in the closed cylindrical system.

Figure 1 .
Figure 1.Horizontal well model in the closed cylindrical system.

Figure 2 .
Figure 2. Fracture porosity definition in a naturally fractured reservoir.

Figure 2 .
Figure 2. Fracture porosity definition in a naturally fractured reservoir.

Figure 3 .
Figure 3. Verification of a horizontal well in a naturally fractured reservoir without stress sensitivity (ε D = 0).
4c are symmetric around the vertical line, 0 D x = , and form a bow shape, and a farther position away from the center of the wellbore corresponds to a lower pressure drop.Figure4dpresents the pressure drop profiles for different D ε at

Figure 4 .
Figure 4. Effect of dimensionless elasticity parameter ( D ε ) on horizontal wellbore pressure-drop profiles at different time.(a)
transient pressure almost overlap together, and the pressure derivative curves almost overlap together, too.As the production time gets longer, especially in the late-flow period, the pressure drop of the model incorporating geomechanics becomes larger than that not incorporating geomechanics, so does the pressure derivative.It can be interpreted by that for the model incorporating geomechanics, the closure of natural fractures gets more and more serious, and more reservoir energy will be consumed to produce the same amount of oil.

Figure 6 .
Figure 6.Flow regime identification of a horizontal well incorporating geomechanics in a naturally fractured reservoir ( 0.05 D ε = ).

Figure 6 .
Figure 6.Flow regime identification of a horizontal well incorporating geomechanics in a naturally fractured reservoir (ε D = 0.05).

Figure 8 .
Figure 8.Effect of dimensionless elasticity parameter ( D ε ) on horizontal well's transient pressure behavior.

Figure 8 .
Figure 8.Effect of dimensionless elasticity parameter ( D ε ) on horizontal well's transient pressure behavior.

Figure 8 .
Figure 8.Effect of dimensionless elasticity parameter (ε D ) on horizontal well's transient pressure behavior.

Figure 9 .
Figure 9.Effect of dimensionless reservoir radius ( eD r ) on a horizontal well's transient pressure behavior.

Figure 9 .
Figure 9.Effect of dimensionless reservoir radius (r eD ) on a horizontal well's transient pressure behavior.

Figure 10 .
Figure 10.Effect of dimensionless reservoir thickness ( D h ) on horizontal well's transient pressure

Figure 11 .
Figure 11.Effect of dimensionless wellbore length ( D L ) on horizontal well's transient pressure behavior.

Figure 11 .
Figure 11.Effect of dimensionless wellbore length ( D L ) on horizontal well's transient pressure behavior.

Figure 11 .
Figure 11.Effect of dimensionless wellbore length (L D ) on horizontal well's transient pressure behavior.

Figure 12 .
Figure 12.Effect of dimensionless wellbore vertical position ( wD z ) on horizontal well's transient pressure behavior.

Figure 12 .
Figure 12.Effect of dimensionless wellbore vertical position (z wD ) on horizontal well's transient pressure behavior.

Figure 13 .
Figure 13.Effect of geomechanical parameters on elasticity parameter ( ε ) in North Truva field.(a) Young's modulus ( E ) ( v = ); (b) Poisson's ratio ( ν ) ( Pa −1 ).6.6.2.Influence of Geomechanical Parameters on the Transient Response of H519We select a horizontal well in KT-I formation, H519, to observe the influence of Young's modulus and Poisson's ratio on this horizontal well's transient pressure behavior.The reservoir and

− 1 ,
the early-and middle-time transient responses of this horizontal well are almost the same, but a slight difference between them occurs in the late-time flow period.Through the enlarged image of the late-time curves, we can see that both the transient pressure and pressure derivative in the pseudo-steady flow regime decrease with the increase of Young's modulus.Similarly, from Figure 15, the early-and middle-time transient response is also not affected by the change of Poisson's ratio, and when Poisson's ratio varies from 0.15 to 0.45, the transient pressure and pressure derivative in the late-time flow period gradually decrease.Overall, the transient response of H519 has a slight negative correlation with Young's modulus and Poisson's ratio in the late-time flow period, especially in the pseudo-steady flow regime.

Figure 14 .
Figure 14.Effect of Young's modulus ( E ) on the transient pressure behavior of H519 in North Truva field.

Figure 15 .
Figure 15.Effect of Poisson's ratio ( ν ) on the transient pressure behavior of H519 in North Truva field.

Figure 14 .
Figure 14.Effect of Young's modulus (E) on the transient pressure behavior of H519 in North Truva field.

Figure 14 .
Figure 14.Effect of Young's modulus ( E ) on the transient pressure behavior of H519 in North Truva field.

Figure 15 .
Figure 15.Effect of Poisson's ratio ( ν ) on the transient pressure behavior of H519 in North Truva field.

Figure 15 .
Figure 15.Effect of Poisson's ratio (ν) on the transient pressure behavior of H519 in North Truva field.As demonstrated in Figure 14, for H519 in North Truva field, when Young's modulus increases from 1.0 × 10 −10 Pa −1 to 4.0 × 10 −10 Pa −1 , and then to 7.0 × 10 −10 Pa −1 , the early-and middle-time transient responses of this horizontal well are almost the same, but a slight difference between them occurs in the late-time flow period.Through the enlarged image of the late-time curves, we can see that both the transient pressure and pressure derivative in the pseudo-steady flow regime decrease with the increase of Young's modulus.Similarly, from Figure 15, the early-and middle-time transient response is also not affected by the change of Poisson's ratio, and when Poisson's ratio varies from 0.15 to 0.45, the transient pressure and pressure derivative in the late-time flow period gradually decrease.Overall, the transient response of H519 has a slight negative correlation with Young's modulus and Poisson's ratio in the late-time flow period, especially in the pseudo-steady flow regime.

( 5 )
For the KT-I formation of the North Truva field, the elasticity parameter decreases with the increase of Young's modulus or Poisson's ratio at a given fracture porosity and this parameter ranges from 8 × 10 −8 Pa −1 to 1.1 × 10 −7 Pa −1 in this formation.Calculation results show that the transient response of H519 has a slight negative correlation with Young's modulus and Poisson's ratio in the late-flow period, especially in the pseudo-steady flow regime.
fracture system, fraction φ m porosity of matrix system, fraction φ f 0 porosity of natural fracture system under the initial condition, fraction φ m0 porosity of matrix system under the initial condition, fraction c 0 oil compressibility, Pa −1 c f fracture compressibility, Pa −1 c m matrix compressibility, Pa

Table 1 .
Definitions of dimensionless variables.

Table 2 .
Basic data used for verification of the model.

Table 3 .
Basic data used in

Table 4 .
Basic data of KT-I formation in North Truva field.

Table 4 .
Basic data of KT-I formation in North Truva field.
inter-porosity flow from matrix to natural fracture, kg/s q point (t) production rate of the point sink, kg/s C wellbore storage coefficient, m 3 /Pa Dimensionless variables: real domain y wD dimensionless distance of a point sink in the y axis, m z wD dimensionless distance of a point sink or a horizontal well in the z axis, m dimensionless pressure in natural fracture system in Laplace domain p mD dimensionless pressure in matrix system in Laplace domain ψ f D0 i-order dimensionless pressure in natural fracture system in Laplace domain Special functions K 0 (x) modified Bessel function (2nd kind, zero order) K 1 (x) modified Bessel function (2nd kind, first order) I 0 (x) modified Bessel function (1st kind, zero order) I 1 (x) modified Bessel function (1st kind, first order) *