# Production Rate Analysis of Fractured Horizontal Well considering Multitransport Mechanisms in Shale Gas Reservoir

Authors: Qi-guo Liu; Wei-hong Wang; Hua Liu; Guangdong Zhang; Long-xin Li; Yu-long Zhao
Journal: Geofluids
 (2018)

---

## Abstract

Shale gas reservoir has been aggressively exploited around the world, which has complex pore structure with multiple transport mechanisms according to the reservoir characteristics. In this paper, a new comprehensive mathematical model is established to analyze the production performance of multiple fractured horizontal well (MFHW) in box-shaped shale gas reservoir considering multiscaled flow mechanisms (ad/desorption and Fick diffusion). In the model, the adsorbed gas is assumed not directly diffused into the natural macrofractures but into the macropores of matrix first and then flows into the natural fractures. The ad/desorption phenomenon of shale gas on the matrix particles is described by a combination of the Langmuir’s isothermal adsorption equation, continuity equation, gas state equation, and the motion equation in matrix system. On the basis of the Green’s function theory, the point source solution is derived under the assumption that gas flow from macropores into natural fractures follows transient interporosity and absorbed gas diffused into macropores from nanopores follows unsteady-state diffusion. The production rate expression of a MFHW producing at constant bottomhole pressure is obtained by using Duhamel’s principle. Moreover, the curves of well production rate and cumulative productionvs. time are plotted by Stehfest numerical inversion algorithm and also the effects of influential factors on well production performance are analyzed. The results derived in this paper have significance to the guidance of shale gas reservoir development.


---

## Body

## 1. Introduction

Shale gas is an unconventional natural gas with the main component of methane existing as adsorbed or free gas in the rich organic shale or its interlayers, the reserve of which is tremendous around the world. With the innovative techniques of horizontal well drilling as well as hydraulic fracturing in North American or other places of the world, the annual production of shale gas in the USA has increased from 564 × 108 m3 in 2007 to 4382 × 108 m3 in 2015. Moreover, well testing and production performance analysis in shale gas reservoir, especially for multiple fractured horizontal well (MFHW), is a hot research topic in recent years [1]. Compared to conventional gas reservoirs, shale gas reservoirs have special features on the accumulation conditions, occurrence manner, percolation mechanism, etc. For instance, shale is a self-sourcing reservoir within which the gas is stored either by compression in the pores and natural fracture as free gas or by adsorption on the surface of the solid material (either organic material or minerals) as adsorbed gas. According to Hill and Nelson [2], the adsorbed gas accounts for up to 85% of total gas, which is dependent on a variety of geologic and geochemical properties [3]. The pore network of shale gas reservoir is complex, which can be further divided into 4 categories: nanoinorganic pore, nanoorganic intragranular pore, macronatural fracture, and large-scaled hydraulic fractures. The first two pore types belong to matrix pores and the latter two belong to fractures [4–7], so the mechanisms of gas flow in shale reservoirs is very complicated. Shale gas reservoirs are generally described by the dual porosity model like Warren-Root model [8], in which the matrix blocks serve as the main storage space of gas and are separated by fractures.Kucuk and Sawyer [9] first analyzed the transient pressure response of wells in shale gas reservoirs by using the dual porosity model (assuming cylindrical and spherical matrix geometries), and the corresponding analytical and numerical results were obtained. However, the analytical model presented in the paper did not take the effects of gas desorption and diffusion into account and the numerical model ignored the effects of gas diffusion.On the basis of isothermal adsorption equation, Bumb and McKee [10] analyzed the transient pressure of a vertical well in coal bed methane (CBM) reservoir by introducing an additional compressibility coefficient in their model. Their method has been popularly used for production performance analysis of unconventional gas. Many researchers followed Bumb and McKee’s [10] approach to consider the ad/desorption effects in their studies; the results of which were proved to be much more realistic than previous studies [11–14].Ertekin et al. [15] utilized a dual-mechanism model for tight gas reservoirs, which considered Darcy flow (driven by pressure gradient) and Fick diffusion (driven by concentration difference). Subsequently, Carlson and Mercer [16] investigated the gas flow behavior in shale gas reservoir by coupling the gas diffusion and desorption with Fick’s law and Langmuir isothermal adsorption equation. And then, using the same method, Chawathe et al. [17] analyzed the pressure and saturation distribution in the reservoir by numerical solution.Recently, many studies related to well pressure and rate transient analysis have been presented [5, 18–33]. El-Banbi [34] utilized a linear flow model to analyze the production performance of a MFHW by ignoring the gas desorption and diffusion. In this model, the reservoir was simplified by continuous medium with slab dual porosity and Warren-Root model [19, 20, 35]. However, most of these studies did not consider desorption or diffusion effects, such as “trilinear” model [24] and other models [36, 37].In fact, a large number of core electron macroscope scanning images show that there are many organic pores existing in the shale matrix, the scale of which ranges from nanometer to macrometer. Previous studies considered the matrix pores as a whole without separating nanometer pores and macropores. Song [38, 39] proposed a “triple porosity” conceptual model to describe the gas storage mechanism, which was equivalent to a dual porosity model combined with gas adsorption. Thereafter, Zhao et al. [40] adopted the model to analyze the pressure transient and rate decline behavior of MFHW in shale gas reservoirs. Recently, a comprehensive mathematical model of a MFHW in a box-shaped reservoir is presented by Zhao et al. [32]. In this model, the reservoir was divided into a dual porosity system, which included natural fractures and nanopores. The gas flow in natural fracture system obeyed the Darcy flow theory, while the Knudsen diffusion and slippage effect were considered for the matrix low permeability system.According to Dehghanpour and Shirdel [41], the dual porosity models have been widely used in naturally fractured shale gas reservoirs with multiple hydraulic fractures. In these models, only nanopores existed in the matrix blocks which were surrounded by a fracture network. Gas flow in the nanopores obeys the Fickian diffusion law and the storage capability of the fracture system is very low, the theoretical productivity of the well from such conceptual models is quite low, and the decline rate of well productivity is very large, which is not consistent with some practical data ([5, 24, 27]).Most physical models proposed by previous investigators did not take into account the nanopores and the macropores in the matrix system [32, 40, 42]. According to the NMR scanning results of the core from YY32 well from the YanChang formation in P.R. China (Figure 1), two peaks are the most important characteristics for many shale gas reservoirs, which means that there are mainly two types of pores in the matrix system: the nanopores and macropores.Figure 1
The pore size distribution map for shale gas reservoir in Ordos Basin (P.R. China).In this paper, a new conceptual model is developed to describe the gas flow in shale gas reservoir. The flow path is presented by into three pore types: natural fractures, macropores, and nanopores. It is assumed that the adsorbed gas first desorbs from the organic particle surface into the nanopore system. It then diffuses into the macropore system. Thereafter, the gas in the macropores is transported into the natural fractures with transient interporosity flow. Based on this physical model, the production performance of a MFHW in rectangular outer-boundary shale gas reservoirs is analyzed.


## 2. Multiple Transport Mechanisms of Shale Gas Reservoirs

### 2.1. Gas Ad/Desorption Mechanism

A large amount of shale gas is adsorbed on the surface of rock particles. As the reservoir pressure depletes, adsorbed gas desorbs from the matrix surface and becomes free gas. This feature differentiates shale gas reservoirs different from conventional reservoirs. The isothermal adsorption curves have the same shape as Langmuir isothermal curve [10], which is frequently used to describe ad/desorption process of shale gas.The occurrence time of gas desorption behavior also depends on the amount of adsorbed gas at initial conditions. As shown in Figure2, if the adsorbed gas is statured under initial conditions, the amount of adsorbed gas is at point B and the gas quickly desorbs from organic surfaces as pressure decreases. However, if the adsorbed gas is unsaturated, the amount of adsorbed gas would be at point C (if reservoir conditions are at point C and there is no free gas in matrix pores) and the gas would not desorb from organic surfaces until reservoir pressure decreases to the pressure at point A, which is called critical desorption pressure. The difference between the critical desorption pressure and the initial reservoir pressure determines the starting time of gas desorption [1].Figure 2
Shale gas isothermal adsorption curve and desorption process diagram.Langmuir isothermal adsorption equation is expressed as(1)Cm=CLpmpL+pm,where Gm is the gas volumetric concentration, m3/m3; GL is the Langmuir adsorption volume of shale gas, representing the limit adsorbence of unit reservoir, m3/m3; pL is the Langmuir pressure, which is the pressure when gas adsorbence reaches 50% of limit adsorbence, Pa; and pm is the pressure in matrix system, Pa.The mass flow rate of desorbed gas from shale reservoir of volumeVb in unit time is
(2)qdes=ρgscVb∂Cm∂t,where qdes is the mass of desorption gas in unit time, kg/s; t is the time, s; Vb is the shale matrix volume, m3; and ρgsc is the shale gas density at standard conditions, kg/m3.If gas desorption happens instantaneously and all desorbed gas transforms into free gas, then (2) can be introduced directly into the continuity equation of the matrix or fracture system, resulting in the steady-state adsorption-desorption model ([10–13]; Civan, 2010; [37]; Wang et al., 2013).


### 2.2. Gas Flow in Natural Fractures and Macropore System

Gas flow in natural and hydraulic fractures and macropores is continuous flow. According to previous analysis (Civan, 2010; [37]), Darcy’s flow equation used for conventional gas reservoir can be adopted to analyze gas flow behavior in fracture and macropore systems in shale reservoirs, that is,
(3)Jl=−plMgZRTklμg∇pl,l=f,mc,where f presents the natural fracture system and the mc presents the macropore system; Z is the gas deviation factor (dimensionless), Z = 1 for ideal gas; R is the gas constant (8.314 J/(mol ⋅ K)); T is the reservoir temperature, K; kl is the permeability in l system, m2; Mg is the gas molar mass, kg/mol; pl is the pressure in l system, Pa.


### 2.3. Gas Diffusion in Nanopores

Gas flow behavior in shale is different from flow in fractures because the pores are of nanoscale. In these nanopores, it is erroneous to assume molecular continuous Darcy flow. Shale gas flow in nanopores includes not only viscous flow but also diffusion. In some ultratight shale reservoirs, the flow of gas molecules occur by only diffusion because the gas in the nanopores exists only as absorbed gas and not free gas.Fick’s law of diffusion is generally used to describe shale gas migration in nanopores and at the surface of kerogen [15, 40, 43–46]. Under the influence of concentration gradient, desorbed gas migrates from high concentration area to low concentration area, and the diffusion stops when gas concentration becomes uniform (as illustrated in Figure 3). Assume that gas molecular diffusion satisfies Fick’s law, and then gas diffusion flux could be expressed as
(4)JF=−MgD⋅∇Cm,where D is the Fick diffusion coefficient, m2/s.Figure 3
Diagram of Fick diffusion.If gas concentration in matrix changed with coordinates, meaning nonequivalent concentration at any time in the matrix, unsteady state diffusion occurs, which can be described by Fick’s second law. Otherwise, the Fick’s first law is used when the concentration is assumed to be the same at all locations in a given time (pseudosteady-state diffusion). It is assumed in this paper that the desorption of the adsorbed gas desorbed from the nanopores into the macropores follows unsteady-state diffusion [15].To simplify the hypothetical shale gas reservoir model, it is assumed that the matrix element is spherical with radiusRm; each matrix element is composed of organics and clay grains, and diffusion is the only gas flow mechanism in the matrix (Figure 4). Therefore, Fick’s second law can be applied to describe gas diffusion from matrix elements into macrofractures or macropores.Figure 4
Conceptual model used in this paper.According to the unsteady-state diffusion theory, gas concentration in the matrixCm is a function of time and space. Therefore, the following mathematical expression can be used to describe gas concentration change in the matrix,
(5)∂Cm∂t=1rm2∂∂rmDrm2∂Cm∂rm.The initial condition of matrix unsteady-state diffusion should satisfy the equation below.(6)Cmt=0,rm=Cipi,where pi is the initial reservoir pressure, Pa.The inner boundary condition can be expressed as(7)∂Cmt,rm=0∂rm=0.Since the outer boundary of matrix is connected with the fracture system or macropore system, pressure at the matrix surface equals that of fracture system, then, the outer boundary can be represented as(8)Cmt,rm=Rm=Cmp,where Rm is the radius of a spherical matrix element, m.Using the above model and relevant boundary conditions, the distribution of gas concentration in the shale matrix can be solved. Then, the following equation can be used to determine gas diffusion flow rate within matrix of volumeVb,
(9)qF=MgVb3DRm∂Cm∂rmrm=Rm.


## 2.1. Gas Ad/Desorption Mechanism

A large amount of shale gas is adsorbed on the surface of rock particles. As the reservoir pressure depletes, adsorbed gas desorbs from the matrix surface and becomes free gas. This feature differentiates shale gas reservoirs different from conventional reservoirs. The isothermal adsorption curves have the same shape as Langmuir isothermal curve [10], which is frequently used to describe ad/desorption process of shale gas.The occurrence time of gas desorption behavior also depends on the amount of adsorbed gas at initial conditions. As shown in Figure2, if the adsorbed gas is statured under initial conditions, the amount of adsorbed gas is at point B and the gas quickly desorbs from organic surfaces as pressure decreases. However, if the adsorbed gas is unsaturated, the amount of adsorbed gas would be at point C (if reservoir conditions are at point C and there is no free gas in matrix pores) and the gas would not desorb from organic surfaces until reservoir pressure decreases to the pressure at point A, which is called critical desorption pressure. The difference between the critical desorption pressure and the initial reservoir pressure determines the starting time of gas desorption [1].Figure 2
Shale gas isothermal adsorption curve and desorption process diagram.Langmuir isothermal adsorption equation is expressed as(1)Cm=CLpmpL+pm,where Gm is the gas volumetric concentration, m3/m3; GL is the Langmuir adsorption volume of shale gas, representing the limit adsorbence of unit reservoir, m3/m3; pL is the Langmuir pressure, which is the pressure when gas adsorbence reaches 50% of limit adsorbence, Pa; and pm is the pressure in matrix system, Pa.The mass flow rate of desorbed gas from shale reservoir of volumeVb in unit time is
(2)qdes=ρgscVb∂Cm∂t,where qdes is the mass of desorption gas in unit time, kg/s; t is the time, s; Vb is the shale matrix volume, m3; and ρgsc is the shale gas density at standard conditions, kg/m3.If gas desorption happens instantaneously and all desorbed gas transforms into free gas, then (2) can be introduced directly into the continuity equation of the matrix or fracture system, resulting in the steady-state adsorption-desorption model ([10–13]; Civan, 2010; [37]; Wang et al., 2013).


## 2.2. Gas Flow in Natural Fractures and Macropore System

Gas flow in natural and hydraulic fractures and macropores is continuous flow. According to previous analysis (Civan, 2010; [37]), Darcy’s flow equation used for conventional gas reservoir can be adopted to analyze gas flow behavior in fracture and macropore systems in shale reservoirs, that is,
(3)Jl=−plMgZRTklμg∇pl,l=f,mc,where f presents the natural fracture system and the mc presents the macropore system; Z is the gas deviation factor (dimensionless), Z = 1 for ideal gas; R is the gas constant (8.314 J/(mol ⋅ K)); T is the reservoir temperature, K; kl is the permeability in l system, m2; Mg is the gas molar mass, kg/mol; pl is the pressure in l system, Pa.


## 2.3. Gas Diffusion in Nanopores

Gas flow behavior in shale is different from flow in fractures because the pores are of nanoscale. In these nanopores, it is erroneous to assume molecular continuous Darcy flow. Shale gas flow in nanopores includes not only viscous flow but also diffusion. In some ultratight shale reservoirs, the flow of gas molecules occur by only diffusion because the gas in the nanopores exists only as absorbed gas and not free gas.Fick’s law of diffusion is generally used to describe shale gas migration in nanopores and at the surface of kerogen [15, 40, 43–46]. Under the influence of concentration gradient, desorbed gas migrates from high concentration area to low concentration area, and the diffusion stops when gas concentration becomes uniform (as illustrated in Figure 3). Assume that gas molecular diffusion satisfies Fick’s law, and then gas diffusion flux could be expressed as
(4)JF=−MgD⋅∇Cm,where D is the Fick diffusion coefficient, m2/s.Figure 3
Diagram of Fick diffusion.If gas concentration in matrix changed with coordinates, meaning nonequivalent concentration at any time in the matrix, unsteady state diffusion occurs, which can be described by Fick’s second law. Otherwise, the Fick’s first law is used when the concentration is assumed to be the same at all locations in a given time (pseudosteady-state diffusion). It is assumed in this paper that the desorption of the adsorbed gas desorbed from the nanopores into the macropores follows unsteady-state diffusion [15].To simplify the hypothetical shale gas reservoir model, it is assumed that the matrix element is spherical with radiusRm; each matrix element is composed of organics and clay grains, and diffusion is the only gas flow mechanism in the matrix (Figure 4). Therefore, Fick’s second law can be applied to describe gas diffusion from matrix elements into macrofractures or macropores.Figure 4
Conceptual model used in this paper.According to the unsteady-state diffusion theory, gas concentration in the matrixCm is a function of time and space. Therefore, the following mathematical expression can be used to describe gas concentration change in the matrix,
(5)∂Cm∂t=1rm2∂∂rmDrm2∂Cm∂rm.The initial condition of matrix unsteady-state diffusion should satisfy the equation below.(6)Cmt=0,rm=Cipi,where pi is the initial reservoir pressure, Pa.The inner boundary condition can be expressed as(7)∂Cmt,rm=0∂rm=0.Since the outer boundary of matrix is connected with the fracture system or macropore system, pressure at the matrix surface equals that of fracture system, then, the outer boundary can be represented as(8)Cmt,rm=Rm=Cmp,where Rm is the radius of a spherical matrix element, m.Using the above model and relevant boundary conditions, the distribution of gas concentration in the shale matrix can be solved. Then, the following equation can be used to determine gas diffusion flow rate within matrix of volumeVb,
(9)qF=MgVb3DRm∂Cm∂rmrm=Rm.


## 3. Solutions of a MFHW in a Rectangular Shale Gas Reservoir

### 3.1. Physical Model

In this section, continuous point source solution for any point in a rectangular shale gas reservoir with closed upper and lower boundaries is derived. The physical model is shown in Figure5. The assumptions for this model are (1) the gas reservoir is homogeneous and anisotropic with horizontal and vertical permeability kh and kz. The width, length, and height of the reservoir are xe, ye, and h; (2) a MFHW with length L is located in the middle of the reservoir and parallel to the upper and lower boundaries. All hydraulic fractures are symmetrical and perpendicular to the wellbore; the half length of the jth fracture is xfi; the fractures are fully penetrating; (3) the well is perforated at the intersections between the fractures and the wellbore, which means that fluid flows from the reservoir to the fractures and then from the fractures to the wellbore. Infinite conductivity fractures are assumed with no pressure drop across the fracture face and along the horizontal wellbore length; (4) each fracture is evenly divided into 2N units with unified flow rate [32, 47].Figure 5
Physical model of a MFHW in a rectangular closed outer-boundary shale gas reservoir.


### 3.2. Continuous Point Source Solution

During the development of shale gas reservoirs, natural fractures are the main paths for gas flowing from the matrix to hydraulic fractures and wellbores (Figure6). Therefore, flow models for these natural fractures are the key to deriving the point source function. The PDE of pseudopressure in cylindrical coordinates of fracture systems for different models are derived in Laplace domain (detailed derivation in the Appendix), which is
(10)1rD∂∂rDrD∂Δm¯f∂rD+∂2Δm¯f∂zD2=fsΔm¯f,where fs is the parameter group of the corresponding model under different percolation mechanism, dimensionless.Figure 6
Diagram of matrix desorbed gas flow through macropores to macrofractures.For unsteady-state diffusion and transient interporosity flow, the expression offs is
(11)fs=ωfs+λmf5gscothgs−1,where gs can be expressed as
(12)gs=15ωmλmf+βλσλθmfscothλθmfs−1s,β=1−ωf−ωm1−ϕf−ϕm.According to the mirror image and superposition principles (the schematic is shown in Figure7) and Newman [48] as well as Ozkan and Raghavan [49], the instantaneous point source solution at the observation point MD′ generated by unit source sink at point MD is [1]
(13)G¯MD,MD′,s=∑k=−∞+∞∑m=−∞+∞∑n=−∞+∞S1,1,1+S2,1,1+S1,2,1+S1,1,2+S2,2,1+S1,2,2+S2,1,2+S2,2,2.Figure 7
Schematic of a point source in rectangular gas reservoir with closed boundaries.In the above equation,Si,j,k represents the pressure response at any point of the gas reservoir caused by each image well due to closed boundary.
(14)Si,j,k=exp−fsx~Di−2kxeD2+y~Dj−2myeD2+z~Dk−2nhD2x~Di−2kxeD2+y~Dj−2myeD2+z~Dk−2nhD2i,j,k=1,2,where
(15)x~D1=xD−xwD,x~D2=xD+xwD,y~D1=yD−ywD,y~D2=yD+ywD,z~D1=zD−zwD,z~D2=zD+zwD.Introducing the following Poisson summation formula, which is [32, 49]
(16)∑n=−∞+∞exp−va2+ξ−2nξe2a2+ξ−2nξe2=1ξeK0av+2∑n=1∞K0av+n2π2ξe2cosnπξξe.Then, the continuous point source solution for rectangular gas reservoirs can be simplified and expressed as [32, 49]
(17)G¯MD,MD′,s=2πhDxeDcoshuy~D1+coshuy~D2usinhuyeD+2∑k=1+∞coskπxDxeDcoskπxwDxeDcoshεky~D1+coshεky~D2εksinhεkyeD+2∑n=1∞cosnπzDhDcosnπzwDhD×coshεnyeD−y~DεnsinhεnyeD+2∑k=1+∞coskπxDxeDcoskπxwDxeDcoshεk,ny~D1+coshεk,ny~D2εk,nsinhεk,nyeD,where
(18)y~D1=yeD−yD−ywDand
(19)y~D2=yeD−yD+ywD.Then, the pressure drop in rectangular reservoir generated by the intensityqscins of the continuous point source is
(20)Δm¯f‐ps=pscTTscqscins2πkfLrefsG¯MD,MD′,s.Introducing (17) into (20) yields
(21)Δm¯f‐ps=pscTTscqscinsπkfLrefhDsπxeDcoshuy~D1coshuy~D2usinhuyeD+2∑k=1+∞coskπxDxeDcoskπxwDxeD×coshεky~D1coshεky~D2εksinhεkyeD+2∑n=1∞cosnπzDhDcosnπzwDhDcoshεnyeD−y~DεnsinhεnyeD+2∑k=1+∞coskπxDxeDcoskπxwDxeDcoshεk,ny~D1+coshεk,ny~D2εk,nsinhεk,nyeD.The above expression is the pressure drop at any point of the reservoir generated by constant intensityqscins in rectangular reservoir. Noted that it can clearly see from the above equation that some infinite summations of eigenfunctions are existed, which generally have slow convergence characteristics. Since the solutions presented in this work are in the Laplace transform domain, we utilized the same approach proposed by Ozkan and Raghavan [49] to improve the computation efficient of source function given in (21). The main routine to overcome this problem is to separate in the infinite-acting and boundary-domained flow contributions in the solution and then recast the slow converging summations and computationally difficult integrals into computationally more efficient forms. In the previous paper [32], the detailed expressions during different time are given, and many authors has also reported them [49, 50].


### 3.3. Solution of the MFHW Producing with Constant Bottomhole Pressure

As is shown in Figure5, a horizontal well with M fully penetrated fractures is drilled in a rectangular outer-boundary reservoir. According to the source function theory, the pressure drop at arbitrary point caused by each discrete element can be obtained by integrating the continuous point source function along the fracture element, which is [40, 51]
(22)Δm¯fixD,yD=∫ΓΔm¯f‐psxD,yD,xwDi,ywDi,sdΓ.Due to the fractures are parallel to thex-axis and all of them are fully penetrated, the above equation can be converted into the following formulation:
(23)Δm¯fixD,yD=∫xmi−ΔLfi/2xmi+ΔLfi/2∫0hΔm¯f‐psxD,yD,xwDi,ywDi,sdzwdxw.Taking (21) into (23) and arranging it yield
(24)Δm¯fi=pscTTscqscLiπkfLrefhDsπxeDΔLficoshuy~D1+coshuy~D2usinhuyeD+2Lref∫−ΔLfDi/2+ΔLfDi/2∑k=1+∞coskπxDxeDcoskπxmDi+αxeDcoshεky~D1+coshεky~D2εksinhεkyeDdα,where qscLi is the continuous line source strength, which has the following relationship with the continuous point source qscins:
(25)qscLi=qscinsh.Defining the following dimensionless production rate and pseudopressure, we have(26)qDi=qscLiΔLfiqscand
(27)mfD=πkfhhTscpscTqscΔmf.Introducing (26) and (27) into (24), the dimensionless pseudopressure drop generated by the ith unit at any reservoir location can be expressed as
(28)m¯fDi=qDiAji,Aji=πsxeDcoshuy~D1+coshuy~D2usinhuyeD+2ΔLfDi∫−ΔLfDi/2ΔLfDi/2∑k=1+∞coskπxDxeDcoskπxmDi+αxeDcoshεky~D1+coshεky~D2εksinhεkyeDdα.Due to the horizontal well is composed of many discrete fracture units, the pressure at any point in the reservoir should be the superposition of the pressure drop caused by all these units at the same point. For the observation pointxDj,yDj, according to the superposition principle, there is [32, 40]
(29)m¯DxDj,yDj=∑i=1M∗2NqDiAjixDi,yDi,xDj,yDj,s.If setting the observation point on each fracture element, then we can obtain2N∗M linear equations and another auxiliary equation of mass conservation for flux can be also obtained; the detailed expressions can be referred to the paper of Zhao et al. [32].When well producing at a constant bottomhole pressure, the dimensionless rate can be defined as follows:(30)qD=qscpscTπkfhTscmfi−mwf.According to previous studies [52], there is the following relation formula between dimensionless pseudopressure and rate,
(31)q¯D=1m¯wD×s2,where the definition of dimensionless pseudopressure in (31) is
(32)mwD=πkfhTscpscTqscΔmf.The above equation is the production rate solution when well producing at a constant pressure.


## 3.1. Physical Model

In this section, continuous point source solution for any point in a rectangular shale gas reservoir with closed upper and lower boundaries is derived. The physical model is shown in Figure5. The assumptions for this model are (1) the gas reservoir is homogeneous and anisotropic with horizontal and vertical permeability kh and kz. The width, length, and height of the reservoir are xe, ye, and h; (2) a MFHW with length L is located in the middle of the reservoir and parallel to the upper and lower boundaries. All hydraulic fractures are symmetrical and perpendicular to the wellbore; the half length of the jth fracture is xfi; the fractures are fully penetrating; (3) the well is perforated at the intersections between the fractures and the wellbore, which means that fluid flows from the reservoir to the fractures and then from the fractures to the wellbore. Infinite conductivity fractures are assumed with no pressure drop across the fracture face and along the horizontal wellbore length; (4) each fracture is evenly divided into 2N units with unified flow rate [32, 47].Figure 5
Physical model of a MFHW in a rectangular closed outer-boundary shale gas reservoir.


## 3.2. Continuous Point Source Solution

During the development of shale gas reservoirs, natural fractures are the main paths for gas flowing from the matrix to hydraulic fractures and wellbores (Figure6). Therefore, flow models for these natural fractures are the key to deriving the point source function. The PDE of pseudopressure in cylindrical coordinates of fracture systems for different models are derived in Laplace domain (detailed derivation in the Appendix), which is
(10)1rD∂∂rDrD∂Δm¯f∂rD+∂2Δm¯f∂zD2=fsΔm¯f,where fs is the parameter group of the corresponding model under different percolation mechanism, dimensionless.Figure 6
Diagram of matrix desorbed gas flow through macropores to macrofractures.For unsteady-state diffusion and transient interporosity flow, the expression offs is
(11)fs=ωfs+λmf5gscothgs−1,where gs can be expressed as
(12)gs=15ωmλmf+βλσλθmfscothλθmfs−1s,β=1−ωf−ωm1−ϕf−ϕm.According to the mirror image and superposition principles (the schematic is shown in Figure7) and Newman [48] as well as Ozkan and Raghavan [49], the instantaneous point source solution at the observation point MD′ generated by unit source sink at point MD is [1]
(13)G¯MD,MD′,s=∑k=−∞+∞∑m=−∞+∞∑n=−∞+∞S1,1,1+S2,1,1+S1,2,1+S1,1,2+S2,2,1+S1,2,2+S2,1,2+S2,2,2.Figure 7
Schematic of a point source in rectangular gas reservoir with closed boundaries.In the above equation,Si,j,k represents the pressure response at any point of the gas reservoir caused by each image well due to closed boundary.
(14)Si,j,k=exp−fsx~Di−2kxeD2+y~Dj−2myeD2+z~Dk−2nhD2x~Di−2kxeD2+y~Dj−2myeD2+z~Dk−2nhD2i,j,k=1,2,where
(15)x~D1=xD−xwD,x~D2=xD+xwD,y~D1=yD−ywD,y~D2=yD+ywD,z~D1=zD−zwD,z~D2=zD+zwD.Introducing the following Poisson summation formula, which is [32, 49]
(16)∑n=−∞+∞exp−va2+ξ−2nξe2a2+ξ−2nξe2=1ξeK0av+2∑n=1∞K0av+n2π2ξe2cosnπξξe.Then, the continuous point source solution for rectangular gas reservoirs can be simplified and expressed as [32, 49]
(17)G¯MD,MD′,s=2πhDxeDcoshuy~D1+coshuy~D2usinhuyeD+2∑k=1+∞coskπxDxeDcoskπxwDxeDcoshεky~D1+coshεky~D2εksinhεkyeD+2∑n=1∞cosnπzDhDcosnπzwDhD×coshεnyeD−y~DεnsinhεnyeD+2∑k=1+∞coskπxDxeDcoskπxwDxeDcoshεk,ny~D1+coshεk,ny~D2εk,nsinhεk,nyeD,where
(18)y~D1=yeD−yD−ywDand
(19)y~D2=yeD−yD+ywD.Then, the pressure drop in rectangular reservoir generated by the intensityqscins of the continuous point source is
(20)Δm¯f‐ps=pscTTscqscins2πkfLrefsG¯MD,MD′,s.Introducing (17) into (20) yields
(21)Δm¯f‐ps=pscTTscqscinsπkfLrefhDsπxeDcoshuy~D1coshuy~D2usinhuyeD+2∑k=1+∞coskπxDxeDcoskπxwDxeD×coshεky~D1coshεky~D2εksinhεkyeD+2∑n=1∞cosnπzDhDcosnπzwDhDcoshεnyeD−y~DεnsinhεnyeD+2∑k=1+∞coskπxDxeDcoskπxwDxeDcoshεk,ny~D1+coshεk,ny~D2εk,nsinhεk,nyeD.The above expression is the pressure drop at any point of the reservoir generated by constant intensityqscins in rectangular reservoir. Noted that it can clearly see from the above equation that some infinite summations of eigenfunctions are existed, which generally have slow convergence characteristics. Since the solutions presented in this work are in the Laplace transform domain, we utilized the same approach proposed by Ozkan and Raghavan [49] to improve the computation efficient of source function given in (21). The main routine to overcome this problem is to separate in the infinite-acting and boundary-domained flow contributions in the solution and then recast the slow converging summations and computationally difficult integrals into computationally more efficient forms. In the previous paper [32], the detailed expressions during different time are given, and many authors has also reported them [49, 50].


## 3.3. Solution of the MFHW Producing with Constant Bottomhole Pressure

As is shown in Figure5, a horizontal well with M fully penetrated fractures is drilled in a rectangular outer-boundary reservoir. According to the source function theory, the pressure drop at arbitrary point caused by each discrete element can be obtained by integrating the continuous point source function along the fracture element, which is [40, 51]
(22)Δm¯fixD,yD=∫ΓΔm¯f‐psxD,yD,xwDi,ywDi,sdΓ.Due to the fractures are parallel to thex-axis and all of them are fully penetrated, the above equation can be converted into the following formulation:
(23)Δm¯fixD,yD=∫xmi−ΔLfi/2xmi+ΔLfi/2∫0hΔm¯f‐psxD,yD,xwDi,ywDi,sdzwdxw.Taking (21) into (23) and arranging it yield
(24)Δm¯fi=pscTTscqscLiπkfLrefhDsπxeDΔLficoshuy~D1+coshuy~D2usinhuyeD+2Lref∫−ΔLfDi/2+ΔLfDi/2∑k=1+∞coskπxDxeDcoskπxmDi+αxeDcoshεky~D1+coshεky~D2εksinhεkyeDdα,where qscLi is the continuous line source strength, which has the following relationship with the continuous point source qscins:
(25)qscLi=qscinsh.Defining the following dimensionless production rate and pseudopressure, we have(26)qDi=qscLiΔLfiqscand
(27)mfD=πkfhhTscpscTqscΔmf.Introducing (26) and (27) into (24), the dimensionless pseudopressure drop generated by the ith unit at any reservoir location can be expressed as
(28)m¯fDi=qDiAji,Aji=πsxeDcoshuy~D1+coshuy~D2usinhuyeD+2ΔLfDi∫−ΔLfDi/2ΔLfDi/2∑k=1+∞coskπxDxeDcoskπxmDi+αxeDcoshεky~D1+coshεky~D2εksinhεkyeDdα.Due to the horizontal well is composed of many discrete fracture units, the pressure at any point in the reservoir should be the superposition of the pressure drop caused by all these units at the same point. For the observation pointxDj,yDj, according to the superposition principle, there is [32, 40]
(29)m¯DxDj,yDj=∑i=1M∗2NqDiAjixDi,yDi,xDj,yDj,s.If setting the observation point on each fracture element, then we can obtain2N∗M linear equations and another auxiliary equation of mass conservation for flux can be also obtained; the detailed expressions can be referred to the paper of Zhao et al. [32].When well producing at a constant bottomhole pressure, the dimensionless rate can be defined as follows:(30)qD=qscpscTπkfhTscmfi−mwf.According to previous studies [52], there is the following relation formula between dimensionless pseudopressure and rate,
(31)q¯D=1m¯wD×s2,where the definition of dimensionless pseudopressure in (31) is
(32)mwD=πkfhTscpscTqscΔmf.The above equation is the production rate solution when well producing at a constant pressure.


## 4. Results Analysis

### 4.1. Model Validation

In this section, we use the solution of conventional dual-porosity reservoirs to prove the correctness of our model. Comparing to the conventional model, the difference between them is that our mode considers the gas desorption and diffusion from the nanopores into macropores. The model presented in this paper can simplify into the conventional dual-porosity model by neglecting the shale gas desorption and diffusion processes. For conventional reservoirs, approximation expressions forfs in (11) can be found in the literatures [8, 9, 53–55]. The expression of fs introduced by Warren and Root for blocked matrix with pseudosteady-state interporosity flow is
(33)fs=ωf1−ωfs+λmf1−ωfs+λmfs.Using de Swaan’s approach [53], the solutions for Warren and Root’s model can be extended to the solutions of transient interporosity flow, which are
(34)fs=ωfs+sλmf1−ωf3stanh31−ωfsλmffor horizontal slab blocks and
(35)fs=ωfs+λmf5151−ωfsλmfcoth151−ωfsλmf−1for spherical blocks.The above three types of matrix block are the most popular idealization model for naturally fractured reservoirs; they have been widely used in the petroleum industry since they are proposed. So, the correctness of our model and method can be validated by comparing the result to the conventional model with spherical blocks.In our model, if the gas diffusion and desorption are neglected, the following relationship must be established, which is(36)ωf+ωm=1.Substituting the above equation into the expression ofgs in (11), which becomes
(37)gs=151‐ωfλmfs.It can be clearly seen that (11) can be simplified into the same solution for spherical blocks described by (35) for naturally fractured reservoirs by taking the above equation into it, which prove the correctness of our solution.


### 4.2. Well Production Performance Analysis

On the basis of the above solution, the semianalytical solution in Laplace domain can be transformed into real space using Stehfest numerical inversion algorithm [56]. Therefore, when the well is produced at a constant bottomhole pressure, the rate decline semilog type curves of production rate and cumulative production vs time can be obtained. The input parameters for simulation sensitivities for the following models are shown in Table 1.Table 1
Synthetic data used in the examples.
Parameters
Value
Parameters
Value
Initial reservoir pressure,pi, MPa
25
Reservoir temperature,T, K
320
Formation thickness,h, m
60
Gas specific gravity,rg, fraction
0.65
Well length,L, m
1200
Constant bottomhole pressure,pwf, MPa
4
Natural fracture system permeability,kf, mD
0.005
Natural fracture system porosity,∅f, fraction
0.01
Macropore system permeability,km, mD
10−6
Macropore system porosity,∅m, fraction
0.05
Langmuir pressure,PL, MPa
4
Langmuir volume,CL (m3/m3)
10
Number of fractures,M (integer)
10
Half fracture length,Lf, m
100
Outer boundary size,xe, m
800
Outer boundary size,ye, m
3000For a well producing at constant bottomhole pressure in a rectangular gas reservoir, the production rate and cumulative production vs. time curves in semilog coordinates with different fracture numbers are shown in Figure8. From the plot, it can be seen that the early time production rate is increasing with the increasing of fracture number, which is mainly because the larger the M is, the bigger the drainage area of the well will be during the early flow period. And when the pressure wave transports to the outer region formation, this effect will be gradually weaken. However, the cumulative production with different M under the same production time in later production period will be the same due to the reserve of them are equivalent.Figure 8
Influence of fracture number on well production performance.Figure9 shows the influence of Langmuir volume CL on the well production performance. Due to the mechanism, physical model in this paper has divided the pores in matrix system into the following two subspace: macropores and nanopores. The adsorbed gas will be desorbed from the particle surface in the matrix and then flow into the macropores by Fick diffusion. When the pressure in the macropore system is lower than the Langmuir pressure, desorption phenomenon will happen. So, when the well producing with the constant bottomhole pressure of 4 MPa, the pressure in the vicinity of the well will be lower than the Langmuir pressure, which will cause the adsorbed gas to desorb from matrix surface. So, we can find from Figure 9 that the production rate is larger with the increasing of Langmuir volume. But this effect is very small. When well is producing for a long time, the cumulative production of well with bigger CL will be larger than the others.Figure 9
Influence of Langmuir volume on well production performance.Figure10 shows the effects of permeability in natural fracture system on the well production performance. Because the natural fracture system is the main flow path of gas flow from the formation into hydraulic fractures, it has a significant influence on well production performance. According to the productivity equation of wells, the well production rate has a proportional relationship with the natural fracture permeability. And due to the constant control volume of our model, the bigger production rate at early flow period will lead to a bigger rate decline in the middle flow period. According to the cumulative production curves on Figure 10, we can clearly see that the larger the kf is, the larger the cumulative production will be in our given time period.Figure 10
Influence of natural fracture permeability on well production performance.Figure11 shows the effects of permeability in macropores on well production rate and cumulative production rate curves. The bigger km represents a stronger supplement of fluid from matrix system into natural fractures. Because the fluid flow from matrix system into natural fractures follows the transient interporosity flow, this supplement will happen immediately after the well opens. The bigger the km is, the larger the well production rate will be, and this trend will be lasting for a long time. So, for shale gas reservoir development, the permeability of the matrix system is a very important parameter to evaluate the well production performance.Figure 11
Influence of permeability in macropore system on well production performance.Figure12 shows the effects of macropores porosity (∅m) on well production performance. The value of ∅m presents the amount of free gas in macropore system of shale matrix system. With the increase of the ∅m, the quantity of free gas flows from macropores into natural fracture system by interporosity will increase. So, when well produces at a constant bottomhole pressure, the higher the porosity is, the larger the well production rate will be. Although the increment of production rate is not quite obvious, for the situation listed in this paper, the increment is still nearly 5000 m3/d. If this increment for conventional gas reservoir could be negligible, it is very considerable for shale gas reservoir due to the long production period.Figure 12
Influence of porosity in macropores on well production performance.Figure13 shows the effects of gas desorption time (τ) on the well production performance. As we know, the gas desorption time presents the gas desorption rate of adsorbed gas desorbed into free gas. The smaller the τ is, the faster the desorption process happens. It can be clearly seen from Figure 12 that the smaller the gas desorption time is, the larger the well production time will be, and the corresponding cumulative production will be larger too.Figure 13
Influence of desorption time on well production performance.


## 4.1. Model Validation

In this section, we use the solution of conventional dual-porosity reservoirs to prove the correctness of our model. Comparing to the conventional model, the difference between them is that our mode considers the gas desorption and diffusion from the nanopores into macropores. The model presented in this paper can simplify into the conventional dual-porosity model by neglecting the shale gas desorption and diffusion processes. For conventional reservoirs, approximation expressions forfs in (11) can be found in the literatures [8, 9, 53–55]. The expression of fs introduced by Warren and Root for blocked matrix with pseudosteady-state interporosity flow is
(33)fs=ωf1−ωfs+λmf1−ωfs+λmfs.Using de Swaan’s approach [53], the solutions for Warren and Root’s model can be extended to the solutions of transient interporosity flow, which are
(34)fs=ωfs+sλmf1−ωf3stanh31−ωfsλmffor horizontal slab blocks and
(35)fs=ωfs+λmf5151−ωfsλmfcoth151−ωfsλmf−1for spherical blocks.The above three types of matrix block are the most popular idealization model for naturally fractured reservoirs; they have been widely used in the petroleum industry since they are proposed. So, the correctness of our model and method can be validated by comparing the result to the conventional model with spherical blocks.In our model, if the gas diffusion and desorption are neglected, the following relationship must be established, which is(36)ωf+ωm=1.Substituting the above equation into the expression ofgs in (11), which becomes
(37)gs=151‐ωfλmfs.It can be clearly seen that (11) can be simplified into the same solution for spherical blocks described by (35) for naturally fractured reservoirs by taking the above equation into it, which prove the correctness of our solution.


## 4.2. Well Production Performance Analysis

On the basis of the above solution, the semianalytical solution in Laplace domain can be transformed into real space using Stehfest numerical inversion algorithm [56]. Therefore, when the well is produced at a constant bottomhole pressure, the rate decline semilog type curves of production rate and cumulative production vs time can be obtained. The input parameters for simulation sensitivities for the following models are shown in Table 1.Table 1
Synthetic data used in the examples.
Parameters
Value
Parameters
Value
Initial reservoir pressure,pi, MPa
25
Reservoir temperature,T, K
320
Formation thickness,h, m
60
Gas specific gravity,rg, fraction
0.65
Well length,L, m
1200
Constant bottomhole pressure,pwf, MPa
4
Natural fracture system permeability,kf, mD
0.005
Natural fracture system porosity,∅f, fraction
0.01
Macropore system permeability,km, mD
10−6
Macropore system porosity,∅m, fraction
0.05
Langmuir pressure,PL, MPa
4
Langmuir volume,CL (m3/m3)
10
Number of fractures,M (integer)
10
Half fracture length,Lf, m
100
Outer boundary size,xe, m
800
Outer boundary size,ye, m
3000For a well producing at constant bottomhole pressure in a rectangular gas reservoir, the production rate and cumulative production vs. time curves in semilog coordinates with different fracture numbers are shown in Figure8. From the plot, it can be seen that the early time production rate is increasing with the increasing of fracture number, which is mainly because the larger the M is, the bigger the drainage area of the well will be during the early flow period. And when the pressure wave transports to the outer region formation, this effect will be gradually weaken. However, the cumulative production with different M under the same production time in later production period will be the same due to the reserve of them are equivalent.Figure 8
Influence of fracture number on well production performance.Figure9 shows the influence of Langmuir volume CL on the well production performance. Due to the mechanism, physical model in this paper has divided the pores in matrix system into the following two subspace: macropores and nanopores. The adsorbed gas will be desorbed from the particle surface in the matrix and then flow into the macropores by Fick diffusion. When the pressure in the macropore system is lower than the Langmuir pressure, desorption phenomenon will happen. So, when the well producing with the constant bottomhole pressure of 4 MPa, the pressure in the vicinity of the well will be lower than the Langmuir pressure, which will cause the adsorbed gas to desorb from matrix surface. So, we can find from Figure 9 that the production rate is larger with the increasing of Langmuir volume. But this effect is very small. When well is producing for a long time, the cumulative production of well with bigger CL will be larger than the others.Figure 9
Influence of Langmuir volume on well production performance.Figure10 shows the effects of permeability in natural fracture system on the well production performance. Because the natural fracture system is the main flow path of gas flow from the formation into hydraulic fractures, it has a significant influence on well production performance. According to the productivity equation of wells, the well production rate has a proportional relationship with the natural fracture permeability. And due to the constant control volume of our model, the bigger production rate at early flow period will lead to a bigger rate decline in the middle flow period. According to the cumulative production curves on Figure 10, we can clearly see that the larger the kf is, the larger the cumulative production will be in our given time period.Figure 10
Influence of natural fracture permeability on well production performance.Figure11 shows the effects of permeability in macropores on well production rate and cumulative production rate curves. The bigger km represents a stronger supplement of fluid from matrix system into natural fractures. Because the fluid flow from matrix system into natural fractures follows the transient interporosity flow, this supplement will happen immediately after the well opens. The bigger the km is, the larger the well production rate will be, and this trend will be lasting for a long time. So, for shale gas reservoir development, the permeability of the matrix system is a very important parameter to evaluate the well production performance.Figure 11
Influence of permeability in macropore system on well production performance.Figure12 shows the effects of macropores porosity (∅m) on well production performance. The value of ∅m presents the amount of free gas in macropore system of shale matrix system. With the increase of the ∅m, the quantity of free gas flows from macropores into natural fracture system by interporosity will increase. So, when well produces at a constant bottomhole pressure, the higher the porosity is, the larger the well production rate will be. Although the increment of production rate is not quite obvious, for the situation listed in this paper, the increment is still nearly 5000 m3/d. If this increment for conventional gas reservoir could be negligible, it is very considerable for shale gas reservoir due to the long production period.Figure 12
Influence of porosity in macropores on well production performance.Figure13 shows the effects of gas desorption time (τ) on the well production performance. As we know, the gas desorption time presents the gas desorption rate of adsorbed gas desorbed into free gas. The smaller the τ is, the faster the desorption process happens. It can be clearly seen from Figure 12 that the smaller the gas desorption time is, the larger the well production time will be, and the corresponding cumulative production will be larger too.Figure 13
Influence of desorption time on well production performance.


## 5. Discussion

This paper established a new mechanism model to describe the gas flow in multiscaled shale gas reservoir, which divided the pores in matrix systems into two subtypes: macropores and nanopores. The gas flow in nanopores obeys the Fick’s second diffusion and in macropores follows the interporosity flow as conventional dual-porosity model. On the basis of the above flow mechanism, the production performance of a MFHW in a rectangular outer-boundary shale gas reservoir is analyzed by semianalytical solution method. Comparing to the solutions of conventional naturally fractured reservoir, the correctness of our model is certified by neglecting the gas desorption and diffusion properties. After that, the effects of corresponding sensitive parameters on well production rate and cumulative production rate curves are analyzed. According to the results, we can conclude that the fracture number only affects the well performance in early flow period; the permeability of natural fractures mainly affects the early performance due to the finite drainage volume in our physical model. Comparing to the permeability of natural fracture system, the effects of permeability in macropore system on early and middle flow periods are obvious. So, for some high rate well in field, the macropore may be an important storage space, and which cannot be neglected or combined with other system to analyze.



---

*Source: 3148298-2018-10-18.xml*