Productivity Model for Shale Gas Reservoir with Comprehensive Consideration of Multi-mechanisms

Multi-stage fracturing horizontal well currently has been proved to be the most effective method to produce shale gas. This method can activate the natural fractures system defined as stimulated reservoir volume (SRV), the remaining region similarly is defined as un-stimulated reservoir volume (USRV). At present, no type curves have been developed for hydraulic fractured shale gas reservoirs in which the SRV zone has triple-porosity dual-depletion flow behavior and the USRV zone has double porosity flow behavior. In this paper, the SRV zone and USRV zone respectively are simplified as cubic triple-porosity and slab dual porosity media. We have established a new productivity model for multifractured horizontal well shale gas with Comprehensive consideration of desorption, diffusion, viscous flow, stress sensitivity and dual-depletion mechanism in matrix. The rate transient responses are inverted into real time space with stehfest numerical inversion algorithm. Type curves are plotted, and different flow regimes in shale gas reservoirs are identified. Effects of relevant parameters are analyzed as well. The whole flow period can be divided into 8 regimes: bilinear flow in SRV; pseudo elliptic flow; dual inter-porosity flow; transitional flow; linear flow in USRV; inter-porosity flow and boundary-dominated flow. The stress sensitivity basically has negative influence on the whole productivity period .The less the value of Langmuir volume and the lager the value of Langmuir pressure, the more lately the inter-porosity flow and boundary-dominated flow occurs. It in concluded that the USRV zone has positive influence on production and could not be ignored.


INTRODUCTION
Shale gas are typical unconventional reservoir due to its ultra-low permeability and porosity.Generally speaking, there is no natural productive capacity for those kinds of reservoirs.Multi-stage fracturing horizontal well currently has been proved to be the most effective way to produce shale gas, and this method can not only create several highconductivity hydraulic fractures, but also activate and connect the existing natural fractures so as to form large spacious network system [1].The zone containing the main high-conductivity hydraulic fractures and large spacious network system both of which can effectively improve the wells performance is defined as SRV (stimulated reservoir volume), and the remaining zone which hardly influenced by the treatment of hydraulic fracturing is similarly defined as USRV (un-stimulated reservoir volume) [2][3][4].
At present, a number of scholars have done large amount of researches about transient rate analysis for shale gas, some analytical and semi-analytical solutions are developed as well.Shale gas reservoir is the classical naturally fractured reservoir (NFR) which contains complex natural fractures and ultra-low permeability.In terms of this kinds of reservoirs, Barenblatt (1960), Warren and Root (1963) originally *Address correspondence to this author at the College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China; Tel: 18810907235; E-mail: 987558984@qq.comproposed the dual-porosity model which assumed pseudo steady state fluid transfer between matrix and fractures [5,6], and then Kazemi (1969), de Swaan (1976) and Ozkan et al. (1987) developed some other dual-porosity models for shale gas reservoirs to enrich the former productivity model, these models assumed unsteady-state (transient) flow condition between matrix and fractures [7][8][9].However, all of these dual-porosity models neglected the diffusion and adsorption phenomenon in shale gas reservoirs.Some scholars investigated amount of field production data and found that these dual-porosity models may not be true in actual reservoirs.An improvement to overcome this drawback is to considerate two different fracture systems with different properties.This system is so-called triple porosity system.Al-Ghamdi and Ershaghi (1996) initiatively proposed the dual fracture triple-porosity model for radial flow [10], and then Liu et al. (2003), Wu et al. (2004) and Dreier (2004) enriched the triple-porosity model [11][12][13], but unfortunately these model still did not considerate the impact of adsorption and diffusion.However, the linear flow stage are apparently identified in some real productivity curves, especially for these fractured reservoirs, therefore, the linear flow models for shale gas are proposed by some scholars.El-Banbi (1998) proposed a linear dual-porosity model in linear fractured reservoirs [14], and originally derived the solutions in Laplace space, but the impact of desorption, diffusion and USRV zone on the production is ignored; Hasan and Al-Ahmadi (2011) proposed a triple-porosity linear flow

Deng et al.
model with consideration of the impact of shale gas desorption and diffusion [15], however, the impact of USRV zone was neglected; Xu et al. (2012) analyzed the effect of USRV zone on shale gas production, at the same time, the impact of desorption diffusion is considered as well [16]; Zhao et al. (2013) proposed triple-porosity spherical flow model for the fractured infinite shale gas reservoirs which considered the impact of diffusion and desorption [17], however, they considered artificial fractures as infinite conductivity.
In terms of these naturally fractured reservoirs, the phenomenon of stress sensitivity is readily observed.Samaniego VF et al. (1980), Raghavan et al. (2004) employed simulation and experiment methods to make some research about the impact of stress sensitivity on the conventional reservoirs [18,19].Pedrosa et al. (1986) first applied the mathematical method to study the stress sensitivity for homogeneous and dual-porosity reservoirs [20]; Wang (2013) develop a dual porosity spherical flowing model with consideration of the stress sensitivity in micro-fractures for shale gas reservoirs [21].
Ezulike Daniel Obinna and Dehghanpour Hassan (2014) first proposed the triple-porosity dual inter-porosity linear flow model [22], that is to say, the gas simultaneously depletes from matrix into micro-fracture and macro-fracture, however, the desorption, diffusion and stress sensitivity in fracture are ignored.
In view of this, shale gas transfer in the reservoir is the result of mutual effects of various percolation mechanisms, and therefore, it is necessary to comprehensively consider the impact of various mechanisms in order to obtain important dynamic parameters for shale gas reservoirs.This paper simplifies the SRV zone and USRV zone as triple-porosity cubic model and dual-porosity slab model respectively, comprehensively taking various mechanisms into account, such as adsorption and diffusion in shale matrix, viscous flow and stress sensitivity in fractures, besides, we assume the gas simultaneously depletes from the matrix into microfractures and macro-fractures in SRV zone and this is the essence of this paper.Laplace transformation, perturbation method are employed to solve this new model.The transient rate responses are inverted into real time space with stehfest numerical inversion algorithm [23].Type curves are plotted, and different flow regimes in shale gas reservoirs are identified.The effects of relevant parameters are analyzed as well.Besides, this model also compares with numerical simulation and exhibits good agreements.

Physical Model
The schematic illustration in Fig. (1a) shows a multistage fracturing horizontal well.Multi-stage fracturing shale gas reservoir is divided into SRV zone and USRV zone, SRV zone and USRV zone is simplified as cubic tripleporosity model and dual-porosity slab model.A horizontal well located in the center of a rectangular closed formation producing at constant wellbore pressure.The other assumptions are as follows: (1) The initial pressure distribution in the reservoir is uniform which equals to P i , the SRV zone contains micro-fractures, macro-fractures and matrix, the USRV zone contains micro-fractures and matrix, the fractures in different zone have different properties.
(2) The macro-fracture is perpendicular to the horizontal well and evenly distributed along the wellbore, the micro-fractures are perpendicular to the macro-fractures as well, the length of reservoir and horizontal wellbore are equivalent, the length of micro-fracture and the width of reservoir respectively equals to y f and y e .
(3) Macro-fractures have finite conductivity and are assumed to be penetrated fully, considering stress sensitivity in macro-fractures.
(4) Flowing is sequential from one medium to another medium.In the SRV zone, only fluid flow from macrofractures to wellbore is considered; The shale gas simultaneously deplete from matrix into micro-fractures and macro-fractures with pseudo-steady state inter-porosity flow; the fluid flow between micro-fractures and macrofractures is unsteady state flowing.In the USRV zone, the fluid flowing from matrix to the fractures is pseudosteady state inter-porosity flow; the connection between SRV zone and USRV zone via the macro-fractures in SRV and micro-fractures in USRV.
(5) Slightly compressible shale gas and compressibility coefficient is constant; (6) Shale gas desorption and diffusion respectively meets the Langmuir isotherm equation and the first law of diffusion; (7) The impact of gravity and capillary pressure is ignored.
This paper considers the simultaneous depletion from matrix into micro-fractures and macro-fractures.To analyze this flow process conveniently, the matrix in SRV zone is artificially divided into two distinct segments which have different permeability and porosity ratio respectively.The schematic illustration in Fig. (1b) shows the depletion process from matrix to micro-fractures and macro-fractures in SRV zone.

Mathematical Model
Based on the mass balance principle, the governing equations respectively in micro-fractures, matrix, and macrofractures in SRV zone and USRV zone with consideration of adsorption, diffusion, viscous flow and stress sensitivity are as follows: Micro-fracture The matrix in SRV zone are divided into two segments, the segment denoted as matrix-1 which permeability and porosity equal to k m1 and m 1 depletes into macro-fractures and another segment denoted as matrix-2 which permeability and porosity equal to k m2 and m 2 depletes into microfractures.Therefore, the governing equations of these two segments are as follows: Matrix-1: Matrix-2: Matrix: Inner boundary condition Inner zone micro-fracture Outer boundary condition To simplify these equations, some dimensionless variables are defined (seen from appendix A) and substituting these dimensionless variables into equation 4-15, the dimensionless governing equations are as follows: SRV Zone Macro-fracture: Inner boundary condition Interface condition Outer boundary condition It's not difficult to find that the macro-fracture equation contains strong nonlinear, the perturbation technology and the Presoda transformation are applied to linearize this equation [20], the formula is as follows: According to the theory conducted by Wang (2013) [24] and Presoda (1986), performing a parameter perturbation in D by defining the following series: Finally, the final governing equation of SRV zone and USRV zone can be changed into the following profiles.

Model Solution
The Laplace transformation is used to solve equation 25-35 the details of this process are listed in appendix B. The final solution of these equations in Laplace space is as follows: However, to analyze the impact of relative parameters and identify the shale gas flowing regimes, the transient rate responses are inverted into real time space with stehfest numerical inversion algorithm.

NUMERICAL SIMULATION VALIDATION
Numerical models were used to validate this new analytical solution proposed in the paper.This validation process is operated in the commercial software-CMG.We assume that the rate of water is ignored and gas is the only phase existing in the reservoir.The basic parameters are listed in the Table 1.The radial coordinate is converted into x-y coordinate, and the value of the permeability of these two directions is different based on the flowing regimes.In SRV zone, the vertical permeability of some grids which are located near the wellbore is assumed to be very high so as to ensure the facts that the gas can flow into the wellbore.It is difficult to simulated the adsorption used CMG, for simplicity, the adsorbed gas is not included.We focus on the tendency of curves.
For gas case, as seen in Fig. (2), the simulation data from CMG are converted into the formation of dimensionless time and dimensionless rate for the convenience of analysis.It is concluded that the curve is in good agreement with the analytical solution at the stage which represents the flowing regime in fractures in SRV zone and stage which represents the flowing regime in fractures in USRV zone, while it deviates a lot at the stage which represents the flowing regime in matrix in SRV zone and stage which represents the flowing regime in matrix in USRV zone.The main reason for the deviation is that adsorption is neglected in simulation model.The deviation also could be caused by using the inappropriate average pressure to calculate the pseudo pressure.At the very early times, the flow is fracturedominated flowing regime in SRV zone, the adsorptive gas has not desorbed from the surface of matrix, so the predictive results are similar with those from the numerical simulation (as seen in stage and ).However, at the middle times and late times, the flow is matrix-dominated and boundary-dominated flowing regime, and the situation is different.The adsorptive gas can desorb from the surface of matrix and depletes into fractures, this process can increase the shale gas production (as seen in stage and ).

TYPE CURVES AND DISCUSSIONS
In this paper, we make some assumptions that SRV zone is cubic triple-porosity model and USRV zone is slab dualporosity model.Some characteristic parameters of these kinds of shale gas reservoirs include: F , f1 , m1 , m2 , f2 , m3 .However, if these parameters are equal to some special value, the new model proposed in this paper can be changed into some other familiar models.For example, f1 is equal to 0, the model of SRV zone can be dual porosity slab model; when the m3 equals to 0, USRV zone can be simplified as homogeneous model; m1 = 0, the gas only depletes from the matrix to micro-fractures in SRV zone (Hasan A. Al-Ahmadi, 2011), in another words, the new model can be similar with the triple-linear model proposed by El-Banbi (1988) as well; m2 = 0, SRV zone can be simplified as dualporosity model (Xu et al, 2013).In short, this new model has universal application to different formations when these characteristic parameters are equal to different values.It is worth mentioning that the stress sensitivity in the macrofractures is considered, this consideration extremely matches the real production condition in fractured reservoirs, especially for the ultra-low permeability shale gas reservoirs .The mechanism of the stress sensitivity will be analyzed in the following sections.
Combined with stehfest numerical inversion algorithm, the type curves concerning dimensionless rate derivative q D ' with respect to t D and dimensionless rate q D with respect to t D are plotted under different condition of formation properties, eight flow regimes can be easily observed by analyzing the following type curves Seen from Fig.

(3):
Regime VIII: Boundary-dominated flow regime.At this time, the boundary has influence on dynamic production of well, the rate and rate derivative curve decrease rapidly.
The paper considers the impact of stress sensitivity in the macro-fractures in SRV zone.Fig. (4) shows that stress sensitivity has great negative influence on horizontal well performance during the whole productivity period.However, with the extension of production time, the negative impact on rate becomes less and less, especially during the boundary-dominated flow stage and the negative impact basically vanishes.The permeability of macro-fractures in SRV zone consecutively decreases with the reduction of pressure, when the pressure wave has reached to the outer boundary, the permeability of the macro-fractures is so small that the influence of stress sensitivity can be ignored.Besides, with the increasing of stress sensitivity coefficient ( D = 0, D = 0.8, D = 1.5), dimensionless rate curve (Fig. 4a) and dimensionless rate derivative curve (Fig. 4b) simultaneously descend.By analyzing the mechanisms, the more the pressure reduces in macro-fractures, the more severely the macro-fractures close, as a result, the permeability of macro-fracture rapidly decreases.
Two most important characteristic parameters are Langmuir pressure and Langmuir volume which reflect the feature of adsorption.P L and V L primarily affect the stage of pseudo-steady state inter-porosity between matrix and fracture in SRV zone, linear flow in USRV zone and outer boundary-dominated flow.The larger the value of P L (P L perspective is equals to 2 MPa, 5 MPa and 10 MPa) or the  smaller the value of V L (V L perspective is equal to 5 sm 3 /m 3 , 10 sm 3 /m 3 , 15 sm 3 /m 3 ), the less the shale gas rate during the whole process of productivity period, and the later the stage of pseudo-steady state inter-porosity between matrix and fractures in SRV zone occurs, and linear flow in USRV zone and outer boundary-dominated flow occurs.
From the analysis of micro-mechanisms, the larger the value of Langmuir volume is and the more the content of adsorbed gas are.When the pressure is smaller than Langmuir pressure, the adsorbed gas starts to desorb from the surface of matrix and spread into fractures so that the decline ratio of pressure and rate slows down.At the same time, due to inter-porosity flowing between matrix and fracture, the negative of outer boundary will delate lately, in another word, the dimensionless rate derivative curve wholly moves right (Fig. 5a, d).In short, due to the gas desorption and diffusion, the decline ratio of pressure becomes slow and the rate is less affected by outer boundary.When shale gas desorption reaches a certain level, the amount of desorption is insufficient to cover the impact of outer boundary, dimensionless rate derivative curve moves downward (Fig. 5b, c).
The paper assumes that gas simultaneously depletes from matrix into micro-fractures and macro-fractures in SRV zone.The matrix is divided into two segments, where, the segment denoted as matrix-1 which permeability and porosity equal to k m1 and m1 depletes into macro-fractures and another segment denoted as matrix2 which permeability and porosity equal to k m2 and m2 depletes into micro-fractures.We introduce the pore volume ratio denoted as which equals to m1 /( m1 + m2 ), and =0, =1, 0< <1 respectively represents: (1) the matrix only depletes into micro-fractures; (2) the matrix only depletes into macro-fractures; (3) matrix simultaneously depletes from matrix into micro-fractures and macro-fractures.The type curves (Fig. 6) under different value of , k m1 and k m2 are plotted.Some conclusions can be obtained via the analysis of curves: (1) k m1 > km 2 , the larger the value of is, the higher the rate of shale gas is (Fig. 6a); (2) k m1 <= km 2 , on the one hand, during the early flowing stage, the larger the value of is, the lower the rate of shale gas is, on the other hand, during the rest of flowing period, the larger the value of is, the higher the rate of shale gas is (Fig. 6b).From the analysis of micro-mechanisms, k m1 <= k m2 < k Fi , as a result, the inter-porosity flowing capacity between matrix-2 and micro-fractures is stronger than that between matrix-1 and macro-fractures, therefore, the larger the pore volume ratio of matrix-1 which depletes into macrofractures is, the less the total amount of shale gas from reservoir during the same productivity period is.In conclusion, as for the actual shale gas reservoir, we can acidize the matrix to increase its permeability so that the shale gas stored in matrix can directly flow into the macro-fractures, this treatment can significantly improve the horizontal well performance.
This paper assumes that USRV zone is dual-porosity slab model including matrix and fracture.The dimensionless rate type curves (Fig. 7) are plotted under different conditions in USRV zone, such as different matrix permeability, different fracture permeability and different width of reservoir.Observing from all of the following four figures, the rate rapidly declines without consideration of USRV zone.Generally speaking, the width of the whole reservoir is much larger than that of SRV zone; therefore, the flow occurs in the USRV zone can last longer than that in SRV zone (Fig. 7a).Fig. (7b) represents dynamic principle of rate under different permeability of matrix, the larger the value of ma-trix permeability is, the more the amount of gas deplete from matrix into fractures in USRV zone is, the later the boundary-dominated occurs.Fig. (7c) represents the dynamic principle of rate under different permeability of fracture in USRV zone, during the early productivity period, the larger the value of fracture permeability is, the more the gas rate is, however, the sooner the pressure wave reaches to the outer boundary, what's worse, the flowing regimes of matrix and micro-fractures in the these two zones are hardly observed (seen from the blue circle in Fig. 7c).Fig. (7d) reflects the impact of the width of reservoir, the wider the reservoir is, the more apparent the flowing regimes we can observe from the type curves.Through the above analysis, we can recognize that it is not reasonable to neglect the impact of USRV zone, especially for these kinds of reservoirs which have great formation properties.

CONCLUSION
New analytical solution for shale gas reservoir with multi-stage fracturing horizontal well has been developed which considers the USRV zone as a dual porosity system and the SRV zone as a triple-porosity system.The solution is more general for type curve analysis both in homogeneous and naturally fractured reservoirs.Numerical simulation model was used to validate the analytical solutions and obtains an excellent agreement.
Analyzing the type curves, shale gas flow is divided into eight regimes: bilinear flow of macro-fractures and microfractures in SRV zone; pseudo-elliptic flow; linear flow of micro-fracture in SRV zone; dual-pseudo steady interporosity flow; transition flow between SRV and USRV zone; linear flow of fracture in outer zone; pseudo -steady state inter-porosity flow; outer boundary-dominated flow.
The impact of adsorption, stress sensitivity in shale gas reservoirs must not be ignored; otherwise, the performance of horizontal well cannot be predicted precisely.It is also concluded that the dual-porosity behavior of USRV zone has a positive effect on production, the larger the value of the permeability of matrix in USRV zone is, the more apparent the positive effect is.The stress sensitivity has negative influence on production during the whole productivity period.This paper initially introduces the triple-porosity dualdepletion model in SRV zone for shale gas.It is concluded that the USRV zone has positive influence on production.Reservoirs could be acidized so as to increase the permeabil-ity of matrix in SRV zone and optimize the performance of horizontal wells for shale gas reservoirs.
Notes: the program can provide if it is necessary.

NOMENCLATURE
Cg gas compressibility, MPa -1 C tF , C tf1 total compressibility of macro-fractures and microfractures in SRV zone, MPa -1 C tm1 ,C tm2 total compressibility of matrix in SRV zone MPa -1 C tf2 total compressibility of micro-fractures in USRV zone, MPa -1 C tm3 total compressibility of matrix in USRV zone, MPa -1 h formation thickness, m L horizontal well length, m y f macro-fracture length, m y e reservoir width, m F , f1 porosity of macro-fractures and micro-fractures in SRV zone m1 , m2 porosity of matrix in SRV zone f2 porosity of micro-fractures in USRV zone m3 porosity of matrix in USRV Zone k F , k f1 permeability of macro-fractures and micro-fractures in SRV zone, D k m1 , k m2 permeability of matrix in SRV zone, D k f2 permeability of micro-fractures in USRV zone, D k m3 permeability of matrix in USRV zone, D P sc pressure at standard condition, MPa q sc well production rate, m 3 /d s Laplace variable t time, h t D dimensionless time T temperature, K T sc temperature at standard condition, K j pseudo pressure, Mpa 2 /cp, j=F, f1, m1,m2,f2,m3 μ viscosity, cp f1-F inter-porosity coefficient from mi-f to ma-f in SRV m2-f1 inter-porosity coefficient from m2 to mi-f in SRV j storability coefficient, j=F, f1, m1, m2, f2, m3, d x, y coordination, m x D , y D dimensionless space q D dimensionless rate pore volume ratio x D , y D dimensionless space A cw the area of wellbore, 2Lh,m 2 m1-F inter-porosity coefficient from m1 to ma-f in SRV subscript D dimensionless superscript Laplace transform , derivative Dimensionless stress sensitivity factor D = ( i wf ) Dimensionless inter-porosity parameter Outer zone j = ( c t ) j

Fig. ( 1
Fig. (1).(a) The illustration of multi-fractured horizontal well (b) The two segments of matrix in SRV zone.
Fig. (3).Type curve and the division of flow regime.

Fig. ( 5 ).
Fig. (5).Type curve under different value of Langmuir parameters.(a) q D vs t D under different Langmuir volume; (b) q' D vs t D under different Langmuir volume; (c) q D vs t D under different Langmuir pressure; (d) q' D vs t D under different Langmuir pressure