Transient flow of a horizontal well with multiple fracture wings in coalbed methane reservoirs

: Horizontal wells with multi-stage fractures have been widely used to improve coalbed methane (CBM) production from coalbed methane reservoirs. The main focus of this work is to establish a new semi-analytical method in the Laplace domain and investigate the transient pressure behavior in coalbed methane reservoirs. With the new semi-analytical method, flow regimes of a multi-fractured horizontal well in coalbed methane reservoirs were identified. In addition, the sensitivities of fracture conductivity, diffusion model, storability ratio, inter-porosity flow coefficient, adsorption index, fracture spacing, fracture asymmetry, non-planar angle, and wellbore storage were studied. Results indicate that six characteristic flow regimes can be identified for multi-fractured horizontal wells in coalbed methane reservoirs, which are bilinear flow, first linear flow, desorption-diffusion flow, first pseudo-radial flow, second linear flow, and second pseudo-radial flow. Furthermore, the sensitivity analysis shows that the early flow is mainly determined by the fracture conductivity, the asymmetry factor, the non-planar angle, and the wellbore storage; while the desorption-diffusion flow regime is mainly influenced by the diffusion model, the storability ratio, the inter-porosity flow coefficient, the adsorption index, and the fracture spacing. Our work can provide a deep insight into the fluid flow mechanism of multi-fractured horizontal wells in coalbed methane reservoirs.


Introduction
Organic fossil energy resources, such as the oil and gas reservoirs, are widely distributed in the world. As a gaseous hydrogen compounds energy source, the development of coalbed methane (CBM) is not only an important supplement to clean and sustainable energy, but also helps to reduce the risks of coal mine development.
With the development of hydraulic fracturing technology, unconventional gases, such as coalbed methane and shale gas, have become major energy sources in recent years [1][2][3][4]. In China, CBM reserves are about 35.81×10 12 m 3 , accounting for 14% of the global CBM reserves [5]. In 2015, 18 billion cubic meters of CBM were produced in China, and most of them were from Shanxi province, Guizhou province, Anhui province, and Henan province [6]. From 2011 to 2015, as many as 11300 new wells were drilled in China for the development of CBM. It is reported that most of the newly drilled wells are vertical wells with a single fracture. However, a few multi-stage fractured horizontal wells have also been drilled in some pilot areas to evaluate the capability of horizontal wells for improving the production of CBM. It was shown that the production rate of horizontal wells with 6-10 fractures was twice that of vertical wells [7][8][9].
Transient pressure analysis is a powerful method, and it is usually used for evaluating fracturing process [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. Cinco-Ley et al. [10] developed the semi-analytical solution to study the transient pressure behavior of a vertical fracture. The bilinear flow and linear flow periods were identified. They also extended this method in dual porosity reservoirs. Also, many researchers have extended the semi-analytical method to a vertical wellbore with multi-wings [11][12][13], horizontal wells with multiple fractures [14][15][16][17][18][19][20][21][22][23][24]. Xing et al. [25] proposed a new model to evaluate the productivity of a well under a pseudo-steady state in an anisotropic reservoir, and found that the optimized well productivity occured when the fractures were perpendicular to the principal permeability. Wang et al. [26] provided a new method to evaluate the production performance of a horizontal well considering the stress sensitivity in reservoir. Their research showed that the stress sensitivity mainly influences the late-flow period. Based on the fracture-wing method and nodal analysis technique, Xing et al. [27] proposed a new semi-analytical model to solve the transient flow problems of horizontal wells with multiple reorientation fractures. However, these models mentioned above focused on transient pressure behaviors of multiply-fractured horizontal wells in conventional reservoirs instead of coalbed reservoirs.
Traditionally, we generally consider oil-gas flow and coalbed methane flow as isothermal flows, which means the temperature in the reservoir is constant. But, in fact, isothermal flow is only a hypothetical approximate treatment method. Some scholars have studied non-isothermal flow and heat conduction. Venart et al. [28] made a comprehensive study of the Knudsen-effect errors of the transient line-source method in fluids, which provides guidance for simultaneous measurement of low-density thermal conductivity and thermal diffusivity values. Alata et al. [29] established a hyperbolic heat conduction model to study the thermal behavior of thermoelectric generators and refrigerators, and presented the transient temperature distributions under different parameters. Khadrawi et al. [30] also developed a hyperbolic heat conduction model to investigate the transient hydrodynamics and thermal behaviors of fluid flow in a micro-channel, which is assumed to be openended and vertical parallel-plate. The effects of Knudsen number and thermal relaxation time on the transient behaviors were discussed in detail. Minea et al. [31] developed a comprehensive and applicable simulation method of nanofluid flow, and presented comparisons of different simulation approaches. The results showed that gravity is a key factor that must be considered in the simulation of nanofluid flow.
The flow in coalbed methane reservoirs is more complex than that in conventional reservoirs, since the coalbed reservoirs have dual porosity consisting of natural fractures and matrix. The free gas is stored in the natural fractures while most gas is stored in the matrix as absorbed state, which can be characterized with the Langmuir adsorption law. The coalbed methane flow mechanism between the matrix system and the natural fractures is always described with the law of diffusion. Two diffusion models for the coalbed methane, the pseudo-steady diffusion and transient diffusion model, have been used in the literature [32][33]. Anbarci and Ertekin [34][35][36] presented a comprehensive study on the transient pressure analysis for the CBM flow in unfractured wells and fractured wells with infinite conductivity.
Guo et al. [37] developed a three-dimensional CBM numerical reservoir simulator to model the flow mechanism of CBM reservoir and predict its production performance. Clarkson et al. [38] proposed new analytical approaches and workflows for investigating the transient pressure of CBM from the vertical, hydraulically-fractured wells and horizontal wells under single and multi-phase flow. Aminian and Ameri [39] presented a quick yet reliable tool for predicting the production performance of coalbed reservoirs. The type curves can be used for parametric studies when the key characteristics are not well established. Nie et al. [40] presented an analytical solution to study the transient flow behavior of a horizontal CBM well. Wang et al. [41] established a novel semi-analytical model to study the influence of asymmetric factor on an asymmetrically fractured well in coalbed methane reservoirs. Zhang et al. [42] proposed a new model of a vertical well with multi-wings in coalbed methane reservoirs, and some important properties of the stimulated reservoir volume were analyzed. Chen et al. [43] developed a new 3D point-sink model to investigate the transient production mechanism of horizontal fractures in a coalbed reservoir, and four flow characteristics for the coalbed methane were identified.
As stated above, these studies mainly focused on the flow model in coalbed methane reservoirs with a single fracture, and most of the models are based on the assumption that the fractures are symmetric and coplanar. With the development of new horizontal well and fracturing technology, more and more horizontal wells with asymmetric and non-planar fractures will be drilled in the future. However, few studies have analyzed the transient pressure behavior of a horizontal well with multiple complex fractures in coalbed methane reservoirs.
In this paper, a semi-analytical method was developed for studying the transient pressure behavior of a multi-fractured horizontal well in a coalbed methane reservoir. The fractures can be modeled by random asymmetry factor and non-planar angle. Furthermore, some complicated fractures can be modeled by combinations of fracture wing with different asymmetry factors and non-planar angles in our study. The effects of fracture conductivity, diffusion model, storability ratio, inter-porosity flow coefficient, adsorption index, fracture spacing, fracture asymmetry, non-planar angle, and wellbore storage on the pseudo-pressure behaviors are discussed.

Mathematical Models
The basic physical assumptions are the same as in Wang et al.'s model [41], except that the coalbed methane reservoir is infinitely large and is developed by a horizontal well instead of a vertical well.
The horizontal well is multi-stage fractured with N asymmetric fractures of finite conductivity. Each fracture is divided into 2 wings by the horizontal well. Therefore, there are 2N fracture wings of different wing lengths and fracture conductivities. The w-th fracture wing is discretized into Nw (w=1, 2, …, 2N) segments (see Figure 1).

Fluid Flow in Coalbed Methane Reservoir
With the definition of pseudo-pressure, the pressure response in Laplace domain of the i-th fracture segment in CBM reservoir satisfies the following equation according to the superposition principle [34][35][36]41], In Equation (2), the term f(s) exhibits different forms with respect to the diffusion model used to describe the flow between the matrix and fractures. The transient state diffusion model can be written as [34,41] ( ) The pseudo-steady state diffusion model can be expressed as

Fluid Flow in the Fracture
As illustrated in Figure 2, the fluids first flow into the fracture and then accumulate into the wellbore. The flow in the fracture is one-dimensional. In Appendix A, we derive a new fracture-wing equation for CBM flow in the fracture The Equation (5) can be furthermore written in discretized form. Furthermore, the pressure for the i-th segment of w-th wing is satisfied with [11] ( )

The Semi-analytical Solution
The pressure and flux along the fracture wing must be satisfied with the following continuity In Laplace domain, the unity condition of the flow rates can be expressed as By solving Equations (1), (5), (7) and (8), we get the wellbore pressure in Laplace domain, which could be converted to the real domain through Stehfest numerical algorithm [44].
As a solution in Laplace domain, the dimensionless wellbore pressure swD ψ with wellbore storage CD can be readily calculated as

Model Validation
To the best of our knowledge, few publications have discussed the pressure behavior of multiple finite-conductivity fractures in CBM reservoirs. In order to validate the new model, we compared our solution with that reported by Zerar [19] for finite conductivity fracture in a homogeneous reservoir (f(s) = s). It is shown that our solutions excellently match with Zerar's solutions before the boundary is felt ( Figure 3).  As shown in Figure. 4, we can divide the flow of fluid into eight stages in CBM reservoirs. Stage 1: In this stage, the pseudo-pressure derivative curve (red line) has a straight line with a slope of 1/4, representing the bilinear flow period. This period can only be observed in low-moderate fracture conductivity.

Flow Characteristics Analysis
Stage 2: A transition period occurs from the bilinear flow to the first linear flow. Stage 3: The first linear flow period, which is characterized by a straight line with a slope of 1/2 on the pseudo-pressure derivative curve.
Stage 4: Because of the effect of desorption and diffusion mechanisms, a typical dip followed by the first linear flow on the pseudo-pressure derivative can be observed. This flow period can continue for a long time and is easy to be observed.
Stage 5: The first pseudo-radial flow period will be observed while the pseudo-pressure derivative curve shows a constant as 1/ (2N). In this period, fluids flow radially from the formation to each fracture individually under the condition that the fracture spacing is sufficiently large compared with the half-length of the fracture.
Stage 6: The second linear flow period, which is characterized by a straight line with a slope of 0.36 on the pseudo-pressure derivative curve, can be observed when the fracture interference is felt.
Stage 7: A transition period occurs from the second linear flow to the second pseudo-radial flow. Stage 8: The second pseudo-radial flow period is distinguished by a horizontal line on the pseudo-pressure derivative curve, which shows a stable value of 0.5.
It should be noted that not all flow periods illustrated in Figure 4 exist for all horizontal wells with multiple fractures in coalbed methane reservoirs. It depends on the combinations of parameters, such as fracture conductivity, diffusion model, storability ratio, inter-porosity flow coefficient, adsorption index, fracture spacing, fracture asymmetry, non-planar angle, and wellbore storage. In the following section, the impacts of these parameters on transient pseudo-pressure responses will be analyzed.

Effects of Parameters on Transient Pseudo-pressure Responses
Fracture conductivity (CfD). Figure 5 presents the pseudo-pressure and pseudo-pressure derivative curves for different dimensionless fracture conductivities. As shown in Figure 5, the fracture conductivity has a major effect on the early pseudo-pressure responses. In the case of low fracture conductivity (CfD = 1), the bilinear flow can be identified clearly and the first linear flow can't be observed. With the increase of the fracture conductivity, the duration of the bilinear flow decreases and the first linear flow gradually dominates the fluid flow in the CBM reservoir. At the moderate fracture conductivity (CfD=10 ~ 50), both the bilinear flow and the first linear flow can be identified. When the fracture conductivity is greater than 100, the first linear flow can be recognized easily, whereas the bilinear flow will be absent.  Figure 6 displays the typical pseudo-pressure curves of TSS and PSS models with other identical parameters. As can be seen in Figure 6, the eight stages presented in the PSS model can also be observed in the TSS model. At a large dimensionless time (tD≥10 2 ), both the TSS and PSS model have the same value on the curves of pseudo-pressure derivative. However, during the intermediate time (10 -2 ≤ tD ≤ 10 2 ), the shapes of the pseudo-pressure derivative curve exhibit huge differences. The characteristic dip on the pseudo-pressure derivative is no longer observed with the TSS model. In the early stages, the diffusion model mainly affects the end of both the bilinear flow and the first linear flow. It can be observed the end of bilinear flow and first linear flow for the PSS model occurs earlier than that of the TSS model. Meanwhile, we notice that in the early flow period (tD < 10 -1 ), the values of pseudopressure and derivative for the PSS model are obviously bigger compared with the TSS model, indicating larger pressure depletions in the PSS model.

Transient state (TSS) diffusion and pseudo-steady state (PSS) diffusion.
Storability ratio (ω). Figure 7 presents the influence of storability ratio ω on the transient pseudo-pressure responses with different fracture conductivities. The storability ratio has a significant effect on the early-stage pseudo-pressure responses.  Table 1 displays the pseudo-pressure and pseudo-pressure derivative data for different storability ratios, and it can be seen that in the early and middle flow periods (tD≤10 2 ), the pseudopressure increases as ω reduces, which indicates that a small storability ratio leads to large pressure depletion. For example, when CfD = 1 and tD = 10 -1 , the pseudo-pressure ΨwD is 0.7861 for the case of ω = 0.05, while they are 0.6951 and 0.4948 for the ω = 0.1 and ω = 0.5, respectively. The storability ratio exerts influence on the flow stage of desorption-diffusion. As shown on the derivative curves, the dip is wider and deeper with the decrease in the storability ratio. The influence of storability ratio depends on the fracture conductivity (Figures 7(a)-(c)). For low fracture conductivity (CfD = 1, Figure 7(a)), the storability ratio influences the bilinear flow period, while it affects the first linear flow period for high conductivity (CfD = 500, Figure 7(c)). For moderate fracture conductivity (CfD = 50, Figure 7(b)), both the bilinear and first linear flow periods will be affected. Figure 8 illustrates the effect of inter-porosity flow coefficient λ on the transient pseudo-pressure responses. It indicates that the inter-porosity flow coefficient mainly affects stage 4, the desorption-diffusion flow period. As the inter-porosity flow coefficient increases, the dip of pseudo-pressure derivative curves moves towards the right. With a big interporosity flow coefficient (λ = 10 5 ), two segments of the second pseudo-radial flow separated by the desorption-diffusion flow can be observed on the derivative curves. Adsorption index (α). Figure 9 demonstrates the impact of adsorption index α on the transient pseudo-pressure responses. It can be seen that the adsorption index mainly impacts the flow periods of stage 4, 5, 6, and 7.

Inter-porosity flow coefficient (λ).
Comparing Figure 9 with Figures 7 and 8, the storability ratio, ω, mainly affects the magnitudes of dip on the pseudo-pressure derivative curve, while the inter-porosity flow coefficient, λ, mainly affects the occurring time of the desorption and diffusion. However, the adsorption index, α, exerts the influence on both the magnitude of dip and the occurring time of desorption and diffusion. In addition, the dips become deeper and wider with the increase of adsorption index. It should be noted that the effects of storability ratio and inter-porosity flow coefficient on the pseudo-pressure responses will gradually disappear after the desorption-diffusion flow period, while the impact of adsorption index will continue to exist in the subsequent flow periods. Fracture spacing (FS). Figure 10 presents the typical pseudo-pressure curves of different fracture spacing with low, moderate and high fracture conductivities. As the fracture spacing increases, the first pseudo-radial flow stage becomes longer, indicating weaker pressure interference between the fractures. As can be seen in Figure 10, the first pseudo-radial flow can be identified clearly when FS is equal to 50, while it may be masked when FS is less than 20. In practice, the dimensionless fracture spacing is usually less than 5, thus, the first pseudo-radial flow may not occur. Another phenomenon shown in Figure 10 is that the dip moves upward as FS reduces. Asymmetry factor (AF). Figure 11 presents the impact of the asymmetry factor on the pseudopressure responses. We assume each fracture has an identical asymmetry factor, AF = 0, 0.4, and 0.8.
Observe that the asymmetry factor mainly influences the early flow periods. The fractures for different conductivities are parallel and coplanar without the wings interference, and the CBM production is mainly supplied by reservoirs far away from the well in the late flow period, the result being that the production performances of asymmetric fractures and symmetric fractures are approximate. For example, when CfD = 50 and tD = 10 3 , the pseudo-pressures for the AF of 0, 0.4 and 0.8 are 2.0077, 2.0098, and 2.0166, respectively.
Under the condition of low fracture conductivity (Figure 11(a)), the bigger the AF is, the earlier the end of bilinear flow occurs. The pseudo-pressure derivative curve also exhibits a deviation from the one fourth slope line responding to the bilinear flow for moderate fracture conductivity ( Figure  11(b)). When the dimensionless fracture conductivity is equal to 500, the first linear flow can also be identified with the impact of the asymmetry factor (Figure 11(c)).  Figure 12 shows the effect of non-planar angle on the pseudo-pressure responses. In our study, the hydraulic fracture system of a horizontal well in coalbed methane reservoirs can be simulated by random combinations of fracture wing, which means that each group of hydraulic fractures can be composed of odd or even fracture wings with different angles or lengths.

Non-planar angle (NPF).
For convenience, the non-planar angles between the two wings of each fracture are considered to be equal in this case. The values of NPF are set to be 30°, 90°, and 180°, respectively. Therefore, the case of NPF = 180° corresponds to the coplanar fracture in our study.
As shown in Figure 12, as the NPF decreases, the pressure interference between the two wings of each fracture reinforces, which eventually leads to a larger pressure drop. It should be noted that the pressure interference caused by the non-planar wings will continue to show a non-negligible effect on the late production response. Such as when CfD = 50 and tD = 10 3 , the pseudo-pressures for the NPF of 30°, 90°, and 180° are 2.1443, 2.0515, and 2.0077, respectively. In the case of low fracture conductivity (Figure 12(a)), a hump caused by wings interference can be observed with a small NPF, such as NPF = 30°, which eventually shortens the duration of the bilinear flow. With the increase in fracture conductivity, the influence of non-planar angle on the pseudo-pressure curves is delayed. Wellbore storage (CD). Figure 13 shows the typical pseudo-pressure curves of a multi-fractured horizontal well with the consideration of wellbore storage CD. The CD is considered to be 10 -7 , 10 -3 , and 10, respectively.
The characteristics of early flow periods depend on the value of CD. With the increase of CD, the duration of the wellbore storage period increases, which is characterized by a straight line with a slope of one. For small value of CD, for example, CD = 10 -7 (Figure 13(a)), the wellbore storage makes the occurrence of flow periods delay and all characteristic flow periods can be identified. However, for big values of CD, like CD = 10 -3 (Figure 13(b)), the wellbore storage can lead to the absence of the early bilinear flow period. For large values of CD, (CD = 10, Figure 13(c)), the wellbore storage will dominate the fluid flow in the whole of the early period, resulting the absence of characteristic flow periods.

Conclusions
Based on our study, several important conclusions can be drawn.
(1) Based on a new fracture wing model, a semi-analytical model has been proposed to obtain the pseudo-pressure responses of horizontal wells with multiple fractures in coalbed methane reservoirs. The effects of fracture conductivity, diffusion model, storability ratio, inter-porosity flow coefficient, adsorption index, fracture spacing, fracture asymmetry, non-planar angle, and wellbore storage on the transient pseudo-pressure responses have been investigated.
(2) Six characteristic flow periods can be identified: bilinear flow, first linear flow, desorptiondiffusion flow, first pseudo-radial flow, second linear flow, and second pseudo-radial flow for a multi-fractured horizontal well in coalbed methane reservoirs.
(3) The fracture conductivity exhibits a major effect on the early pseudo-pressure responses. In the case of low fracture conductivity, the bilinear flow can be recognized clearly, and the first linear flow can't be observed. With an increase in the fracture conductivity, the duration of the bilinear flow decreases and the first linear flow gradually dominates the fluid flow in the coalbed methane reservoir. At moderate fracture conductivity, both the bilinear flow period and the first linear flow period can be identified. When the dimensionless fracture conductivity is greater than 100, the first linear flow can be recognized easily, whereas the bilinear flow will not occur.
(4) Results show that several parameters affect the pseudo-pressure responses in coalbed methane reservoirs. As the parameter value increases, both the storability ratio and the inter-porosity flow coefficient make the lowest point in the dip move in the top right orientation; the adsorption index leads to the move of the lowest point towards the bottom left corner and the fracture spacing makes the lowest point drop vertically. Sensitivity analysis also indicates that the asymmetry factor mainly impacts the bilinear flow periods. A small non-planar angle will cause strong pressure interference between the wings, eventually making the end of the bilinear flow or the first linear flow occur earlier. Small wellbore storage merely delays the occurrence of flow periods, while large wellbore storage may conceal the bilinear flow and the first linear flow period.  For the fracture wing shown in Figure 2, the fluid flow in the wing is one-dimension in the Cartesian coordinate. The length of the wing is Lf with height h and width wf, and the fracture permeability is kf.
The dimensionless fracture pseudo-pressure is given as ( )