A coupling flow model for fractured horizontal well and anisotropic tight gas reservoir

The development of multistage fracturing technology in horizontal wells is a great impulsion to the successful development of unconventional resources. The hydraulic fractures distribute regularly along the horizontal wellbore, forming a seepage channel for fluids in tight gas reservoir and greatly improving the productivity of horizontal wells. Based on Green function and Neumann product principle, we establish a flow model of fractured horizontal well coupled with anisotropic tight gas reservoir under both unsteady state and pseudo-steady state and propose a method to solve this model. The calculation results show that flow rate of horizontal well under the early unsteady state is larger than that under the pseudo-steady state. There is no interference among fractures in the early unsteady state, and flow rate is in direct proportion to fracture numbers. Affected by frictional and acceleration pressure drop, flow rate of the end fractures is obviously larger than other fractures in pseudo-steady state. The permeabilities in different directions have great influence on well flow rate distribution. With the increasing Kx, the interference between the fractures is reduced, and the flow distribution is more balanced. When Ky becomes larger, the interference between fractures are stronger, and the “U” shape distribution of the wellbore flow is more significant.


Introduction
In recent years, with the development of hydraulic fracturing technologies, multistage hydraulic fracturing has gradually become the key technology for tight gas reservoir development. After hydraulic fracturing, large fracture networks are generated by connection of the induced high conductivity fractures with the existing natural fractures. The hydraulic fractured tight gas reservoir not only has high preliminary production, but also is more favorable for long-term stable production (Gale et al. 2007;Zeng et al. 2013;Pang et al. 2016;Yang et al. 2017;Ren et al. 2018). The flow in formation involves the coupling of fracture, horizontal well, and reservoir. Due to the complexity of fractures in tight gas reservoir, it is difficult to accurately describe the fractured horizontal well and reservoir flow mechanism (Ozkan 2011;Xu et al. 2015;Paul et al. 2016;He et al. 2017). Medeiros et al. (2008) discussed the performance and productivity of fractured horizontal wells in heterogeneous, tight gas formations. They found that if hydraulic fracturing affected stress distribution to create or rejuvenate natural fractures around the well, the productivity of the system would be significantly increased. Based on the fluid properties, Wang et al. (2014) established the formation characteristics and development measures of Chang 8 formation tight oil reservoir in Ordos Basin through the methods of reservoir engineering and numerical simulation. Hyeonseog and Seil (2015) carried out various rate transient analysis techniques (flowing material balance, type curve matching, square root time plot, and analytical model analysis) to estimate ultimate recovery (EUR), stimulated reservoir volume (SRV) and relevant information. Egberts and peters (2015) developed a simulation tool to calculate the production performance of a radial well, which provided an efficient way to assess the production potential of a stimulated well for many different well designs. Jimenez et al. (2017) proposed a method based on the gas production analysis technique in conjunction with pressure transient analysis for diagnostics specialized of linear flow, which improves the confidence in reserve estimations, specifically in early time production. Qin et al. (2018) allowed each hydraulic fracture of multifractured horizontal well consists of two individual segments with their own properties and presented a novel approach to evaluate effective fracture properties.
To date, there is few analytical or semi-analytical models to calculate the production of fractured horizontal well in tight gas reservoir, and the coupling model of conventional gas reservoir and fractured horizontal well is not suitable for tight gas reservoir (Larsen and Hegre 1994;Lian et al. 2010;Zhao et al. 2014;Ma et al. 2015;Jia et al. 2016;Wang et al. 2017;Tang et al. 2018). Therefore, it is necessary to further study the coupled flow of tight gas reservoir and fractured horizontal well to clarify the flow characteristics near the fracture and wellbore. In this paper, by using the Green function and Newman product principle, the coupling model of the anisotropy tight gas reservoir and fractured horizontal well is established, and a solving method is also given. Taking a tight gas reservoir data as an example, the sensitivity parameters of the model are analyzed.

Assumptions
The schematic diagram of tight gas model is shown in Fig. 1, and the main assumptions are as follows: (1) There is a tight gas reservoir, with no-flow boundaries in all directions. The length, width, and height of the drainage volume are a, b, and h, respectively.
(2) The permeability of the matrix is so low that it can be ignored. The permeability of natural fractures in x, y, and z directions is expressed as K x , K y , and K z . The equivalent porosity of natural fractures and matrix is φ.
(3) A horizontal well placed in the tight gas reservoir with radius r w and length L. The well is placed from (x 0 , y 1 , z 0 ) to (x 0 , y 2 , z 0 ) and is parallel to the y-axis.
(4) The horizontal segment has a number of N f segments of hydraulic fractures. The fractures are perpendicular to the wellbore, equidistant distribution and pass through the thickness of the whole gas reservoir, the fractures have a length x f , and width y d .
(5) Initially, at time t = 0, the pressure is uniform throughout the reservoir and is equal to p 0 .

Single fracture flow model
Defining the pseudo-pressure function as follows (Li et al. 2009): By using instantaneous point sink function (Green function) and Neumann product principle (Gringarten and Ramey 1973), the pressure drops at any location ( 0 ≤ x ≤ a, 0 ≤ y ≤ b, 0 ≤ z ≤ h ) of a single fracture vertical to wellbore at any time (t > 0) can be expressed as follows: cos n x 0 a cos n x a factors in x, y, and z directions, respectively; q is the mass flow rate under pseudo-pressure; S 1 , S 2 and S 3 are instantaneous point sink functions located at (x 0 , y 0 , z 0 ). According to equations of state of the actual gas, the mass flow rate of the fractured horizontal well is transformed into volumetric flow rate under the standard condition. It can be written as: According to Eq. (1) and (5), the following equation is obtained: where q is the volumetric flow rate under the standard condition, m 3 /s; p sc is the pressure under the standard condition; T is the reservoir temperature, K; T sc is the temperature under the standard condition, K.
In Eq. (7), after integration to instantaneous point sink functions, time t is the only unknown variable, so the total pressure drop from the fracture is given by (Penmatcha and Aziz 1998):

Coupling model of fractured horizontal well and tight gas reservoir
In order to establish the fractured horizontal well and tight gas coupling model, we firstly divide the horizontal well into N f segments, in which segment 1 is assumed to lie closest to the heel of the well, while segment N f lies closest to the toe of the well. There is a fracture in the middle of each segment, and the fractures are uniformly distributed. There is no fluid inflow in the horizontal section, so the total production of horizontal well is equal to the flow rates all of fractures. As shown in Fig. 2, there are 6·N f unknown variables in the model, respectively: (1) Tight gas reservoir nodal pressure: p 1 , p 2 , ..., p N f ; (2) Right nodal pressure hydraulic fracture: Flow rate from the tight gas reservoir node: q 1 , q 2 , ..., q N f ; (6) Flow rate between hydraulic fractures: q f ,1 , q f ,2 , ..., q f ,N f .

Mass conservation equations
According to mass conservation principles, the wellbore nodal flow rate equals to the tight gas flow rate from the reservoir and fractures. For simplicity, we assume that the tight gas density is constant, and we can get N f equations:

Relationship between nodal pressures of hydraulic fracture
Because of gas flow from fracture to wellbore, the pressure at the right side of the fracture is not equal to the pressure at the left side. Therefore, we adopt the average pressure of two sides of fracture as the central nodal pressure of the hydraulic fracture (Lian et al. 2011):

Pressure drop in hydraulic fracture
Assume the i th central fracture nodal coordinate is (x 0 , y i , z 0 ) , the tight gas reservoir nodal coordinate is ( x i , y i , z i ), if the fracture is anisotropic, the following equations can be used to transform the fracture into isotropic: Because the cross-sectional area of fractures is larger than that of wellbore, the flow from the edge of the fractures to the wellbore can be considered as radial flow ignoring gravity. The area has thickness y d , permeability K 0 , well skin factor S, flow radius r f , and boundary pressure p i .
In the n th time step, the flow equation can be written as follows:

Continuity equations
The time step is assumed as Δt , and t = nΔt , the spatial and temporal superposition of pressure drop caused by all N f fracture production is calculated, and we can get the continuity equation of the i th fracture at the n th time step: where F ij is the effect of fracture j on fracture i.

Constraint equations
If the flow rate constraint is used, the equation would be: where Q max is the maximum flow rate. The flow of tight gas into the wellbore presents an elliptical flow, the major axis is a = x f , and minor axis is b = h. In order to facilitate the simulation, the elliptic flow is transformed into a circular flow with equal radius, as shown in Fig. 3. The seepage radius can be calculated: If the bottom-hole pressure constraint is used, the equation would be: where p min,wf is the minimum bottom-hole pressure. Figure 4 shows skeleton diagram of the gas flow model in fracture and horizontal well. In the horizontal well section between two fractures, there is no gas flow from reservoir to wellbore, and the frictional pressure drop mainly exists in the pipe. Because of gas flow from fracture to wellbore, there is acceleration of pressure drop.

Wellbore friction pressure drop
The frictional pressure drop of horizontal well section between two fractures could be written as follows: where M air is molecular mass of air, g/mol; γ g is the relative density of gas, dimensionless; Δl j is the length of the j th section; f j is the friction factor of the j th section; R is ideal gas constant, 8.314 J⋅mol −1 ⋅K −1 ; D is the diameter of wellbore, m.
The friction factor f is different in laminar and turbulent flow. The key parameter to distinguish turbulent flow and laminar flow is Reynolds number, which is expressed as follows: Rf ,j−1 − P 2 Lf ,j = 1.0 × 10 −6 f j M air g R 16p 2 sc q 2 f ,j TZ 2 D 5 T 2 sc Δl j j = 2, 3, ..., N f According to the classification criteria, the state of fluid is laminar flow when the Reynolds number is less than 2100, the state of fluid is turbulent flow when the Reynolds number is greater than 2100. The Reynolds number is related to the seepage velocity, and the flow velocity depends on the flow rate.
According to the study of Ouyang et al. (1996), the formulas for calculating friction factor in laminar and turbulent flow are as follows:

Acceleration pressure drop
As shown in Fig. 4, it is assumed that the flow velocity at the right and left sides of the fracture i is v 1,i and v 2,i , respectively: The acceleration pressure drop is: Substituting (23) and (24) into (25), the following equation can be obtained:  Eqs.9,10,14,15,16 or 17,18,25), the model can be solved. In order to simplify the solution process, the flow from the tight gas reservoir node and central nodal pressure of hydraulic fracture is eliminated by using Eqs. 9 and 10. After simplification, there are 4·N f unknown quantities, which are shown, respectively, as follows: (1) Tight gas reservoir nodal pressure: p 1 , p 2 , ..., p N f ; (2) Right nodal pressure hydraulic fracture: (3) Left node pressure of hydraulic fracture:  (17)) to form the following equation group: The equation group (26) has nonlinear characteristics, so the Newton iterative method is used to solve them. At the right hand, only flow rates are unknown variables, we can build Jacobi matrix by calculating the coefficients. The initial values of a set of solutions are given firstly, and then, cyclic iteration is performed. When the error between certain iteration and the next iteration is less than 10 -4 , the iteration is stopped and the results are adopted. We use VC ++ to compile the procedure, and the calculation flow diagram is shown in Fig. 5. By adjusting the input parameters and convergence criteria, the calculation process can be converged.

Calculation examples and sensitive parameter analysis
Taking the tight gas reservoir data as an example, the sensitive parameters of the model are analyzed by using the proposed model. The parameters of tight gas reservoir and fractured horizontal well are shown in Table 1. Figure 6 shows the relationship between horizontal well production rate and time with different number of hydraulic fractures. The minimum bottom-hole pressure is set as 5 MPa. Because of the low permeability of tight gas reservoir, when the fractured horizontal well begins to produce, only the reservoir near the fracture flows. At this time, there is no interference between the fractures, and the flow state is the early unsteady state. In Fig. 6, it can be seen that the production of fractured horizontal wells is proportional to the number of fractures. As time increases, the pressure wave gradually propagates outward. When the pressure waves of adjacent fractures overlap, interference occurs between the fractures, and the production of each fracture begins to decrease. The more fractures there is, the stronger is the interference, and the larger is the decline of single fracture production.

Production analysis of fractured horizontal well
When the pressure wave propagates to the outer boundary, the flow enters pseudo-steady state, and the flow rate is gradually stabilized. Because tight gas reservoirs are mostly in the pseudo-steady state during the development process, there is an optimized number of fractures. In this example, when the number of fractures is larger than 10, the production of the pseudo-steady state changes little, so it is enough to use 10 fractures. Figures 7 and 8 show the relationship between the total flow rate and the length of the fractured horizontal well at 0.01 and 10 days, respectively. It can be seen that in the unsteady state, with the increasing wellbore length, the total flow rate of horizontal well decreases. Because there is no interference between the fractures in the unsteady state, when the fracture parameters are the same, the total flow is affected by the wellbore friction. With the increase in the horizontal well length, the friction pressure drops of the wellbore increases, leading to the lower production of the toe and heel fractures, which results in the decrease of the total flow rate. This is obviously different from the pseudo-steady state. In the pseudo-steady state, with the increase of the horizontal well length, the interference between fractures decreases, and the horizontal well production increases gradually.  Flow rate and pressure distribution along the fractured horizontal well Figure 9 shows the distribution of fracture flow rate in unsteady state (0.01 days) and pseudo-steady state (10 days). The number of hydraulic fractures is 10, and the maximum flow rate of fractured horizontal wells is 100000 m 3 /d. It can be seen from Fig. 9 that in the unsteady stage, the flow rate of the heel fracture is obviously larger than that of the other position, and the closer the toe is, the smaller the flow rate of the fracture is. Because of the existence of wellbore pressure drop, the pressure drop of the heel is greater than that of the toe, resulting in uneven distribution of flow. In the pseudosteady state, the pressure in the wellbore is balanced, and the flow distribution becomes relatively balanced. Figure 10 shows the pressure drop of the heel of the horizontal well variation with time, considering the impact of the friction and acceleration. The impact of the friction on pressure drop is small, but the influence of the acceleration on pressure drop is even smaller. Only when the length of the horizontal well is very short and the flow rate is very high, the effect of acceleration would be great. Figure 10 also shows that the unsteady flow phase lasts about 7 days, and then the pseudo-steady state begins, and the pressure difference between infinite conductivity model and finite conductivity model in the pseudo-steady state is about 35 kPa. According to the research results of this study, it is concluded that the impact of friction on low permeability oil and gas reservoirs is much smaller than that of high permeability reservoirs. Figure 11 shows the pressure distribution along the wellbore of fractured horizontal wells with 10 fractures at 10 days. The pressure drops at the heel and toe are 1,044 kPa and 997 kPa, respectively. Therefore, the pressure drop in the wellbore is about 46 kPa. Under this pressure drop, the fluid in the wellbore can flow from the toe to the heel. In the wellbore, the flow rate at the heel is larger than that at the toe, and the friction and acceleration pressure drop are larger too. Figure 12 shows the influence of the changing K x , K y and K z on wellbore flow distribution. Due to the thin thickness of tight gas reservoir, and the fracture runs through the thickness of the gas reservoir, the fluid is mainly from X and Y directions, and K x and K y have great influence on well flow rate distribution. X direction is parallel to the fracture, with the increase of K x , the interference between the fractures decreases and the flow rate distribution is more balanced. Y direction is perpendicular to the fracture. When K y is small, the energy is mainly obtained from X direction, and the distribution of flow rate is relatively balanced. When K y becomes large, the interference between fractures is stronger, more gas flows into the wellbore from fractures at both ends, and the "U" shape distribution of the wellbore flow is more significant.

Conclusions
The flow mechanics of fractured horizontal well in anisotropic tight gas reservoir is studied in detail. The flow characteristics under both unsteady state and pseudo-steady state are also analyzed. The distribution of pressure and flow rate along the wellbore are compared. The main conclusions are as follows.
(1) The flow pattern of the horizontal well is unstable during the early stage of production, and the decline rate The distance from the heel of the well, m Fig.11 Pressure along the well at 10 days in pseudo-steady state of production is large. But after a short time, the decline rate will become stable, and the flow enters the pseudosteady state. In the early stage of unsteady flow, the flow rate is proportional to the number of the fractures. In the pseudo-steady state, the fractures are interfered and the flow rate of each fracture is reduced.
(2) When the number of fractures is constant, the total flow rate of horizontal well decreases slightly with the increasing wellbore length in the unsteady stage, while the total flow rate of horizontal well increases slightly with the increasing wellbore length in the pseudo-steady stage.
(3) The influence of anisotropy on tight gas reservoir is studied, with the increasing K x , the interference between the fractures decreases, and the flow distribution is more balanced. When K y becomes larger, the interference between fractures are stronger, and the "U" shape distribution of the wellbore flow is more significant.

Compliance with ethical standards
Conflict of interest The authors declared that they have no conflicts of interest in this work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.