Performance evaluation of multiple fractured horizontal wells in shale gas reservoirs

With the rapid development of hydraulic fracturing technology, shale gas has become a significant natural gas resource in China. At present, SRV (stimulated reservoir volume) fracturing technology is widely used to exploit shale gas reservoirs. To evaluate the shale gas reserves and design a development plan, a simple and rapid method is needed to predict the production and performance of the multistage fractured horizontal wells with stimulated reservoir volume in shale gas reservoirs. Network fractures with different lengths and widths are not symmetrical to the horizontal well. To account for the complexity of the fracture system and desorption characteristics of shale gas, a semianalytical multiregion seepage model has been established to analyze the performance of multiple fractured horizontal wells. In this paper, a triple‐porosity dual‐permeability model is developed for simulating the flow mechanism of shale gas. Different lengths and widths of network fractures, different fracture intervals, and asymmetry distribution are considered in this model. By employing the principle of Poisson's superposition and Bessel transformation, the fundamental solution of instantaneous point source function has been derived to calculate the pressure response of each hydraulic fracture. Then, the mathematical model for unsteady performance evaluation in multiple fractured horizontal wells has been obtained to simultaneously solve the pressure response of every hydraulic fracture. The sensitivity analysis in terms of the lengths and widths of network fractures, fracture intervals, and the asymmetry distribution of a complex fracture system is presented. The research results show that the SRV morphology will exert a great influence on the performance of multiple fractured horizontal wells.

complicated flow characteristics in shale reservoirs, horizontal well and volume fracturing technology are used to improve gas well productivity. In general, there are differences in fluid flow ability between shale matrix and fractures, and these differences were taken into consideration to analyze the fluid flow characteristics in shale gas reservoirs. [5][6][7] However, it is a challenge to quantify the flow characteristics in a complex network fracture for evaluating the production.
In the process of hydraulic fracturing, bi-wing fractures, complex fractures, or network fractures would be formed. The conventional analytical or empirical method is no longer applicable to calculate the performance of multistage fractured horizontal wells. In the conventional analytical model, hydraulic fractures are often simplified as bi-wing for the convenience of calculation. That all the lengths of the fractures are equal and perpendicular to the horizontal direction is one of the key assumptions. In the actual situation, both the simple planar fractures and network fractures are identified by the microseismic monitoring. To describe network fractures, discrete fracture models have been widely used to simplify the complex fracture network, including pre-existing natural fractures and hydraulic fractures. 8,9 Jiang et al 10,11 developed a novel hybrid fracture model (UDFM-MINC) that effectively integrates the discrete fracture-matrix model (DFM) for simulating the SRV of shale formations. The discrete fractures model (DFM) can better solve the noncoupling problem of gas flow between matrix and fractures. Jia and Cheng et al 12,13 built a numerical discrete fracture model for simulating the gas flow behavior in multiscale fracture networks. In this model, large-scale fractures have complex geometry and high conductivity, and a small-scale fracturematrix system provides the gas source.
The conventional assumption of fluid flow is the gas seepage from the matrix to microfractures, and then to hydraulic fractures. However, gas flows from the matrix to both microfractures and hydraulic fractures simultaneously. Therefore, a quadrilinear flow model was proposed to improve the assumption. 14,15 Considering fracture interferences between primary fractures and secondary fractures, Chen and Liao et al 16,17 proposed a new semianalytical model to analyze transient pressure of SRV areas in horizontal wells. For describing the heterogeneity of shale reservoirs, the fractal theory has been introduced to account for the complexity of fracture networks. 18,19 Wang and Su et al 20 set up a semianalytical fractal model coupled with a tri-linear diffusive equation to diagnose the flow regime and the size of SRV. Ren and Guo established a novel model with fractional Darcy's law to calculate the anomalous diffusion in stochastic natural fractures. 21 Due to the accessibility and accuracy of multilinear flow models, plenty of improved models have been introduced. 22 Sureshjani and Clarkson broke the flow model of SRV into five parts for analyzing the complex fracture geometry and the straight-line behavior of transient linear flow, which might interfere with the production of SRV. 23 Zhang and Guo et al 24,25 extended a tri-linear flow model and quadruple-porosity model to simulate well production performance of fracture network systems in shale gas reservoirs. With considering different fracture intensity of SRV, Yuan and Su divided the multilinear flow of horizontal wells into seven contiguous-flow regions, which completed within SRV in tight gas reservoirs. 26 Deng et al 27,28 built an enhanced five-region model considering of size of SRV, permeability of SRV, fracture number and fracture intervals. They divided the SRV into different permeability regions, which is in better accord with the morphology of fracture networks. Jia et al 29,30 proposed an improved model to describe the tri-linear flow model considering fracture-production interaction and property heterogeneity of SRV.
Because of the influence of induced stress and reservoir properties, the lengths of the simple planar fractures are different, and part of the fractures are not symmetrical about the horizontal well. In order to describe pressure behaviors of multiple fractured horizontal wells (MFHW) in shale gas reservoirs with more accuracy, lots of semianalytical models have been set up to study the influence of finite-conductivity fractures, different fracture lengths, different fracture intervals, different angles between the fractures, and the horizontal well. 31,32 By dividing the area of multiple fractured horizontal wells into hydraulic fractures, SRV, and unstimulated reservoir, a composite mode is presented to describe the gas flow characteristics of shale reservoirs. 33,34 According to the arrival time of pressure to the SRV outer boundary, the size of SRV can be determined. 35,36 For studying the performance of multiple fractures, different from the familiar point source function method, a novel slab-source function has been applied to describe the flow characteristics. 37, 38 Zhao and Pu et al 39,40 created a new production performance model of complex fractures in fractured tight sandstone gas reservoirs and studied the influence of fracture angles and asymmetric hydraulic fractures in SRV on production. Zhu et al 41 proposed a numerical simulation of thermal recovery based on hydraulic fracture heating technology in shale gas reservoir. Behmanesh et al 42 established an analytical modeling to simulate the linear flow in single-phase tight oil and tight gas reservoirs.
In this paper, we proposed a novel semianalytical multiregion fluid flow model which considered the triple-porosity and dual-permeability characteristics in shale gas reservoirs. In addition to the effects of different morphology of SRV, such as the length and width of network fractures, network fracture number, fracture intervals and the asymmetry distribution of SRV on pressure response have also been figured out. Comparing with the computation complex and time-consuming numerical simulation method, reservoir parameters could be fit and production performance could be analyzed in a simple and accurate way according to the proposed semianalytical model and SRV parameters from microseismic interpretation. Finally, a case study was proposed to validate the accuracy of the research method established in this paper.

MATHEMATICAL MODELS OF SHALE GAS
Shale gas is a type of unconventional natural gas whose source and reservoir are in the same formation. According to the reservoir geological characteristics and the flow mechanisms of shale gas reservoirs, a triple-porosity/dual-permeability shale gas flow physical model is developed (Figure 1). Because of this triple-porosity medium, the flow process of gas in shale reservoirs can be divided into three stages as desorption, diffusion, and seepage. Gas desorbs from the matrix and then diffuses in organic holes, and finally the desorbed gas and free gas flow through microfractures to the bedding.
The higher the organic matter content, the more development of organic matter there is and, the larger the porosity. Although shale matrix permeability is very low, as a result of the existence of microfractures, bedding, and nature fractures, there is still a larger production rate in shale gas reservoirs. After shale gas is desorbed from organic hole surfaces, desorbed gas, and free gas flow to the microfractures and natural fractures in a steady stream. Then, shale gas flows out from microcracks into the wellbore through hydraulic fractures.
The actual complex pore structure was simplified into an ideal capillary model, and the assumptions of the mathematical model are as follows: 1. The desorption process of shale gas follows the rules of Langmuir isotherm adsorption, and the initial state on the isothermal adsorption curve. 2. There is no temperature change in the process of desorption, diffusion, and seepage.
3. Shale gas diffuses directly from the matrix to the microfractures, and the diffusion status is unsteady diffusion. 4. The flow of gas in microfractures is radial laminar flow that follows Darcy's seepage law. 5. Shale reservoirs are homogeneous and incompressible; temperature, porosity, and permeability are constant and isotropic. 6. Only single-phase gas flow is considered in shale reservoirs, regardless of the influence of gravity and capillary force. 7. The production of a well is a constant.
Similar to the derivation process of the traditional dual medium seepage model, the triple-porosity/dual-permeability model in shale gas reservoirs is presented in Appendix A.
By dimensionless transformation, the shale gas flow mathematical model of dimensionless transformation is as follows: Gas diffusion in matrix Gas flow in microfracture Gas flow in hydraulic fracture

Initial conditions
Inner boundary conditions Organic pore, microfracture, and bedding in shale

Outer boundary conditions
where r D is the dimensionless radius defined by r D = r ; storage ratio in microfracture mf = mf c mf ∕ ; storage ratio in hydraulic fracture hf = hf c hf ∕ ; dimensionless pressure:

| Different width and length of stimulated reservoir volume
Unlike regular fracturing wells, the number of artificial fractures is greater and the spread of fractures is more complex in shale gas reservoirs. When gas starts to flow from the shale reservoir into the horizontal well through stimulated reservoir volume, the flow stage and characteristic is more complex. Different from the conventional single cluster perforation multistage fracturing in sandstone reservoirs, multicluster perforation multistage fracturing is widely used in shale reservoirs. Because of the influence of induced stress and reservoir properties, the length and width of stimulated reservoir volume (SRV) are also not the same (Figure 2). In order to calculate more conveniently, the mirror reflection method should be used. With the mirror reflection method, closed stimulated reservoir volume could be regarded as the superposition sum of infinite point source in the same type on the longitudinal direction.
The basic solution of instantaneous point source function in closed boundary has been derived as following, The gas diffusion in matrix is described in Equation (1), and then, the dimensionless function is as follows Initial conditions Inner boundary condition

Outer boundary conditions
As M = r iD c D , then function (7) can be transformed into By Laplace transformation, function (11) can be transformed into The general solution of function (12) is where sinh is hyperbolic sine function, sinh (x) = e x −e −x 2 ; cosh is the hyperbolic cosine function, cosh (x) = e x +e −x 2 . Based on the initial and boundary conditions, the coefficients A and B are, So function (13) (26) is the exponential equation and the calculation of infinite series. The common mathematic method can be hard to solve it; therefore, this function must be simplified through some mathematic methods. Poisson's superposition formula is similar with the basic solution of instantaneous point source function in closed boundary, In the left of Poisson's formula, by multiplying t  Assuming production is constant, the center of the fracture is (0, 0, Z e ∕2), the length of SRV is 2L f , the width of SRV is h. Integral point source basic solution for the X direction, Laplace solution of bottom hole pressure response function of vertical fracture is Integral Equation (7) for the X direction, Laplace solution of bottom hole pressure response function of SRV is where (38) Microseismic diagram and symmetry distribution of stimulated reservoir volume

(B) SRV (A) Microseismic
The mathematical model of multistage fractured horizontal wells with stimulated reservoir volume can be obtained by applying the superposition principle (as shown in Figure 3). Figure 3 reflects that the first fracture is affected by the second fracture, the third fracture, the fourth fracture, and until to the n'th fracture. In addition, the first fracture will affect the second fracture, the third fracture, the fourth fracture, and until to the n'th fracture. By that analogy, all the fractures interact with each other. The pressure and productivity can be obtained by the following expression.
These functions can be written into equation, which is a matrix form. where If the production of different stages has been determined, the bottom hole pressure response model of SRV can be calculated.
Duhamel's principle can be used to solve these models considering the wellbore storage and skin effect because the limit of the point source function. Assume that well production is q D (t D ), the corresponding wellbore pressure response function P WD (t D ) can be obtained through convolution method.
where P wD is dimensionless wellbore pressure, q D (t D ) is dimensionless rate, S is skin factor. According to the definition of wellbore storage factor, q D (t D ) can be written as follows: Laplace transformation for function (44) and (45): The expression between dimensionless wellbore pressure and the dimensionless wellbore pressure affected by the wellbore storage and skin effect can be obtained: where u is lagrangian, C D is dimensionless wellbore storage coefficient.

| Asymmetry distribution of stimulated reservoir volume
Due to the influence of anisotropy reservoir properties and the in situ stress field, some of the SRV are not symmetrical about the horizontal well ( Figure 4).
In the actual hydraulic fraction process, these test data are disorder. For the evaluation of fracturing results, microseism analysis can be used (as shown in Figure 4A). Each cluster fractures have its microseismic signals, which is helpful to analyze the fractures' characteristics. By classifying different microseismic signals, the asymmetry distribution of simulated reservoir volume model can be simplified (as shown in Figure 4B).

| Gas flow characteristics
Flow pattern of a multistage fractured horizontal well with SRV can be divided into eight flow stages (as shown in Figure 5): (a) wellbore storage stage, which the curve of pressure and pressure derivative curve is superposition and the slope of them is 1; (b) transition flow stage, which reflects the degree of pollution near the bottom; (c) linear flow in SRV, in which the slope of pressure and pressure derivative curve is 0.5; (d) the first radial flow stage, which gas flows to SRV in an ellipse shape; (e) the second radial flow stage in which the slope of pressure and pressure derivative curve is 0.5; (f) gas flow in the microfractures, which the slope of pressure derivative curve is relatively smooth; (g) the transition stage when gas desorbs and diffuses from matrix to microfractures; and (h) comprehensive radial flow stage, which reflects that the radial flow of whole reservoir pressure reaches balance, and the pressure derivative curve appears as a flat line.
To validate the accuracy of the mathematical model proposed in this paper, we established an embedded discrete fracture model (EDFM), which coupled the fluid flow in matrix and fractures. In this paper, the MRST open source program and MATLAB programming were used to simulate the production of shale gas reservoirs. A shale gas reservoir model was established, and these basic parameters are listed in Table 1. The distribution diagram of the complex fracture network is shown in Figure 6.
After the simulation is completed, the pressure-time curve and pressure-pressure derivative curve were plot, as shown in Figures 7 and 8.
In Figure 11, we can get that the shape of the pressure & pressure derivative-time curve is coincided to the flow pattern of multistage fractured horizontal well with SRV, which indicated the mathematical model proposed in this paper is correct. Figure 9 shows the influence of space between SRV on pressure characteristics. The greater the spacing is, the more obviously the transition of linear flow in SRV and first radial flow is. With the increase of SRV spacing, the dimensionless pressure would be smaller, which reflects the lesser pressure drop loss in the process of seepage. Figure 10 shows the influence of width of SRV on pressure characteristics. The width of SRV mainly affects the early fracture linear flow stage. The bigger the width of SRV is, the higher the dimensionless pressure is, and the closer the slope of the pressure derivative is to 0.5. Figure 11 shows the influence of s length of SRV on pressure characteristics. The bigger the length of SRV is, the longer the duration time of fracture linear flow is. At the same time, dimensionless pressure reduces and the duration time of the dimensionless pressure derivative curve increases gradually.

| CASE STUDY
Wei 201-H1 is the first shale gas horizontal well in China, which is located in the Longmaxi formation. To illustrate the SRV method is useful, in this paper, we chose this well's fracturing data to conduct the case analysis. The well test interpretation results were shown in Figures 12 and 13. The well test interpretation results indicated the shale gas reservoir is homogeneous, the fracture type is vertical and the number of vertical fractures is 11, and the skin factor is 0.0279.

| CONCLUSIONS
In order to evaluate the shale gas reserves and design a development plan, a simple and accurate method is needed to predict the production and performance of the multistage fractured horizontal wells. In this paper, a triple-porosity/ dual-permeability model has been developed for simulating F I G U R E 1 2 Log-log curve of pressure F I G U R E 1 3 Pressure-time curve the flow mechanism of shale gas. The conclusions are as follows: 1. The basic solution of point source function of a gas flow model with complex fractures system can be solved by applying a combination of mirror reflection theory, the integral of Bessel function and Poisson's superimposed formula. 2. Type curves of dimensionless pressure, pressure derivative, and dimensionless production are calculated and plotted as log-log type curves, which could reflect the characteristics of gas flow in a multistage fractured horizontal well with stimulated reservoir volume. 3. Flow pattern of a multistage fractured horizontal well with SRV can be divided into eight flow stages. Influence factors on gas flow characteristics (length and width of a fracture network, intervals between two fracture clusters, asymmetry distribution) and production declination of shale gas have been analyzed and discussed.
Based on these research achievements, in the later study, we plan to study the gas flow characteristics by using numerical simulation method. According to an actual shale gas reservoir, we will establish an embedded discrete fracture model (EDFM) to simulate the actual production process in shale gas reservoir. By changing fracture shapes with different levels of complexity, the daily gas production and tired gas production will be compared on the basis of MRST open source program and MATLAB programming. These simulated results will be different under different fracture shapes, which will provide a better guidance for the optimum selection for well fracturing design programs.