A Model for a Multistage Fractured Horizontal Well with Rectangular SRV in a Shale Gas Reservoir

Although great success has been achieved in the shale gas industry, accurate production dynamic analyses is still a challenging task. Long horizontal wells coupling with mass hydraulic fracturing has become a necessary technique to extract shale gas efficiently. In this paper, a comprehensive mathematical model of a multiple fractured horizontal well (MFHW) in a rectangular drainage area with a rectangular stimulated reservoir volume (SRV) has been established, based on the conceptual model of “tri-pores” in shale gas reservoirs. Dimensionless treatment and Laplace transformation were employed in the modeling process, while the boundary element method was used to solve the mathematical model. The Stehfest numerical inversion method and computer programing techniques were employed to obtain dimensionless type curves, production rate, and cumulative production. Results suggest that 9 flow stages can be observed from the pseudopressure derivative type curve when the reservoir and the SRV are large enough. The number of fractures, SRV permeability, and reservoir permeability have no effect on the total production when the well is abandoned. As SRV and reservoir permeability increases, the production rate is much higher in the middle production stage. Although the SRV scale and its permeability are very important for early and intermediate production rates, the key factors restricting the shale gas production rate are the properties of the shale itself, such as adsorbed gas content, natural fractures, and organic content. The proposed model is useful for analyzing production dynamics with stimulated horizontal wells in shale gas reservoirs.


Introduction
Shale gas, together with coalbed methane, which was once the source rock of oil and gas. Shale gas is a typical unconventional gas stored with adsorbed, free gas ( [1]; Li Yong et al., 2019). Shale gas consists mainly of methane, and there are many kinds of interlayers in shale strata, such as siltstone, silty mudstone, and politic siltstone. Currently, shale gas extraction is developing rapidly and has played an important role in the US energy industry since 2000. It has spread quickly to China and Australia as a potential energy source [2]. For example, the annual production of shale gas in 2007 in the US was 56.4 billion cubic meters, while it reached 447.0 billion cubic meters in 2016; its proportion to the world annual production is 12% [3]. With the success of shale gas development, shale gas has the potential to become the primary energy source in the future.
Due to the special formation characteristics and pore types of a shale gas reservoir, exploiting shale gas requires hydraulic fracturing. Specifically, i.e., horizontal wells with massive hydraulic fracturing [4][5][6]. Massive hydraulic fracturing of horizontal wells in shale gas reservoirs creates a fracture network around the well called stimulated reservoir volume (SRV). The SRV, which is the main method of increasing shale gas production, consists of hydraulic fractures and a fracture network. Therefore, it is important to devise an effective and accurate model to describe the SRV and to help understand the flow characteristics of shale gas underground.
Kucuk and Sawyer [7] proposed a model to analyze the adsorption and desorption characteristics in fractured reservoirs and derived a numerical solution to the model considering the effects of desorption and slippage. Lancaster and Gatens III [8] analyzed postfractured well testing data of eastern Devonian shale, which was described as a dualporosity reservoir, and those data were used to estimate fracture length and fracture conductivity. Devonian shale wells showed that fracture half-lengths calculated from conventional postfracture pressure-buildup tests were shorter than designed, based on which Johnston and Lee [9] proposed a new method to determine net pay thickness and presented results suggesting that conventional tests of those wells are inadequate.
Brown et al. [10] and Ozkan et al. [11] assumed that the main flow type in fractured horizontal wells is linear, and established a trilinear flow model for a multistage fractured horizontal well in a shale gas reservoir using the model presented by Lee and Brockenbrough (Lee and Brockenbrough 1986). Based on the linear-flow theory, a plate dualporosity model, and the Warren-Root block dual-porosity model, Bello and Watenbargen [12] and Al-Ahmadi et al. [13] proposed a model to analyze the production dynamics and pressures of MFHWs in shale gas reservoirs. Nobakht et al. (Nobakht et al., 2012) first proposed the conceptual model of different combinations of horizontal wells, porous media types, and hydraulic fractures in shale gas reservoirs, and then analyzed the different flow stages during the development process. Assuming SRV exists only in a limited region around the fractures and fracture half-lengths equal to the width of the reservoir, Stalgorova and Mattar [14] derived a modified trilinear flow model and obtained the solution of this model in Laplace space. Xu et al. [15] proposed an unsteady seepage model of a MFHW at constant rate production in a shale gas reservoir. The model assumed that the flow types in a fractured horizontal well are alternately linear and transitional, and the fracture network exists in a limited region around the hydraulic fractures.
To summarize, there are two types of models to analyze well production performance of MFHWs with SRVs in shale gas reservoirs. One is the linear flow model, which is based on the assumption that all fractures have the same length and conductivity and are spaced uniformly along the horizontal well ( [11][12][13][14][15][16][17][18][19][20][21]; Zhao et al., 2016b; [22][23][24][25]), and the main related schematics of linear models are listed in Table 1. Although it is common field practice to design equally spaced hydraulic fractures, the requirements of equal properties and fracture lengths of shale gas reservoirs are hard to satisfy because of the development of natural fractures and their anisotropy and heterogeneity. The other type of flow model is the numerical simulation model based on finite-element, finite-difference, or other numerical methods [26][27][28][29][30][31][32][33][34][35][36].
Normal shale gas reservoirs have a dual-porosity nature, with macropores and micropores in matrix, and microfractures made by hydraulic fracturing (Ge & Zhang, 2012; [38][39][40]). The output process of shale gas from the reservoir can be subdivided into 4 stages: (1) the gas desorbed from the matrix particle surfaces, (2) diffusion from the matrix to the macropores, (3) flows from the macropores to the microfractures, and (4) Darcy flows through the fracture network to the bottom hole.
Well testing analysis is an important dynamic monitoring tool for petroleum engineers to accurately determine the parameters of oil and gas reservoir engineering, to evalu-ate the stimulation in the shale gas reservoir, and to provide more effective plans for the oil field. In addition, a suitable mathematical model for horizontal wells with SRVs in shale gas reservoirs is vital to well testing.
As we can see from the above discussion, both analytical methods and linear flow models have their own deficiencies in obtaining production dynamics in shale gas reservoirs. On one side, analytical methods, such as point source and Green's function method, can only deal with seepage problems of shale gas in regular-shaped reservoirs, such as reservoirs with circular or infinite outer boundaries. However, the practical boundary of shale gas reservoirs is quite complex considering the effects of SRV, which in consequence cannot be handled by this traditional analytical method. On the other, linear flow models usually put the whole reservoir as a rectangular area, and subdivide the reservoir region into some small rectangular areas, as shown in Table 1. The subdivision principle of the small rectangular sections is the gas flow capacities in formations, i.e., a bigger flow capacity in the SRV region and a weaker flow capacity as the region. In linear flow models, a quarter of the area around a hydraulic fracture is usually studied for simplicity, which brings into another assumption that all hydraulic fractures are uniformly spaced and with the same length. These assumptions are also not in accordance with practical conditions of shale gas reservoirs. Besides, other flow patterns, such as redial flow and elliptical flow, and interferences among fractures cannot be reflected from the linear flow models. Thus, a more comprehensive and accurate model considering both the rectangular composite nature of a shale gas reservoir and the fracture heterogeneities need to be established in an effort to analyze production dynamics better in shale gas reservoirs. In this paper, we present a percolation model for a multistage horizontal well in a rectangular shale gas reservoir with an SRV considering unsteady transport theory, which cannot be tackled by analytical methods or linear flow conceptions. The boundary element method is employed to solve the proposed mathematical model, with parameter sensitivity analysis conducted for production dynamics. Results show that the model we propose should be useful for analyzing stimulated horizontal wells in shale gas reservoirs.

Transport Mechanism of Shale Gas
2.1. Characteristics of Shale Gas Reservoir. The characteristics of a shale gas reservoir are very different from those of a conventional gas well, in both the structure and the flow mechanism. The typical classification method is that of Loucks et al. [41]. They classified the pore types of the shale matrix and the natural microfractures; their classification is comprised of the following 4 types: intergranular pores, intragranular pores, organic pores, and natural microfractures ( Figure 1 shows the main pore types of a shale gas reservoir). The first 2 pore types are related to the mineral grains. The difference between them is that intergranular pores develop between grains and intragranular pores develop within grains. Organic pores are intrapores of organic matter. In shale, the main pores are micropores, nanopores, and microfractures, so approximately 20% to 80% of the gas is absorbed by the 2 Geofluids 3 Geofluids particles in the matrix, and the rest is free in the microfractures. Gas flow behaviors in shale nanopores are quite different from that of gas flow in conventional reservoirs [42][43][44][45]. The percolation mechanism includes diffusion, adsorption/desorption, and slip flow [46][47][48].

Adsorption
Models of Shale Gas. The physical adsorption defined by Brunauer et al. [49] occurs when gas contacts a solid surface: some gas molecules are captured by the solid because of the intermolecular force. This means that if the volume of gas is constant, gas pressure drops; if the pressure is constant, the volume of gas decreases. The gas molecules captured by the solid are said to be adsorbed, and the solid is called an adsorbent.
Because shale gas adsorption models are numerous, it is necessary to choose the best one before using it. Here, we use a nonlinear fitting method to analyze the degree of fit between adsorption models and some adsorption data from subcritical and supercritical states.
In the original state of a shale gas reservoir, adsorbed and desorbed gasses are in dynamic pressure equilibrium. When the equilibrium is disturbed after production, the gas adsorbed in the organic matter surfaces begins to desorb, and then the desorbed gas becomes free and is stored in the  Congruent with the previous analysis, the Langmuir isotherm is one of the most popular models used to describe the gas adsorption/desorption process. The amount of gas adsorbed by a solid surface is given by the Langmuir equation below, characterizing the desorption processes as a function of pressure at constant temperature [33]: where V is the volume of the adsorbed gas (cm 3 /g), P m is the gas pressure (MPa), V L is the Langmuir volume (cm 3 /g), P L is the Langmuir pressure (1/MPa). In addition, the expression of desorbed gas mass flux that desorbs from a reservoir whose volume is V b during a unit of time is where q des is the mass flow rate of gas desorption from the reservoir of volume V b (kg/s), ρ gsc is the density of shale gas under standard conditions (kg/m 3 ). Incorporating equation (1) into equation (2) results in Introducing pseudopressure into equation (3), it can be rewritten as [50] It is assumed that the desorbed process of shale gas is instantaneous and all the desorbed gas will become free gas, so equation (4) can be used as a continuity equation of a fracture network or matrix.

Model for Gas Flow in a Fractured Shale
Gas Reservoir Macropores exist in a practical shale matrix, which may lead to an error if gas storage and flow are ignored in these macropores. Therefore, researchers have devised a "tri-pore" conceptual model, where the shale gas reservoir consists of matrix pores, macropores, and microfractures. When the adsorbed gas desorbs from the surface of the matrix particles, it enters the macropores first and then flows into the microfractures [39,[50][51][52]. Figure 2 shows the physical model. Song [39] first proposed a physical model concerning transfusion in the macropores of the matrix. But it was only a physical seepage-mechanism model, which was not modeled mathematically. In this paper, a mathematical model for this physical conception is proposed and solved. The interporosity flow of free gas between macropores and microfractures can be divided into 2 types: (1) a pseudosteady interporosity flow model and (2) For unsteady interporosity and spherical matrix particles, the continuity equation in the matrix can be expressed as follows: Because the initial pressure is equal throughout the reservoir, the initial condition in the matrix equals to the initial reservoir pressure p i :

Geofluids
Because seepage of the gas in the macropores is symmetrical, the inner boundary condition is Because the outer boundary of the spherical matrix particles is connected to the fracture system, the outer boundary pressure of the matrix is equal to the pressure of the fracture system: When the flow between the macropores and the microfractures is unsteady, the interporosity flow rate is By the Langmuir isotherm equation, the desorption rate of shale gas q des can be written as By introducing dimensionless variables into the above mathematical models and then reformatting them by several simple transformations (a detailed derivation process is described in Appendix A, while the dimensionless variables are shown in Appendix B), the dimensionless seepage equation of the fracture system can be obtained as where the expression of f ðsÞ is 3.2.2. Pseudosteady Interporosity Flow Model. For pseudosteady gas flow from the macropores to the fractures, incorporating pseudosteady interporosity flow rate and pseudopressure into the seepage equation of the fracture system, the seepage equation of fracture system can be written as For pseudosteady gas flow from the matrix to the fractures, the seepage equation of the macropores is where c d is an additional compressibility that considers the effect of desorption. The expression of c d is Incorporating the above dimensionless variables into the seepage equation of the matrix and the fracture system, the dimensionless seepage equations of the matrix and the fracture system are Using the Laplace transform method, equations (17) and (18) can be rewritten as According to equations (19) and (20), the seepage equation of the fracture system is where 4. MFHWs in a Rectangular Composite Shale Gas Reservoir with an SRV 4.1. Physical Model. In a shale gas or tight-gas reservoir, hydraulic fracturing is widely used to improve the productivity of the reservoir. A MFHW can make the cracks contact the matrix completely, form flow channels, and ensure efficient extraction from an unconventional reservoir [52,[55][56][57][58].
In 2008, Mayerhofer et al. [57] studied the change of cracks in Barnett shale by microseismic technology and proposed the SRV concept. The SRV refers to an effective reservoir that is made by hydraulic fracturing, with a composite 6 Geofluids fracture network formed as a consequence. SRV can maximize the contact surface between matrix and fractures so that the seepage distance of oil and gas in any direction from matrix to cracks will be shortened, which improves the effective permeability of the entire reservoir. According to the microseismic monitoring map for an MFHW in Barnett shale reported by Fisher et al. [56], the physical model ( Figure 3: left) can be approximated by the rectangular composite shale gas reservoir conceptual model (Figure 3: right) with enough accuracy, which is better than the circular composite conceptual model for long horizontal well cases.
To simplify the analysis, we use the following basic assumptions for this model: the MFHW is located at the center of a closed rectangular reservoir; the reservoir length, width, and height are x e , y e , and h, respectively; the fractures are perpendicular to the horizontal well and randomly distributed along it; the effective well length is L (from the leftmost to the rightmost fractures); each fracture is symmetrical or asymmetrical, and the fracture half-lengths are L fui and L fli ; the length and width of the SRV region are x m and y m , respectively; the number of fractures is M, and every fracture has infinite conductivity; no flow from the matrix to the wellbore and the pressure drop caused by the gas flow through the wellbore are ignored.

Boundary Element Model of a Composite Rectangular Gas
Reservoir considering the SRV. According to the physical model described above, the governing equation for gas flow in an inner fracture network system is ( [59][60][61][62], 2019) For the outer region The cohesive conditions on the surface between the inner and outer regions are Because the outer boundary of our model is closed, the outer boundary condition is Incorporating the dimensionless variables defined in Appendix B into equations (23) and (24), the mathematical model in the Laplace space becomes

Geofluids
When using the boundary element method to solve the model above, the fundamental solutions of differential equations are needed. It is assumed that the fundamental solutions of the inner and outer regions in the shale gas reservoir are E 1 ðR D , R D ′ , sÞ and E 2 ðR D , R D ′ , sÞ, respectively, which satisfy Since the fractures are assumed completely open in this paper, the expressions of fundamental solutions E 1 and E 2 are Both sides of equations (31) and (32) are multiplied by m f1D and m f2D , respectively, and the equations can be obtained by Both sides of equations (23) and (28) are multiplied by E 1 and E 2 , respectively, and we have Using equations (34) and (35) minus equations (36) and (37), respectively, the equations can be transformed as Integrating equations (38) and (39) on the inner region Ω 1 and outer region Ω 2 , respectively, we have 2  3  5  6  7  8  9  1  4 Inner SRV region Outer reservoir region Figure 4: Grid of the inner and outer boundaries and the sequence of the discrete elements in a rectangular composite gas reservoir. 8 Geofluids We use Green's second identity as follows [61,62]: where SðS = ∑ i S i Þ is the entire boundary of region Ω. Based on the characteristics of equation (42) and the delta function, the differential equation of the inner and outer regions can be changed into an integral equation on the boundary, and the boundary integral equations of equations (40) and (41) are There are two kind of boundaries (the outer and the inner) in the gas reservoir. It is assumed that the outer boundary can be discretized into NO elements and the inner boundary can be discretized into NI elements, so the total number of the elements for the boundaries is NO + NI. The serial number of inner boundary elements is clockwise, while it is counterclockwise for outer boundary elements, as shown in Figure 4.
For the integral equation on the boundary, R D is an arbitrary point on the boundary, while R D ′ is an arbitrary point in the reservoir region. R S1D and R S2D are the dimensionless coordinates of source points in the inner and outer regions, respectively. If R D ′ is not only in the reservoir region but also on the boundaries, equations (43) and equations (44) can be rewritten as where θ is a constant related to the geometric shape of R D ′ , and the expression of θ is [62] where α is the angle of the boundary tangent line at R D ′ . Incorporating the fundamental solutions E 1 and E 2 into equations (45) and (46), a series of linear equation sets can be obtained by solving the reservoir points R D ′ along the boundary points and fracture points, after which the dimensionless pseudopressure response curves of the fractured well can be obtained by solving the equation sets above. Due to the complexity of the process, its details are not shown in this article; interested readers may consult these papers for details: [52,59,60,[63][64][65][66][67][68][69][70][71].

Pressure and Production
Type Curve Analysis. The boundary element method is employed to acquire type curves of dimensionless bottomhole pressures at constant rate, and the production rate and cumulative production at a constant pressure are based on the parameters in Table 2. Figure 5 shows the well-test type curves of reservoirs with different sizes and SRVs. The following flow stages can be observed when the reservoir and SRV sizes are large enough: (a) Stage 1: early wellbore storage (b) Stage 2: transition flow between wellbore storage and early linear flow of a fracture system. A hump is in the pseudopressure derivative curve for this stage,    Figure 5: Effect of reservoir size on type curves.  Figure 7: Effects of number of fractures on production rate and cumulative production curves.
t D /C D When the sizes of the reservoir and the SRV are reasonable, the pseudopressure and pseudopressure derivative curves are shown in Figures 5-11 as the plots with yellow pellets. After the early linear flow of the fracture system, the pressure reaches the region close to the reservoir boundary, so the curves turn upward, and then the pseudopressure and pseudopressure derivative curves are almost parallel to each other as a straight line, but the slope is not 0.5. When the pressure reaches the reservoir boundary, the 2 curves overlap. Figures 6 and 7 show the effects of the fracture number on well-test type curves (dimensionless pseudopressure and its derivative) and production curves (production rate and cumulative production), respectively. As shown in Figure 6, for a constant reservoir size and SRV, the number of hydraulic fractures affects mainly the early flow periods following K f1 = 0.05 K f1 = 0.25 K f1 = 1.25 Figure 9: Effects of SRV permeability on production rate and cumulative production curves.
t D /C D Figure 10: Effect of reservoir permeability on well-test type curves.
the wellbore storage flow period. When the pressure wave transports to the outer region, the effects will vanish and all pressure and pressure derivative curves will coincide. When a well produces at a constant rate, the higher the fracture number, the higher the early production; this effect lasts for just a few years. Figures 8 and 9 show the effects of different flow region permeabilities on well-test type curves and production curves, respectively. The effects occur in the middle flow stage. The position of the well-test type curves is lower when the SRV permeability is higher, but the difference between the well-test type curves becomes smaller and smaller as the pressure reaches the reservoir boundary. Figure 9 shows that when SRV permeability increases, the early production rate is higher than the previous one, the production rate is much higher in the middle stage, and the SRV permeability has no effect on the final cumulative production. Figures 10 and 11 show the effects of reservoir permeability on type curves and production curves, respectively. Because the reference permeability in our model is the reservoir permeability, the reservoir permeability has the opposite effect on type curves compared to that of SRV. The higher the reservoir permeability, the higher the loca-tion of the type curves (as shown in Figure 10), which does not mean that the drawdown pressure is higher. When the well is producing at constant pressure, the gas supply is timely due to the high permeability of the reservoir, and gas production rate increases with reservoir permeability increasing (as shown in Figure 11). Even in the later production period when the well has been producing for 1000 d and 2000 d, the well production rate for a reservoir with a higher permeability is still much higher, as shown in Table 3. Therefore, a desert area is very important for the exploitation and development of unconventional shale gas reservoirs.

Conclusions
We established a comprehensive mathematical model to describe a multistage fractured horizontal well in a rectangular composite shale gas reservoir with an SRV, and the solution was also obtained by a boundary element method. The main conclusions we derived are as follows: (1) A composite rectangular model of a fractured horizontal well using the "tri-pores" conceptual model is proposed in this paper, and the complex model was solved by BEM, which has not been comprehensively reported before (2) The number of fractures affects the well production rate at primarily the early flow period; the effect on the later flow period is weak (3) The permeability in the SRV region has a substantial effect on the early well production rate; the higher the permeability, the greater the rate will be  Figure 11: Effects of reservoir permeability on production rate and cumulative production curves.