An Analysis for the Influences of Fracture Network System on Multi-Stage Fractured Horizontal Well Productivity in Shale Gas Reservoirs

Deliang Zhang 1,2,*, Yu Dai 1, Xinhua Ma 1, Liehui Zhang 2, Bing Zhong 1, Jianfa Wu 1 and Zhengwu Tao 3 1 PetroChina Southwest Oil & Gas Field Company, Chengdu 610000, Sichuan, China; daiyucup@petrochina.com.cn (Y.D.); xinhuam@petrochina.com.cn (X.M.); zhongbing@petrochina.com.cn (B.Z.); wu_jianfa@petrochina.com.cn (J.W.) 2 Petroleum Engineering College, Southwest Petroleum University, Chengdu 610500, Sichuan, China; zhangliehui@vip.163.com 3 PetroChina Tarim Oilfield Company, Korla 841000, Xinjiang, China; taozw-tlm@petrochina.com.cn * Correspondence: maxwellzdl@gmail.com; Tel.: +86-028-86015496


Introduction
Shale gas reservoir is a typical unconventional resource, which has played an irreplaceable role in world energy consumption since the beginning of this century.The reservoir is characterized by self-sealing and self-sourcing, multi-scale porosity and extremely tight matrix with ultra-low permeability.Moreover, proven by huge amount of field practices, fractured horizontal well could be the only way to obtain commercial production from shale gas reservoir [1,2].In recent years, with the consideration of complex boundary constraints, dual-storage and multi-scale seepage mechanisms, flow dynamic models of fractured horizontal wells in shale gas reservoir are presented by many scholars.However, as a significant indicator of production rate, fracture network system still need comprehensive discussions.In this paper, based on the previous studies of multi-scale storage and seepage mechanisms and considered representative fracture network models, the partial differential flow equations in shale gas reservoirs are established.Then, with the combination of control-volume finite element method (CVFEM) numerical solution [3] and discrete fracture network algorithm, the production decline curves and pressure transient distributions are obtained.After that, the influences of diffusion, adsorbed gas, Langmuir volume and net-fracture system on production performances are analyzed.In particular, this paper focuses on the discussions about fracture network system.Different from conventional reservoirs, the pore diameter distribution in shale is characterized by a bimodal curve occuring in the 10-1000 nm range [4].Due to the brittleness of shale, after stimulation operation, fracture network system would become an indispensable part of the porosity.Obviously, the multimodal pore space may lead to the multi-mechanisms and complicated state equations of the gas flow in shale, such as ad-desorption, diffusion, slippage and stress sensitivity.
As reported, in shale formations, natural gas is generally trapped as both free gas and adsorbed gas [5].The existence of kerogen and higher contents of clay minerals in shale has undoubtedly caused the larger surface of nanopores [6,7].In contrast to conventional reservoirs, the content of gas adsorbed on the surfaces might be up to 85%, while the other 15% is held as free gas in macropore spaces [8].More importantly, according to Mengal's research [9], during the production period, the adsorbed gas can be extracted through the pressure depletion, which is an important contribution to keep the steady productivity in later period.In recent works, scholars [10][11][12] hold out that the ad-desorption theory in some shale should be described by the Brunauer-Emmett-Teller (BET) isotherm.However, in this paper, the ad-desorption law conforms to the Langmuir formula, which is proven by most observations [13][14][15][16].
Since the nanopore volume contributes the majority of shale porosity, gas flow in shale matrix should be described by multi-mechanisms (Figure 1).Ertekin et al. (1986) establish a dual mechanism model to describe the fluid flow in tight porous media [17].Comparing the Fick diffusion with the Klinkenberg equation, the gas flow in micropores is calculated by a dynamic apparent permeability factor.Many studies adopted a similar idea to model shale gas seepage mechanisms by modified apparent permeability coefficient.According to the definition, Roy and Raju (2003) obtained the Knudsen number by calculating the ratio of the gas molecular mean free path and shale pore radius [18].The results show that the slip flow and transitional flow are the typical flow regime in shale, which is different from the linear flow equation of conventional gas reservoirs.Ozkan et al. (2010) studied the flow mechanisms in shale matrix and micro fracture system [19], and developed a dual-mechanisms model by coupling diffusion with Darcy flow to describe the gas flow in shale matrix.After that, the influences of diffusion, adsorbed gas, Langmuir volume and net-fracture system on production performances are analyzed.In particular, this paper focuses on the discussions about fracture network system.Different from conventional reservoirs, the pore diameter distribution in shale is characterized by a bimodal curve occuring in the 10-1000 nm range [4].Due to the brittleness of shale, after stimulation operation, fracture network system would become an indispensable part of the porosity.Obviously, the multimodal pore space may lead to the multi-mechanisms and complicated state equations of the gas flow in shale, such as ad-desorption, diffusion, slippage and stress sensitivity.
As reported, in shale formations, natural gas is generally trapped as both free gas and adsorbed gas [5].The existence of kerogen and higher contents of clay minerals in shale has undoubtedly caused the larger surface of nanopores [6,7].In contrast to conventional reservoirs, the content of gas adsorbed on the surfaces might be up to 85%, while the other 15% is held as free gas in macropore spaces [8].More importantly, according to Mengal's research [9], during the production period, the adsorbed gas can be extracted through the pressure depletion, which is an important contribution to keep the steady productivity in later period.In recent works, scholars [10][11][12] hold out that the ad-desorption theory in some shale should be described by the Brunauer-Emmett-Teller (BET) isotherm.However, in this paper, the ad-desorption law conforms to the Langmuir formula, which is proven by most observations [13][14][15][16].
Since the nanopore volume contributes the majority of shale porosity, gas flow in shale matrix should be described by multi-mechanisms (Figure 1).Ertekin et al. (1986) establish a dual mechanism model to describe the fluid flow in tight porous media [17].Comparing the Fick diffusion with the Klinkenberg equation, the gas flow in micropores is calculated by a dynamic apparent permeability factor.Many studies adopted a similar idea to model shale gas seepage mechanisms by modified apparent permeability coefficient.According to the definition, Roy and Raju (2003) obtained the Knudsen number by calculating the ratio of the gas molecular mean free path and shale pore radius [18].The results show that the slip flow and transitional flow are the typical flow regime in shale, which is different from the linear flow equation of conventional gas reservoirs.Ozkan et al. (2010) studied the flow mechanisms in shale matrix and micro fracture system [19], and developed a dual-mechanisms model by coupling diffusion with Darcy flow to describe the gas flow in shale matrix.Javadpour et al. (2007) put forward an improved method for the Knudsen number calculation [20].Then, the flow regimes are classified through the improved formula.Considering shale pore structures with the calculated results, Javadpour (2009Javadpour ( , 2013) presented a gas apparent permeability to describe the seepage mechanism in nano-pore [21,22], which is a combination of viscous flows and Knudsen diffusion.Civan (2010Civan ( , 2011Civan ( , 2013) used a simplified second-order slip equation and several empirical parameters to model a seepage permeability formula in tight porous media, and by incorporating the Knudsen diffusion into the slippage equation [23,24], an apparent permeability formula for shale gas reservoir was proposed [25].According to Javadpour's results, Swami and Javadpour et al. (2007) put forward an improved method for the Knudsen number calculation [20].Then, the flow regimes are classified through the improved formula.Considering shale pore structures with the calculated results, Javadpour (2009Javadpour ( , 2013) presented a gas apparent permeability to describe the seepage mechanism in nano-pore [21,22], which is a combination of viscous flows and Knudsen diffusion.Civan (2010Civan ( , 2011Civan ( , 2013) used a simplified second-order slip equation and several empirical parameters to model a seepage permeability formula in tight porous media, and by incorporating the Knudsen diffusion into the slippage equation [23,24], an apparent permeability formula for shale gas reservoir was proposed [25].According to Javadpour's results, Swami and Settari (2012) coupled Darcy flow, Knudsen diffusion and slip flow to establish the nano-scale seepage mechanism in shale [26].As shown by the calculations, the apparent permeability of shale gas reservoir increased with the pressure depletion, which indicates the multi-mechanism has important influence on the later production period.Based on molecular dynamic theories, Deng et al. (2014) proposed a multi-mechanism model considering the gas diffusion, slippage effect and adsorption desorption by solving the Beskjok-Kaeniadakis equation [27].The permeability is treated as a function of the adsorption layer, which is not easy to determine for field practices.According to historical works, gas diffusion, slip flow, and Darcy flow, as the main flow regimes of gas flow in shale, are considered in this paper.Moreover, to facilitate the calculation, the apparent permeability is derived as a function of the Knudsen number.
Shale formations have a high content of lime, mica, calcareous and other brittle minerals.Take the Longmaxi formation, a representative shale formation in the Sichuan basin, as an example, the average brittle mineral content is up to 60%.After hydraulic stimulation, the structural fractures, shrinkage fractures and overpressure fractures can be connected by the artificial fractures.The widely distributed fracture network system with high conductivity is defined as the stimulated reservoir volume (SRV), which plays a significant role for commercial production in unconventional reservoirs.Thus far, different methods are used by scholars to study the influences of the SRV in shale.Wu (2009) established a triple-porosity model to describe the gas flow in shale [28].In this model, the stimulated reservoir is divided into micro and macro fracture system.Through detailed analysis, the effects of SRV system on economic benefits in unconventional reservoirs are revealed; however, un-stimulated zone and characteristics flow mechanisms of shale fail to be discussed.Based on the classical Warren-Root dual-porosity model, an equivalent Multi-stage fractured horizontal (MFH) model in shale [29] is presented by Bello and Watenbargen (2010).Although the flow regimes and production dynamics is discussed in detail in the work, the diffusion and adsorption effects in shale are neglected.Moreover, the model assumes that linear flow is the only regime in MFH well, which is not accurate and could lead to an unconvinced result in radial or elliptic flow periods.Since the linear flow assumption is convenient for derivation and solution, the linear flow model is widely used for research on unconventional reservoir [30][31][32][33].In fact, because of the early flow around fractures is linear, the linear flow results are in good agreement with the accurate solutions in the early period of time.Nevertheless, the existence of fracture with high conductivity makes the duration of the linear flow not assumable.Therefore, the liner flow model is more suitable to describe the flow in homogeneous and unconventional reservoirs without fracture network system.To overcome this limitation, scholars try to study the gas flow through numerical simulation.With the consideration of multi-scale flow mechanisms and ad-desorption, Jiang (2015) used dual permeability model to characterize the production dynamic in shale [34].The results show that the desorption gas is an important production supplements in low pressure period.Wang (2015) used finite difference method to solve a MFH well flow model in shale [35].The model combines the slippage effect, viscous flow and Knudsen diffusion mechanism of shale gas reservoirs.The calculation indicates that the fracture network system could form a low pressure area in shale reservoirs which can stimulate the gas desorption.Although the literature presents several significant and useful suggestions to shale reservoirs development, limited by the orthogonal grids, conventional bi-wing fracture and explicit algorithm, the contribution of fracture network system for gas extraction is still lacking from the discussion.
In this paper, to study the influences of characteristic mechanisms on flow dynamic in shale, the models give a comprehensive consideration of multi-scale flow mechanisms, such as diffusion flow, desorption and slippage flow.More significantly, both discrete fracture model and dual-porosity model are presented to elaborate the stimulation of fracture network for the gas flow rate (Figure 2).Then, with the assistances of unstructured grids, district fracture network (DFN) model and control volume finite element method (CVFEM), the flow equations are solved implicitly.Ultimately, the influences of characteristic parameters on production rate are discussed by sensitivity analyses, which provide some insights into the study of flow dynamic in shale.

Model Assumption
To describe the gas flow in shale, some idealizations and assumptions are made (Figure 3): (1) The model is derived for single-phase isothermal gas flow in shale.
(2) Both natural fractures and hydraulic fractures are main flow paths, and filled with the free gas.
Additionally, shale matrix is filled with part of the free gas and almost all of the adsorption gas.

Model Assumption
To describe the gas flow in shale, some idealizations and assumptions are made (Figure 3): (1) The model is derived for single-phase isothermal gas flow in shale.
(2) Both natural fractures and hydraulic fractures are main flow paths, and filled with the free gas.
Additionally, shale matrix is filled with part of the free gas and almost all of the adsorption gas.

Adsorption/Desorption
According to the similarity of coalbed methane, scholars have long noticed that the adsorption gas is a significant feature in shale.As Hill reported [8], approximately 30-60% of the gas occurs by adsorption mechanism.Based on the analysis of experiments, the majority of current research uses Langmuir's isotherm to describe the adsorbed gas [36][37][38][39].Moreover, the theory is also adopted in

Adsorption/Desorption
According to the similarity of coalbed methane, scholars have long noticed that the adsorption gas is a significant feature in shale.As Hill reported [8], approximately 30-60% of the gas occurs by adsorption mechanism.Based on the analysis of experiments, the majority of current research uses Langmuir's isotherm to describe the adsorbed gas [36][37][38][39].Moreover, the theory is also adopted in many commercial simulators [40].As the Langmuir equation shows, under constant temperature condition, the relationship between gas pressure and adsorbed volume can be expressed as: where V L and P L are the Langmuir volume and Langmuir pressure, respectively, which represent the adsorption capability and reflect the characteristic of the curve.

Flow Regimes and Mechanisms in Micro-Pores
Historic research of shale pore structures reveals that the pore distribution is bimodal.Figure 4 shows representative pore distributions curves of Longmaxi formation Sichuan basin, where the diameter of micropores range 10-100 nm, which is extremely tight compared to conventional formations.According to core experiments, the matrix permeability of shale ranges from 10 −5 to 10 −4 mD.

Adsorption/Desorption
According to the similarity of coalbed methane, scholars have long noticed that the adsorption gas is a significant feature in shale.As Hill reported [8], approximately 30-60% of the gas occurs by adsorption mechanism.Based on the analysis of experiments, the majority of current research uses Langmuir's isotherm to describe the adsorbed gas [36][37][38][39].Moreover, the theory is also adopted in many commercial simulators [40].As the Langmuir equation shows, under constant temperature condition, the relationship between gas pressure and adsorbed volume can be expressed as: where VL and PL are the Langmuir volume and Langmuir pressure, respectively, which represent the adsorption capability and reflect the characteristic of the curve.

Flow Regimes and Mechanisms in Micro-Pores
Historic research of shale pore structures reveals that the pore distribution is bimodal.Figure 4 shows representative pore distributions curves of Longmaxi formation Sichuan basin, where the diameter of micropores range 10-100 nm, which is extremely tight compared to conventional formations.According to core experiments, the matrix permeability of shale ranges from 10 −5 to 10 −4 mD.
inter-porosity flow flow in fractures adsorbed gas desorbed gas free gas The ratio of the mean free path to the characteristic length is defined as the Knudsen number, which is known as the key parameter to recognize flow regimes, and can be calculated by where L is effective pore radius, and λ is the molecular mean free path, which is given by where k B is the Boltzmann constant, and k B = 1.3806 × By substituting the real gas state equation into Equation (4), we obtain the molecular mean free path of real gas by λ = 3.16 × 10 9 µZ p πRT 2M g (5) According to the calculation results of Knudsen number, the regimes of gas flow in porous medium can be classified as free-molecular flow, transition flow, slip flow and continuum flow (Figure 5).Because the pore diameter range of conventional reservoirs is larger than 1 µm, the flow regimes of these formation is governed by continuum flow, which can be described by Darcy's law.Moreover, taking pore distributions of shale matrix into account, the corresponding regimes in shale micro-pores are slip flow and transition flow (the yellow region in Figure 5).Thus, the comprehensive flow model in shale must be a coupling mechanism of viscous flow, slip flow and Knudsen flow.

π p δ
where kB is the Boltzmann constant, and kB = 1.3806 × 10 −23 J/K; δ is the molecule collision distance in nm.Because it is difficult to calculate the collision distance, Civan (2010) put forward Equation (4) to calculate molecular mean free path.9 g 3.16 10 2 By substituting the real gas state equation into Equation (4), we obtain the molecular mean free path of real gas by 9 g 3.16 10 2 According to the calculation results of Knudsen number, the regimes of gas flow in porous medium can be classified as free-molecular flow, transition flow, slip flow and continuum flow (Figure 5).Because the pore diameter range of conventional reservoirs is larger than 1 μm, the flow regimes of these formation is governed by continuum flow, which can be described by Darcy's law.Moreover, taking pore distributions of shale matrix into account, the corresponding regimes in shale micro-pores are slip flow and transition flow (the yellow region in Figure 5).Thus, the comprehensive flow model in shale must be a coupling mechanism of viscous flow, slip flow and Knudsen flow.

Slip Flow
When Knudsen number is between 0.001 and 0.1, the flow caused by the collision between gas molecules is more obvious than the flow caused by the collision between gas molecules and pore surfaces.The flow regime is dominated by the slip flow, which can be expressed as

Slip Flow
When Knudsen number is between 0.001 and 0.1, the flow caused by the collision between gas molecules is more obvious than the flow caused by the collision between gas molecules and pore surfaces.The flow regime is dominated by the slip flow, which can be expressed as

Knudsen Flow
When Knudsen number is more than 0.1, the probability of gas molecules collision is equal to the collision between gas molecules and pore surfaces.The gas transport is defined as Knudsen diffusion.Previous works have proven that Knudsen diffusion contributes a significant part of gas transport flux in nanopores.Assuming that the molecule rebounds in accordance with the basic law of diffuse reflection, the Knudsen diffusion for real gas can be calculate by Energies 2018, 11, 414 7 of 19

Total Gas Flux
As shown in Figure 5, gas transport in shale micropores is governed by slip flow and transition flow regimes.Therefore, the flux equation is a coupling mechanism including viscous flow, slip flow and Knudsen flow, which makes the total flux in shale much more complex than conventional formation.Based on the mechanism research of the slip flow and Knudsen diffusion, and with the assumption that the interaction between the two mechanisms is negligible, the total flux equation can be obtained by combining Equations ( 6) and (7).

Flow-Governing Equations
According to the motion equation and mass balance law, the gas flow in shale is governed by the following equations.
For the shale matrix system For the fracture system For the hydraulic (discrete) fracture system where Q is the source intensity at wellbore, q * fF is the inter-porosity flux between fracture system and hydraulic fracture system, and q * mf is the inter-porosity flux between matrix and fracture system, which can be calculated by The flux between fracture and matrix system should be calculated accurately since the pressure difference is obvious between the two systems.Nevertheless, because of the low dimension of hydraulic fracture (2D fracture patches in 3D or Pseudo-3D reservoirs, 1D fracture lines in 2D reservoirs), and a better flow capacity in both fracture region and hydraulic fracture system, the pressure at discrete hydraulic fracture meshes can be considered equal to the adjacent fracture blocks.Based on the coupling condition, for fracture blocks containing discrete hydraulic fractures, the equations can be solved as a pseudo system.
q des is the desorption volume, which can be calculated form isothermal ad-desorption curves.Once the Langmuir adsorption isotherm is used to fit ad-desorption experiment results, q des can be written as k ζe (ζ = m, f, F) is the apparent permeability for each system.According to the flow mechanisms study in Section 2.3, and considering the stress sensitivity of fractures, it can be expressed by

Numerical Solution
Because of the accurate description for complex boundary conditions, unstructured grid system is preferred to mesh the complex flow region of shale reservoirs, rather than Cartesian grid system.In this section, the control volume finite element method (CVFEM) based on the unstructured triangle mesh is used to obtain numerical solutions (Figure 6).
kζe (ζ = m, f, F) is the apparent permeability for each system.According to the flow mechanisms study in Section 2.3, and considering the stress sensitivity of fractures, it can be expressed by

Numerical Solution
Because of the accurate description for complex boundary conditions, unstructured grid system is preferred to mesh the complex flow region of shale reservoirs, rather than Cartesian grid system.In this section, the control volume finite element method (CVFEM) based on the unstructured triangle mesh is used to obtain numerical solutions (Figure 6).Because of the different flow and desorption mechanisms, the governing Equations ( 9)-( 11) are distinguished.However, because the nonlinear items would be iterated by implicit solution, the matrix, fracture and hydraulic fracture systems can be expressed in a basic form.

∇ •
Equation ( 15) can be written in an equivalent integral form: Similar to the finite difference method, the partial differential items can be written as discrete equations using CVFEM: For the typical triangle element (Figure 7), transmissibility of CVFEM is where A is the area of the triangle element, and Energies 2018, 11, 414  Based on discrete fracture model (Figure 8) [41][42][43], the equivalent equation for the control volume region containing discrete fractures is, Through Equation ( 25), the element discretization equation of fracture and hydraulic fracture system can be written,

Results and Discussion
In Since the partial differential of time is independent of the space domain, the discrete form of Equation ( 16) is The physical meaning of Equation ( 19) is the material balance in the area controlled by node i of the triangle element.For each triangle element, the balance equations of three nodes can be obtained For the desorption item, it can be considered as V E is the adsorbed gas volume per unit rock framework volume (m 3 /m 3 ).
For the inter-porosity flux item, it can be written as According to the preceding derivation, the balance equations of the systems can be expressed as: W (e) is the inter-porosity flow matrix of element, is the element adsorption matrix, V Based on discrete fracture model (Figure 8) [41][42][43], the equivalent equation for the control volume region containing discrete fractures is, Energies 2018, 11, x FOR PEER REVIEW 11 of 19 Based on discrete fracture model (Figure 8) [41][42][43], the equivalent equation for the control volume region containing discrete fractures is, Through Equation ( 25), the element discretization equation of fracture and hydraulic fracture system can be written,

Results and Discussion
In  Due to the difficulty of precisely describing the branch fractures (either geological recognition or reservoir engineering calculation), the dual-porosity composite model can be used to represent a multi-stage fractured horizontal well (MFHW) in shale gas reservoir (Figure 10).Due to the difficulty of precisely describing the branch fractures (either geological recognition or reservoir engineering calculation), the dual-porosity composite model can be used to represent a multi-stage fractured horizontal well (MFHW) in shale gas reservoir (Figure 10).Due to the difficulty of precisely describing the branch fractures (either geological recognition or reservoir engineering calculation), the dual-porosity composite model can be used to represent a multi-stage fractured horizontal well (MFHW) in shale gas reservoir (Figure 10). Figure 11 compares the production rate curves of unconventional reservoirs and conventional homogeneous reservoir.At the beginning of production, high conductivity fracture network means more flow paths from reservoir to wellbore, which results in a higher production.At the latter part of production, it can be seen from the pressure distribution figures (Figure 12) that either dual-porosity model or discrete fracture model can form a large low pressure region, which improves the desorption gas production as well as reduces the flow resistance to the well.

SRV area (dual-porosity medium)
SRV area (dual-porosity medium) Figure 11 compares the production rate curves of unconventional reservoirs and conventional homogeneous reservoir.At the beginning of production, high conductivity fracture network means more flow paths from reservoir to wellbore, which results in a higher production.At the latter part of production, it can be seen from the pressure distribution figures (Figure 12) that either dual-porosity model or discrete fracture model can form a large low pressure region, which improves the desorption gas production as well as reduces the flow resistance to the well.According to the results discussion, both the dual-porosity model and the discrete fracture model can reflect the characteristic performance of SRV in shale.The contribution of fracture networks is summarized as (1) They improve the free gas flow capability in the early period.
(2) They increase the flow interface between the fractured and un-fractured zone.
(3) They stimulate the desorption of the adsorbed gas; Notably, once the fracture distribution can be obtained quantitatively, the DFN model may present a more accurate description of the gas flow dynamic and the inter-porosity flux, especially for the early period.In the mid-later period, because the low pressure region formed by the fracture network system may hide the discreteness of fractures, the dual-porosity model would obtain the satisfied results.Although the DFN model needs more complicated treatments and calculations than the dual-porosity model, limited by uncertainty of the fracture network, the dual-porosity model is still a useful and recommended method for numerical simulation in shale.In view of this, the  According to the results discussion, both the dual-porosity model and the discrete fracture model can reflect the characteristic performance of SRV in shale.The contribution of fracture networks is summarized as (1) They improve the free gas flow capability in the early period.
(2) They increase the flow interface between the fractured and un-fractured zone.
(3) They stimulate the desorption of the adsorbed gas; Notably, once the fracture distribution can be obtained quantitatively, the DFN model may present a more accurate description of the gas flow dynamic and the inter-porosity flux, especially for the early period.In the mid-later period, because the low pressure region formed by the fracture network system may hide the discreteness of fractures, the dual-porosity model would obtain the satisfied results.Although the DFN model needs more complicated treatments and calculations than the dual-porosity model, limited by uncertainty of the fracture network, the dual-porosity model is According to the results discussion, both the dual-porosity model and the discrete fracture model can reflect the characteristic performance of SRV in shale.The contribution of fracture networks is summarized as (1) They improve the free gas flow capability in the early period.
(2) They increase the flow interface between the fractured and un-fractured zone.
(3) They stimulate the desorption of the adsorbed gas; Notably, once the fracture distribution can be obtained quantitatively, the DFN model may present a more accurate description of the gas flow dynamic and the inter-porosity flux, especially for the early period.In the mid-later period, because the low pressure region formed by the fracture network system may hide the discreteness of fractures, the dual-porosity model would obtain the satisfied results.Although the DFN model needs more complicated treatments and calculations than the dual-porosity model, limited by uncertainty of the fracture network, the dual-porosity model is still a useful and recommended method for numerical simulation in shale.In view of this, the dual-porosity model is used to analyze the sensitivity in the next section.

Sensitivity Analyses
Different from conventional reservoirs, gas flow in shale is affected by many characteristic parameters, such as hydraulic fracture numbers, hydraulic fracture length, hydraulic fracture height, conductivity Langmuir volume and stress sensitivity factor.

Hydraulic Fracture Numbers (N)
Figure 13 shows the influence of hydraulic fracture number on production performance of the horizontal well.With the increase of hydraulic fracture numbers, gas production increases; such effect becomes smaller with increase of the fracture numbers.For the same fracture and horizontal well length, there is an optimal number for hydraulic fractures.At the late production period, flow will form an approximately rectangular area of low pressure, which can be considered as a potential body.Thus, the effect of number of fractures on the production becomes very small.After 500 days, the production rates of horizontal wells with different numbers of fractures are almost equal, and, according to the material balance principle, cumulative production of each well is close to the same after a long time.
Energies 2018, 11, x FOR PEER REVIEW 14 of 19 Figure 13 shows the influence of hydraulic fracture number on production performance of the horizontal well.With the increase of hydraulic fracture numbers, gas production increases; such effect becomes smaller with increase of the fracture numbers.For the same fracture and horizontal well length, there is an optimal number for hydraulic fractures.At the late production period, flow will form an approximately rectangular area of low pressure, which can be considered as a potential body.Thus, the effect of number of fractures on the production becomes very small.After 500 days, the production rates of horizontal wells with different numbers of fractures are almost equal, and, according to the material balance principle, cumulative production of each well is close to the same after a long time.

Hydraulic Fracture Length (xF)
Figure 14 shows the effect of the fracture length on the production of a MFH well in shale reservoir.As the fracture length increases, production rate and cumulative production increase.This is because the fracture half-length increase will significantly improve the single well controlling area with the same horizontal well length.The figure shows the production performance of a MFH well with fracture half-length of 50 m, 100 m and 150 m, corresponding, respectively, to model regions 1200 m × 300 m, 1200 m × 400 m and 1200 m × 500 m.Therefore, in the development of shale gas reservoir, extension of the fracture length can obtain higher production.3.2.2.Hydraulic Fracture Length (x F ) Figure 14 shows the effect of the fracture length on the production of a MFH well in shale reservoir.As the fracture length increases, production rate and cumulative production increase.This is because the fracture half-length increase will significantly improve the single well controlling area with the same horizontal well length.The figure shows the production performance of a MFH well with fracture half-length of 50 m, 100 m and 150 m, corresponding, respectively, to model regions 1200 m × 300 m, 1200 m × 400 m and 1200 m × 500 m.Therefore, in the development of shale gas reservoir, extension of the fracture length can obtain higher production.
reservoir.As the fracture length increases, production rate and cumulative production increase.This is because the fracture half-length increase will significantly improve the single well controlling area with the same horizontal well length.The figure shows the production performance of a MFH well with fracture half-length of 50 m, 100 m and 150 m, corresponding, respectively, to model regions 1200 m × 300 m, 1200 m × 400 m and 1200 m × 500 m.Therefore, in the development of shale gas reservoir, extension of the fracture length can obtain higher production.

Hydraulic Fracture Height (hF)
The fracture penetration ratio β, which is the ratio of hydraulic fracture height to reservoir thickness, is defined as

Hydraulic Fracture Height (h F )
The fracture penetration ratio β, which is the ratio of hydraulic fracture height to reservoir thickness, is defined as Figure 15 shows the influence of hydraulic fracturing degree on production of a MFH well in shale gas reservoir.Fracture penetration ratio directly affects the early production.After some days, an approximate rectangular low pressure region is created surrounding the multistage fractured horizontal well.At this time, the effect of fracture penetration ratio is not obvious.
Figure 15 shows the influence of hydraulic fracturing degree on production of a MFH well in shale gas reservoir.Fracture penetration ratio directly affects the early production.After some days, an approximate rectangular low pressure region is created surrounding the multistage fractured horizontal well.At this time, the effect of fracture penetration ratio is not obvious.

Hydraulic Fracture Conductivity (FcD)
The effect of hydraulic fracture conductivity on production is shown in Figure 16.It can be seen that production rate and cumulative production of the fractured horizontal well increase with increase of fracture conductivity.This effect is more obvious during the early-time flow period and diminishes with production time.Note that the increase of early production with fracture conductivity tends to maximize at certain level.As shown in the plot, for hydraulic fracture conductivity of 20 mD•m, the production rate and cumulative production are close to those of a horizontal well with infinite conductive fractures.
According to previous analysis, there is an optimal value of fracture conductivity for hydraulic fracturing design.For fracture conductivity of 20 mD•m, the calculated fracture permeability is 20D for fracture width of 0.001 m.Obviously, this is easy to achieve during normal fracturing operation.Therefore, it is not worth trying to improve well production through increasing hydraulic fracture conductivity.

Hydraulic Fracture Conductivity (F cD )
The effect of hydraulic fracture conductivity on production is shown in Figure 16.It can be seen that production rate and cumulative production of the fractured horizontal well increase with increase of fracture conductivity.This effect is more obvious during the early-time flow period and diminishes with production time.Note that the increase of early production with fracture conductivity tends to maximize at certain level.As shown in the plot, for hydraulic fracture conductivity of 20 mD•m, the production rate and cumulative production are close to those of a horizontal well with infinite conductive fractures.
According to previous analysis, there is an optimal value of fracture conductivity for hydraulic fracturing design.For fracture conductivity of 20 mD•m, the calculated fracture permeability is 20D for fracture width of 0.001 m.Obviously, this is easy to achieve during normal fracturing operation.Therefore, it is not worth trying to improve well production through increasing hydraulic fracture conductivity.
increase of fracture conductivity.This effect is more obvious during the early-time flow period and diminishes with production time.Note that the increase of early production with fracture conductivity tends to maximize at certain level.As shown in the plot, for hydraulic fracture conductivity of 20 mD•m, the production rate and cumulative production are close to those of a horizontal well with infinite conductive fractures.
According to previous analysis, there is an optimal value of fracture conductivity for hydraulic fracturing design.For fracture conductivity of 20 mD•m, the calculated fracture permeability is 20D for fracture width of 0.001 m.Obviously, this is easy to achieve during normal fracturing operation.Therefore, it is not worth trying to improve well production through increasing hydraulic fracture conductivity.

Langmuir Volume (V m )
The influences of Langmuir volume on production performance is shown in Figure 17.It is clearly observed that VL mainly influences the curves in the mid-later flow period, which is the inter-porosity flow between matrix and natural fracture system.The larger VL is, the greater production rate for a constant BHP; after the quick decline period, production curve becomes flat.This is because a bigger value of Langmuir volume represents a larger amount of adsorbed gas in shale matrix which decreases the formation pressure drop and supplements the gas capacity.The influences of Langmuir volume on production performance is shown in Figure 17.It is clearly observed that VL mainly influences the curves in the mid-later flow period, which is the inter-porosity flow between matrix and natural fracture system.The larger VL is, the greater production rate for a constant BHP; after the quick decline period, production curve becomes flat.This is because a bigger value of Langmuir volume represents a larger amount of adsorbed gas in shale matrix which decreases the formation pressure drop and supplements the gas capacity.

Stress Sensitivity Coefficient (α)
As Figure 18 shows, the production rate and accumulative production both decrease when stress sensitivity increases.According to Equation ( 14), when pressure declines, the higher stress sensitivity leads to higher decrease of natural fracture permeability.If we assume the stress sensitivity coefficient is 0.01, 0.05 and 0.1, after 3000 days, the recovery is 52%, 49% and 47%, respectively.

Stress Sensitivity Coefficient (α)
As Figure 18 shows, the production rate and accumulative production both decrease when stress sensitivity increases.According to Equation ( 14), when pressure declines, the higher stress sensitivity leads to higher decrease of natural fracture permeability.If we assume the stress sensitivity coefficient is 0.01, 0.05 and 0.1, after 3000 days, the recovery is 52%, 49% and 47%, respectively.
As Figure 18 shows, the production rate and accumulative production both decrease when stress sensitivity increases.According to Equation ( 14), when pressure declines, the higher stress sensitivity leads to higher decrease of natural fracture permeability.If we assume the stress sensitivity coefficient is 0.01, 0.05 and 0.1, after 3000 days, the recovery is 52%, 49% and 47%, respectively.

Conclusions
This paper focuses on the flow rate performance of a MFH well in shale.Based on the discrete fracture model and dual-porosity model, the influences of fracture network system is analyzed.Moreover, with the consideration of gas adsorption/desorption, stress dependence and multi-scale flow mechanisms in the equations, the sensitivities of the characteristic parameters in shale are discussed.The main conclusions of this paper can be summarized as follows: (1) The contribution of fracture networks on productivities can be described as:

Conclusions
This paper focuses on the flow rate performance of a MFH well in shale.Based on the discrete fracture model and dual-porosity model, the influences of fracture network system is analyzed.Moreover, with the consideration of gas adsorption/desorption, stress dependence and multi-scale flow mechanisms in the equations, the sensitivities of the characteristic parameters in shale are discussed.The main conclusions of this paper can be summarized as follows: (1) The contribution of fracture networks on productivities can be described as: (2) Once the fracture distribution can be obtained quantitatively, the DFN model may present a more accurate description of the gas flow dynamic and the inter-porosity flux in the early flow period.Nevertheless, limited by uncertainty of the fracture network, the dual-porosity model is still a useful and recommended method for numerical simulation in shale.(3) Sensitivity analyses reveal that the increases of hydraulic fracture number, length and Langmuir volume represent a higher flow rate, and the marginal effect of fractures number should be noticed in fracturing engineering.

Figure 1 .
Figure 1.The flow in nan-pores of shale reservoirs.

Figure 1 .
Figure 1.The flow in nan-pores of shale reservoirs.

Figure 2 .
Figure 2. Hydraulic multiple fractured horizontal well with the stimulated reservoir volume (SRV) and its simplified models.

( 3 )
The inter-porosity flux between shale matrix and fracture system is described by the pseudo-steady mechanism.(4) With the consideration of the stress sensitivity, a dual-porosity composite model and a DFN model are established to reflect the effects of fracture networks.(5)The main hydraulic fractures are assumed to be finite-conductivity.

Figure 2 .
Figure 2. Hydraulic multiple fractured horizontal well with the stimulated reservoir volume (SRV) and its simplified models.

Figure 3 .
Figure 3. Gas storage and flow mechanisms in shale gas reservoirs.

Figure 3 .
Figure 3. Gas storage and flow mechanisms in shale gas reservoirs.

Figure 4 .
Figure 4.The pore distributions curves of Longmaxi formation Sichuan basin.

Figure 4 .
Figure 4.The pore distributions curves of Longmaxi formation Sichuan basin.

Figure 5 .
Figure 5.The Knudsen number of different pore diameters.

Figure 5 .
Figure 5.The Knudsen number of different pore diameters.

Figure 6 .
Figure 6.The triangle meshes and control volume.Figure 6.The triangle meshes and control volume.

Figure 6 .
Figure 6.The triangle meshes and control volume.Figure 6.The triangle meshes and control volume.

Figure 8 .
Figure 8. Concise description of the discrete fracture model network (DFN) model.
is the cumulative matrix of element, N

Figure 8 .
Figure 8. Concise description of the discrete fracture model network (DFN) model.

Figure 8 .
Figure 8. Concise description of the discrete fracture model network (DFN) model.

Figure 9 .
Figure 9. DFN model mesh of a fracture network system in shale.

Figure 9 .
Figure 9. DFN model mesh of a fracture network system in shale.

Figure 9 .
Figure 9. DFN model mesh of a fracture network system in shale.

Figure 10 .
Figure 10.Dual-porosity model of a fracture network system in shale.

Figure 10 .
Figure 10.Dual-porosity model of a fracture network system in shale.

Figure 11 .
Figure 11.Curves of production rate and accumulative production of a multi-stage fractured horizontal well (MFHW) in shale.

Figure 12 .
Figure 12.The pressure distributions of different SRV models (3000 days).

Figure 11 . 19 Figure 11 .
Figure 11.Curves of production rate and accumulative production of a multi-stage fractured horizontal well (MFHW) in shale.

Figure 12 .
Figure 12.The pressure distributions of different SRV models (3000 days).

Figure 13 .
Figure 13.Effect of hydraulic fracture numbers increase on production of a MFHW in shale.

Figure 13 .
Figure 13.Effect of hydraulic fracture numbers increase on production of a MFHW in shale.

Figure 14 .
Figure 14.Effect of the fracture length on the production of shale gas reservoirs.

Figure 14 .
Figure 14.Effect of the fracture length on the production of shale gas reservoirs.

Figure 15 .
Figure 15.Effect of Hydraulic Fracturing Degree on Production of a MFHW in Shale.

Figure 15 .
Figure 15.Effect of Hydraulic Fracturing Degree on Production of a MFHW in Shale.

Figure 16 .Figure 16 .
Figure 16.Effect of Hydraulic Fracture Conductivity on Production of a MFHW in Shale.

Figure 17 .
Figure 17.Effect of Langmuir volume on production performance of a MFHW in shale.

Figure 17 .
Figure 17.Effect of Langmuir volume on production performance of a MFHW in shale.

Figure 18 .
Figure 18.Effect of stress sensitivity coefficient on production performance of a MFHW in shale.

Figure 18 .
Figure 18.Effect of stress sensitivity coefficient on production performance of a MFHW in shale.
free gas flow capability in the early period.(b)They increase the flow interface between the fractured and un-fractured zone.(c) They stimulate the desorption of the adsorbed gas.
Rgas constant, 8.314 J/(mol•K) T temperature, K V L Langmuir volume, sm 3 /m 3 V E adsorbed gas volume per unit rock volume at a constant pressure, sm 3 /m 3 Z the ratio of real gas volume and ideal gas volume, m 3 /m 3