A semi-analytical two-phase flow model of fractured horizontal well with complex fracture networks in natural fractured reservoirs

Complex fracture networks are easily developed along the horizontal wellbore during hydraulic fracturing. The water phase increases the seepage resistance of oil in natural fractured reservoir. The flow regimes become more intricate due to the complex fractures and the occurrence of two-phase flow. Therefore, a semi-analytical two-phase flow model is developed based on the assumption of orthogonal fracture networks to describe the complicate flow regimes. The natural micro-fractures are treated as a dual-porosity system and the hydraulic fracture with complex fracture networks are characterized explicitly by discretizing the fracture networks into multiple fracture segments. The model is solved according to Laplace transformation and Duhamel superposition principle. Results show that seven possible flow regimes are described according to the typical curves. The major difference between the vertical fractures and the fracture networks along the horizontal wellbore is the fluid “feed flow” behavior from the secondary fracture to the main fracture. A natural fracture pseudo-radial flow stage is added in the proposed model comparing with the conventional dual-porosity model. The water content has a major effect on the fluid total mobility and flow capacity in dual-porosity system and complex fracture networks. With the increase of the main fracture number, the interference of the fractures increases and the linear flow characteristics in the fracture become more obvious. The secondary fracture number has major influence on the fluid feed capacity from the secondary fracture to the main fracture. The elastic storativity ratio mainly influences the fracture flow period and inter-porosity flow period in the dual-porosity system. The inter-porosity flow coefficient corresponds to the inter-porosity flow period of the pressure curves. This work is significantly important for the hydraulic fracture characterization and performance prediction of the fractured horizontal well with complex fracture networks in natural fractured reservoirs.


Introduction
Hydraulic fracturing is an important stimulation technique to increase the productivity of horizontal wells. Due to the large quantity of micro-fractures developed in natural fractured tight oil reservoir, complex fracture systems are easily generated along the horizontal wellbore after hydraulic fracturing treatments. The intricate arrangement of fractures in fracture networks leads to great changes of formation seepage in tight oil reservoir. Many researchers at home and abroad developed a large number of mathematical models to predict the pressure transient behavior and evaluating the fracture network characteristics of horizontal wells with complex fracture networks. Some researchers (Brohi et al., 2011;Brown et al., 2009;Lee and Brockenbrough, 1986;Olarewaju and Lee, 1989;Ozkan et al., 2011;Wooden and Coble, 1990) used a tri-linear flow model to analyze the pressure transient behavior of the fractured horizontal well with finite conductivity fractures. Horne and Terneng (1995) et al. considered the pressure interference in the fractures and proposed the well test model of fractured horizontal well with Duhamel superposition principle. Later, many researchers (Chen et al., 2015;Imam, 2013;Wan and Aziz, 1999;Zarzer and Bettam, 2003;Zeng et al., 2016) also used the analytical/semi-analytical method to obtain the type curves of fractured horizontal well in consideration of the pressure interference in the fractures.
The hydraulic fractures in aforementioned models are usually bi-wing fractures orthogonal to horizontal wellbore. However, micro-seismic monitoring and core tests (Cipolla et al., 2011;Mayerhofer et al., 2008;Weng, 2014) indicate there are many complex fracture networks created by hydraulic fracturing treatments. It will result in errors using these models to describe the complex fracture networks. Due to the intricate flow behavior in complex fracture networks, numerical methods (Cipolla and Wallace, 2014;Harikesavanallur et al., 2010;Moinfar et al., 2013;Rasdi and Chu, 2012) are used to simulate the production performance of the fractured horizontal well. Jones et al. (2013) analyzed the pressure transient behavior of horizontal well with fracture complexity by simulating the fracture with refined high-permeability grid cells. Mirzaei and Cipolla (2012) proposed a new method to automatically generate unstructured grids making simulation of fracture complexity more accurate. However, a large quantity of grid points will cause numerical dispersion in iterative calculation and make the computation much more time consuming. Therefore, semi-analytical methods (Chen and Raghavan, 1997;Dehghanpour, 2014a, 2014b;Guo and Evans, 1993) are introduced to study the transient pressure in complex fracture networks. Lin et al. (2012) and Hwang et al. (2013) predicted the performance of the fractured horizontal well with fracture networks by semi-analytical method. However, the finite fracture conductivity of the horizontal well is not taken into consideration in the model, leading to unreliable results. Xu et al. (2015) proposed a semi-analytical model to estimate the performance of the fractured well with fracture networks. The secondary fractures are described by dualporosity media, which is not in line with the actual physical model. Jia et al. (2016) and Chen et al. (2017) developed a semi-analytical model by taking the fracture networks as orthogonal main fractures and secondary fractures. The pressure transient behavior of the fractured horizontal well is analyzed and the pressure interference in fractures is considered in the model. These models assumed single-phase flow conditions for fractured horizontal well with secondary fractures, which is inaccurate to describe the oil/water two phase flow conditions for the fractured horizontal well.
Dealing with the two-phase flow conditions, Perrine (1956) and Martin (1959) (P-M) presented and verified an empirical method by substitution the single-phase compressibility and mobility into the sum of total mobility and compressibility of the multiphase system. Matthews and Russell (1967) proposed the fluid parameter calculation method of the P-M theory. Al-khalifah et al. (1989) and Hatzignatiou and Reynolds (1996) presented the pressure square method for the multiphase flow system. It is indicated that the pressure square method is more suitable for multiphase system with gas phase, while the P-M method is appropriate for oil/water two-phase system. Ezulike et al. (2013), Adefidipe et al. (2014), andXu et al. (2016) introduced a dynamic relative permeability function, calculated by the cumulative production data of water and hydrocarbon, to describe the multiphase flow behavior. Based on the material balance theory, Yang et al. (2016) and Chen et al. (2018) proposed semi-analytical method to analyze the two-phase flow preformation of the horizontal well with complex fracture networks. However, they only considered the two-phase flow in the fracture networks and the flow in the matrix was single phase.
In this work, we developed a semi-analytical method to analyze the oil/water two-phase flow pressure behavior in natural fractured tight oil reservoirs. The proposed model considered the twophase flow both in the fracture networks and dual-porosity system in fractured reservoir, which is more accurate to describe the flow and pressure behavior for the fractured horizontal well with complex fracture networks. Moreover, the pressure interferences and the finite conductivity of fractures are also considered in the model. The natural fracture is treated as a dual-porosity system and the hydraulic fracture with complex fracture networks are characterized explicitly by discretize the fracture networks into multiple fracture segments. The combined type curves are developed to investigate the properties of water content and hydraulic fracture parameters on transient pressure behaviors in detail. This work provides guidance for fracture parameter optimization and productivity estimation of fractured horizontal well with complex fracture networks in natural fractured tight oil reservoirs.

Physical model
The reservoir is assumed as a dual-porosity and horizontal infinite reservoir with constant temperature and uniform initial pressure. A fractured horizontal well producing at a constant rate is located in the center of the reservoir. The hydraulic fractures are supposed as complex fracture networks in which the main fractures is orthogonal to the secondary fractures, as shown in Figure 1. The natural fractured reservoir is described by the Warren-Root model, which is two overlapping continuum with uniform and orthogonal fractures and matrix. The fractures provide the conductive pathway for fluid flow, whereas the matrix is treated as a source term to fractures. The mass exchange between matrix and fracture system is assumed to be pseudo-steady. The fluid is slightly compressible and obeys Darcy's law. The hydraulic fractures are finite conductivity vertical fractures and completely penetrate the formation. The fluid flow at the end of the fracture wing and the flow from reservoir into the horizontal wellbore is ignored. The fluid in the secondary fractures first flows into the main fractures, and then flows from the main fractures into the horizontal wellbore. The gravity effect of the two-dimensional system is neglected.

Mathematical model
According to the oil/water two-phase flow theory and Duhamel superposition principle, the fluid flow in the complex fracture networks was coupled with fluid flow in dual-porosity reservoir system. Then the semi-analytical two-phase flow model can be developed and solved through fracture discretization and coupled matrix solution.

Analytical solution of fluid flow in dual-porosity system
According to the dual-porosity model proposed by Warren and Root, the oil/water two-phase flow diffusion equation in natural fractured reservoir was proposed as follows: Fracture system Matrix system Combining equations (1) to (4), we obtain where M tf denotes the two-phase total mobility in fracture system and M tm is the total mobility in the matrix system. The initial condition The inner boundary condition The outer boundary condition In order to simplify calculation, a set of dimensionless variables were introduced as follows: The dimensionless form of equations (5) to (11) can be written as where f w1 is the water content for oil/water two-phase flow in dual-porosity system, The initial condition The inner boundary condition The outer boundary condition Introducing Laplace transform, equations (12-17) can be expressed as The initial condition The inner boundary condition The outer boundary condition where u denotes the Laplace variable with time, p fD is the dimensionless pressure, and q D is the dimensionless production rate in dual-porosity system of Laplace domain. Substitute equation (14) into equations (18) and (19), we can obtain Equation (24) can be rewritten as Substitute equation (25) into equation (23), we can obtain where, Equation (26) can be written in the form of zero-order virtual argument Bessel equation, The general solution form of Equation (28) can be written as By substituting equations (21) and (22) into equation (29), it is concluded The point source solution of the dimensionless bottom-hole pressure in dual-porosity system can be expressed as Discretizing the main fracture into several fracture segments, the line source solution for a main fracture segment can be expressed in integral form Discretizing the secondary fracture into several fracture segments, the line source solution for a secondary fracture segment in integral form yields Analytical solution of fluid flow in complex fracture networks According to the method, the continuity equation of oil/water two-phase flow in a fracture wing can be described by where ρ o , ρ w is the density of oil and water, kg/m 3 ; v oF , v wF is the flow velocity of oil and water in the fracture, cm/s; q oF , q wF is the flow rate of oil and water, cm 3 /s; w F is fracture width, m; ∅ F is porosity of the fracture, %.
Combining equations (34) and (35), and ignoring the fracture compressibility, the continuity equation can be rewritten as where M tF denotes the total mobility of oil and water in the fracture, The initial condition The inner boundary condition The outer boundary condition Assuming the fluid flow in the fracture is one-dimensional steady-state flow, the dimensionless form for equations (36) to (40) can be expressed as where, The initial condition The inner boundary condition The outer boundary condition The solution of equations (41) to (45) can be solved by integral method as The pressure drawdown in the main fracture wing can be expressed in the Laplace domain The pressure drawdown in the secondary fracture wing in the Laplace domain can be written as where p FD is the dimensionless pressure in fracture; q FD is the dimensionless flow rate of per-unit fracture length in fracture; q FwD is the flow rate from the fracture wing into the wellbore. The subscripts m and s denote the main fracture and secondary fracture, respectively. Figure 2 shows the fracture discretization for a fracture network unit as shown in the dashed box in Figure 1. The fracture network can be divided into multiple fracture segments and each segment can be approximately regarded as a uniform flux segment. Assuming the number of main fracture is n m , and the number of secondary fracture which is orthogonal to the main fracture is 2n s . Each main fracture is divided into Mm segments, and each secondary fracture wing is discretized into M s segments. Therefore, the total number of main fractures is n m × 2M m , and the total number of secondary fractures is 2n s × 2M s × n m . Based on Duhamel superposition principle, the pressure response for the main fracture segment can be expressed as p αD1ij,lm = (1 − f w1 ) q D,lm y wD,lm +ΔL fDm,lm /2 ywD,lm−ΔL fDm,lm/2 k 0 (x D,ij − x wD,lm ) 2 + (y D,ij − α) 2 (1 − f w1 )uf (u) dα (49) Figure 2. Fracture discretization of a fracture network unit in complex fracture networks.

Model discretization and coupled solution
The pressure response for each discrete fracture segment of the secondary fracture can be written as x wD,lm +ΔL fDs,lm /2 x wD,lm −ΔL fDs,lm /2 k 0 (x D,ij − α) 2 + (y D,ij − y wD,lm ) 2 (1 − f w1 )uf (u) dα (50) Thus, the transient pressure at each discrete fracture segment with superposition principle can be expressed as where the subscripts i, j indicate the jth fracture segment in the ith fracture wing. The subscripts l, m are the mth fracture segments in the lth fracture wing. There is Fredholm integral equation including in equations (47) and (48), therefore, numerical discretization method is used to solve the equations. By discretizing each main fracture wing into M m fracture segments, equation (47) can be rewritten as When the fluid flows at the intersection of the main fracture and the secondary fracture, the flow rate for the cross fracture segment also includes the fluid in the secondary fracture wing which is intersected orthogonally to the main fracture. The Q FD,ij . satisfies Q FDm,ij = q FDm,ij + q FDs ·L fDs L fDm for cross fracture segment q FDm,ij for independent segment By discretizing each secondary fracture wing into M s fracture segments, equation (48) can be rewritten as where, l = 1，2，… 4n s × n m ; j = 1, 2, … M s ; Δx Dl = L FDs M s , x D,lm = (i − 1 / 2)Δx Dl . Considering the continuity of the pressure and flow rate along the fracture face, the continuity condition for the fracture segment is According to the fracture flow normalization conditions, the assumption of constant production rate satisfies Combining the equations from equations (49) to (59), we can obtain a 2n m × M m + 4n s × n m × M s + 1 order matrix for reservoir system and a 2n m × M m + 4n s × n m × M s + 1 order matrix for complex fracture network system. By solving the coupled matrix, we can get the dimensionless bottom-hole pressure in time domain by using the Stehfest numerical algorithm.

Model validation
To ensure the accuracy and reliability of the proposed model, the numerical simulation software (ECLIPSE 300) is used to simulate the two-phase flow behavior in the horizontal well with complex fracture networks. The basic parameters for model validation are listed in Table 1. Figure 3 shows the pressure curves of the numerical model comparing with the results in the proposed model. From the figure, it is found that the results of our semi-analytical model match well with the numerical model.

Model features analysis
The typical curves of fractured horizontal well with complex fracture networks under two-phase flow conditions were plotted based on the proposed model. Different influencing factors on the behaviors of the pressure curves were also analyzed in this section.  Figure 4 illustrates the two-phase type curves of fractured horizontal well with complex fracture networks in fractured reservoir. According to the feature of the curves, the formation flow can be divided into seven stages. The first stage is the fracture linear flow period. The fluid flow linearly and vertically to the fracture surface. The pressure derivative curve at this stage is a straight line with a slope of 1/2. The second stage is the fluid "feed flow" period from the secondary fractures into the main fractures.

Features of typical curves
There is a dip on the pressure derivative curve, indicating the fluid feed capacity from the secondary fractures to the main fractures. The third stage is the pressure interference period between the main fractures and the secondary fractures. The fourth stage is the natural fracture linear flow period. The fluid in the natural fracture of the dual-porosity system flows linearly and parallel to the main fracture surface. The fifth stage is the natural fracture pseudo-radial flow stage. The fluids in the natural fracture of the dual-porosity system flow pseudo-radially around the complex fracture networks. The sixth stage is the inter-porosity flow stage between the fracture and the matrix in dual-porosity system. The seventh stage is the formation pseudo-radial flow stage. The pressure wave has spread to the area out of the fractured horizontal well at this stage. The fluids in the dual-porosity system flow pseudo-radially around the complex fracture networks of the fractured horizontal well.

Sensitivity analysis
1. Effect of main fracture conductivity and secondary fracture conductivity The influences of dimensionless fracture conductivity on the pressure curves of the fractured horizontal well are shown in Figure 5. As seen from Figure 5(a), main fracture conductivity has major influence on the fracture linear flow period and the feed flow period from the secondary fracture into the main fracture. As the main fracture conductivity gets greater, the linear flow capacity of the fluids from the dual-porosity system to the main fracture increases, resulting in the decrease of the flow rate from the dual-porosity system into the secondary fracture. Therefore, the dip of the fluid feed period on the pressure derivative curve gets smaller. From Figure 5(b), the secondary fracture conductivity mainly influences the fluid feed flow period between the main fractures and secondary fractures. With the increase of the secondary fracture conductivity, the earlier the fluid flows from the secondary fracture to the main fracture, and the dip on the pressure derivative curves becomes larger. 2. Effect of water content in the dual-porosity system or complex fracture networks Figure 6(a) is the influences of water content in the dual-porosity system on the pressure curves of the fractured horizontal well. By comparison of the curves, it is found that water content in dual-porosity system mainly has an effect on natural fracture linear flow and pseudo-radial flow period, inter-porosity flow period and the formation pseudo-radial flow period. As the water content in dual-porosity system increases, the total mobility of the formation  fluid gets greater, the pressure and pressure derivative curves at the formation flow stage (stage 4 to stage 7) move downward to bottom of the figure. Figure 6(b) describes the pressure behavior of the fractured horizontal well in complex fracture networks. As seen from the figure, the water content in complex fracture networks has a major influence on the fluid feed flow period from the secondary fracture into the main fractures. The larger the water content in the complex fracture system, the deeper the dip on the pressure curves become and the lower the curve is offset. 3. Effect of elastic storativity ratio and inter-porosity flow coefficient Figure 7(a) shows the influences of fracture elastic storativity ratios on the pressure curves. From the figure, the elastic storativity ratio has a major effect on the inter-porosity flow period between the fracture and matrix in the dual-porosity system. Moreover, it also has influence on the fracture flow period of the complex fracture networks. When the value of elastic storativity ratio gets greater, the inter-porosity flow period occurs later and the dip on the pressure curves of the inter-porosity flow period moves more upward to the right of the curves. The pressure curves of the fracture flow period move more downward and the natural fracture linear flow period occurs later. Figure 7(b) illustrates the pressure behavior of the fractured horizontal well with different the inter-porosity flow coefficients. From the figure, the inter-porosity flow coefficient mainly influences the inter-porosity flow period between the fracture and matrix in the dual-porosity system. As the value of inter-porosity flow coefficient becomes higher, the easier the fluid exchange between the fracture and matrix and the earlier the formation pseudo-radial flow stage occurs. Thus, the dip moves closer to the left with the increase of the inter-porosity flow coefficient. 4. Effect of the main fracture number and secondary fracture number The influences of main fracture number on the pressure curves are shown in Figure 8(a). By comparison and analysis, the main fracture number has major influence on the fracture linear flow period, the fluid feed flow period and the pressure interference period. When the main fracture number gets greater, the linear flow characteristics in the fracture become more obvious. The pressure curves moves downward from the original position. The interference of the fractures increases due to the decrease of the main fracture distance. Thus, the fluid feed flow period and the pressure interference period becomes short and the formation linear flow period occurs earlier. Figure 8(b) shows the influence of secondary fracture number on the pressure behavior of the fractured horizontal well. From the figure, the secondary fracture mainly has an effect on the fluid feed flow period of the pressure curves. As the secondary fracture number gets larger, the fluid feed capacity from the secondary fracture to the main fracture increases. Therefore, the fluid feed period occurs earlier and the dip on the pressure curves moves downward when the fracture number increases.

Conclusions
In this work, a semi-analytical two-phase flow model was developed to analyze the intricate flow behavior of the fractured horizontal well. The pressure transient behavior of the complex fracture networks and the natural fracture system was analyzed by the proposed model under two-phase flow conditions. Based on the investigation, the following conclusions are drawn:  fluid "feed flow" period, (3) pressure interference period, (4) natural fracture linear flow period, (5) natural fracture pseudo-radial flow stage, (6) inter-porosity flow stage, (7) the formation pseudo-radial flow stage. 2. The significant difference for the flow regimes between the vertical fractures and complex fracture networks is the fluid "feed flow" stage from the secondary fractures to the main fractures. The pressure behavior of the developed model adds a natural fracture pseudo-radial flow stage compared to the conventional dual-porosity model. 3. The main fracture conductivity mainly influences the fluid linear flow capacity from the dualporosity system to the main fractures. The secondary fracture conductivity has a major effect on the fluid feed flow stage of the pressure curves. 4. The variation of water content will have an impact on the fluid total mobility and flow capacity of the formation fluid. The water content in dual-porosity system mainly influences the formation flow period of the pressure curves. The water content in complex fracture networks has major influence on the fluid "feed flow" stage of the pressure curves. 5. The elastic storativity ratio has a major effect on the fracture flow period and inter-porosity flow period in the dual-porosity system. The inter-porosity flow coefficient mainly influences the inter-porosity flow period of the dual-porosity system. When fracture number gets greater, the interference of the fractures increases and the linear flow characteristics become more obvious.