A Mathematical Pressure Transient Analysis Model for Multiple Fractured Horizontal Wells in Shale Gas Reservoirs

Multistage fractured horizontal wells (MFHWs) have become the main technology for shale gas exploration. However, the existing models have neglected the percolation mechanism in nanopores of organic matter and failed to consider the differences among the reservoir properties in different areas. On that account, in this study, a modified apparent permeability model was proposed describing gas flow in shale gas reservoirs by integrating bulk gas flow in nanopores and gas desorption from nanopores. The apparent permeability was introduced into the macroseepage model to establish a dynamic pressure analysis model for MFHWs dual-porosity formations. The Laplace transformation and the regular perturbation method were used to obtain an analytical solution. The influences of fracture half-length, fracture permeability, Langmuir volume, matrix radius, matrix permeability, and induced fracture permeability on pressure and production were discussed. Results show that fracture half-length, fracture permeability, and induced fracture permeability exert a significant influence on production. A larger Langmuir volume results in a smaller pressure and pressure derivative. An increase in matrix permeability increases the production rate. Besides, this model fits the actual field data relatively well. It has a reliable theoretical foundation and can preferably describe the dynamic changes of pressure in the exploration process.


Introduction
Shale gas is known as a key resource to meet the increasing world energy demand because of its rich reserves and extensive distribution [1].To be more specific, many pieces of research have shown that shale gas reservoirs are characterized by extremely low permeability values.Consequently, hydraulic fracturing simulations and horizontal drilling are developed for commercial exploitation.In shale gas reservoirs, the complex pore structure, macro/microfracture network distribution, and microscale reservoir properties make the gas flow characteristics become more complicated than those in conventional reservoirs.Gas flow in organic matter, natural fractures, secondary fractures, and hydraulic fractures is controlled by different mechanisms [2].To be more specific, 20%-80% of shale gas exists in the adsorbed state on the surface of organic particulates.Desorption of adsorbed gas in organic matter affects gas transport mechanisms significantly [3,4].Understanding the pressure-and rate-transient behaviors of multistage fractured horizontal wells in shale gas reservoirs is of great importance for production forecasting, well placement, and configuration optimization [5].
In recent studies [6][7][8], gas transport mechanisms in organic matter can be described by continuum flow, slip flow viscous flow, the Knudsen diffusion, and transition flow.Wang and Li [9] made a comparison of collision frequency in gas-gas and gas-solid interactions.They concluded that the collisions between gas molecules and the solid wall account for a large percentage of the total collisions.Therefore, the Knudsen diffusion should be considered when studying gas transport in nanopores.Darabi et al. [10] also reported that, under typical shale gas reservoir conditions, the Knudsen diffusion dominates gas transport and the contribution to cumulative production can reach 20%.Adsorbed gas desorption from the surface of organic matter also plays a significant role in gas transport [11,12].Swami et al. [13] developed a modified model to investigate the contribution of adsorbed gas to the final recovery.They found that desorption of adsorbed gas can increase the pore diameter, reduce tortuosity, and cause extra slippage at the solid boundary.Yu and Sepehrnoori [14] analyzed the shale gas development history for a period of 30 years in North America and concluded that different shale samples demonstrate quite different adsorption capacities.
Mayerhofer et al. [15] studied the difference between conventional double wings symmetric fractures and hydraulic fractures of Barnet shale and were the first to propose the simulated reservoir volume (SRV) conception.In SRV, multistage hydraulic fracturing connects existing natural fractures, which generates a large fracture network [16].Kucuk and Sawyer [17] firstly developed an analytical model to investigate transient pressure in shale gas reservoirs.However, the model ignored diffusion flow and the desorption effect.Lee and Brockenbrough [18] analyzed production data coming from vertical wells and developed an analytical solution to describe transient pressure.Ozkan et al. [19] proposed a dual-mechanism, that is to say, a dual-porosity horizontal well model for shale gas reservoirs, which includes hydraulic fractures, as well as the inner region and the outer region.The model can adopt the diffusive flow mechanism but fail to take adsorbed gas into consideration.Stalgorova and Mattar [20][21][22] extended the trilinear flow model to a five-region model for MFHWs in homogeneous shale gas reservoirs.In order to investigate the effect of natural fractures on matrix permeability, Apaydin et al. [23] combined composite blocks into an analytical trilinear flow model.Zhao et al. [24,25] analyzed the pressure response and production performance of MFHWs in shale gas reservoirs and they introduced the source function theory to characterize SRV.Liu et al. [26] discussed the desorption and diffusion effects in the analytical model and studied the productivity-decline law.Wang [27] considered the stress sensitive effect of natural fractures and hydraulic fracture angles in the semianalytical model.His research implied that the stress sensitive effect has a significant impact on bottom hole flowing pressure.
In recent years, many general semianalytical models, considering microscopic seepage mechanisms, have been proposed with the rapid development of computer technology [24,25,[28][29][30].These models can describe predominant flow regimes for shale well but fail to consider all of the typical characteristics and possible situations in shale gas reservoirs.In the present study, a novel apparent permeability model was put forward by taking into consideration viscous flow, slip flow, transition flow, gas adsorption/desorption [31], and the poromechanical response [32,33].The novelty of this study concentrates on the improvement of the semianalytical model to conduct pressure-and rate-transient analysis in SRV, by considering microscale percolation mechanisms and heterogeneity.The rest of this study can be described in the following way.In Section 2, the apparent permeability model for gas transport in shale nanopores is introduced.The semianalytical model for MFHWs in a heterogeneous shale gas reservoir, by considering microscale percolation mechanisms in a dual-porosity formation, is introduced in Section 3. Log-log dimension pressure, dimension pressure derivative, and production type curves are plotted firstly, and then a sensitivity analysis of formation parameters is conducted in Section 4.An actual field case is studied in Section 5. Finally, some conclusions are drawn in Section 6.

Apparent Permeability Model for
Gas Transport in Shale Nanopores Viscous flow is strongly prominent when the short molecule free path is far less than the pore diameter.In other words, the transport mechanism is governed by viscous flow in macropores (where the pore diameter > 50 nm).The Hagen-Poiseuille equation can describe the molar flux [34] as stated below: where   is the viscous flow flux, mol/(m 2 * s);  is porosity, dimensionless;  is tortuosity, dimensionless;  is the pore radius, m;  is viscosity, Pa * s;  is the universal gas constant, J/(mol * K);  is temperature, K;  is pressure, Pa; and  is the shale gas transport distance, m.Gas viscosity gradually deviates from the traditional viscosity definition for high Knudsen numbers.Karniadakis et al. [35] modified the viscosity definition in the Knudsen layer by considering the rarefaction effect.
where  eff is the viscosity considering rarefaction effect, Pa * s;  is the rarefaction coefficient, dimensionless; and Kn is the Knudsen number, dimensionless.
According to the Darcy equation, the apparent permeability of viscous flow can be expressed as [36]: where  V is the apparent permeability, taking into account the effect of rarefaction, m 2 ;  std is the gas mole volume at standard conditions, 22.414 * 10 −3 m 3 /mol; and  ∞ is the permeability of viscous flow, m 2 .

The Knudsen Diffusion.
When there is shale gas transport through micropores under low pressure, the molecule free path is long and equal to the pore diameter.The Knudsen diffusion is prominent because collisions between molecules and the micropores wall are more frequent than the intermolecular collisions.The flux of the Knudsen diffusion can be expressed as [34] where  K is the Knudsen diffusion flux, mol/(m 2 * s);  K is the Knudsen diffusion coefficient, m 2 /s; and  is the gas concentration, m 3 /mol.The Knudsen diffusion coefficient is obtained as where  is the gas molar mass, kg/mol.The rough wall has a significant effect on the Knudsen diffusion.Darabi et al. [10] proposed a function to describe the roughness effect as follows: where  eff is the Knudsen diffusion coefficient which takes into consideration the roughness effect, m 2 /s;  is the mean free path of gas molecules, m; and   is the fractal dimension of the rough wall, dimensionless.By combining ( 4) and ( 6), the apparent permeability of the Knudsen diffusion can be obtained and written as follows: 2.3.Weight Factor.In the actual shale gas reservoirs, different transport mechanisms exist at the same time.Therefore, the total shale gas flux should be a weighted summation of the viscous flow flux and the Knudsen diffusion flux, based on their different contributions.The weight factors of the viscous flow and the Knudsen diffusion flow are defined as the ratios of intermolecular collisions and molecular/porewall collisions to total collisions, respectively.Based on the definitions for a gas molecule free path [11] and collision numbers [37], the weight factors of viscous flow and Knudsen diffusion can be approximated as follows: Therefore, the total apparent permeability in nanopores is With a decrease in pressure, the gas molecule free path increases immediately.The weight factor of the Knudsen diffusion also increases.The piece of research carried out by Wu et al. [38] showed that, under the 10 nm radius condition, viscous flow dominates when the pressure is larger than 4.02 * 10 5 Pa; otherwise, the Knudsen diffusion dominates.

Sorption-Induced Swelling
Response.Organic matter is described by weak strength and strong sensitivity to stress change.With an increase in reservoir pressure, adsorbed gas begins to desorb.Gas desorption results in shrinkage of the organic matrix and to an increased effective hydraulic diameter [39].According to the solid deformation theory [40] and the Langmuir isotherm equation [41], the relationship between the degree of solid deformation and the reservoir pressure can be written as follows: where Δ is the solid deformation degree, dimensionless;   is the Langmuir volume, m 3 /kg;   is the density of the shale matrix, kg/m 3 ;  is the shale matrix Young modulus, Pa;   is the Langmuir pressure, Pa;  ini is the initial pressure, Pa; and   is the rock compressibility, Pa −1 .Meanwhile, based on Seidle's model [42], the approximation relation between the porosity and the solid deformation degree is given by Assume that the pore volumes are proportional to the gas flow channels.In accordance with the capillary model, the effective hydraulic diameter can be obtained as 2.5.Gas Adsorption/Desorption.Because of the large surface area and oil-wet characteristic, nanopores have a strong adsorbed gas capacity.The mass balance equation, which considers the adsorbed gas, can be expressed as where  is the shale gas density, kg/m 3 ;   is the gas compressibility, Pa −1 ; and  is the mass of adsorbed gas per unit volume of shale, kg/m 3 .The apparent permeability, which considers desorption of absorbed gas, can be expressed as Cui et al. [12] defined the effective adsorption porosity as According to the Langmuir isotherm equation, the effective adsorption porosity can be expressed as Horizontal well Finally, the apparent permeability model for gas transport in shale nanopores can be expressed as

Productivity Prediction Model for MFHWs
3.1.Physical Model. Figure 1 shows the diagram of MFHWs in a shale gas reservoir.The reservoir is divided into seven contiguous regions: two upper-reservoir regions (regions 6 and 7), two outer-reservoir regions (regions 4 and 5), two inner-reservoir regions (regions 2 and 3), and a hydraulic region (region 1).Not all reservoir zones between primary hydraulic fractures are stimulated by induced fractures.1D linear flow is assumed within unstimulated region and the flow direction depends on the location.As shown in Figure 1, a region of higher permeability around each fracture, the socalled SRV area, is introduced to represent fracture branching (region 2).The SRV region occupies part of the space between fractures and the flow in the other part (region 3) towards region 2, parallel to the wellbore.Regions 4 and 5 are in the outer reservoir beyond the tips of the hydraulic fractures and they connect with regions 2 and 3 in  direction.Because of the higher permeability in region 2, the pressure decreases in this area first.Therefore, the flow direction in regions 4 and 5 is towards regions 2 and 3, perpendicular to the wellbore.Although the main extending direction of the hydraulic fractures is horizontal ( direction in Figure 1), we also need to consider fractures height in vertical direction ( direction in Figure 1).Based on this consideration, regions 6 and 7 are defined in the upper reservoir beyond the tips of the hydraulic (connecting with regions 2 and 3 in  direction) and the flow direction is vertical to the wellbore, towards the lower reservoir.
The main assumptions considered for this model are listed below: (1) The shale reservoir has been fractured and has a constant reservoir thickness.All hydraulic fractures are symmetric and arranged uniformly along the horizontal well.Length and conductivity of the hydraulic fractures are the same.
(2) Natural fracture permeability is stress-dependent and shale gas flow in the matrix is driven by concentration difference.
(3) The initial pressure throughout the reservoir is uniform.
(4) Hydraulic fractures penetrate the entire pay formation.Fracture interference is ignored.After SRV, large fracture networks can be generated in region 2 and it has higher natural fracture permeability compared to the other regions.
(5) The shale gas flow in the reservoir is considered to be an isothermal flow and the effect of gravity and capillary pressure are negligible.
(6) Nonsimulated areas are described by a dual-porosity system consisting of a shale matrix and natural fractures.
(7) The only route for shale gas to reach the horizontal wellbore is through hydraulic fractures.The flow rate of MFHW is the sum of the flow rate from every hydraulic fracture.

Mathematical Model.
According to the mass conservation equation, the motion equation, and the state equation, the partial differential formulas for shale gas flow in a matrix system and a fracture system can be obtained as follows [19,27]: where   is a characteristic dimensionless parameter of stress sensitivity and relates to the permeability modulus . is defined as  = (/2), where   =   exp(−Δ  ).
Pedrosa's substitution [43] is used to weaken the nonlinearity of (19).By setting To solve (18) and ( 20), these two equations were converted to the Laplace domain.And then the definite boundary continuity is introduced to develop an integral seepage flow differential equation.The final equation can be written as where A series of variables and dimensionless variables are defined in Table 1.

Region 7.
In this upper-reservoir region, the governing diffusivity equation for shale gas flow in the -direction becomes where At the top of the reservoir ( =  2 ), the boundary condition is no-flow: On the interfaces ( =  1 ) of regions 7 and 3 and regions 7 and 5, pressure continuity condition can be obtained as Equation ( 22) is solved simultaneously with ( 24) and ( 25).The final solution for (22) can be given as 3.2.2.Region 6.Consistent with region 7, the governing diffusivity equation for shale gas flow in the -direction becomes where ] .
The boundary condition at the top of the reservoir The boundary conditions on the interfaces of regions 6 and 2 and regions 6 and 4 are The solution for (27)  (31)

Region 5.
In this region, shale gas flow takes place in the -direction and the -direction.So the governing diffusivity equation for region 5 is where Along the horizontal well ( = 0), there is a no-flow boundary: Flux across the interface of regions 5 and 7 ( =  1 ) is continuous; therefore, Substituting (34) and ( 35) into (32), each of its terms from 0 to  1 is integrated with respect to   .The governing diffusivity can be simplified as In this region, shale gas flow in the -direction and with a closed boundary condition ( =   ) is expressed as Pressure continuity condition between regions 5 and 3 (at  =  1 ) is expressed as By solving (36), the solution can be obtained as where 3.2.4.Region 4. In a similar way, shale gas flow in this region takes place in the -direction and the -direction.The governing diffusivity equation for region 4 is where The same method is also used in region 5, where the governing diffusivity equation can be simplified as No-flow condition at the outer-reservoir boundary ( =   ) is given as On the interface of regions 4 and 2, the boundary condition is The solution for ( 41) is given as where 3.2.5.Region 3.This region is adjacent to regions 5 and 7, so shale gas flow direction takes place along the -direction, direction, and -direction.The governing diffusivity equation becomes where For the -direction, the boundary condition used is a no-flow condition for a horizontal well ( = 0) and a flux continuity condition is used on the interface of regions 5 and 3 ( =  1 ): For the -direction, (51) can be obtained: By substituting (50) and ( 51) into (48), the final governing diffusivity equation becomes Between two fractures, there is a no-flow boundary (at  =  2 ): On the interface of regions 3 and 2, the pressure continuity condition is written as The solution for (48) is given as where (56) 3.2.6.Region 2. Similarly, as region 2 is adjacent to hydraulic fracture, the governing diffusivity equation is Consistent with region 3, the simplified equation can be written as where At the interface of regions 2 and 3, the boundary condition is characterized by flux continuity: And between the hydraulic fracture level and region 2, the boundary condition is Solving ( 58) to (61), where (63)

Region 1.
In this region, the diffusivity equation is given by where   =     / 2 .The boundary condition can be stated by the expressions: The solution can be obtained: where Assuming that pressure distribution along the horizontal well is uniform, the dimensionless down-hole pressure can be solved by setting   = 0 [44]: Shale gas flow in hydraulic fractures is not completely linear.The nonlinear streamlines result in an additional pressure drop.Mukherjee and Economides [45] presented a formulation to take this additional pressure drop into consideration.

25
×10 −3 ×10 −3 r = 6 nm apparent permeability (adsorbed gas) r = 8 nm apparent permeability (adsorbed gas) r = 6 nm apparent permeability (bulk gas) r = 8 nm apparent permeability (bulk gas) Based on the Duhamel principle and the superposition principle, the dimensionless flowing bottom hole pressure (FHBP), which takes into consideration the wellbore storage effect and the skin effect, can be expressed as [30] The FHBP for MFHWs, when considering the stress sensitivity effect, can be obtained as The relationship between the dimensionless rate under the condition of constant pressure and dimensionless pressure caused by the constant rate is described as [46]

Discussion and Analysis
4.1.Apparent Permeability in the Shale Matrix.Figure 2 shows the variation of the apparent permeabilities against the average pressure.The input parameters are listed in Table 2.
Based on Figure 2, one can observe that, in the early period, with a decrease of the average pressure, the bulk permeability  app,bulk decreases dramatically, but the apparent permeability caused by gas desorption  app,ad increases greatly.In the later period,  app,bulk increases slightly with an increase of the pressure. app,ad also increases in a similar way with an increase in pressure.The explanation is that, at high pressure, the Knudsen diffusion increases significantly.As can be seen in Figure 2, viscous flow dominates at high pressure and the Knudsen diffusion dominates at low pressure.Figure 2 also indicates that the pore diameter has a great impact on the apparent permeability.With a decrease of the pressure, nanopores become larger and larger under the influence of matrix shrinkage, so the apparent permeability increases gradually.However, there is an inflection point between the apparent permeability and the average pressure.Furthermore, pore diameter is bigger and the time of the inflection point for apparent permeability occurring appears later.The apparent permeability decreases faster after the inflection point.The results presented in Figure 2 agree well with those presented by Wu et al. [38].

Effects of Hydraulic
Fracture Half-Length.Figures 3 and  4 show the effect of hydraulic fracture half-length on pseudopressure, pseudo-pressure derivative, and dimensionless production.The input parameters are listed in Table 3.According to the results, fracture half-length affects the whole production process.The longer the hydraulic fracture half-length, the more significant the pressure drop and the higher the production rate.The contact area between the hydraulic  fracture and the shale gas reservoir is larger when the fracture half-length is longer.Therefore, when gas flow resistance is reduced, more shale gas can flow to the fracture.Furthermore, larger flow rate results in significant pressure drop and a larger production.5 and  6 present the curve of the relationships of dimensionless time and pseudo-pressure and dimensionless production rate under different hydraulic fracture permeability.The figures show that hydraulic fracture permeability dominates all of the production period except for when we have early pure wellbore storage regime and boundary flow regime.The results are quite different from those observed in a single porosity tight reservoir where hydraulic fracture permeability only exerts an influence on early time pressure behavior [47].

Effect of Hydraulic Fracture Permeability. Figures
The figures also indicate that, with an increase in fracture permeability, the pseudo-pressure drop and the production rate growth tend to slow down.Therefore, there is an optimal fracture permeability for a certain reservoir.

Effect of the Langmuir Volume.
The effects of the Langmuir volume on pressure and production are plotted in Figures 7 and 8.The Langmuir volume is defined as the maximum volume of shale gas that can be absorbed to a unit weight organic matrix at infinite pressure.The Langmuir volume is a significant parameter for gas adsorption and desorption.As illustrated in Figures 7 and 8, the larger the Langmuir volume is, the lower the pressure drop and the production will be.This is because a larger Langmuir volume results in a greater amount of absorbed gas which provides an   extra contribution for well reserves.However, the Langmuir volume does not affect the linear flow stage as only shale gas in hydraulic fractures contributes to production at that period.It can also be seen that, with an increase of the Langmuir volume, the dimensionless rate decreases slowly.This results in longer production time.

Effect of Matrix Radius. Figures 9 and 10
show the effect of matrix block radius on pseudo-pressure, pseudo-pressure derivative, and dimensionless production rate.The shape and size of the matrix blocks exert a very significant effect on the interporosity flow stage.In this model, the value of the matrix block radius determines the magnitudes of the interporosity flow coefficient and the storativity ratio.According to the results, the matrix block radius mainly affects the transitional stage and the stage of bilinear flow.That is to say, the larger the matrix block radius, the more significant the pressure drop of the transitional stage.Consequently, the stage of bilinear flow occurs later.When the matrix block radius is large enough, bilinear flow disappears.Meanwhile, the magnitude of the matrix block radius also affects the production in the early and middle stages.With an increasing matrix block radius, the production rate gradually decreases, and the production decline rate decreases as well.In the late stage of production, the production rates under different matrix block radii are approximately the same [30].fracture system.Under ultra-low matrix permeability conditions, the pressure variation is more sensitive.

4.7.
Effect of Induced Fracture Permeability.Figures 13 and 14 represent the effect of induced fracture permeability on pressure and production.The induced fracture permeability of shale reservoirs is an extremely important reservoir parameter, as shale gas transport to hydraulic fracture should take place through induced fracture systems.Therefore, the induced fracture permeability determines the scale of gas well productivity.The results show that the higher the induced fracture permeability, the more significant the pressure drop during the transitional period.With a continuous increase in induced fracture permeability, the duration of the stage of bilinear flow becomes shorter and shorter until it completely disappears.Figure 13 indicates that the induced fracture permeability has a more significant influence on pressure than matrix radius and matrix permeability.Meanwhile,   the induced fracture permeability also affects the scale of production rate; that is, the higher the permeability, the larger the scale.However, as far as production rate is concerned, the secondary fracture permeability only affects its scale, but not its decline rate.

Case Study
The MFHW is located in the transitional zone between the Sichuan Basin and the Yunnan-Guizhou Plateau (400-1,200 m in altitude) and mainly developed a marine and continental sedimentary facies.The specific stratum parameters are listed in Table 4.
As seen in Figure 15, the production data show that there is a relatively high degree of correlation between pressure and production; the production data present an obvious response to the pressure fluctuations caused by the test carried out in the early stage; therefore, the data have a certain degree of reliability.As shown by a comparative analysis on the actual production data and production curve, the measured data are deviated from the typical curve by a relatively significant extent.This is mainly because of the skin effect and wellbore storage effect on the early stage seepage.The fitting curve shows that the data points are concentrated in the stage of linear flow of matrix (see Figure 16), mainly because the extremely low permeability of matrix-fractures has extended the duration of the stage of nonsteady state linear flow.

Conclusion
In this study, an analytical model was presented to simulate shale gas flow through MFHW in a shale reservoir.Based on microscale percolation mechanisms in nanopores, a seven-linear-flow model for dual-porosity formations was established.This model takes viscous flow, the Knudsen diffusion, the sorption-induced swelling response, the

2. 1 .
Viscous Flow.The total shale gas flux is composed of a viscous flow flux and a Knudsen diffusion flux component.

Figure 1 :
Figure 1: Diagram of MFHWs and of the seven-flow-region model.

Figure 3 :
Figure 3: Effect of hydraulic fracture half-length on pressure.

Figure 9 :Figure 10 :
Figure 9: Effect of matrix radius on pressure.

Figure 13 :
Figure 13: Effect of induced fracture permeability on pressure.

Figure 14 :
Figure 14: Effect of induced fracture permeability on production.

Table 4 :
Parameters of the typical well.