Research on transient flow theory of a multiple fractured horizontal well in a composite shale gas reservoir based on the finite-element method

Journal of Natural Gas Science and Engineering 33 (2016) 587e598

Contents lists available at ScienceDirect

Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse

Research on transient flow theory of a multiple fractured horizontal well in a composite shale gas reservoir based on the finite-element method Rui-han Zhang a, b, *, Lie-hui Zhang a, Rui-he Wang b, Yu-long Zhao a, b, **, De-liang Zhang a a b

State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, China China National Oil and Gas Exploration and Development Corporation, Beijing 00724, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 9 December 2015 Received in revised form 23 May 2016 Accepted 25 May 2016 Available online 28 May 2016

A shale gas reservoir, as a special type of gas reservoir, is significantly different from a conventional gas reservoir. Gas flow in shale involves a complex process of multiple flow mechanisms (including adsorption, diffusion, slippage, and viscous flow in multi-scaled systems of nano-to macro-porosity) that substantially deviates from Darcy flow. The multiple hydraulic fracturing procedure not only creates the stimulated reservoir volume (SRV) to improve production but also causes the flow in the shale to become more complex. In this work, we first propose a dual porosity continuum and discrete fracture model to describe the flow of MFH well in a rectangular composite shale gas reservoir. The finite element method based on the unstructured 3D tetrahedral meshes is used to obtain a numerical solution. Next, we assume a transient well index to correct the early calculation error caused by the mesh precision and obtain the complete pressure transient and production decline curves. Sensitivity analyses show that the Langmuir volume, equivalent apparent permeability, hydraulic fracture distribution and size of SRV have great impacts on the transition flow regime and production performance. The research results obtained in this paper can provide theoretical guidance for efficient and large-scale development of shale gas reservoirs. © 2016 Elsevier B.V. All rights reserved.

Keywords: Composite shale gas reservoir Multiple flow mechanisms Finite element method Multiple fractured horizontal well Stimulated reservoir volume (SRV)

1. Introduction Shale gas is a form of self-sealing and self-sourcing, unconventional gas that is mainly stored in clay shale with clusters of discontinuous micro-fractures and an extremely tight matrix (Zhang et al., 2015a). To reduce the amount of pollution entering the environment, as a form of clean energy, demand for natural gas is increasing year by year. Shale gas is gradually becoming a key component of the world’s gas energy supply (Yu et al., 2014; Guo et al., 2015). Taking China as an example, its natural gas demand was 219 billion cubic meters in 2015 and is expected to grow to 411

* Corresponding author. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, China. ** Corresponding author. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, China. E-mail addresses: [email protected] (R.-h. Zhang), [email protected] (Y.-l. Zhao). http://dx.doi.org/10.1016/j.jngse.2016.05.065 1875-5100/© 2016 Elsevier B.V. All rights reserved.

billion cubic meters by the end of the year 2020. However, the conventional natural gas production in China is expected to be only 210 billion cubic meters by 2020. This situation reflects the importance of developing unconventional energy sources such as shale gas. (Source: http://www.zgnyb.com/html/trq/2015/0924/270.html) The United States is very advanced regarding shale gas development, According to the EIA (Energy Information Administration) statistics, the U.S. shale gas production reached 110 billion cubic meters in 2011 (Sun, 2014). However, the engineering of the efficient development of shale gas, especially regarding the special seepage law and its complex mechanism of gas flow, lag far behind in field practices. As a special type of gas reservoir, according to Wang et al. (2009), Javadpour (2009) and Yan et al. (2015), shale gas has complex multimodal pores, such as pores in organic matter (mainly kerogen), in tight inorganic matter, in micro natural fractures and in hydraulic fractures, causing shale gas to be transported by multiple flow mechanisms. Klinkenberg (1941) first discovered the existence of the slip flow of gas in porous media. When gas flows in nano-scale shale pores, the dimensions of these pores are

588

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

within the range of the mean free path of gas molecules, which causes the gas to slip on the rock surface (Florence et al., 2007; Rushing et al., 2007). According to the Knudsen number, many scholars (Ertekin et al., 1986; Roy et al., 2003) believe that diffusion flows occur in shale matrix nano-pores; these flows can be divided into Knudsen diffusion (Kn
10), Fick diffusion (Kn  0.1) and transition diffusion (0.1 < Kn < 10). These nano- and micro-scale pores are very important for shale gas because they can adsorb and store much shale gas (Aguilera, 2010). As reported, more than 85% of shale gas is adsorbed in these pores (Hill and Nelson, 2000). Thus, many researchers (Brunauer et al., 1940; Hatman, 2009; Zhang et al., 2014) used experimental results to developed theories to interpret the adsorption-desorption behavior of shale gas, with the most widely used theory obviously being the classical Langmuir isothermal adsorption theory (Langmuir, 1916). To date, researchers have attempted to describe gas flow in shale via a combination of the multiple flow mechanisms. Kucuk and Sawyer (1980) first studied the pressure transient behavior of a vertical well in a shale gas reservoir based on the dual-porosity media. Subsequently, some researchers (Bumb and McKee, 1988; Carlson and Mercer, 1991; Gao et al., 1994) represented by Bumb and McKee (1988) assumed an additional adsorption coefficient in Langmuir isotherm theory to study the effects of adsorption and desorption on the transient behavior. However, these works all neglected the diffusive flow in nano- and micro-scale pores. Carlson and Mercer (1991) studied the pressure behavior of a vertical well in shale by introducing the effects of desorption and diffusion into the classical dual-porosity model. They used the average pressure method to obtain the solution and ignored the slippage. Swami (2012, Swami et al., 2013) formulated a flow model that considers the Knudsen diffusion, slippage, adsorption, gas dissolved in kerogen bulk in nano-scale pores and Darcy flow in nature fractures. This model, which is based on dual-porosity media, was validated by laboratory data. Kuuskraa et al. (1992) proposed that the shale gas reservoir should be regarded as “tri-porosity” media containing three interconnected systems: conventional micropores system, nano-pore system, and natural fracture system. Based on the “tri-porosity” media model, Schepers et al. (2009); Dehghanpour and Shirdel (2011) coupled diffusion and adsorption flows into the matrix system. However, the flow in “triporosity” media was all described by Darcy flow. Owing to the extremely low permeability of a gas reservoir, field practice indicated that the key enabling technologies to obtain the industrial production from a shale gas reservoir have been long horizontal well and multiple fracturing stages (Mcclure and Horne, 2013). Scholars have presented many models to describe the flow dynamics and mechanisms of fractured horizontal wells in shale gas reservoirs. Al-Ahmadi et al. (2010); Aboaba and Cheng (2010) used the linear flow model to describe the typical curves and production decline of fractured horizontal wells in a shale gas reservoir, without considering the effect of adsorption and diffusion. Wang (2013) combined the point source function and the perturbation method to propose a well testing model for a MFH well in a shale gas reservoir that takes into account multiple flow mechanisms, stress-sensitivity in natural fractures, and hydraulic fracture orientation. However, the models proposed above all ignored the effect of the stimulated reservoir volume (SRV) on the shale gas reservoir. To concisely describe the fracture network (natural or induced) around the hydraulic fractures, Garipov et al. (2016) established the complex fracture network around the hydraulic fractures due to the stress. However, determining how to improve the accuracy of fractures distribution is a difficult task, and the large amount of short discrete fractures make the efficiency of calculation low. Xu et al. (2015a,b) and Zhang et al. (2015b) all assumed the MFH

well to be in a circular composite reservoir, where the fracture network is simplified by a dual-porosity model to describe the performance of the well in a shale gas reservoir. However, as Fig. 1 shows, the rectangular drainage area better conforms to field practice. Some scholars (Brohi et al., 2011; Ozken et al., 2011; Stalgorova and Mattar, 2012; Zhang et al., 2015a) represented by Zhang et al. (2015a) presented the effective trilinear flow or extended five region flow model, where the SRV is the inner region to describe the performance of the MFH well in a rectangular composite shale reservoir. These models take into account the multiple flow mechanisms of shale gas and are easy to solve. However, these models can only calculate the performance of the linear regime when the hydraulic fractures are uniform; they cannot reflect the complete flow regime response, such as the pseudo-radial flow regime and interference between fractures. To solve this complex MFH well in a rectangular composite shale reservoir, Zhao (2015) established a semi-analytical model and used the finite-difference method based on structured meshes. However, they did not take the hydraulic fracture distribution into account. Some researchers (Moinfar et al., 2013; Jiang and Younis, 2015) proposed a coupled dual continuum for composite reservoir and discrete fractures for the hydraulic fracture model for the simulation of an unconventional reservoir. Fan (2013) and Sun (2014) used this coupled model to study the pressure behavior and production performance of a MFH well in a shale gas reservoir. They used the finite-element method based on unstructured meshes to obtain the numerical solution without considering the dual-porosity media of a natural reservoir or the diffusion flow mechanism in shale. In view of the above-described situation, this paper presents a proposed dual-porosity continuum and discrete fracture model to describe the flow of a MFH well in a rectangular composite shale gas reservoir. This model takes into account the multiple flow mechanisms of a shale gas reservoir such as adsorption, diffusion, and viscous flow in the fracture network (natural or induced) and the stimulated reservoir volume. The finite element method based on the unstructured 3D tetrahedral meshes is used to obtain a numerical solution. We assume a transient well index to describe the early period of the MFH well pressure response. The main flow regimes are identified as follows: early-linear in hydraulic fracture; inter-porosity flow in the SRV and outer region; desorption and diffusion flow; and the later pseudo-radial flow. Sensitivity analyses show that the Langmuir volume, diffusion, slippage, hydraulic fracture distribution and size of SRV have great impacts on the transition flow regime and production performance. 2. Mathematical model 2.1. Model assumption Fig. 2 shows a schematic diagram of the volume fractured horizontal well in a composite shale gas reservoir. The x-y-axes represent the horizontal directions, and the z-axis represents the vertical direction. The inner SRV and outer region have distinct properties. To derive a practical model, the following idealizations and simplifying assumptions are made: 1. The model is for isothermal single-phase gas flow through a pseudo-pressure transformation. 2. A rectangular composite matrix block with dual-porosity idealization is used to simulate the stimulated zone, and the un-stimulated area is also assumed to be a dual-porosity porous medium that may have distinct properties. 3. Flow from the reservoir into the horizontal well is only through hydraulic fractures. Natural fractures are the main flow paths

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

589

Fig. 1. MFH well with SRV in a composite shale gas reservoir.

z

x

y

(a) 3D tetrahedral meshes for a rectangular composite reservoir (The red region is the SRV area; the blue plane represents the hydraulic fractures, which are dispersed by Delaunay triangulations; the horizontal well is simplified by a 1D solid that neglects flow resistance in the wellbore) z

1

Fracture angle

i

Hydraulic fractures N

y x h

z

(b) Diagram of a MFH well in shale gas

Horizontal well

Partially fractured

Fig. 2. Schematic diagram of the MFH well in a composite shale gas reservoir. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

and are full of free gas. The absorbed gas is trapped on the surfaces of the organic materials. 4. The inter-porosity transfer from the matrix to the net fracture system can be described by Fick’s first law, which corresponds to the pseudo-steady diffusion mechanism. 5. Infinite-conductivity hydraulic fractures are assumed to be distributed along the horizontal well. Moreover, the stress distribution and reservoir anisotropy may affect the fracture extension direction and open degree. In this paper, hydraulic fractures are

assumed to be at an arbitrary azimuth angle or partially fractured, which is more practical than the existing flow models. 2.2. Mathematical model and numerical solution A shale gas reservoir has a complex pore structure. The multiscale pore sizes and the unique storage properties cause shale gas to migrate by multiple flow mechanisms, including desorption from the organic surfaces, diffusion and slippage flow in nano-scale pores, and viscous flow in natural fractures (as shown in Fig. 3).

590

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

Fig. 3. Micro flow mechanism in a shale reservoir.

There are a large number of micro fractures and a certain number of macro pores in the shale gas reservoir. Moreover, under the action of hydraulic fracturing, hydraulic fracture networks are formed around the well. The scale of this system is relativity large. We can use Darcy’s law to describe flow in this system:

v¼

k

mg

Vp

kapp v¼ Vp

mg

(2)

kapp ¼ Fkm þ Dg cg mg Klinkenberg (1941), who first discovered the existence of the slip flow of gas in porous media, described slip flow velocity as follows:




8pRT F ¼1þ M

0:5

mg pavg rm




 2 1 f

S0 Mg vq ¼ NA vt

VL pL

RTsc

ðpL þ pmc Þ2

8
 > q dn 8RT 0:5 > > > DK ¼ Kn
10 > > t 3 pM > > > < Dg ¼ DF ¼ q km T Kn  0:1 > t 3pmg dg > > > > > 1  > > > : D1 þ D1 0:1 < Kn < 10 K

(3)

(7)

(4)

h i vm l2 qkam mmD2  mfD2 ¼ ut2 mD2 vtD

(8)

h i vm V2 mfD2 þ l2 qkam mmD2  mfD2 ¼ u2 uf2 fD2 vtD

(9)

Flow in the SRV region For the matrix system

(10)

For the fracture system

h i vm M12 V2 mfD1 þ M12 l1 qkam mmD1  mfD1 þ q* ¼ u1 uf 1 fD1 vtD (11)

K

(5)

According to the Langmuir adsorption isotherm (Langmuir, 1916), the following equation can be obtained:

V p ¼ VL pL þ p

Based on the assumption of the physical model, the seepage equations for the shale gas reservoir can be derived by combining the mass conservation equation, the equation of state and motion (Wang, 2013). Thus, we can obtain the dimensionless form of flow equations for a multiple fractured horizontal well in a composite reservoir as given by Eqs. 8e11. In view of the length of the paper, the specific derivation and the solution process are shown in Appendices A, B and C. Flow in the outer region For the matrix system

h i vm M12 l1 qkam mmD1  mfD1 ¼ ut1 mD1 vtD

To describe desorption mechanisms from the matrix cores, the fractional coverage at equilibrium is given by

kads P kdes þ kads P

vpmc vt

For the fracture system

The diffusivity constant, Dg, can be calculated using the following correlation given by Li et al. (2014):

q¼

q*des ¼

  Mg psc 1  fmc  ff rs

(1)

According to the primary work of Javadpour (2009), the concept of equivalent apparent permeability based on micro-continuum theory is used to combine the gas slippage in nanopores and diffusion.

q¼

consistently justified by the change of fractional coverage (Shabro et al., 2011). We can obtain:

(6)

The effective gas desorption flux per volume must be

To solve Eqs. 8e11, the inner and outer boundary conditions should be considered. The boundary condition assumed here is the rectangular sealed boundary of the reservoir, which can be substituted by 0 vectors according to the finite element methods (Chen et al., 2003). The process of solving the inner boundary condition of an infinite-conductivity fractured horizontal well using the numerical well model method is shown in Appendix C. 3. Model verification In this paper, the model of a rectangular composite reservoir

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

4. Results and discussion 4.1. Pressure profiles simulation Fig. 6 shows the pressure profiles for different times at a constant production rate. In this section, we assume that the 1200-m length of the MFH well is surrounded by a 1600  400 m SRV region. The pressure initially declines around the MFH well and then waves outward to the reservoir boundary. In detail, pressure waves spread along the hydraulic fracture direction; the velocity of the pressure wave in the SRV region is higher than in the unstimulated outer region.

1E+2

mwD,m'wD*tD

1E+1

CD=0;S=0;Lf=30m;L=1200m;Nf=5 xe=ye=50000m; L=300m

1E+0

1E+1

Moinfar method 1E+0

FEM KAPPA

m'wD*tD

with infinite-conductivity fractured horizontal wells is established for a shale gas reservoir. If the multiple flow mechanisms and the induced fracture networks around the hydraulic fractures (SRV) region are neglected, then the model can be applied for a MFHW in conventional single-porosity gas reservoirs. To verify the accuracy of the finite element method, comparisons are made with the solution using the semi-analytical method (Zhao, 2015; Zhang et al., 2015b). As shown in Fig. 4, the type curves and the derivative type curves for a constant production rate from our method are in good agreement with those from semi-analytical methods. However, as the green dashed lines shows, at the early linear flow period, the numerical method would produce a calculation error caused by the mesh precision. We should assume a transient well index (TWI) to correct the early calculation accuracy following the approach of (Al-Mohannadi et al., 2004; Aguilar et al., 2007). Moreover, in this section, we compare the predictions of FEM with the EDFM approach proposed by Moinfar et al. (2013). Both simulations are applied for a 2000-ft-length horizontal well with three equally-spaced 200-ft-length transverse fractures in a bounded rectangular conventional reservoir. As shown in Fig. 5, the dimensionless pressure derivative curve of the FEM in this paper is lower than the curve of the EDFM in the early period and midperiod; this difference occurs because Moinfar et al. (2013) assumed that the hydraulic fracture has an aperture of 0.03 ft and a permeability of 10,000 md, whereas we assume that the fracture is of infinite conductivity. The larger the fracture conductivity is, the lower the pressure and pressure derivate curve will be. In addition, a dimensionless pressure derivative curve of infinite conductivity MFHW in the same model is obtained using commercial software KAPPA. We conclude that our method can match the results of the commercial software for a conventional reservoir and that our simulator can more comprehensively simulate the production and pressure performance of MFHW in a shale gas reservoir considering multiple flow mechanisms and SRV regions.

591

1E-1

1E-2 1E-2

1E-1

1E+0

tD

1E+1

1E+2

1E+3

Fig. 5. Comparison of the dimensionless WBP derivate versus the dimensionless time calculated using FEM, EDFM and commercial software KAPPA.

4.2. Type curves for a MFH well in a composite shale reservoir In this section, we will analyze the pressure response obtained using the finite element method in this paper. The basic parameters used are from Table 1. Fig. 7 shows the type curves of pseudopressure and the derivative for a MFH well with the induced fracture network in a rectangular composite shale gas reservoir. If we assume that the supply boundary and the volume of stimulated area are large enough, then the curves may be cataloged into the following flow periods. Stage 1: early linear flow exists in the hydraulic fracture system and is characterized by a slope of 1/2 in both the pseudo-pressure and pseudo-pressure derivative curves on a log-log scale; the position of the straight line is controlled by the permeability of the SRV, and the duration time is related to the permeability, fracture length and space. Stage 2: an early pseudoradial flow period in the hydraulic fracture system that is characterized by a horizontal line with the value of 1/(2*NF*M12); it is sometimes hidden if the width of the SRV is similar to the hydraulic fracture length. Stage 3: a transition period and the pseudo-radial flow period in the SRV region; this period is extremely sensitive to the permeability of the SRV and is characterized by a horizontal line with the value of 1/(2*M12), i.e., the higher is the permeability, the lower the linear position. Stage 4: a transition period occurs from the pseudo-radial flow period in the SRV region to the later pseudo-radial flow period in the outer region. Stage 5: a later pseudo radial flow period during which the fluid exchange between the matrix and the fracture system reaches a dynamically balanced state; thus, the value of the horizontal line is 0.5. Stage 6: boundary response period, which is characterized by a slope of 1 in both the pseudo-pressure and pseudo-pressure derivative curves. As shown by the green line in Fig. 7, if we assume that the sizes of the volume in the SRV and the shale gas reservoir are within a reasonable range, after the early linear flow in the fracture system, the pressure wave transmits to the near boundary quickly and the curve goes upwards. The curves of pseudo-pressure and the pseudo-pressure derivative become parallel lines, and the horizontal line is subject to the combined action of inter-porosity transfer and boundary response.

1E-1 Semi-analytical solution Numerical solution TWI model

1E-2 1E-3 1E-4

1E-2

1E+0

1E+2 tD

1E+4

1E+6

Fig. 4. Comparison of the semi-analytical solution and the numerical solution.

4.3. Sensitivity analyses In this section, we will discuss the effects of some key parameters of the pressure response and production behavior in shale gas. The basic parameters are given in Table 1. To clearly display the differences, all production decline curves are drawn with dimensional variables in the real domain.

592

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

Fig. 6. Pressure profiles after different times simulated using the FE method.

Table 1 Data used in the analysis of the examples. Initial reservoir pressure, pi, MPa Formation thickness, h, m Gas specific gravity, rg, fraction Slip-flow coefficient, F Langmuir pressure, PL, MPa Initial viscosity, mi, cp SRV region Fracture permeability, kf1, mD Fracture porosity, :f1 Matrix permeability, km1, mD Matrix porosity, :m1 Length, m Width, m

1E+2 5

0.5 Horizontal line

2 1

1/(2*M12)

1E-2

1/(2NF*M12)

1E-3 1E-3

1E+1

1E-1 xm=1600m,ym=400m,VL=0 xm=1600m,ym=400m,VL=3 xm=3200m,ym=800m,VL=3

1E-2

0.5 1E-2

0.01 0.02 0.0001 0.12 6000 1500

1E+0

3

1E-1

330 30 1 10e6 3 0.01

6

4

1E+0

mwD,m'wD*tD

0.1 0.1 0.0001 0.12 1600 400

xe=6000m,ye=1500m;xm=1600m,ym=400m xe=ye=50000m;xm=ym=5000m

1E+1

Reservoir temperature, T, K Half fracture length, xf, m Well production rate, qsc, 104 m3/d Knudsen diffusion coefficient, Dg, m2/s Langmuir volume, VL, cm3/g Initial compressibility of gas, cgi, MPa1 Outer region Fracture permeability, kf2, mD Fracture porosity, :f2 Matrix permeability, km2, mD Matrix porosity, :m2 Reservoir length, m Reservoir width, m

mwD,m'wD*tD

1E+2

15 50 0.65 1.5 4 0.015

1E-1 1E+0 1E+1 1E+2 1E+3 1E+4 1E+5 1E+6 1E+7 tD

1E-3 1E-3

1E-2

1E-1

1E+0

1E+1 tD

1E+2

1E+3

1E+4

1E+5

Fig. 7. Type curves for a MFH well in a composite reservoir.

Fig. 8. Effect of Langmuir volume (VL) on the type curves.

4.3.1. Effect of Langmuir volume Figs. 8 and 9 show the effect of the Langmuir volume (VL) on the well testing type curves and the production decline performance of a MFH well in a rectangular composite reservoir. Clearly, VL mainly influences the curves in the mid-later flow period, which corresponds to the inter-porosity flow between matrix and the natural fracture system. The larger the value of VL is, the greater the production rate for a constant WBP. This behavior is observed because a higher value of Langmuir volume represents a larger amount of gas adsorbed on the surface of the shale matrix. As the pressure declines, gas desorbs from the surfaces, thereby delaying the velocity of the pressure wave and supplementing the gas capacity. Moreover, for an extremely tight reservoir, the SRV region has an important impact on the pressure and production performance.

The larger SRV represents a more effective stimulation, which definitely leads to a greater production. Fig. 9 shows in detail the production rate decline when considering the VL influence. Regarding the production of the MFH well in a shale gas reservoir, the production rate is observed to decline quickly after 1e2 years. In addition, the larger VL and SRV could result in a greater production rate. After 1000 days of production, the production rate of VL ¼ 3 cm3/g and a larger SRV is 1.46  104 m3/d; however, it is only 0.61  104 m3/d for the case of VL ¼ 0 cm3/g and smaller SRV. If we set the economic limit production as 104 m3/d, the production duration will last approximately 1800 days for the situation of VL ¼ 3 cm3/g and a larger SRV, which is 3.5 times the duration of the situation of VL ¼ 0 cm3/g and smaller SRV.

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

593

1E+2

1E+2

θ=3

θ=5

Production Rate (104m3/d)

Productin Rate(104m3/d)

θ=1 1E+1

1E+1

1E+0

1E+0 xm=1600m,ym=400m,VL=0 xm=1600m,ym=400m,VL=3 xm=3200m,ym=800m,VL=3 1E+1

1E+2 Time(d)

1E+3

1E+4

1E-1 1E-1

1E+5

Fig. 9. Effect of Langmuir volume (VL) on the production performance curves.

4.3.2. Effect of the apparent permeability The analysis above suggests that gas flow in multiple pore types of a shale matrix leads to serious deviation from Darcy’s flow. An integrated method of equivalent apparent permeability is used to describe the slippage and diffusion flow in the shale matrix system. As shown in Fig. 10, to study the equivalent apparent permeability on shale gas production performance, we assume some different values of q, which represents the ratio of the apparent permeability to the primitive shale matrix permeability. The apparent permeability is observed to mainly influence the curves in the mid-later flow periods; the larger the value of q is, the higher is the gas production rate. At early times, the majority of production is the free gas is compressed in fractures; subsequently, the gas that was compressed and adsorbed on the surface of the shale matrix began diffusing into the fracture system. Coefficient q is a part of the interporosity coefficient l; as the pressure declines, a larger q indicates that more gas diffuses into the fracture system, ultimately leading to a higher and more durable production rate. After 1000 days, the production rate of q ¼ 5 is almost 2.37 times that without considering slippage. The diffusion effect is shown as the blue line in Fig. 11. 4.3.3. Effect of the number of hydraulic fractures Natural shale formation consists of clusters of discontinuous micro-fractures and an extremely tight matrix, and as a result, the natural production is extremely low. Hydraulic fracturing technology is the main means to improve the production in shale gas development. Figs. 12 and 13 show the effect of the hydraulic fractures number on the well testing type curves and the production decline performance. More specifically, the number of fractures mainly influences the curves in the early linear-flow period.

1E+0

1E+1

1E+2 Time (d)

1E+3

1E+4

1E+5

Fig. 11. Effect of apparent permeability (q) on the production performance curves.

1E+2 NF=3

NF=5

NF=7

1E+1 1E+0

mwD,m'wD*tD

1E+0

1E-1 1E-2 1E-3 1E-3

1E-2

1E-1

1E+0

tD

1E+1

1E+2

1E+3

1E+4

Fig. 12. Effect of fracture number (NF) on the type curves.

1000 NF=3 Production Rate(104m3/d)

1E-1 1E-1

NF=5

NF=7

100

10

1

0.1 1E-2

1E-1

1E+0

1E+1 1E+2 Time (d)

1E+3

1E+4

1E+5

Fig. 13. Effect of fracture number (NF) on the production performance curves.

1E+2 θ=1

θ=3

As the number increases from 3 to 7, the production rate increases, but the range of increase decreases gradually. For a horizontal well with a fixed length, an optimal value of fracture number exists, considering the costs and benefits.

θ=5

1E+1

mwD,m'wD*tD

1E+0 1E-1 1E-2 1E-3 1E-3

1E-2

1E-1

1E+0

tD

1E+1

1E+2

1E+3

Fig. 10. Effect of apparent permeability (q) on the type curves.

1E+4

4.3.4. Effect of distribution of hydraulic fractures As shown in Fig. 14(a)e(d), in this section, we assume some different situations of the distribution of hydraulic fractures to study the respective changes in the production decline rate. From the results shown in Fig. 15, we can conclude that the longer the fracture length, as shown in situation (c), the higher the production rate will be. When the length is the same, the distribution of situation (a) has a higher production rate value compared to the partial fracturing situation (b) and the dip angle fractures situation (d)

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

z

1E+2

z

y

(a) Vertical and fully fractured x y

y

(b) Vertical and partially fractured x y

(c) Longer fracture length than (a)

(d) Fractures with arbitrary angles

Fig. 14. Diagram of the fracture distribution.

situation(a) Situation(b) Situation(c) situation(d) D-value(a vs d)

1E+2

4

2 1E+0

1E-1

1E+0 1E+1 Time (d)

1E+2

0 1E+3

Fig. 15. Effect of fracture distribution on the production performance curves.

because the uniform spacing and vertical fractures reduce the influence between fractures and fully implemented hydraulic fracturing results in a larger SRV region. These conclusions can provide a theoretical basis for the hydraulic fracturing operation in the field practice. 4.4. Comparison of traditional and rectangular composite models Figs. 16 and 17 compare the type curves and production performance of a fractured horizontal well in a shale gas reservoir without SRV model and in a rectangular composite reservoir with SRV model. In general, the existence of the SRV can obviously improve the flow capacity, leading to a lower pressure drop and a higher production rate. For a constant rate condition, the curves of

1E+2

Reservoir with SRV Reservoir without SRV

mwD,m'wD*tD

1E+1 1E+0 1E-1 1E-2 1E-3 1E-3

1E-1

1E+1 tD

1E+3

1E+0 Reservoir with SRV Reservoir without SRV 1E-1 1E-1

1E+0

1E+1 Time(d)

1E+2

1E+3

6

1E+1

1E-1 1E-2

1E+1

Fig. 17. Comparison of the production performances in a reservoir with and without SRV.

Production Rate Difference (104m3/d)

Production Rate (104m3/d)

1E+3

Productin Rate(104m3/d)

594

1E+5

Fig. 16. Comparison of the type curves in reservoirs with and without SRV.

pseudo-pressure and derivative for the MFH well in the composite reservoir model are lower than those for the model without SRV. After 100 days of production, with the consideration of SRV effects, the production rate under the constant BHP is 4.007  104 m3/d; however, it is just 1.440  104 m3/d for a reservoir without SRV.

5. Conclusions Considering the multiple flow mechanisms of a shale gas reservoir, we investigated a comprehensive transport model for MFH wells in a rectangular composite shale gas reservoir and analyzed the transient pressure response and production performance. The main conclusions of this paper are as follows: (1) A dual porosity and discrete fracture model considering the multiple flow mechanisms to describe the flow of a MFH well in the composite shale gas reservoir is established, and the solution is obtained using the FE method based on the 3D tetrahedral meshes. The numerical solution was verified by comparison with the analytical method, revealing a perfect agreement. Moreover, the simulator in this paper can easily describe the production performance of a horizontal well with asymmetric hydraulic fractures (each fracture has a different length, different angle, and different fracturing efficiency) without establishing the complex Green function calculation formula as an analytical method. (2) The fracture network induced by hydraulic fracturing has a great effect on the well production rate. The production rate for a fractured well considering the SRV region is much higher than the model without SRV. The larger the SRV is, the higher the production rate will be, and also the higher the fracturing costs will be. We should find a balance between better benefits and lower costs. (3) The distribution characteristics of the hydraulic main fractures are always ignored by the existing model; however, they may contribute significantly to the well production performance. The longer the fracture length is, the higher the production rate will be. For a horizontal well with fixed length, an optimal fracture number exists. The vertical and fully hydraulic fractures obtain a higher production rate than the dipping or partial hydraulic fractures. (4) Some key parameters have great effects on typical curves. A large Langmuir volume (VL) indicates that more gas will be desorbed with the pressure declines, resulting in a higher production rate at mid-later periods. The slippage and diffusion flow manly influence the inter-porosity flow

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

between the matrix and natural fractures; the higher the slippage and diffusion coefficients, the deeper the trough and the higher the production rate will be.

595

inter-porosity coefficient, 1/m2 number of sites available for gas adsorption per surface area Avogadro’s number, dimensionless

a S0 NA

Acknowledgments This work was supported by the National Natural Science Foundation of China (Key Program) (Grant No. 51534006), the National Science Fund for Distinguished Young Scholars of China (Grant No. 51125019). The authors appreciate the feedback of the reviewers and editors, whose critical comments were very helpful in preparing this article. Nomenclature cd cg Dg F h hD k kapp Lref m Mg NF p pL qdes qm qsct R rm T t ta u vkm VL vm x xD y yD z z zD Subscript D f F g m sc t

additional desorption compressibility, MPa1 gas compressibility, MPa1 diffusivity coefficient, m2/s gas slippage factor, dimensionless effective thickness, m dimensionless effective thickness, m permeability, mD apparent permeability, mD reference length, m pseudo-pressure, MPa2/(mPa$s) gas molecular weight, g/mol number of net nodes include hydraulic fractures, dimensionless pressure, MPa Langmuir pressure, MPa diffusion rate from matrix particle to micro-pore, m3/s inter-porosity flux rate from matrix to natural fracture system, m3/s total well production rate, m3/s gas constant, 8.314 J/(mol$K) radius of matrix core, m reservoir temperature, K time, s pseudo-time, dimensionless slippage velocity in matrix, m/s Knudsen diffusion velocity in matrix, m/s Langmuir volume, cm3/g flow velocity in matrix, m/s distance on x direction, m dimensionless distance on x direction, m distance on y direction, m dimensionless distance on y direction, m distance on z direction, m gas compressibility factor, fraction dimensionless distance on z direction, m

Appendix A. Derivation of the solution for a MFH well in a composite shale gas reservoir Considering multiple flow mechanism processes and interporosity transfer described by Fick’s law, the mass conservation equation is established as follows: Outer region For the matrix system

qm ¼

  v fm rgm vt

þ qdes

(A-1)

where qm is the inter-porosity flux rate given by Fick’s first law:

qm ¼ a


 kapp pm Mg p Mg pm  f pf mg zRT zRT

(A-2)

In summary, the right side of the equation is changed to:

  v fm rgm vt

Mg pm vpm cgm RT z vt   Mg psc 1  fm  ff rs

þ q*des ¼ fm þ

RTsc

VL pL vpm ðpL þ pm Þ2 vt (A-3)

where cgm is the gas compressibility coefficient without considering the adsorption-desorption process. The second term is the effective gas desorption flux per volume, as described by Eq. (7). To obtain Eqs. 8e11, the pseudo-pressure is defined by:

Zx   p’ dp’ m px ¼ 2 mg z p

(A-4)

p0

Consequently, the differential governing equations of matrix flow are:

2a

 i f m cmg vmðp Þ kapp h m g m mðpm Þ  m pf ¼ vt kf kf fm mg cd vmðpm Þ þ vt kf

(A-5)

where dimensionless natural fracture system hydraulic fracture gas matrix system standard state total

Geek symbols mg viscosity, mPa$s rg gas density, g/cm3 rs shale rock density, g/cm3 f porosity, fraction



cd ¼

rs T psc 1  fm  ff fm

Tsc

 VL pL z ðpL þ pm Þ2 pm

(A-6)

Following the approach proposed by Escobar et al. (2013), a pseudo-time function to account for the time dependence of the gas viscosity and the total system compressibility can be defined as

Zt ta ¼ t0

mi cgi dt mðtÞcg ðtÞ

(A-7)

Ultimately, Eq. (A-3) can be linearized into the following forms:

596

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

2a

 i f m c kapp h m gi mgi vmðpm2 Þ mðpm2 Þ  m pf2 ¼ vta kf2 kf2 fm mgi cd vmðpm2 Þ þ vta kf2

(A-8)

For a natural fracture system Seepage equation in the fracture network could be expressed as:



 v ff rf  V rf vf þ vt



    fm cmgi þff cfgi ff cfgi 1   1 ; uf1 ¼  u1 ¼  fm cmgi þff cfgi fm cmgi þff cfgi 1þ2 1 ð1 ¼ SRV;2     fm cmgi þff cfgi ff cfgi 2   2 ; uf2 ¼  u2 ¼  fm cmgi þff cfgi fm cmgi þff cfgi 1þ2

¼ qm

(B-3)

(A-9) Transfer coefficient:

After Eq. (A-9) is linearized as above, we obtain:

     i f m c vm pf kapp h f gi fi 2 mðpm2 Þ  m pf2 ¼ V m pf 2 þ 2a vta kf2 kf2

l1 ¼ a

km 2 L ; kf 1 ref

SRV region Seepage equations in the SRV region for the matrix and the natural fracture system can be linearized as: For the matrix system

2a

kapp kf1

 i f m c m gi mgi vmðpm1 Þ mðpm1 Þ  m pf1 ¼ vta kf1 fm mgi cd vmðpm1 Þ þ vta kf1

km 2 L kf2 ref

(B-4)

fm1 cd1 fm cmgi þ ff cfgi

 ;

ad2 ¼ 

1

fm2 cd2 fm cmgi þ ff cfgi



(B-5) 2

Mobility ratio:

(A-11)

. kf 1 mgi . ¼ kf 2 mgi

(B-6)

Permeability coefficient:

qkam ¼

kapp km

(B-7)

Dimensionless distances:

(A-12)

where q is the point source from the infinite-conductivity hydraulic fractures; moreover, flow from the reservoir into the horizontal well is only through the hydraulic fractures, and the production directly from the surface of the horizontal well is negligible. q can be expressed as:

2psc T q¼ qsc dðx  x0 Þdðy  y0 Þdðz  z0 Þ Tsc kf 1

ad1 ¼ 

M12

For the natural system

   i kapp h mðpm1 Þ  m pf1 þ q V2 m pf 1 þ 2a kf1   vm pf 1 ff1 mgi cfi ¼ vta kf1

l2 ¼ a

Adsorption coefficient:

(A-10)

h

2

¼ outerÞ

xD ¼

x y z ; yD ¼ ; zD ¼ Lref Lref Lref

(B-8)

In particular, point sources from the infinite-conductivity hydraulic fracture can be defined as:

q¼

2psc TL2ref Tsc kf1 L3ref

¼ 2p (A-13)

pk hTsc qsci dðxD  x0D ÞdðyD  y0D ÞdðzD  z0D Þ   f2 jqsct jpsc T

qsci hD dðxD  x0D ÞdðyD  y0D ÞdðzD  z0D Þ jqsct j M12 (B-9)

where i is the point source or element number containing the hydraulic fractures of the horizontal well. Appendix B. Dimensionless terms For convenience, we derive the flow solution in terms of dimensionless variables, considering: Dimensionless pseudo-pressure:

mD ¼

pkf2 hTsc qsct psc T

Dmx ðx ¼ m; f Þ  

Dmf ¼ mðpi Þ  m pf ;

(B-1)

kf 2 fm cmgi þ ff cfgi



Storage coefficient:

1þ2

C.1 Discretization We choose a tetrahedral element as an example. According to the Galerkin method, we assume the following shape function or basic function:

Nm ¼

Dmm ¼ mðpi Þ  mðpm Þ

1 ðam þ bm x þ cm y þ dm zÞ; 6V

m ¼ i; j; k; l

(C-1)

And the natural coordinate expressions satisfy the following conditions:

Dimensionless time:

tD ¼ 

Appendix C. Numerical solution for pressure analysis model

mgi L2ref

ta

(B-2)

∭ La1 Lb2 Lc3 Ld4 dV ¼ V

where

a!b!c!d! 6V ða þ b þ c þ d þ 3Þ!

(C-2)

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

 1  1  1 V¼  61 1

xi xj xk xl

yi yj yk yl

2

 zi  zj  zk  zl 

The displacement function is:

mD ¼

l X

Nv mDv

(C-3)

6 6 6 6 6 6 6 6 B¼6 6 6 6 6 6 6 6 4

v¼i

where v is the number of nodes in a tetrahedral, and mDv is the dimensionless pressure at the node v. C.2 Matrix Formula and Solution We extract a grid element (e) for analysis, substituting Eq. (C-3) into Eq. (8) ~ Eq. (11), after transforming the pressure equation into weak form by integrating over the volume of the element; as a result, we obtain: For the matrix system:

    vmmDx dV ¼ 0 ∭ Nv Mx lx qkam mmDx  mfDx þ utx vt

(C-4)

D

Ve

For the natural fracture system:

 ∭ Nv Mx V2 mfDx þ Mx lx qkam



Ve

 vmfDx dV mmDx  mfDx  ux uf x vtD 

¼0 (C-5)

x ¼ SRV; outer

vNj vxD

vNj vyD

vNk vxD

vNk vyD

vNl vxD

vNl vyD

3 vNi vzD 7 7 7 7 32 3T 2 vNj 7 7 Ni Ni 7 vzD 7 76 7 6 7; Tm ¼ NNT ¼ 6 Nj 76 Nj 7 7 4 Nk 54 Nk 5 vNk 7 7 Nl Nl vzD 7 7 7 7 vNl 5 vzD

Finally, we achieve the whole domain pressure matrix by superimposing each element matrix according to Eq. (C-8):

Kmnþ1 ¼F fD

(C-8)

where m(nþ1) is the pressure value of each node for the natural fD fracture system in the next step; K is the stiffness matrix; and F is the load vector. However, F should be modified to add sources qi* when considering some elements that include source nodes from the hydraulic fractures. Because the flow capacity in hydraulic fractures and horizontal wellbores is much higher than in the formation, which is extremely obvious in shale reservoirs, we can assume that they have infinite conductivity. Therefore,

mnþ1 ¼ mnþ1 fD fDx¼1;2…N

(C-9)

And the dimensionless production from each source node qi* satisfies the condition: NF X

2phD q*i ¼ 1

(C-10)

1 M12

By combining Eq. (C-9) and Eq. (C-10) with the whole domain pressure matrix described in Eq. (C-8), we obtain an (N þ NF)  (N þ NF) linear equation system, which can solve the N þ NF unknowns of mfD1,mfD2, … mfDN,q1*, … qNF*. The matrix for this system can be expressed as

x ¼ outer x ¼ SRV

To solve the coupled equation of the matrix and fracture system, we first calculate the pressure of the natural fracture system at the next time step explicitly and then substitute the results into Eq. (C7) to obtain the next time step pressure of the matrix system. A simple form will be:

Ke mnþ1 ¼ Fe fD

(C-6)

b nþ1 c n mnþ1 wD ¼ amfD þ amwD

(C-7)

where

utx ;a ¼ b þc b ¼ Mx lx qkam ; c ¼ Dt 
bc ux uf x þ Tm Ke ¼ Mx ∭ BT BdVm þ a Dt V Fe ¼

vNi vyD

i¼1

where

Mx ¼

vNi vxD

597

ux uf x bc Tm mnmD þ T mn a Dt m fD

Km ¼ F

(C-11)

where

2

k11 6 $$ 6 6 kNP 1 6 6 kN1 K¼6 6 $$ 6 6 kðNþN Þ1 P 6 4 kðNþN þ1Þ1 P $$ 0 B B B B B m¼B B B B B @

mfD1 mfD2 « mfDN q*1 q*2 « q*NF

$ $ $ $ $ $ $ $

k1NP $$ kNP NP $$ $$ 1 $$ $$

1nþ1 C C C C C C C C C C A

$ $ $ $ $ $ 1 $ 0

;

B B B B B F¼B B B B @

k1N $$ $$ kNN $$ $$ $$ $$ F1 F2 « FN 0 0 « 1

$ $ $ $ $ $

k1ðNþNP Þ $$ 2phD $$ $$ $$ $$ 1

$ $ $ $ $ $ $ 1

3 k1ðNþNF Þ 7 $ 7 7 $ 7 7 $ 7 7 $ 7 7 $ 7 5 $ kðNþNF ÞðNþNF Þ

1 C C C C C C C C C A

We can obtain the well bore pressure through mfDw (the average pressure of well point grid) and the transient well index as Aguilar

598

R.-h. Zhang et al. / Journal of Natural Gas Science and Engineering 33 (2016) 587e598

et al. (2007). The discrete numerical Laplace transform method is used to convert the real well bore pressure into Laplace space (Xu et al., 2015a,b). Thus, the dimensionless total production rate of a MFH well with a constant pressure drop can be obtained by

qD ¼

1 s2 mwD

(C-12)

where

qD ¼

qsct psc T kf 2 hTsc Dm

References Aboaba, A., Cheng, Y., 2010. Estimation of fracture properties for a horizontal well with multiple hydraulic fractures in gas shale. In: Paper SPE 138524 Presented at the SPE Eastern Regional Meeting, Morgantown, West Virginia, USA. Aguilar, C., Ozkan, E., Kazemi, H., 2007. Transient behavior of multilateral wells in numerical models: a hybrid analytical-numerical approach. In: Paper SPE 104581 Presented at the 15th SPE Middle East Oil & Gas Show and Conference. Kingdom of Bahrain. Aguilera, R., 2010. Flow units: from conventional to tight gas to shale gas reservoirs. In: Paper SPE 132845 Presented at the SPE Trinidad and Tobago Energy Resources Conference, Trinidad, Spain. Al-Ahmadi, H.A., Almarzooq, A.M., Watenbargen, R.A., 2010. Application of linear flow analysis to shale gas wells-field cases. In: Paper SPE 130370 Presented at the SPE Unconventional Gas Conference, Pittsburgh, Pennsylvania, USA. Al-Mohannadi, N., Ozkan, E., Kazemi, H., 2004. Grid system requirements in numerical modeling of pressure-transient tests in horizontal wells. In: Paper SPE 92041 Presented at the SPE International Petroleum Conference. Peubla, Mexico. Brohi, I., Pooladi-Darvish, M., Aguilera, R., 2011. Modeling fractured horizontal wells as dual porosity composite reservoirs-application to tight gas, shale gas and tight oil cases. In: Paper SPE 144057 Presented at the SPE Western North American Region Meeting, Anchorage, Alaska, USA. Brunauer, S., Deming, L.S., Edwards, W., 1940. On a theory of the der Waals adsorption of gases. J. Am. Chem. Soc. 62 (7), 1723e1732. Bumb, A.C., McKee, C.R., 1988. Gas-well testing in the presence of desorption for coalbed methane and Devonian shale. SPE. Form. Eval. 3 (1), 179e185. Carlson, E.S., Mercer, J.C., 1991. Devonian shale gas production: mechanisms and simple models. J. Pet. Technol. 43 (4), 476e482. Chen, Z.X., Huan, G.R., Li, B.Y., 2003. Modeling 2D and 3D horizontal wells using CVFA. Commun. Math. Sci. 1 (1), 30e43. Dehghanpour, H., Shirdel, M., 2011. A triple porosity model for shale gas reservoirs. In: Paper SPE 149501 Presented at the Canadian Unconventional Resources Conference, Calgary, Alberta, Canada. Ertekin, T., King, G.R., Schwerer, F.C., 1986. Dynamic gas slippage: a unique dualmechanism approach to the flow of gas in tight formations. SPE Form. Eval. 1 (1), 43e52. Escobar, F.H., Martinez, J.A., Matilde, M.M., 2013. Pressure transient analysis for a reservoir with a finite-conductivity fault. CT&F - Ciencia. Tecnol. Futuro. 5 (2), 5e18. Fan, D.Y., 2013. Well Test Theory and Interpretation Method of Multi-fractured Horizontal Wells Based on the Discrete Fracture Model. China University of Petroleum, Dongying, China. Florence, F.A., Rushing, J.A., Newsham, K.E., et al., 2007. Improved permeability prediction relations for low-permeability sands. In: Paper SPE 107954 Presented at SPE Rocky Mountain Oil&Gas Technology Symposium, Denver, Colorado, USA. Gao, C., Lee, W.J., Spivey, J.P., et al., 1994. Modeling multilayer gas reservoirs including sorption effects. In: Paper SPE 29173 Presented at the SPE Eastern Regional Meeting. Charleston, West Virginia. Guo, C., Xu, J., Wu, K., et al., 2015. Study on gas flow through nano pores of shale gas reservoirs. J. Fuel 143, 107e117. Garipov, T.T., Karimi-Fard, M., Tchelepi, H.A., 2016. Discrete fracture model for coupled flow and geomechanics. Comput. Geosci. 1e12. Hatman, C., 2009. Shale Gas Core Analyses Required for Gas Reserve Estimates. Hill, David G., Nelson, C.R., 2000. Gas productive fractured shales: an overview and update. Gas. TIPS 6, 4e13. Javadpour, F., 2009. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). J. Can. Petrol. Technol. 48 (8), 16e21. Jiang, J.M., Younis, R.M., 2015. A generic physics-based numerical platform with hybrid fracture modelling techniques for simulating unconventional gas

reservoirs. In: Paper SPE 173318 Presented at the SPE Simulation Symposium, Houston, Texas, USA. Klinkenberg, L.J., 1941. The permeability of porous media to liquids and gases. In: API Drilling and Production Practice. Kucuk, Fikri, Sawyer, Walter K., 1980. Transient flow in naturally fractured reservoirs and its application to Devonian gas shales. Paper SPE 9397 presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas. Kuuskraa, V.A., Wicks, D.E., Thruber, J.L., et al. Geologic and reservoir mechanisms controlling gas recovery from the Antrim shale. Paper SPE 24883 Presented at the SPE Annual Technical Conference and Exhibition, Washington D.C, USA. Langmuir, I., 1916. The constitution of fundamental properties of solids and liquids. J. Am. Chem. Soc. 38, 2221e2295. Li, D.M., Jiang, H.Q., L, J.J., 2014. The impact of diffusion type on multiscale discrete fracture model numerical simulation for shale gas. Nat. Gas. Sci. Eng. 20, 74e81. Mcclure, M.W., Horne, R.N., 2013. Discrete Fracture Network Modeling of Hydraulic Stimulation. Springer Berlin. Moinfar, A., Varavei, A., Sepehrnoori, K., et al., 2013. Development of a coupled dual continuum and discrete fracture model for the simulation of unconventional reservoirs. In: Paper SPE 163647 Presented at the SPE Reservoir Simulation Symposium, Woodlands, Texas, USA. Ozkan, E., Brown, M., Raghavan, R., et al., 2011. Comparison of fractured-horizontalwell performance in tight sand and shale reservoirs. SPE Reserv. Eval. Eng. 14 (2), 248e259. Roy, S., Raju, R., Chuang, H.F., et al., 2003. Modeling gas flow through microchannels and nanopores. J. Appl. Phys. 93 (8), 4870e4879. Rushing, J.A., Perego, A.D., Sullivan, R., et al., 2007. Estimating reserves in tight gas sands at HP/HT reservoir conditions: use and misuse of an Arps decline curve methodology. In: Paper SPE 109625 Presented at the SPE Annual Technical Conference and Exhibition. Anaheim, California, USA. Schepers, K.C., Gonzalez, R.J., Koperna, G.J., et al., 2009. Reservoir modeling in support of shale gas exploration. In: Paper SPE 123057 presented at the Latin American and Caribbean Petroleum Engineering Conference, Cartagena de Indias, Colombia. Shabro, V., Torres-verdín, C., Javadpour, F., 2011. Numerical simulation of shale-gas production: from pore-scale modeling of slip-flow Knudsen diffusion and Langmuir desorption to reservoir modeling of compressible fluid. In: Paper SPE 144355 Presented at the SPE North American Unconventional Gas Conference and Exhibition, Woodlands, Texas, USA. Stalgorova, E., Mattar, L., 2012. Practical analytical model to simulate production of horizontal wells with branch fractures. In: Paper SPE 162515 Presented at the SPE Canadian Unconventional Resources Conference, Calgray, Alberta, Canada. Sun, H., 2014. Multi-scale Simulation Theory and Method of Gas Transport in Shale Gas Reservoirs. China University of Petroleum, Dongying, China. Swami, V., 2012. Shale gas reservoir modeling: from nanopores to laboratory. In: Paper SPE 163065 Presented at the SPE International Student Paper Contest at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA. Swami, V., Settari, A.T., Javadpour, F., 2013. A numerical model for multi-mechanism flow in shale gas reservoirs with application to laboratory scale testing. In: Paper SPE 164840 Presented at the EAGE Annual Conference & Exhibition Incorporating SPE Europec, London, United Kingdom. Wang, F.P., Reed, R.M., John, A., et al., 2009. Pore networks and fluid flow in gas shales. In: Paper SPE 124253 Presented at the 2009 SPE Annual Technical Conference and Exhibition. New Orleans, Louisiana, USA. Wang, H.T., 2013. Performance of multiple fractured horizontal wells in shale gas reservoirs with consideration of multiple mechanisms. J. Hydrol. 510, 299e312. Xu, J.C., Guo, C.H., Wei, M.Z., et al., 2015a. Production performance analysis for composite shale gas reservoir considering multiple transport mechanisms. Nat. Gas. Sci. Eng. 26, 382e395. Xu, J.C., Jiang, R.Z., Teng, W.C., 2015b. Nonlinear flow characteristics and horizontal well pressure transient analysis for low-permeability offshore reservoirs. Mat. Pro. Eng. 1e13. Yan, B.C., Wang, Y., Killough, J.E., 2015. Beyond dual-porosity modeling for the simulation of complex flow mechanisms in shale reservoirs. Comput. Geosci. 1e23. Yu, W., Zhang, T., Du, S., et al., 2014. Numerical study of the effect of uneven proppant distribution between multiple fractures on shale gas well performance. Fuel 142, 189e198. Zhang, D.L., Zhang, L.H., Guo, J.J., et al., 2015a. Research on the production performance of multistage fractured horizontal well in shale gas reservoir. Nat. Gas. Sci. Eng. 26, 279e289. Zhang, D.L., Zhang, L.H., Zhao, Y.L., et al., 2015b. A composite model to analyze the decline performance of a multiple fractured horizontal well in shale reservoirs. Nat. Gas. Sci. Eng. 26, 999e1010. Zhang, L.H., Li, J.C., Tang, H.M., et al., 2014. Fractal pore structure model and multilayer fractal adsorption in shale. Fractals-complex Geomet. Patterns Scaling Nat. Soc. 22 (3), 321e329. Zhao, Y.L., 2015. Research on Transient Seepage Theory of Fractured Wells with Complex Percolation Mechanisms in Multiscale Shale Gas Reservoir. Southwest Petroleum University, Chengdu, China.