Pressure transient analysis for multi-wing fractured vertical well in coalbed methane reservoir based on fractal geometry and fractional calculus

In this paper, a semi-analytical mathematical model of pressure transient analysis (PTA) for multi-wing fractured vertical well (MWFV) in coalbed methane (CBM) reservoir is proposed, which considers the complexity of porous media by fractal geometry, the anomalous diffusion based on fractional calculus and the stress sensitivity represented by the exponential expression. Then through line source theory, dimensionless transformation, Pedrosa transformation, and other methods, the solution of the bottom hole pressure is obtained. At last, the PTA curve is presented by the Stehfest inversion method, and seven flowing regimes are identified. It can be observed that after introducing fractal geometry and fractional calculus, the PTA curve is quite different from the traditional curve in the slope of the derivative curve. The influences of related parameters are discussed, including mass fractal dimension, anomalous diffusion coefficient, number of hydraulic fractures, fracture angle, and stress sensitivity factor. Relevant results can provide better guidance to understand the CBM production performance in MWFV.


Introduction
It was reported in BP ENERGY OUTLOOK in 2016 that to maintain the stability and increment of the global economic (Wei et al., 2021a), the energy demand will be in dramatic growth and fossil fuel will account for 80% of the whole energy supplement on a global scale (Wei et al., 2021b). For the aim to alleviate the burden of traditional petroleum energy supply, the unconventional natural gas resource is gradually getting people's attention (Nie et al., 2012). More importantly, the unconventional natural gas resource has an abundant distribution worldwide. For example, in China, a large number of unconventional natural gas resources such as tight sandstone gas (Sulige Gas Field) and coalbed methane (Taiyuan Formation) are developed in the Ordos Basin. Also, in the United States, numerous tight gas resources are stored in the western part of the Rocky Mountains, and in the San Juan Basin, there is also a large amount of coalbed methane (Miskimins and Jennifer, 2009). Coalbed methane (CBM), which is a member of the family of unconventional natural gas resources, is beginning to get much exploitation due to its high thermal efficiency and low pollution (Dong et al., 2017).
The permeability of CBM reservoirs is usually low and vertical fracturing technology is widely introduced to the exploitation of CBM to achieve high productivity . Meanwhile, during the fracturing period, multiple fractures are formed ( Figure 1) in the effect of stress instead of ideal single planner fracture which is known as multi-wing fractured vertical well (MWFV).
Pressure transient analysis (PTA) is an efficient method for the monitoring of the performance of production well, and many researchers studied the PTA curve of MWFV. Choo and Wu (1987) researched MWFV in the gas reservoir. Through numerical method, they obtained the PTA curve and noticed that there is an obvious difference in the PTA curve between multi-wing fractures and single fracture when the dimensionless conductivity of fractures is small than 10. But the numerical method is complicated in representation and time-consuming in the calculation, a semianalytical method was proposed later according to the research achievement of Cinco-Ley H (Cinco et al., 1978;Cinco and Meng, 1988). Craig and Blasingame (2006) constructed a vertical fractured well model with two fractures (one is secondary fracture) and proposed the solution of bottom hole pressure at the constant production rate. They analyzed the impacts of the length of two fractures, Figure 1. Multiple fractures connected to the vertical wellbore (Chen et al., 2016). the conductivity of fractured, and the angle between two fractures. Luo and Tang (2015) studied the PTA curve of MWFV in detail based on the semi-analytical method. They found that there is a slight 'hump' in the pressure derivative curve which is called 'interference regime'. Research on the PTA curve of MWFV in CBM reservoirs is also within the scope of researchers. Chen et al. (2016) established the corresponding model for MWFV in CBM reservoirs which took absorption, diffusion, and stress sensitivity into consideration. They also drew the corresponding PTA curve and analyzed the influence of various factors in detail. Zhang et al. (2018) proposed a model for MWFV in CBM reservoirs which coupled the stimulated region through the composite model. A transition regime in the PTA curve between the stimulated region and unstimulated region was put forward and the radius of the stimulated region was analyzed.
However, all research mentioned above is based on a uniform research background: the Euclidean geometric space, in which the reservoir is a homogeneous porous media and natural fractures (or cleats in CBM reservoirs) are distributed uniformly. While this is opposite to the fact that the distribution of cleats is extremely random in reservoirs (Wang et al., 2016). That is to say, the complexity of reservoirs is ignored in the above research. In order to modify this problem, the fractal theory is applied to revise. Based on fractal theory, Chang and Yortsos (1990) derived the formula of permeability and porosity (fractal permeability and porosity) which are the powerlaw function of reference distance. This method can better reflect the complexity of porous media and the reservoirs described is closer to reality and is widely used in much research. Through fractal permeability and porosity, Beier (1994) conducted research on the PTA curve of a vertical fractured well and found this method can match the field data better. Flamenco-Lopez and Camacho-Velazquez (2001) investigated the PTA curve of naturally fractured reservoirs which took fractal characteristics into consideration. Cossio et al. (2012) utilized the tri-linear model (Lee and Brockenbrough, 1986) to study the pressure performance of the fractal vertical fractured well. Also, Wang et al. (2015aWang et al. ( , 2015b introduced fractal theory into the tri-linear model to represent the stimulated region in multi-staged horizontal wells and analyzed the properties of the stimulated region on the PTA curve. Wang et al. (2019) coupled the fractal theory with the fivelinear model (Stalgorova and Mattar, 2013) to study the imbibition effect on production performance of the multi-staged horizontal well.
Meanwhile, it is approved that due to the extremely random distribution of natural fractures, the characteristic of transport in porous media is transformed into anomalous diffusion (Metzler et al., 1994;Park et al., 2000). Aimed at anomalous diffusion, fractional calculus is proposed to describe this phenomenon, that is to say, modifying the darcy equation by fractional calculus. This method is widely used in the study of PTA considering anomalous diffusion (Raghavan, 2012;Raghavan and Chen, 2013;Razminia et al., 2015). In recent years, Ozcan et al. (2014) combined anomalous diffusion with the tri-linear model to obtain an analytical solution of a fractured horizontal well in tight formation and suggested that the new tri-linear anomalous-diffusion model is useful to make a transient analysis. Aimed at shale gas reservoirs, Ren and Guo (2015) studied the pressure response of multiple fractured horizontal wells considering anomalous diffusion, but the conductivity of hydraulic fractures was not accounted. Gu et al. (2019) derived a FFDM model (fractally fractional diffusion model) to study the PTA curve of multiple fractured horizontal wells (infinite conductivity) with the stimulated region. Jiang et al. (2021) established the PTA curve for multi-fractured horizontal wells in CBM reservoirs based on fractal and fractional calculus, but they also assumed that hydraulic fractures are infinite conductivity.
Through extensive literature research, lots of work had been done on the pressure transient of MWFV and fractal reservoir, but there is no study about the PTA research of MWFV that both consider the fractal geometry, anomalous diffusion, and stress sensitivity in CBM reservoir. In this paper, we derive a corresponding mathematical model based on fractal theory and fractional calculus and conduct some relevant analysis.
The following of this paper is arranged as: 1. The physical model of MWFV in CBM reservoirs is introduced and some assumptions are made to simplify this problem. 2. Based on the physical model, the mathematical model is established, and then the mathematical model is solved to obtain the pressure response of MWFV. 3. The verification of the results of this paper is conducted and a detailed analysis of the PTA curve is performed.

Physical model
The physical model for MWFV in coalbed methane is illuminated in Figure 2. The basic hypotheses are: 1. During the fracturing period, M hydraulic fractures centered on the wellbore are formed, which penetrate the reservoir and each of which has the same length (L F ), permeability (K F ), and with (W F ). The pressure drop along the fracture is considered. The gas firstly flows into the fractures, then flows into the wellbore through fractures. 2. The MWFV produces at the constant rate (q sc ), the whole reservoir shares the initial pressure (p i ). 3. The cleat system which is embedded in the Euclidean matrix system is represented as fractal geometry to consider the complexity of porous media, and the flow in the cleat system obeys anomalous diffusion. The gas flows from the matrix system to the cleat system by pseudo diffusion. 4. Ignoring the temperature change, capillary force, and gravity. Only a single gas phase is considered.

Establishment and solution of mathematical model
Model establishment The whole system consists of two elements: reservoir and hydraulic fractures, therefore, the model establishment also includes two parts.
Model of reservoir. According to Appendix B, the flowing equation of the gas in the cleat system considering fractal geometry, anomalous diffusion, and stress sensitivity can be formulated as: Based on the dimensionless variable listed in Appendix A, the dimensionless form of equation (1) is: The gas from the matrix system to the cleat system is described by Fick's first law (Xu et al., 2015): The corresponding initial condition, outer boundary condition, and inner boundary are: Model of hydraulic fractures. There are M hydraulic fractures in total and have the same properties, so it is reasonable to write the flowing equation for the i-th fracture only, and the other fractures' flowing equations are the same. According to Appendices A and B, the dimensionless flowing equation can be described as: The relevant boundary conditions are: Model solution Solution for reservoir model. Equation (2) is strongly nonlinear due to the consideration of the stress sensitivity, therefore, Pedrosa transformation and Perturbation transformation are adopted to handle this problem. Pedrosa transformation is shown as: Through equation (10), equations (2)-(6) is transformed into: Then, according to Perturbation transformation: Due to the tiny value of γ pD , equation (16) can be also formulated as: So, equations (11)-(15) can be furtherly written as: ( At last, taking Laplace transformation (in regard to t D ) to equations (18)-(22) and combing equations (18) and (19): Where, Making n = d − θ−1, m = d − 1, equation (23) is reshaped to: The general solution of equation (27) is: Taking the outer boundary equation (24) and inner boundary equation (25) into equation (29), the unknown coefficient a and b are: Where, Finally, Equation (32) is the line-source solution for the reservoir model.
Solution for hydraulic fractures model. The method to solve hydraulic fractures model is analogous to solution for reservoir model, after Pedrosa transformation, Perturbation transformation, and Laplace transformation, Equations (7)-(9) are transformed into: Substituting equation (35) into equation (33) and taking double integral: The double integral in the left side of equation (36) is: Then, equation (36) can be written as: Assuming that the pressure is equal at the interface between fractures and reservoir, equation (38) can also be expressed as:

Pressure response of MWFV
To obtain the pressure response of MWFV, discretization and superposition principles are utilized to cope with this problem. Each hydraulic fracture is split into N segments, as shown in Figure 3. Therefore, the length of each segment can be calculated by the length of the fracture, and the midpoint coordinate of the j-th segment in the i-th fracture can also be calculated. Applying the superposition principle, the total pressure response at any point r D or (x D , y D ) can be obtained by superimposing pressure responses of all segments at that point: By the relation between radial coordinates and cartesian coordinates, equation (40) can be written as: Where, Point (x D , y D ) can be placed at the midpoint of the segment, and equation (40) is further written as: Combining partial integration and numerical integration, equation (43) is finally transformed into: By taking all midpoints into equation (44) And through Duhamel's principle (Van Everdingen and Hurst, 1949), the wellbore storage and skin effect are taken into consideration: Then by Stehfest numerical inversion (Stehfest, 1970), we can get the solution in the real-time domain.

Verification
In order to verify the accuracy of the fractal solution of this paper, a comparison between the fractal solution (solved by the semi-analytical method in this paper) and the numerical solution (solved by COMSOL). As shown in Figures 4 and 5, the red solid line has good consistent with the green hollow dotted line and the relative error is acceptable, which approve that the fractal solution is reliable. Also, making d = 2, θ = 0, the model in this paper can degenerate to the model analogous to Ren's model. As shown in Figure 6, there is a good agreement between these two results, which is explaining the correctness of the model in this paper.

Typical PTA curve
According to different features of the PTA curve, seven flowing regimes are split: I: Wellbore storage and skin effect regime. II: Bi-linear flowing regime. The slope of the derivative curve is bigger than 1/4 compared to the traditional model. III: Interference regime. A slight hump in the derivative curve is a feature in this regime. IV: Linear flowing regime. The derivative curve tends to have a slope greater than 1/2, which is also different from the traditional model. V: Transition regime. VI: Diffusion regime. For pseudo steady diffusion, there is an obvious dip in the derivative curve, which illustrates that CBM is flowing from matrix to cleat under the action of concentration difference. But for transient diffusion, there is also a dip which is less obvious than that if the diffusion is pseudo steady. VII: Later radial flowing regime. After considering fractal geometry and anomalous diffusion, the derivative curve is no more a horizontal line. (Figure 7) Discussions Mass fractal dimension and anomalous diffusion index. As shown in Figures 8 and 9. Mass fractal dimension and anomalous diffusion index have an opposite impact on the PTA curve. When the mass fractal dimension becomes smaller, the position of curves in the Bi-linear flowing regime, Interference regime, and Linear flowing regime become up, while the position of curves in the Diffusion regime and Later radial flowing regime becomes down. Instead, the PTA curve shows an inverse phenomenon as the anomalous diffusion index increase. In addition, the increment of anomalous diffusion index also makes flowing regimes after flowing regime I delay.
Number of hydraulic fractures. In Figure 10, the impact of the number of hydraulic fractures is mainly in Bi-linear flowing, Interference, and Linear flowing regime. First, the position of the PTA curve in these regimes goes down as the number of fractures increases. Additionally, the hump of the derivative curve in the Interference regime becomes more obvious, which means more fractures represent worse interference.
Hydraulic fracture length. With the increment of hydraulic fracture length, the stimulated region becomes larger, which means the energy wastage required for gas flowing becomes smaller. So, it can be observed from Figure 11 that a bigger fracture length makes the position of the PTA curve go down in Bi-linear flowing, Interference, and Linear flowing regimes. Meanwhile, the duration of Bi-linear flowing is longer, and the duration of the Linear flowing regime is the opposite.  Angle between adjacent hydraulic fracture. When the angle between adjacent hydraulic fractures decreases, the Interference regime appears early and the hump in the derivative curve is more obvious (Figure 12). A smaller angle represents that the adjacent fractures are very close to each other, which means the interference between adjacent fractures is easier to appear and the greater energy loss in the whole production. Therefore, the relevant phenomenon is formed.
Hydraulic fracture conductivity. In Figure 13, it can be seen that with the increment of fracture conductivity, the duration of the Bi-linear flowing regime becomes shorter, but the duration of the Linear flowing regime becomes longer. If the value of fracture conductivity is large enough, which means the fracture has infinite conductivity, the Bi-linear flowing regime will be covered up, leaving the Linear flowing regime only.
Stress sensitivity factor. With the increment of stress sensitivity factor, the PTA curve in Later radial flowing regime upwarps (Figure 14). This is because the stress sensitivity effect will be reflected  when the pressure difference increases to a certain extent, which leads to the increment of flowing resistance and energy wastage and is not conducive to production. Therefore, it shows the above phenomenon.

Field application
In this part, one field application is conducted based on the model above. Well X3 is an MWFV of CBM reservoir in China. Through the micro-seismogram result, there are five fracture wings formed during the fracturing period. The matching of the PTA curve is shown in Figure 15 and relevant interpretation results are listed in Table 1. It should be noted that due to the short time of shut-in, the whole flowing regimes are not complete. Also, the interpreted adsorption coefficient, transfer coefficient, and storage coefficient are not necessarily accurate because the Diffusion regime is not identified.

A new mathematical model for multi-wing fractured vertical well in coalbed methane reservoir
is established which takes the complexity of porous media, anomalous diffusion, and stress sensitivity into consideration. Then the solution is obtained by a semi-analytical method. 2. The corresponding PTA curve is drawn and divided into seven different flowing regimes. The main difference compared to the traditional one is reflected in the slope of the derivative curve. The slope is bigger than 1/4 in the Bi-linear flowing regime or 1/2 in the Linear flowing regime. And the derivative is no more a horizontal line in the Later radial flowing regime. 3. Sensitivity analysis reveals that mass fractal dimension and anomalous diffusion index have a great influence on the whole PTA curve, also the increasing anomalous diffusion index leads to the delay of flowing regimes after regime I. Parameters related to the geometry of hydraulic fractures mainly affect the middle flowing regime. Stress sensitivity factor cause the upward of later flowing regimes.      (3) 4. Field application shows that based on the proposed model, the PTA curve has a good fitting result which represents the practicality of this paper's model.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Key S&amp;T Special Projects, The National Natural Science Foundation of China, (grant number 2017ZX05008-004-004-001, 41972129).

ORCID iD
Haoyuan flux density of the j-th segment in the i-th hydraulic fracture, m 3 /(s·m) q sc production rate of the multi-wing fractured vertical well, m 3 /s r radial distance, m S skin factor, dimensionless t time, s T temperature, K s Laplace transform variable, dimensionless x, y Cartesian coordinates, m x w , y w Cartesian coordinates of the line-sink, m V L Langmuir volume, m3/m3 Z compressibility factor, fraction β adsorption coefficient for matrix, dimensionless γ permeability modulus, Pa-1 ω storage coefficient, dimensionless λ transfer coefficient, dimensionless θ anomalous diffusion index, dimensionless α the order of fractional derivative, dimensionless μ viscosity, Pa·s ϕ porosity, fraction σ perturbation deformation function r L ref β = p sc q sc T πk cri hT sc ψ L ψ i (ψ L + ψ c )

Reservoir model
In the radial coordinates, the gas flowing in the cleat system succumbs to the continuity equation: Here, the fractional calculus is introduced to describe the anomalous diffusion (Ren and Guo, 2015): The pseudo-pressure is: Taking equations (49) and (50), the continuity equation is: The complexity of the porous media is represented by the fractal permeability and porosity proposed by Chang and Yortsos (1990): And permeability modulus based on the pseudo-pressure is introduced to describe the stress sensitivity: