The three-zone composite productivity model for a multi-fractured horizontal shale gas well

Due to the nano-micro pore structures and the massive multi-stage multi-cluster hydraulic fracturing in shale gas reservoirs, the multi-scale seepage flows are much more complicated than in most other conventional reservoirs, and are crucial for the economic development of shale gas. In this study, a new multi-scale non-linear flow model was established and simplified, based on different diffusion and slip correction coefficients. Due to the fact that different flow laws existed between the fracture network and matrix zone, a three-zone composite model was proposed. Then, according to the conformal transformation combined with the law of equivalent percolation resistance, the productivity equation of a horizontal fractured well, with consideration given to diffusion, slip, desorption, and absorption, was built. Also, an analytic solution was derived, and the interference of the multi-cluster fractures was analyzed. The results indicated that the diffusion of the shale gas was mainly in the transition and Fick diffusion regions. The matrix permeability was found to be influenced by slippage and diffusion, which was determined by the pore pressure and diameter according to the Knudsen number. It was determined that, with the increased half-lengths of the fracture clusters, flow conductivity of the fractures, and permeability of the fracture network, the productivity of the fractured well also increased. Meanwhile, with the increased number of fractures, the distance between the fractures decreased, and the productivity slowly increased due to the mutual interfere of the fractures.

changed with the Knudsen number were not considered in all of the above-mentioned flow models. Mi et al. [10] (2014) presented a multi-scale flow model based on the correction of the different diffusion and slip coefficients. However, it was more complex to calculate. Therefore, there remained an urgent need to develop a new multi-scale non-linear flow model for shale gas reservoirs which could be applied to a productivity model. The permeability of a tight shale reservoir is extremely low. A massive multi-stage multi-cluster hydraulic fracturing stimulation technique is the major means by which to improve gas production in shale gas reservoirs, and has the ability to form multi-fractures and complex fracture networks. The half-length of the hydraulic fracture clusters, flow conductivity of the fractures, and number of fractures, all have great effects on the production [11]. At the present time, multi-porosity medium models are widely used for the productivities of horizontal wells [12][13]. Rasheed and Robert [14] (2010) developed a linear dual porosity model to approximate the shale gas completions with a horizontal well and multiple hydraulic fractures, and five flow regimes are identified from their study. Ozkan and Brown [15] (2011) presented an analytical tri-linear flow solution to simulate the pressure transient and production behavior of fractured horizontal wells in unconventional reservoirs. A new flow theory should be built to predict the potential productions of multi-stage hydraulic fracturing wells in shale gas reservoirs.
In this study, a multi-scale non-linear flow model of a shale gas reservoir was adopted, with consideration given to the slippage effects under the different diffusion coefficients, which was combined with the different weights according to the Knudsen number. Then, based on a three-zone composite model, the productivity equation of a horizontal fractured well, with consideration given to the diffusion, slip, desorption, and absorption, was built. An analytic solution was derived, and the interference of the multi-cluster fractures was analyzed, in order to provide a theoretical basis for the production optimization.

Simplified non-linear model in nano-micro pores
The Beskok-Karniadakis model was applied to the continuous, slip, transition, and free molecular flows, and was introduced as follows: The Knudsen number ( n K ), which is used to distinguish the flow regimes, was defined by Knudsen (1934) as follows: The rarefication coefficient α , which was the experimental correction of the bulk viscosity μ, was as follows: Where b is the slip coefficient, b = 0 when the slip boundary was in the first-order slip conditions, and b = -1 at second-order slip conditions; λ is the mean free-path of a molecule covered between the molecular collisions, m; and r is pore radius, m. Kn < 0.1 in the continuous and slip flow regime, where the second and higher order terms can be ignored. Therefore, the Beskok-Karniadakis model was divided from the Taylor expansion, and used to perform the first two items.
Meanwhile, for Kn ≥ 0.1, the polynomial correction coefficient ξ was introduced to modify the Klinkenberg equation, which could then guarantee a higher accuracy of the simplified binomial in calculation as follows: The approximate linear function was respectively piecewise fitted using a least square method as follows [14][15] The polynomial correction coefficient ξ was piecewise fitted using a least square method as shown in Table 1, and then plotted using MATLAB as seen in Figure 1.
Figure1 shows the comparison between the modified and Beskok-Karniadakis models. The curves of the simplified model were fitted using a least square method, and were much smoother. Also, the fitting errors were found to be much fewer in number with the higher precision.

Multi-scale motion with consideration given to the effects of the diffusion and slippage
2.2.1. Diffusion types. At the present time, most scholars believe that only the Knudsen diffusion exists in shale gas reservoirs. However, as shown in Figure 2, with the increased pore radius, the Knudsen number decreases, which leads to more molecular collisions with the pore walls, and this type of diffusion transforms the Knudsen diffusion to a Fick diffusion. The diffusion mechanisms have been suggested to be Fick diffusions in the macrospores, and Knudsen diffusions in microspores. When the pore radius increased to 50 nm, the characteristics of a Fick diffusion were observed.
As shown in Figure 3, the Fick diffusion coefficient was irrelevant to the pores' radius, while the Knudsen diffusion coefficient increased with the increasing radius. Therefore, the actual diffusion coefficient must be calculated in correspondence to the range of the Knudsen number.

Slippage effect.
The results of this study indicated that the shale gas slippage effect was combined with the gas molecule slip on the capillary walls, and the molecular diffusion of the gas inside the capillary tubes. It was also a result of a combination of the concentrations and pressure fields. The flow driven by pressure was simulated using Darcy's law, and driven by the concentrations associated with the slippage, was also simulated by the diffusion law [16][17][18][19][20].
(1) Flow regime as a continuum slippage flow when Kn＜ 0.1 A Fick diffusion (molecule diffusion) is caused by a gas molecule collision when the gas phase molecular free path is smaller than the pore size of a porous medium. The nano-micro pores, which have complex geometric structures of gas molecular diffusion channels, were crisscrossed due to twisting. Therefore, the diffusion coefficient was amended to introduce porosity and tortuosity, with consideration given to the openness and curvature degree of the pore. The Fick diffusion coefficient [19 The gas phase molecular mean free path ： The flow velocity ： Where D F is the Fick diffusivity, m 2 /s; K 0 is the permeability of the matrix, mD; K B is the Boltzmann constant, K B =1.38× 10 -23 J/K; T represents the reservoir temperature, K; µ is the gas viscosity, mPa·s; δ is the gas molecular collision diameter, m; p represents the reservoir pressure, Pa; r g is the radius of the gas molecules, m; φ represents the porosity;and τ is the tortuosity.
(2)Flow regime as a Knudsen flow when Kn＞ 10 A Knudsen diffusion is caused by a molecule collision with the pore walls when the gas phase molecular mean free path is equal to or greater than the pore size of a porous medium. Porosity and tortuosity were also introduced to modify the diffusion coefficient as follows: The Knudsen diffusion coefficient [20]: The gas phase molecular mean free path ： = 2 w zRT M p π µ λ (10) The flow velocity ： 0 0 3 16 Where D K is the Knudsen diffusivity, m 2 /s; Z is the gas compressibility factor, Z = 0.89, f; R is the universal gas constant, J• (mol•K) -1 ; and M w represents the molecular weight, f. (3)Flow regime as a transition flow when 0.1＜Kn＜10 A transition diffusion is caused by a collision between a molecule and the pore wall when the gas phase molecular mean free path is at the same scale as the pore size of a porous medium, which is a combination of the Fick and Knudsen diffusions.
The flow velocity ： The gas flow in a shale gas reservoir consists of a Darcy flow driven by the pressure, and the slippage flow influenced by the diffusion. For the different scales and different flow regimes, and in accordance to the range of the Knudsen number, Eq. (8), Eq. (11), and Eq. (13) were established by the corresponding diffusion coefficient. Then, combined with different weights ζ β , a multi-scale nonlinear flow model of a shale gas reservoir was adopted, with consideration given to the slippage effect under the different diffusion coefficient as follows: Song Fuquan et al. [21] measured the different diffusion coefficients under different pressures and pore scales. In this study, based on the flow velocity shown above, Figure 4 was constructed according to the parameters obtained in the experiment. The flow driven by the Knudsen diffusion was found to be much higher than the others. Therefore, it was necessary to simulate the different diffusion laws according to the different Knudsen numbers. Furthermore, the multi-scale flow model established in this study provided a convenient process to further the efficient development of theories for shale gas reservoirs.

Productivity equation of a horizontal fractured well
In this study, a composite fracture network system was proposed due to the permeability differences between the fracture, fracture network, and matrix zones. A three zones coupled composite elliptical model was constructed to describe the fractured reservoir, as seen in Figure 5.
In Figure 5: Then, based on the conformal transformation combined with the law of equivalent percolation resistance, the productivity equation of a horizontal fractured well's interference by multi-cluster fractures, with consideration given to the diffusion, slip, desorption, and absorption, was built as shown below:

Zone I (clusters of the major fractures) flow equation
The zone of the clusters of the major fractures was simplified to a rectangular flow zone. When the gas flowed from the boundary of the fractured rectangular flow zone to the wellbore, due to the much higher half-length of the major fractures when compared to the radius of the horizontal wellbore, it was observed as a radial flow near the wellbore area, where the flow radius was h/2, and the boundary pressure was p 1 . Therefore, the flow in the zone of the clusters of the major fractures could be divided into two parts: a linear flow along the clusters of major fractures, and a radial flow near the well bore.
The productivity equation of the linear flow was: The productivity equation of the radial flow was: The percolation resistance in the clusters of the major fractures was expressed as ： ( ) 1 sc sc sc sc Where p sc is the standard pressure, MPa; T sc represents the standard temperature, K; Z sc is the standard gas compressibility factor, f; H represents the formation thickness, m; T is the formation temperature, K; r w is the radius of the wellbore, m; p w represents the pressure at the well's bottom, MPa; w f is the width of the clusters of the major fractures, m; x f is the half-length of the clusters of the major fractures, m; K f is the permeability of the cluster of the major fractures, mD.

Zone II (fracture network) ellipse flow equation
The percolation resistance in the fracture network was expressed as ： The volume flow rate of the gas desorption was expressed as： Where K N represents the permeability of the fracture network, mD; r N is the radius of the fracture network, m; p L is the Langmuir pressure, MPa; and ρ c represents the density of the rock in the shale gas reservoir, kg/m 3 .

Zone III (shale matrix) ellipse flow equation
The productivity equation of the matrix, with consideration given to the different gas diffusion and slippage effects, was as follows： ( )

Productivity equation of the multi-scale fractured horizontal wells
The volume flow rates of the three zones connected in series, according to an equivalent flow resistance method. In this study, it was assumed that the productivity and pressures in the interface were equal, or q sc11 + q sc12 = q sc2 + q d2 = q sc3 + q d3 . The productivity formula of the multi-scale fractured horizontal wells, with consideration given to diffusion, slip, desorption, and absorption, was solved as follows:

Productivity equation of a multi-stage multi-cluster fractured horizontal well
The interference between the adjacent two elliptic flow discharged areas was found to be equivalent to the reduction of percolation resistance in accordance with an equivalent percolation resistance method. It was determined that only the interfered fractures influenced the productivity as follows: If both two cases existed, then the productivity equation of the multi-stage multi-cluster fractured horizontal well was a combination of the above two equations.

Results and discussion
In this study, based on the productivity of the fractured horizontal well (Eqs. 26 to 27), the marine shale gas reservoir parameters of the Lower Silurian Longmaxi Formation in southern China is displayed in Table 2. A MATLAB programming calculation was used to analyse the influence of each of the parameters on the productivity of the fractured horizontal well, with consideration given to the diffusion, slip, desorption, and absorption factors, such as the half-length of the fracture cluster, flow conductivity of the fractures, and permeability of the fracture network. The interference of the multicluster fractures caused by the increasing fracture clusters was then analysed. Also, the contributions made by the free and desorbed gases to the production were also clarified. The influence of the different half-length of the fracture clusters on the production of the fractured horizontal well can be seen in Figure 6. Figure7 shows the influences of the different flow conductivities of the fracture clusters on the production of the fractured horizontal well. The gas production was found to be increased with the increased half-lengths, flow conductivity of the fractures clusters, and permeability of the fracture network. Meanwhile, the increments of production gradually decreased as a result of the lowered percolation resistance. Figure 8 shows that the different cluster distances of the fractures influenced the production of the fractured horizontal well. The gas production was found to be increased with the increased number of fracture clusters. For instance, by assuming that the length of the horizontal wellbore was 900 m, when the number of fracture clusters increased to eight, the fracture interval decreased, and the fractures began to reciprocally interfere. This resulted in the productivity of a single cluster fracture becoming decreased, and the increasing rate of the production slow down. Figure 9 shows the influence of the different number of fracture clusters on the production of the fractured horizontal well. The existence of the fractures reduced the percolation resistance, and expanded the fluid flow region. When the distance between two fracture clusters was 100 m, there was no interference of flow discharge observed between the two clusters. However, the interference became more serious with the decreased intervals between the clusters, and as a result, the production of the gas well became lower.

Conclusions
(1) The different diffusion and slip correction coefficients had a significant effect on the multi-scale flow mechanism. The diffusion of the shale gas was found to be mainly in the transition diffusion and Fick diffusion regions. The flow capability was determined by the porous pressure and radius, according to the Knudsen number. When the pore radius increased to 50 nm, it was determined that the type of diffusion transformed from a Knudsen diffusion to a Fick diffusion. Therefore, a multiscale non-linear flow model should be incorporated into shale gas reservoir simulators. (2) In this study, according to the nonlinear flow characteristics in a multi-scale fractured reservoir, which were not coincident with Darcy's law, a multi-scale flow non-linear model was introduced. A three zone elliptical flow model was used in the calculation of the productivity of a horizontal fractured well, with consideration given to the diffusion, slip, desorption, and absorption, and was found to be much more applicable for shale gas recovery processes.
(3) The results of this study indicated that, with the increased half-lengths of the fracture clusters, flow conductivity of the fractures, and permeability of the fracture network, the productivity of the fractured well increased. The aforementioned were found to play significant roles during the shale gas production. In this study, it was assumed that the length of the horizontal wellbore was 900 m, and when the number of the fractures increased to eight, the distance between the fracture clusters was found to decrease to100, which resulted in the productivity being slowly increased due to the interference of the mutual fractures.