Transient Pressure Analysis of a Multiple Fractured Well in a Stress-Sensitive Coal Seam Gas Reservoir

This paper investigates the bottom-hole pressure (BHP) performance of a fractured well with multiple radial fracture wings in a coalbed methane (CBM) reservoir with consideration of stress sensitivity. The fluid flow in the matrix simultaneously considers adsorption–desorption and diffusion, whereas fluid flow in the natural fracture system and the induced fracture network obeys Darcy’s law. The continuous line-source function in the CBM reservoir associated with the discretization method is employed in the Laplace domain. With the aid of Stehfest numerical inversion technology and Gauss elimination, the transient BHP responses are determined and analyzed. It is found that the main flow regimes for the proposed model in the CBM reservoir are as follows: linear flow between adjacent radial fracture wings, pseudo-radial flow in the inner region or Stimulated Reservoir Volume (SRV), and radial flow in outer region (un-stimulated region). The effects of permeability modulus, radius of SRV, ratio of permeability in SRV to that in un-stimulated region, properties of radial fracture wings, storativity ratio of the un-stimulated region, inter-porosity flow parameter, and adsorption–desorption constant on the transient BHP responses are discussed. The results obtained in this study will be of great significance for the quantitative analyzing of the transient performances of the wells with multiple radial fractures in CBM reservoirs.

fractured wells with multiple radial artificial fractures in stress-sensitive CBM reservoirs. Recently, Wei et al. [10,11] and Yuan et al. [12] numerically investigated the flow mechanism for fractured horizontal wells in shale gas reservoirs, and there are many similar features between CBM reservoirs and shale gas reservoirs, such as the gas adsorption-desorption, gas diffusion, and stress sensitivity of the reservoir permeability, etc.
Numerous mathematical models have been adopted to investigate the fluid flow in CBM reservoirs with varied assumptions. King et al. [13] numerically simulated the gas-water flow in micropores of CBM reservoirs. Later, Anbarci and Ertekin [14,15] provided a novel well model in a CBM reservoir. In their model, two different flow regimes, steady state and pseudo-steady state, were considered under the effect of varied inner and outer boundary conditions. Then, Engler and Rajtar [16] investigated the BHP responses for the horizontal well in a CBM reservoir utilizing Fourier and Laplace transformation technologies. Subsequently, Clarkson et al. [17] analyzed the fractured well performances in terms of BHP and production rate in a CBM reservoir. In their model, the hydraulic fractures possess finite-conductivity. Recently, Nie et al. [18] established a semi-analytical model for a horizontal well in a coal seam, in which adsorption-desorption, diffusion, and Darcy flow were taken into account. More recently, Zhao et al. [19] obtained an analytical solution for the transient BHP response of a fractured well in a CBM reservoir. In their model, the induced fracture network in SRV were treated as an inner region with more desirable petro-physical properties.
As stimulation technologies develop, high-energy gas fracturing technology, a method employed to develop unconventional hydrocarbon reservoirs, is able to increase the number of artificial fractures and further enlarge the range of linear flow. More specifically, multi-wing radial artificial fractures can be obtained along the wellbore after the stimulation [20][21][22], which have been verified based on core analysis and microseimic image results [9,23,24] (Figure 1). Some analytical and numerical models have been proposed to analyze the transient BHP and rate performances with consideration of multi-wing artificial fractures. For example, Choo and Wu [25] derived a new numerical solution for multiple fractured vertical wells to investigate the BHP response. Later, Tiab [26] analyzed the transient BHP response of the model, taking asymmetrically-distributed hydraulic fractures with finite conductivity into account by employing the Tiab's Direct Synthesis Technique. Recently, Zhang et al. [20] proposed a well-testing model considering the fractured well with multiple radial hydraulic fractures in a composite CBM reservoir in order to simulate the transient pressure and rate by means of continuous line-source function. In addition, refracturing technology is also able to generate multiple radial artificial fractures. For example, Hou et al. [27] sketched an analytical solution to calculate the angle between adjacent hydraulic fractures initiated after refracturing. To verify the model, the authors compared the data obtained from oil wells in practical fields with the newly developed solution. However, most models proposed above are not able to consider the effect of SRV, which is a key element in enhancing gas recovery.
Usually, stimulation treatment serves two main purposes, generating fracture network and artificial fractures near the wellbore. The generated fracture network close to the wellbore in unconventional hydrocarbon reservoirs caused by stimulation treatment are generally termed as SRV. SRV in this Energies 2020, 13, 3849 3 of 20 study is considered as an inner region near the wellbore, which includes the multi-wing artificial fractures and the induced fracture network. To differentiate the SRV and un-stimulated region in CBM reservoirs, the CBM reservoir in this study is divided into two regions, each can be described by the dual-porosity system. Specifically, the inner region (SRV) possesses more desirable petro-physical properties, such as larger porosity and permeability.
Unlike conventional natural gas, unconventional natural gas such as CBM and shale gas are generally absorbed on the mineral particle surface in coal seam matrix and the CBM flow in reservoirs is subject to multiple transport mechanisms [28][29][30]. More specifically, as reservoir pressure decreases, the absorbed CBM molecules are able to desorb from the surface of mineral particles. Then, the gas molecules can be driven towards fracture under the effect of concentration difference (diffusion). Finally, due to pressure difference, the flow of CBM from coal seam matrix can be observed in fractures, which can be characterized by Darcy law, see Figure 2.
Energies 2020, 13, x FOR PEER REVIEW 3 of 21 Usually, stimulation treatment serves two main purposes, generating fracture network and artificial fractures near the wellbore. The generated fracture network close to the wellbore in unconventional hydrocarbon reservoirs caused by stimulation treatment are generally termed as SRV. SRV in this study is considered as an inner region near the wellbore, which includes the multi-wing artificial fractures and the induced fracture network. To differentiate the SRV and un-stimulated region in CBM reservoirs, the CBM reservoir in this study is divided into two regions, each can be described by the dual-porosity system. Specifically, the inner region (SRV) possesses more desirable petro-physical properties, such as larger porosity and permeability.
Unlike conventional natural gas, unconventional natural gas such as CBM and shale gas are generally absorbed on the mineral particle surface in coal seam matrix and the CBM flow in reservoirs is subject to multiple transport mechanisms [28][29][30]. More specifically, as reservoir pressure decreases, the absorbed CBM molecules are able to desorb from the surface of mineral particles. Then, the gas molecules can be driven towards fracture under the effect of concentration difference (diffusion). Finally, due to pressure difference, the flow of CBM from coal seam matrix can be observed in fractures, which can be characterized by Darcy law, see Figure 2. As is well known, transient pressure analysis (TPA) is suitable to determine key reservoir parameters and to monitor the transient BHP performance of gas wells. The research interest of this study is to propose a well-testing model, considering multiple transport mechanisms of CBM in the reservoir, to investigate the transient BHP response of a fractured well with multiple artificial fractures in a stress-sensitive CBM reservoir. The main characteristics of type curves obtained in this work are discussed. The model proposed in the study can be useful in well testing interpretations and production transient analyses of unconventional gas reservoirs.

Model Description
After stimulation treatment, both multiple artificial fractures and micro-fracture network can be initialized near the vertical wellbore, thus leading to the more favorable properties (i.e., higher permeability) in this stimulated region than that in the un-stimulated region. Therefore, the whole reservoir can be simplified as a composite reservoir system. As shown in Figure 3, a multiple fractured vertical well (MFVW) in a composite CBM reservoir is proposed in this study. As is well known, transient pressure analysis (TPA) is suitable to determine key reservoir parameters and to monitor the transient BHP performance of gas wells. The research interest of this study is to propose a well-testing model, considering multiple transport mechanisms of CBM in the reservoir, to investigate the transient BHP response of a fractured well with multiple artificial fractures in a stress-sensitive CBM reservoir. The main characteristics of type curves obtained in this work are discussed. The model proposed in the study can be useful in well testing interpretations and production transient analyses of unconventional gas reservoirs.

Model Description
After stimulation treatment, both multiple artificial fractures and micro-fracture network can be initialized near the vertical wellbore, thus leading to the more favorable properties (i.e., higher permeability) in this stimulated region than that in the un-stimulated region. Therefore, the whole reservoir can be simplified as a composite reservoir system. As shown in Figure 3, a multiple fractured vertical well (MFVW) in a composite CBM reservoir is proposed in this study.
Additionally, for the sake of deriving a semi-analytical solution for this model, other assumptions are as follows: (1) The CBM reservoir can be radially divided into two regions, region 1 and region 2. Region 1, the inner region, includes the micro-fracture network. Whereas region 2 represents the un-stimulated region. (2) Due to the influences of stimulation treatment, a multiple radial fractures model is adopted to veritably describe the dynamic flow process of CBM. (3) Infinite-conductivity artificial fractures are considered. (4) The constant production rate of MFVW is defined as q sc , however, the production rates at different locations of a certain fracture are unique.
Energies 2020, 13, 3849 4 of 20 (5) The permeability of the micro-fracture network in SRV is considered stress-sensitive. (6) The pseudo-steady state gas diffusion and unsteady state gas diffusion from matrix to fracture network both are considered in this model.  Additionally, for the sake of deriving a semi-analytical solution for this model, other assumptions are as follows: 1) The CBM reservoir can be radially divided into two regions, region 1 and region 2. Region 1, the inner region, includes the micro-fracture network. Whereas region 2 represents the un-stimulated region. 2) Due to the influences of stimulation treatment, a multiple radial fractures model is adopted to veritably describe the dynamic flow process of CBM. 3) Infinite-conductivity artificial fractures are considered. 4) The constant production rate of MFVW is defined as qsc, however, the production rates at different locations of a certain fracture are unique. 5) The permeability of the micro-fracture network in SRV is considered stress-sensitive. 6) The pseudo-steady state gas diffusion and unsteady state gas diffusion from matrix to fracture network both are considered in this model.

7)
The outer boundary is infinite and the bottom boundary and top boundary are both impermeable. 8) Ignoring capillary pressure and gravity.

Continuous Line-Source Solution
In this study, the number of the radial fracture wings is set to M. Each fracture wing is composed of N segments. For each radial fracture wing, the corresponding serial number and angle with respect to X-axis positive direction are defined in a clockwise direction for ease of later programming, as can be seen in Figure 4. As the purpose of this portion is to obtain a continuous line-source solution for a MFVW in a CBM reservoir, the seepage differential equations for two regions should be developed respectively. Then, the continuous line-source function in this composite CBM reservoir can be derived by coupling the governing equation in SRV with that in the un-stimulated region. Zhang

Continuous Line-Source Solution
In this study, the number of the radial fracture wings is set to M. Each fracture wing is composed of N segments. For each radial fracture wing, the corresponding serial number and angle with respect to X-axis positive direction are defined in a clockwise direction for ease of later programming, as can be seen in Figure 4. As the purpose of this portion is to obtain a continuous line-source solution for a MFVW in a CBM reservoir, the seepage differential equations for two regions should be developed respectively. Then, the continuous line-source function in this composite CBM reservoir can be derived by coupling the governing equation in SRV with that in the un-stimulated region. Zhang et al. [20] proposed a novel line-source function for a CBM reservoir considering Darcy flow and Knudson diffusion. In this study, an improved continuous line-source function for a MFVW in a composite CBM reservoir with consideration of stress sensitivity is presented based on Zhang et al.'s work, see Equation (1). The relevant dimensionless parameter definitions can be found in Appendix A and the derivations of continuous line-source functions under the unsteady state diffusion and pseudo-steady state diffusion in a composite CBM reservoir are shown in Appendices B and C, respectively.

Mathematical Model of MFVW
Since it is difficult to obtain the transient pressure response by direct integration of Equation (1), the discretization and superposition methods can be applied in this study. For example, each fracture is composed of N radial small segments and the flux density inside each segment can be considered uniform. Therefore, the transient pseudo pressure response at an arbitrary position in the composite CBM reservoir, generated due to the production of the ith segment of the jth radial fracture, can be obtained by the integration of Equation (1) along the small segment: the jth radial fracture, respectively, which can be expressed as follows: where j  is the angle of the jth radial fracture, Dj r  is the dimensionless radial length of each segment in the jth radial fracture.

Mathematical Model of MFVW
Since it is difficult to obtain the transient pressure response by direct integration of Equation (1), the discretization and superposition methods can be applied in this study. For example, each fracture is composed of N radial small segments and the flux density inside each segment can be considered uniform. Therefore, the transient pseudo pressure response at an arbitrary position in the composite CBM reservoir, generated due to the production of the ith segment of the jth radial fracture, can be obtained by the integration of Equation (1) along the small segment: whereq fD(i,j) is the dimensionless flux density of the ith segment of the jth radial fracture and x D(i,j) and y D(i,j) are the horizontal axis and vertical axis of the middle point of the ith segment of the jth radial fracture, respectively, which can be expressed as follows: where θ j is the angle of the jth radial fracture, ∆r D j is the dimensionless radial length of each segment in the jth radial fracture. Substituting Equations (5) and (6) into Equation (2) results in: where v = r/L ref .
Energies 2020, 13, 3849 6 of 20 Equation (7) can be rewritten as follows based on the variable substitution method: Because it is time-consuming to directly program Equation (9), variable substitution can be applied again (integration variable v can be displaced by integration variable σ), thus Equation (9) can be further transferred into the following equation: where r D(i,j) is the dimensionless radial distance from wellbore to the middle point of the ith segment in the jth radial fracture. Therefore, the transient pseudo pressure response at the kth segment of the mth radial fracture caused by the ith segment of the jth radial fracture can be provided as: Since the transient pressure response at the kth segment of the mth radial fracture is caused by the production of all the segments (including the kth segment of the mth radial fracture itself), it can be determined by using the following superposition principle: where ξ 1fD0(ij,km) is the transient pressure response at the kth segment of the mth radial fracture caused by the ith segment of the jth radial fracture. In addition, due to the infinite-conductivity radial fractures, we can obtain: According to Equation (19), for M × N segments, we can obtain M × N algebraic equations with (M × N + 1) unknowns, ξ wD ,q fD(1,1) ,q fD(1,2) ,q fD(N,M) . As a result, one more linear algebraic equation is required to obtain the final solution.
where Equation (20) represents the total production rate of the well being equal to the summation of the production rates of all M × N segments in the Laplace domain. Additionally, the Duhamel's principle can be employed to incorporate both wellbore storage and skin effects. In this study, the expression presented by Kucuk and Ayestaran [31] and Gringarten et al. [32] is applied to consider the influences of the wellbore storage and skin factor: where ξ wDv is the dimensionless pseudo BHP with consideration of wellbore storage and skin factor in the Laplace domain.
In addition, the BHP response of the well in the time domain can be determined by applying the Stehfest numerical inversion method [33][34][35][36].
Finally, the term ξ wDv (t D ) in Equation (22) is the zero-order perturbation solution for the dimensionless pseudo BHP. Therefore, the real dimensionless pseudo BHP in this study, which takes stress sensitivity in the micro-fracture network into account, can be determined by employing Equation (23) [4,5]:

Model Validation
In order to validate the developed model, two comparisons in terms of pseudo pressure response (PPR) and pseudo pressure derivative (PPD), one between the developed model and commercial software and the other one between the developed model and the model proposed by relevant researchers, are conducted. Because it is hard to consider multi-wing fractured wells in software, the'conventional fractured well with two-wing fractures is adopted, where the angles of right and left radial fracture are set to 0 • and 180 • , respectively. The effect of SRV is ignored, suggesting identical petro-physical properties of the inner region and outer region. In addition, because the model in Saphir considers that the desorbed CBM molecules directly enter the natural fracture network, we utilized f (u) to replace f 1 (u) and f 2 (u) in the Laplace domain. Figure 5a shows the satisfactory comparison results, indicating the developed model is reliable.
The second comparison between the developed model and the model proposed by Zhang et al. [20] focuses on the effect of multi-wing fractures on the PPR and PPD. To be more specific, desorption and diffusion of CBM molecules are ignored and the petro-physical properties of the inner and outer region are identical. Based on the results from Figure 5b, the developed model is in good agreement with the proposed model, again, denoting the reliability of the developed model.

Type Curves for the Proposed Model
According the description in Section 3, since there are (M×N+1) algebraic equations with (M×N+1) unknowns, the combination of Equation (18) to Equation (23) and the utilization of computer programming could result in the time-dependent PPR and PPD, see Figure 6. Based on Figure 6, the type curves resulting from the proposed model can be divided into the following eight flow regimes: Stage 1: Wellbore storage. The PPR and PPD both are straight lines with unit slope during this flow regime. Stage 4: Pseudo-radial flow in SRV. As the pressure wave propagates, the pseudo-radial flow in SRV can be observed during this regime and the curve of the PPR exhibits a horizontal line. It is worth noting that this flow regime is subject to the radial fracture length-SRV radius ratio.
Stage 5: Short-time transition flow regime from pseudo-radial flow in the inner region to radial fluid flow in the outer region.
Stage 6: Radial flow regime in the natural fracture system of the outer region. As the pressure wave propagates farther, the SRV-centered radial flow occurs in the outer region. This flow regime is characterized as a flat trend in the PPD curve during this period.
Stage 7: Diffusive flow regime (matrix-dominated flow regime). The gas concentration difference between the natural fracture system and matrix can be expected with the production of CBM residing in the natural fracture system in the outer region, CBM molecules start to desorb from the surfaces of mineral particles in the matrix and diffuse into the natural fracture system under the effect of concentration difference between the natural fracture system and coal seam matrix. For the pseudo-steady state diffusion model, an obvious "dip" appears in the PPD curve during this flow regime, however, for the unsteady state diffusion model, a less obvious "dip" can be found in the PPD curve during this flow regime.
Stage 8: Pseudo-radial flow in the unstimulated regime. A dynamic balance is achieved for the gas transfer between the natural fracture system and the coal seam matrix in the un-stimulated region and the PPD curve exhibits a horizontal line, whose vertical-axis value is 0.5.

Type Curves for the Proposed Model
According the description in Section 3, since there are (M × N + 1) algebraic equations with (M × N + 1) unknowns, the combination of Equation (18) to Equation (23) and the utilization of computer programming could result in the time-dependent PPR and PPD, see Figure 6. Based on Figure 6, the type curves resulting from the proposed model can be divided into the following eight flow regimes: Energies 2020, 13, x FOR PEER REVIEW 9 of 21 It is worth noting that the upward trend both in the PPR and PPD at late flow periods can be observed if the stress-sensitive effect of the fracture network in SRV is considered, representing that more pressure depletion is required for the production of CBM in stress-sensitive reservoirs, see Figure 6.
Additionally, because the un-steady state diffusion model is more practical in most cases, the un-steady state diffusion model is applied here to perform the following PPR/PPD-sensitive analyses. Figure 7 illustrates the three main flow regimes resulting from the proposed model.  Stage 4: Pseudo-radial flow in SRV. As the pressure wave propagates, the pseudo-radial flow in SRV can be observed during this regime and the curve of the PPR exhibits a horizontal line. It is worth noting that this flow regime is subject to the radial fracture length-SRV radius ratio.
Stage 5: Short-time transition flow regime from pseudo-radial flow in the inner region to radial fluid flow in the outer region.
Stage 6: Radial flow regime in the natural fracture system of the outer region. As the pressure wave propagates farther, the SRV-centered radial flow occurs in the outer region. This flow regime is characterized as a flat trend in the PPD curve during this period.
Energies 2020, 13, 3849 9 of 20 Stage 7: Diffusive flow regime (matrix-dominated flow regime). The gas concentration difference between the natural fracture system and matrix can be expected with the production of CBM residing in the natural fracture system in the outer region, CBM molecules start to desorb from the surfaces of mineral particles in the matrix and diffuse into the natural fracture system under the effect of concentration difference between the natural fracture system and coal seam matrix. For the pseudo-steady state diffusion model, an obvious "dip" appears in the PPD curve during this flow regime, however, for the unsteady state diffusion model, a less obvious "dip" can be found in the PPD curve during this flow regime.
Stage 8: Pseudo-radial flow in the unstimulated regime. A dynamic balance is achieved for the gas transfer between the natural fracture system and the coal seam matrix in the un-stimulated region and the PPD curve exhibits a horizontal line, whose vertical-axis value is 0.5.
It is worth noting that the upward trend both in the PPR and PPD at late flow periods can be observed if the stress-sensitive effect of the fracture network in SRV is considered, representing that more pressure depletion is required for the production of CBM in stress-sensitive reservoirs, see Figure 6.
Additionally, because the un-steady state diffusion model is more practical in most cases, the un-steady state diffusion model is applied here to perform the following PPR/PPD-sensitive analyses. Figure 7 illustrates the three main flow regimes resulting from the proposed model.  It is worth noting that the upward trend both in the PPR and PPD at late flow periods can be observed if the stress-sensitive effect of the fracture network in SRV is considered, representing that more pressure depletion is required for the production of CBM in stress-sensitive reservoirs, see Figure 6.
Additionally, because the un-steady state diffusion model is more practical in most cases, the un-steady state diffusion model is applied here to perform the following PPR/PPD-sensitive analyses. Figure 7 illustrates the three main flow regimes resulting from the proposed model.

Effect of the Ratio of Permeability in the Inner Region on that in the Outer Region
In the following subsections, sensitivity analyses are performed in terms of several key parameters for the development of CBM based on the developed model and programming. Usually, SRV can be generated after stimulation treatment for the cost-effective development of unconventional hydrocarbon. As a result, the determination of permeability of SRV is of great significance [37]. Figure 8 demonstrates the impacts of the ratio of the permeability in the SRV to that in the un-stimulated region, M12, on the PPR/PPD for a rate-constant production. The PPR/PPD decreases with increasing M12 in both linear and pseudo-radial flow regimes, denoting that the higher permeability of SRV is able to decrease energy consumption for the rate-constant production. Therefore, the initialization of SRV is of great importance for the effective development of CBM.

Effect of the Ratio of Permeability in the Inner Region on that in the Outer Region
In the following subsections, sensitivity analyses are performed in terms of several key parameters for the development of CBM based on the developed model and programming. Usually, SRV can be generated after stimulation treatment for the cost-effective development of unconventional hydrocarbon. As a result, the determination of permeability of SRV is of great significance [37]. Figure 8 demonstrates the impacts of the ratio of the permeability in the SRV to that in the un-stimulated region, M 12 , on he PPR/PPD for a rate-constant production. The PPR/PPD decreases with increasing M 12 in both linear and pseudo-radial flow regimes, denoting that the higher permeability of SRV is able to decrease energy consumption for the rate-constant production. Therefore, the initialization of SRV is of great importance for the effective development of CBM.
Energies 2020, 13, x FOR PEER REVIEW 10 of 21 Figure 8. Effect of the ratio of the permeability in the SRV to that in the un-stimulated region on the pseudo pressure response (PPR) and pseudo pressure derivative (PPD). Figure 9 illustrates the influence of the varied SRV radius (r1) on the PPR/PPD while the rate production keeps constant. The parameter r1D is the dimensionless radius of SRV in the composite CBM reservoirs, including the induced micro-fracture network and multiple fracture wings. As can be seen in Figure 9, the SRV radius dramatically affects the duration time of the pseudo-radial flow in SRV. More specifically, larger dimensionless radius of SRV corresponds to the longer duration time of the pseudo-radial flow regime in SRV, denoting that more pressure (or energy) is required for production.

Effect of Permeability Modulus
Since stress-sensitive permeability of the micro-fracture network in SRV is considered in this study, Figure 10 illustrates the impact of the permeability modulus (γ D ) on the PPR/PPD for the rate-constant production. Based on Figure 10, a smaller permeability modulus corresponds to the lower PPR/PPD (less obvious upward trend) in late flow regimes, which represents that the existence of stress sensitivity leads to larger pressure depletion in reservoirs compared with the no-stress-sensitivity case.  Figure 9 illustrates the influence of the varied SRV radius (r 1 ) on the PPR/PPD while the rate production keeps constant. The parameter r 1D is the dimensionless radius of SRV in the composite CBM reservoirs, including the induced micro-fracture network and multiple fracture wings. As can be seen in Figure 9, the SRV radius dramatically affects the duration time of the pseudo-radial flow in SRV. More specifically, larger dimensionless radius of SRV corresponds to the longer duration time of the pseudo-radial flow regime in SRV, denoting that more pressure (or energy) is required for production.

Effect of the Radius of SRV Region
Energies 2020, 13, x FOR PEER REVIEW 10 of 21 Figure 8. Effect of the ratio of the permeability in the SRV to that in the un-stimulated region on the pseudo pressure response (PPR) and pseudo pressure derivative (PPD). Figure 9 illustrates the influence of the varied SRV radius (r1) on the PPR/PPD while the rate production keeps constant. The parameter r1D is the dimensionless radius of SRV in the composite CBM reservoirs, including the induced micro-fracture network and multiple fracture wings. As can be seen in Figure 9, the SRV radius dramatically affects the duration time of the pseudo-radial flow in SRV. More specifically, larger dimensionless radius of SRV corresponds to the longer duration time of the pseudo-radial flow regime in SRV, denoting that more pressure (or energy) is required for production.

Effect of Permeability Modulus
Since stress-sensitive permeability of the micro-fracture network in SRV is considered in this study, Figure 10 illustrates the impact of the permeability modulus (γ D ) on the PPR/PPD for the rate-constant production. Based on Figure 10, a smaller permeability modulus corresponds to the lower PPR/PPD (less obvious upward trend) in late flow regimes, which represents that the existence of stress sensitivity leads to larger pressure depletion in reservoirs compared with the no-stress-sensitivity case.

Effect of Permeability Modulus
Since stress-sensitive permeability of the micro-fracture network in SRV is considered in this study, Figure 10 illustrates the impact of the permeability modulus (γ D ) on the PPR/PPD for the rate-constant production. Based on Figure 10, a smaller permeability modulus corresponds to the lower PPR/PPD (less obvious upward trend) in late flow regimes, which represents that the existence of stress sensitivity leads to larger pressure depletion in reservoirs compared with the no-stress-sensitivity case.  Figure 11 illustrates the impact of radial fracture wing angle symmetry on the PPR/PPD. As shown in Figure 11, the PPR/PPD increases with the increasing θ3 and θ6 and with the decreasing θ1 and θ4 in the later period of linear flow. This can be explained by the fact that the fracture interference becomes more severe as the angle between adjacent radial fractures decreases. That is to say, uniformly distributed radial fracture wings are able to weaken the fracture interference, and, as a result, reduce the energy consumption.  Figure 12 presents the effect of the number of radial fracture wings on the PPR/PPD. Based on Figure 12, the number of radial fracture wings impacts the PPR/PPD drastically. The PPR/PPD decreases with the increasing number of radial fracture wings from two to six. Therefore, the generation of multiple radial fracture wings is able to reduce the pressure depletion (energy consumption) required for production. The results obtained here are of great importance for the practical stimulation design.  Figure 11 illustrates the impact of radial fracture wing angle symmetry on the PPR/PPD. As shown in Figure 11, the PPR/PPD increases with the increasing θ 3 and θ 6 and with the decreasing θ 1 and θ 4 in the later period of linear flow. This can be explained by the fact that the fracture interference becomes more severe as the angle between adjacent radial fractures decreases. That is to say, uniformly distributed radial fracture wings are able to weaken the fracture interference, and, as a result, reduce the energy consumption.

Effect of Radial Fracture Angle Symmetry
Energies 2020, 13, x FOR PEER REVIEW 11 of 21 Figure 10. Effect of permeability modulus on the PPR and PPD. Figure 11 illustrates the impact of radial fracture wing angle symmetry on the PPR/PPD. As shown in Figure 11, the PPR/PPD increases with the increasing θ3 and θ6 and with the decreasing θ1 and θ4 in the later period of linear flow. This can be explained by the fact that the fracture interference becomes more severe as the angle between adjacent radial fractures decreases. That is to say, uniformly distributed radial fracture wings are able to weaken the fracture interference, and, as a result, reduce the energy consumption.  Figure 12 presents the effect of the number of radial fracture wings on the PPR/PPD. Based on Figure 12, the number of radial fracture wings impacts the PPR/PPD drastically. The PPR/PPD decreases with the increasing number of radial fracture wings from two to six. Therefore, the generation of multiple radial fracture wings is able to reduce the pressure depletion (energy consumption) required for production. The results obtained here are of great importance for the practical stimulation design.  Figure 12 presents the effect of the number of radial fracture wings on the PPR/PPD. Based on Figure 12, the number of radial fracture wings impacts the PPR/PPD drastically. The PPR/PPD decreases with the increasing number of radial fracture wings from two to six. Therefore, the generation of multiple radial fracture wings is able to reduce the pressure depletion (energy consumption) required for production. The results obtained here are of great importance for the practical stimulation design.  Figure 13 illustrates the impact of the length of radial fracture wings on the PPR/PPD for a rate-constant production. The linear flow regime and pseudo-radial flow regime in SRV can be affected by the length of radial fracture wings for a constant SRV radius, see Figure 13. The PPR/PPD decreases with the increasing length of radial fracture wings from 10m to 30m and the duration time of the pseudo-radial flow in SRV becomes shorter as the length of radial fracture wing increases. More specifically, the radial fracture wings with smaller length increase the pressure depletion during production.

Effect of the Storativity Ratio of the Outer Region
Since the double-porosity system is considered in this study, Figure 14 shows the effect of the storativity ratio of the outer region on the PPR/PPD. The storativity ratio of the outer region mainly affects the diffusive flow regime in the outer region, see Figure 14. A lower storativity ratio of the outer region corresponds to a wider and deeper concave during this flow regime.  Figure 13 illustrates the impact of the length of radial fracture wings on the PPR/PPD for a rate-constant production. The linear flow regime and pseudo-radial flow regime in SRV can be affected by the length of radial fracture wings for a constant SRV radius, see Figure 13. The PPR/PPD decreases with the increasing length of radial fracture wings from 10m to 30m and the duration time of the pseudo-radial flow in SRV becomes shorter as the length of radial fracture wing increases. More specifically, the radial fracture wings with smaller length increase the pressure depletion during production.  Figure 13 illustrates the impact of the length of radial fracture wings on the PPR/PPD for a rate-constant production. The linear flow regime and pseudo-radial flow regime in SRV can be affected by the length of radial fracture wings for a constant SRV radius, see Figure 13. The PPR/PPD decreases with the increasing length of radial fracture wings from 10m to 30m and the duration time of the pseudo-radial flow in SRV becomes shorter as the length of radial fracture wing increases. More specifically, the radial fracture wings with smaller length increase the pressure depletion during production.

Effect of the Storativity Ratio of the Outer Region
Since the double-porosity system is considered in this study, Figure 14 shows the effect of the storativity ratio of the outer region on the PPR/PPD. The storativity ratio of the outer region mainly affects the diffusive flow regime in the outer region, see Figure 14. A lower storativity ratio of the outer region corresponds to a wider and deeper concave during this flow regime.

Effect of the Storativity Ratio of the Outer Region
Since the double-porosity system is considered in this study, Figure 14 shows the effect of the storativity ratio of the outer region on the PPR/PPD. The storativity ratio of the outer region mainly affects the diffusive flow regime in the outer region, see Figure 14. A lower storativity ratio of the outer region corresponds to a wider and deeper concave during this flow regime.

Effect of the Inter-Porosity Flow Parameter
The other important parameter resulting from the double-porosity system is the inter-porosity flow parameter. Figure 15 demonstrates the effect of the inter-porosity flow parameter on the PPR/PPD for the rate-constant production. Based on the results provided by Figure 15, the inter-porosity flow parameter affects the diffusion flow regime drastically: the diffusion flow regime occurs later as the inter-porosity flow parameter decreases.  Figure 16 shows the effect of the adsorption-desorption constant on the PPR/PPD. The adsorption-desorption constant mainly affects the diffusion flow regime, as can be seen in Figure 16. More specifically, a higher adsorption-desorption constant corresponds to a deeper and wider concave, which represents the diffusion regime. The adsorption-desorption constant is adopted in this study to represent the amount of gas adsorbed at the surface of mineral surfaces and a higher adsorption-desorption constant denotes more adsorbed gas existing in the coal seam matrix. Therefore, more adsorbed gas is able to desorb and diffuse into the fractures for a larger adsorptiondesorption constant during production.

Effect of the Inter-Porosity Flow Parameter
The other important parameter resulting from the double-porosity system is the inter-porosity flow parameter. Figure 15 demonstrates the effect of the inter-porosity flow parameter on the PPR/PPD for the rate-constant production. Based on the results provided by Figure 15, the inter-porosity flow parameter affects the diffusion flow regime drastically: the diffusion flow regime occurs later as the inter-porosity flow parameter decreases.
Energies 2020, 13, x FOR PEER REVIEW 13 of 21 Figure 14. Effect of the storativity ratio of outer region on the PPR and PPD.

Effect of the Inter-Porosity Flow Parameter
The other important parameter resulting from the double-porosity system is the inter-porosity flow parameter. Figure 15 demonstrates the effect of the inter-porosity flow parameter on the PPR/PPD for the rate-constant production. Based on the results provided by Figure 15, the inter-porosity flow parameter affects the diffusion flow regime drastically: the diffusion flow regime occurs later as the inter-porosity flow parameter decreases.  Figure 16 shows the effect of the adsorption-desorption constant on the PPR/PPD. The adsorption-desorption constant mainly affects the diffusion flow regime, as can be seen in Figure 16. More specifically, a higher adsorption-desorption constant corresponds to a deeper and wider concave, which represents the diffusion regime. The adsorption-desorption constant is adopted in this study to represent the amount of gas adsorbed at the surface of mineral surfaces and a higher adsorption-desorption constant denotes more adsorbed gas existing in the coal seam matrix. Therefore, more adsorbed gas is able to desorb and diffuse into the fractures for a larger adsorptiondesorption constant during production.  Figure 16 shows the effect of the adsorption-desorption constant on the PPR/PPD. The adsorption-desorption constant mainly affects the diffusion flow regime, as can be seen in Figure 16. More specifically, a higher adsorption-desorption constant corresponds to a deeper and wider concave, which represents the diffusion regime. The adsorption-desorption constant is adopted in this study to represent the amount of gas adsorbed at the surface of mineral surfaces and a higher adsorption-desorption constant denotes more adsorbed gas existing in the coal seam matrix. Therefore, more adsorbed gas is able to desorb and diffuse into the fractures for a larger adsorption-desorption constant during production.

Conclusions
This work presents a semi-analytical model for a vertical well with multiple radial fracture wings in a stress-sensitive CBM reservoir. The corresponding PPR and PPD in the well bottom are determined and discussed. The conclusions obtained in this work are shown as follows: 1) The linear flow between adjacent radial fracture wings, the radial flow in the SRV, and the radial flow in the outer region are the three main flow regimes for the proposed model in CBM reservoirs.
2) The impact of stress sensitivity on the PPR/PPD is obvious. The existence of stress sensitivity in the micro-fracture network results in larger pressure depletion in later flow regimes.
3) The SRV, which includes the micro-fracture network and radial fracture wings, is able to reduce the pressure depletion. As the size of the SRV becomes smaller, the transition flow regime from pseudo-radial flow in the SRV to radial flow in the outer region occurs earlier. 4) The properties and distribution of multiple radial fracture wings affect the PPR/PPD in CBM reservoirs drastically. An increase in the number of radial fracture wings leads to a decrease in the pressure depletion when producing. The well model with an un-uniform fracture wing distribution requires more energy consumption compared to that with a uniform fracture wing distribution when producing at the same rate. 5) The storativity ratio and the inter-porosity flow parameter of the outer region, two parameters in a double-porosity system, mainly affect the diffusion flow regime, where the storativity ratio of the outer region represents the capacity of gas supply for the matrix and the inter-porosity flow parameter denotes the diffusion occurrence time. In addition, the adsorption-desorption constant characterizes the amount of adsorbed gas in the matrix.

Acknowledgments:
The authors would also like to thank the reviewers and editors whose critical comments were very helpful in preparing this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
This work presents a semi-analytical model for a vertical well with multiple radial fracture wings in a stress-sensitive CBM reservoir. The corresponding PPR and PPD in the well bottom are determined and discussed. The conclusions obtained in this work are shown as follows: (1) The linear flow between adjacent radial fracture wings, the radial flow in the SRV, and the radial flow in the outer region are the three main flow regimes for the proposed model in CBM reservoirs. (2) The impact of stress sensitivity on the PPR/PPD is obvious. The existence of stress sensitivity in the micro-fracture network results in larger pressure depletion in later flow regimes. (3) The SRV, which includes the micro-fracture network and radial fracture wings, is able to reduce the pressure depletion. As the size of the SRV becomes smaller, the transition flow regime from pseudo-radial flow in the SRV to radial flow in the outer region occurs earlier. (4) The properties and distribution of multiple radial fracture wings affect the PPR/PPD in CBM reservoirs drastically. An increase in the number of radial fracture wings leads to a decrease in the pressure depletion when producing. The well model with an un-uniform fracture wing distribution requires more energy consumption compared to that with a uniform fracture wing distribution when producing at the same rate. (5) The storativity ratio and the inter-porosity flow parameter of the outer region, two parameters in a double-porosity system, mainly affect the diffusion flow regime, where the storativity ratio of the outer region represents the capacity of gas supply for the matrix and the inter-porosity flow parameter denotes the diffusion occurrence time. In addition, the adsorption-desorption constant characterizes the amount of adsorbed gas in the matrix.
where L ref is reference length, which can be replaced by the length of radial fracture, r f , in this study. Some dimensionless parameters regarding the dual-porosity system are defined as follows: The dimensionless gas concentration difference in matrix can be defined as: The ratio of the permeability in SRV to that in un-stimulated region is given as follows: The dimensionless wellbore storage coefficient is defined as: In addition, the following dimensionless variables can be obtained: