Investigation of the Effect of Natural Fractures on Multiple Shale-Gas Well Performance Using Non-Intrusive EDFM Technology

The influence of complex natural fractures on multiple shale-gas well performance with varying well spacing is poorly understood. It is difficult to apply the traditional local grid refinement with structured or unstructured gridding techniques to accurately and efficiently handle complex natural fractures. In this study, we introduced a powerful non-intrusive embedded discrete fracture model (EDFM) technology to overcome the limitations of exiting methods. Through this unique technology, complex fracture configurations can be easily and explicitly embedded into structured matrix blocks. We set up a field-scale two-phase reservoir model to history match field production data and predict long-term recovery from Marcellus. The effective fracture properties were determined thorough history matching. In addition, we extended the single-well model to include two horizontal wells with and without including natural fractures. The effects of different numbers of natural fractures on two-well performance with varying well spacing of 200 m, 300 m, and 400 m were examined. The simulation results illustrate that gas productivity almost linearly increases with the number of two-set natural fractures. Furthermore, the difference of well performance between different well spacing increases with an increase in natural fracture density. A larger well spacing is preferred for economically developing the shale-gas reservoirs with a larger natural fracture density. The findings of this study provide key insights into understanding the effect of natural fractures on well performance and well spacing optimization.


Introduction
Typical shale gas reservoirs consist of free gas and adsorbed gas.Economic shale gas production has been enabled by the advanced technologies of multiple horizontal wells and multi-stage hydraulic fracturing.The U.S. Energy Information Administration (EIA) estimates that dry shale gas production (0.48 trillion cubic meters) accounts for around 62% of the total U.S. dry natural gas production in 2017 [1].The Marcellus shale formation is the most productive gas field in the United States.It has been reported that natural fractures are commonly observed in many shale gas formations based on outcrops, cores and image logs [2].Engelder et al. [3] observed that there are two sets of natural fractures or two regional joint sets (J 1 and J 2 sets) in the Marcellus shale formation, which significantly affect the successful hydraulic fracturing treatment.A better understanding of natural fracture effects on well performance plays an important role in well spacing optimization of shale gas reservoirs.
The reactivation of pre-existing natural fractures during the hydraulic fracturing process is an important factor to create complex fracture networks, which have been observed and predicted by many examples of microseismic event patterns [4,5] and complex fracture propagation models [6][7][8][9][10][11]. Complex fracture networks will create a large productive fracture surface area, which is necessary to maximize shale-gas production.Based on the history matching results with an actual shale-gas well, Yu et al. [12] compared the long-term well productivity between simple and complex fracture geometries and found that complex fracture networks can produce 36.4% more gas recovery after 30 years than the simple fractures.Yu et al. [13] built a synthetic single shale-gas well model including 11 planar hydraulic fractures and 200 natural fractures.The authors found that the well performance of 200 two-set natural fractures is much better than that of 200 one-set natural fractures due to the formation of a much more complex connected fracture network.An increase of gas recovery of about 23.2% was achieved by the 200 two-set natural fractures when compared to the base model without natural fractures.Although there are many existing reservoir simulation studies for well spacing optimization in shale gas reservoirs [14][15][16][17][18][19], the influence of natural fractures on multiple shale-gas well performance with varying well spacing has not been well examined and understood.
Although there are many analytical or semi-analytical models to simulate shale-gas reservoirs [20][21][22][23], numerical reservoir simulation is needed to accurately model multiple shale-gas well production due to complexity of natural and hydraulic fracture configurations and complex two-phase flow physics [24][25][26][27].The traditional local grid refinement (LGR) method with structured grids is difficult to model complex fractures explicitly [28,29].Although the LGR method with unstructured grids has the capability to handle complex fractures, the computational efficiency is a big issue.In addition, advanced parallel computing power is generally needed.Furthermore, when performing sensitivity studies and history matching with varying fracture geometries such as fracture number, length, and height, re-gridding of matrix blocks containing fractures is required.In our previous work, we developed an innovative non-intrusive embedded discrete fracture model (EDFM) technology in conjunction with any third-party reservoir simulators with structured grids to accurately and efficiently handle any complex hydraulic and natural fractures [30,31], which provides a unique solution to overcome the above limitations of existing methods.The non-intrusive feature means that we do not need to get access to the source code of commercial reservoir simulators and just modifying their input files.
In this study, we introduced the non-intrusive EDFM technology in conjunction with a commercial reservoir simulator of CMG-GEM [32] to simulate shale gas production considering complex natural fractures and two-phase flow (gas and water).Based on an actual shale-gas well with available gas and water flow rates from the Marcellus shale formation, we build a field-scale reservoir model to perform history matching with gas and water flow rates under the flowing bottomhole pressure (BHP) constraint.After history matching, we predict a 30-year production forecasting.Subsequently, we extend the history-matched reservoir model with effective fracture properties to include two horizontal wells with three different well spacing values of 200, 300 and 400 m.In addition, the impacts of different numbers of natural fractures such as 100, 1000, 5000, and 100,000 on multiple shale-gas well performance are investigated.The effect of natural fracture density on well spacing optimization and pressure distribution is discussed.

Methodology
The non-intrusive EDFM technology in conjunction with the commercial reservoir simulator was originally developed by Xu et al. [30,31].Here, we only briefly introduced this powerful technology.Based on the non-intrusive EDFM technology, both 3D complex hydraulic and natural fractures in the physical domain can be directly and explicitly embedded into the simple structured matrix grids, as demonstrated in Figure 1a.Based on the intersections between complex fractures and matrix grids, a number of extra fracture grids will be generated, as shown in Figure 1b [30].In Figure 1, two 3D fractures are explicitly embedded into three matrix grids in the physical domain.Correspondingly, an additional four fracture grids are created in the computational domain.Multi-phase flow such as gas and water between matrix and fracture grids can be conveniently simulated through calculating transmissibility between these non-neighboring connections (NNCs) grids, as shown in different colors of arrows of Figure 1b, which is a general feature for commercial reservoir simulators to handle faults before.A non-intrusive EDFM preprocessor has been developed to automatically check the complex intersections of fractures and matrix grids and calculate the transmissibility and other physical properties of the additional fracture grids such as porosity and depth, which are needed to input for the commercial reservoir simulator.two 3D fractures are explicitly embedded into three matrix grids in the physical domain.
Correspondingly, an additional four fracture grids are created in the computational domain.Multiphase flow such as gas and water between matrix and fracture grids can be conveniently simulated through calculating transmissibility between these non-neighboring connections (NNCs) grids, as shown in different colors of arrows of Figure 1b, which is a general feature for commercial reservoir simulators to handle faults before.A non-intrusive EDFM preprocessor has been developed to automatically check the complex intersections of fractures and matrix grids and calculate the transmissibility and other physical properties of the additional fracture grids such as porosity and depth, which are needed to input for the commercial reservoir simulator.There are three types of NNCs including NNC type I, which is the connection between matrix and fracture, NNC type II, which is the connection between fracture segments in a single fracture, and NNC type III, which is the connection between intersecting fracture segments [30].The flow rate ( l q , m 3 /s) of phase l between NNC grids is calculated by the following equation: where λl represents the relative mobility of phase l (cp −1 ), Δpl represents the pressure difference between NNC grids (Pa), TNNC represents the transmissibility factor of NNC grids (mD-m), which can be calculated by: There are three types of NNCs including NNC type I, which is the connection between matrix and fracture, NNC type II, which is the connection between fracture segments in a single fracture, and NNC type III, which is the connection between intersecting fracture segments [30].The flow rate (q l , m 3 /s) of phase l between NNC grids is calculated by the following equation: where λ l represents the relative mobility of phase l (cp −1 ), ∆p l represents the pressure difference between NNC grids (Pa), T NNC represents the transmissibility factor of NNC grids (mD-m), which can be calculated by: where k NNC represents the matrix grid permeability for NNC type I and average fracture permeability for NNC types II and III (mD); A NNC represents the contact area between NNC grids (m 2 ); d NNC represents the connection distance between different NNC grids (m).It should be mentioned that the non-intrusive EDFM technology only deals with the transmissibility factor and does not involve relative phase mobility calculation.Hence, it can be applied in both black oil simulation and compositional simulation.
For the connection between fractures and wellbore, the following effective wellbore index (W I, mD-m) will be calculated [30]: where k f is the fracture permeability (mD), w f is fracture width (m), L is fracture segment length (m), W is fracture segment height (m), and r w is the wellbore radius (m).The validation of this methodology against the traditional LGR method can be found in our previous work [30,31].The efficiency of the EDFM method is much higher than the LGR method, especially when dealing with multiple wells with a large number of fractures [33].It has been widely applied to model well interference due to complex fracture hits [34,35], automatic history matching for shale reservoirs [36][37][38], gas injection for enhanced unconventional oil recovery [39][40][41], and naturally fractured reservoir simulation [42].

History Matching
An actual Marcellus shale-gas well with available 486-day gas and water production data was selected to perform history matching.The well has 11 perforations stages with a cluster spacing of 15.24 m and was hydraulically fractured using 16,770 m 3 water and 2,895,956 kg proppant.The lateral length of horizontal wellbore is about 1382 m.There are 10 stages with nine clusters per stage and one stage with seven clusters, so 97 effective hydraulic fractures were assumed in the following simulation study.
We set up a field-scale shale-gas reservoir model with two-phase flow (gas and water) using a compositional numerical reservoir simulator [32], as shown in Figure 2. The gas type is methane.The dimension of the reservoir model without embedding hydraulic fractures is 1585 m × 622 m × 27.6 m and it is assumed that the reservoir thickness was fully penetrated by hydraulic fractures.Hence, the effective hydraulic fracture height is 27.6 m.The matrix grid size is 6 m × 12 m × 27.6 m in x, y, and z directions.
Energies 2019, 12, x FOR PEER REVIEW 4 of 16 where kNNC represents the matrix grid permeability for NNC type I and average fracture permeability for NNC types II and III (mD); ANNC represents the contact area between NNC grids (m 2 ); dNNC represents the connection distance between different NNC grids (m).It should be mentioned that the non-intrusive EDFM technology only deals with the transmissibility factor and does not involve relative phase mobility calculation.Hence, it can be applied in both black oil simulation and compositional simulation.
For the connection between fractures and wellbore, the following effective wellbore index (WI , mD-m) will be calculated [30]: where kf is the fracture permeability (mD), wf is fracture width (m), L is fracture segment length (m), W is fracture segment height (m), and rw is the wellbore radius (m).The validation of this methodology against the traditional LGR method can be found in our previous work [30,31].The efficiency of the EDFM method is much higher than the LGR method, especially when dealing with multiple wells with a large number of fractures [33].It has been widely applied to model well interference due to complex fracture hits [34,35], automatic history matching for shale reservoirs [36][37][38], gas injection for enhanced unconventional oil recovery [39][40][41], and naturally fractured reservoir simulation [42].

History Matching
An actual Marcellus shale-gas well with available 486-day gas and water production data was selected to perform history matching.The well has 11 perforations stages with a cluster spacing of 15.24 m and was hydraulically fractured using 16,770 m 3 water and 2,895,956 kg proppant.The lateral length of horizontal wellbore is about 1382 m.There are 10 stages with nine clusters per stage and one stage with seven clusters, so 97 effective hydraulic fractures were assumed in the following simulation study.
We set up a field-scale shale-gas reservoir model with two-phase flow (gas and water) using a compositional numerical reservoir simulator [32], as shown in Figure 2. The gas type is methane.The dimension of the reservoir model without embedding hydraulic fractures is 1585 m × 622 m × 27.6 m and it is assumed that the reservoir thickness was fully penetrated by hydraulic fractures.Hence, the effective hydraulic fracture height is 27.6 m.The matrix grid size is 6 m × 12 m × 27.6 m in x, y, and z directions.The number of matrix grids is 13,260.Hydraulic fractures were modeled using the non-intrusive EDFM method.After embedding 97 hydraulic fractures with fracture half-length of 127 m into matrix grids, an extra number of fracture grids of 2040 was generated along the x direction.The fracture grid size remains the same as the matrix grid size.The model dimension with hydraulic fractures becomes 1828.8 m × 622 m × 27.6 m.It is assumed that both matrix and fracture grids have the same gas-water relative permeability curve, which is from the work by Yu et al. [38].In addition, the experimental measurements of gas desorption are considered in the simulation model [38].The other reservoir and fracture parameters are listed in Table 1.During the history-matching process, we applied the measured flowing BHP for the reservoir simulation constraint, as shown in Figure 3a.In order to achieve good match results with measured gas and water flow rates, we mainly tune matrix permeability, fracture half-length, fracture conductivity, and fracture water saturation because these parameters are very sensitive for tuning and performing history matching.
The number of matrix grids is 13,260.Hydraulic fractures were modeled using the non-intrusive EDFM method.After embedding 97 hydraulic fractures with fracture half-length of 127 m into matrix grids, an extra number of fracture grids of 2040 was generated along the x direction.The fracture grid size remains the same as the matrix grid size.The model dimension with hydraulic fractures becomes 1828.8 m × 622 m × 27.6 m.It is assumed that both matrix and fracture grids have the same gas-water relative permeability curve, which is from the work by Yu et al. [38].In addition, the experimental measurements of gas desorption are considered in the simulation model [38].The other reservoir and fracture parameters are listed in Table 1.During the history-matching process, we applied the measured flowing BHP for the reservoir simulation constraint, as shown in Figure 3a.In order to achieve good match results with measured gas and water flow rates, we mainly tune matrix permeability, fracture half-length, fracture conductivity, and fracture water saturation because these parameters are very sensitive for tuning and performing history matching.It should be mentioned that the fracture water saturation is higher than the matrix water saturation because water is more difficult to transport into the matrix than into the fracture due to lower matrix permeability.Figure 3b,c present the comparison of gas and water flow rates between filed data and simulation results.The relative error of history-matching results for gas flow rate is about 10%.It can be illustrated that reasonable agreements were obtained.The final history-matching tuning parameters were determined as matrix permeability of 0.000525 mD, fracture half-length of 127 m, fracture conductivity of 0.73 mD-m, and fracture water saturation of 37.4%.It should be noted that fracture half-length has an important impact on multiple shale-gas well placement.In general, longer fracture half-length will result in larger well spacing.

Production Forecasting
After the history-matching process, we performed long-term production forecasting for 30 years.Since there are no available actual well control information for prediction, a constant flowing BHP of 3.45 MPa after the history-matching period was assumed for the simulation constraint, as illustrated in Figure 4a. Figure 4b,c show the prediction of gas flow rate and cumulative gas production incorporating a short-term field data, respectively.As shown, the estimated ultimate recovery (EUR) after 30 years for this shale-gas well is about 413 × 10 6 m 3 .Figure 4d,e show the prediction of water flow rate and cumulative water production incorporating the short-term field data, respectively.It can be seen that the water flow rate becomes negligible after around three years.The cumulative water production after 30 years for this shale-gas well is about 747 m 3 .Figure 5 compares the pressure distribution after the history-matching period and 30 years of production, clearly illustrating different drainage area and production interference intensity between hydraulic fractures.It can be observed that a stronger production interference between hydraulic fractures occures after 30 years of production when compared to that after the history-matching period.In addition, a larger effective drainage area can be generated after 30 years of production.It should be mentioned that the fracture water saturation is higher than the matrix water saturation because water is more difficult to transport into the matrix than into the fracture due to lower matrix permeability.Figures 3b,c present the comparison of gas and water flow rates between filed data and simulation results.The relative error of history-matching results for gas flow rate is about 10%.It can be illustrated that reasonable agreements were obtained.The final history-matching tuning parameters were determined as matrix permeability of 0.000525 mD, fracture half-length of 127 m, fracture conductivity of 0.73 mD-m, and fracture water saturation of 37.4%.It should be noted that fracture half-length has an important impact on multiple shale-gas well placement.In general, longer fracture half-length will result in larger well spacing.

Production Forecasting
After the history-matching process, we performed long-term production forecasting for 30 years.Since there are no available actual well control information for prediction, a constant flowing BHP of 3.45 MPa after the history-matching period was assumed for the simulation constraint, as illustrated in Figure 4a.Figures 4b,c show the prediction of gas flow rate and cumulative gas production incorporating a short-term field data, respectively.As shown, the estimated ultimate recovery (EUR) after 30 years for this shale-gas well is about 413 × 10 6 m 3 .Figures 4d,e show the prediction of water flow rate and cumulative water production incorporating the short-term field data, respectively.It can be seen that the water flow rate becomes negligible after around three years.The cumulative water production after 30 years for this shale-gas well is about 747 m 3 .Figure 5 compares the pressure distribution after the history-matching period and 30 years of production, clearly illustrating different drainage area and production interference intensity between hydraulic fractures.It can be observed that a stronger production interference between hydraulic fractures occures after 30 years of production when compared to that after the history-matching period.In addition, a larger effective drainage area can be generated after 30 years of production.

Varying Well Spacing without Natural Fractures
Based on the history-matching results from the field-scale shale-gas model, we extended the reservoir model to include two horizontal shale-gas wells with varying well spacing and each well has 97 hydraulic fractures.The effect of natural fractures was not considered in this basic model.The model dimension without embedding hydraulic fractures is 1585 m × 975 m × 27.6 m and fractures fully penetrate the reservoir thickness.The matrix grid size is 6 m × 12 m × 27.6 m in x, y, and z directions.The number of matrix grids is 20,800.Three different well spacing values such as 200 m, 300 m, and 400 m were investigated, as depicted in Figure 6.The placement of hydraulic fractures from two horizontal wells has a staggered pattern.As shown in Figure 6a, there is about 21% overlap of fracture half-length.The distance between neighboring fracture tips of two wells is 47.5 m and 148 m for well spacing of 300 m and 400 m, respectively, as shown in Figures 6b,c.After embedding 194 hydraulic fractures into matrix grids, an extra number of fracture grids of 4320 were generated along

Varying Well Spacing without Natural Fractures
Based on the history-matching results from the field-scale shale-gas model, we extended the reservoir model to include two horizontal shale-gas wells with varying well spacing and each well has 97 hydraulic fractures.The effect of natural fractures was not considered in this basic model.The model dimension without embedding hydraulic fractures is 1585 m × 975 m × 27.6 m and fractures fully penetrate the reservoir thickness.The matrix grid size is 6 m × 12 m × 27.6 m in x, y, and z directions.The number of matrix grids is 20,800.Three different well spacing values such as 200 m, 300 m, and 400 m were investigated, as depicted in Figure 6.The placement of hydraulic fractures from two horizontal wells has a staggered pattern.As shown in Figure 6a, there is about 21% overlap of fracture half-length.The distance between neighboring fracture tips of two wells is 47.5 m and 148 m for well spacing of 300 m and 400 m, respectively, as shown in Figures 6b,c.After embedding 194 hydraulic fractures into matrix grids, an extra number of fracture grids of 4320 were generated along

Varying Well Spacing without Natural Fractures
Based on the history-matching results from the field-scale shale-gas model, we extended the reservoir model to include two horizontal shale-gas wells with varying well spacing and each well has 97 hydraulic fractures.The effect of natural fractures was not considered in this basic model.The model dimension without embedding hydraulic fractures is 1585 m × 975 m × 27.6 m and fractures fully penetrate the reservoir thickness.The matrix grid size is 6 m × 12 m × 27.6 m in x, y, and z directions.The number of matrix grids is 20,800.Three different well spacing values such as 200 m, 300 m, and 400 m were investigated, as depicted in Figure 6.The placement of hydraulic fractures from two horizontal wells has a staggered pattern.As shown in Figure 6a, there is about 21% overlap of fracture half-length.The distance between neighboring fracture tips of two wells is 47.5 m and 148 m for well spacing of 300 m and 400 m, respectively, as shown in Figure 6b,c

Varying Well Spacing with Natural Fractures
Next, the effect of natural fractures with different numbers was considered in the basic reservoir model with varying well spacing.Four different numbers of natural fractures such as 100, 1000, 5000, and 10,000 were investigated in this study.It should be mentioned that hydraulic fractures are nonplanar if considering the interaction between pre-existing natural fracture with hydraulic fracture using fracture propagation models [6][7][8].However, hydraulic fractures are assumed to be planar in this study.We did not apply the fracture propagation model to predict complex fracture network.Figure 7 displays the natural fracture distribution of four different cases using an example of the well spacing of 300 m.We applied a statistical method to generate the locations of natural fractures, which are normally distributed in the reservoir model.For each case, there are two sets of natural fractures and each set has the same number of natural fractures.One set of natural fractures has an orientation of 45 degrees with respect to the x axis.Another set has an orientation of 135 degrees along the x axis.The natural fractures have a range of length from 30 m to 91 m.The natural fracture conductivity for each case is assumed to be 0.03 mD-m with natural fracture aperture of 0.003 m.The natural fractures fully penetrate the reservoir thickness.It should be mentioned that the impacts of natural fracture conductivity, length, and orientation are not examined in this study.After embedding natural and hydraulic fractures together into matrix grids with well spacing of 300 m, an extra number of fracture grids of 5520, 15,840, 62,000, and 119,440 were generated for the number of natural fractures of 100, 1000, 5000, 10,000, respectively, along the x direction.As can be seen that, the extra number of fractures grids is larger than the number of matrix grids of 20,800 for the 5000 and 10,000 natural fractures.With the increasing number of natural fractures, there are more interactions between

Varying Well Spacing with Natural Fractures
Next, the effect of natural fractures with different numbers was considered in the basic reservoir model with varying well spacing.Four different numbers of natural fractures such as 100, 1000, 5000, and 10,000 were investigated in this study.It should be mentioned that hydraulic fractures are non-planar if considering the interaction between pre-existing natural fracture with hydraulic fracture using fracture propagation models [6][7][8].However, hydraulic fractures are assumed to be planar in this study.We did not apply the fracture propagation model to predict complex fracture network.Figure 7 displays the natural fracture distribution of four different cases using an example of the well spacing of 300 m.We applied a statistical method to generate the locations of natural fractures, which are normally distributed in the reservoir model.For each case, there are two sets of natural fractures and each set has the same number of natural fractures.One set of natural fractures has an orientation of 45 degrees with respect to the x axis.Another set has an orientation of 135 degrees along the x axis.The natural fractures have a range of length from 30 m to 91 m.The natural fracture conductivity for each case is assumed to be 0.03 mD-m with natural fracture aperture of 0.003 m.The natural fractures fully penetrate the reservoir thickness.It should be mentioned that the impacts of natural fracture conductivity, length, and orientation are not examined in this study.After embedding natural and hydraulic fractures together into matrix grids with well spacing of 300 m, an extra number of fracture grids of 5520, 15,840, 62,000, and 119,440 were generated for the number of natural fractures of 100, 1000, 5000, 10,000, respectively, along the x direction.As can be seen that, the extra number of fractures grids is larger than the number of matrix grids of 20,800 for the 5000 and 10,000 natural fractures.
With the increasing number of natural fractures, there are more interactions between natural and hydraulic fractures, resulting in more well interference due to complex fracture connections.

Effect of Natural Fractures on Well Performance
The effects of different numbers of natural fractures on cumulative gas production are compared in Figures 8-10 for well spacings of 200, 300 and 400 m, respectively.It can be clearly observed that a larger number of natural fractures contributes to more production.The incremental EUR after 30 years compared to the case with well spacing of 200 m without natural fractures is about 1.5%, 13.8%, 69%, and 130.7% for the fracture number of 100, 1000, 5000, and 10,000, respectively.When compared to the case with well spacing of 300 m, the incremental EUR after 30 years of production is about 1.4%, 13.9%, 66.7%, and 124.3% for the fracture number of 100, 1000, 5000, and 10,000, respectively.When compared to the case with well spacing of 400 m, the incremental EUR after 30 years of production is about 1.6%, 15.3%, 72.6%, and 131% for the fracture number of 100, 1000, 5000, and 10,000, respectively.For 100 natural fractures, there is only a smaller number of natural fractures connecting with hydraulic fractures, leading to a slight increase of cumulative gas production.However, for 10,000 natural fractures, there is a much more complex connected fracture network between natural and hydraulic fractures, resulting in a big increase of contact area between fracture and matrix and a large increase of cumulative gas production.The CPU time for the model with well spacing of 300 m is about 5, 7 min 17 and 30 min, corresponding to the number of natural fractures of 100, 1000, 5000, and 10,000, respectively.Hence, the natural fractures play an important role in long-term multiple shale-gas well performance, which can be easily and efficiently handled by the non-intrusive EDFM technology.

Effect of Natural Fractures on Well Performance
The effects of different numbers of natural fractures on cumulative gas production are compared in Figures 8-10 for well spacings of 200, 300 and 400 m, respectively.It can be clearly observed that a larger number of natural fractures contributes to more production.The incremental EUR after 30 years compared to the case with well spacing of 200 m without natural fractures is about 1.5%, 13.8%, 69%, and 130.7% for the fracture number of 100, 1000, 5000, and 10,000, respectively.When compared to the case with well spacing of 300 m, the incremental EUR after 30 years of production is about 1.4%, 13.9%, 66.7%, and 124.3% for the fracture number of 100, 1000, 5000, and 10,000, respectively.When compared to the case with well spacing of 400 m, the incremental EUR after 30 years of production is about 1.6%, 15.3%, 72.6%, and 131% for the fracture number of 100, 1000, 5000, and 10,000, respectively.For 100 natural fractures, there is only a smaller number of natural fractures connecting with hydraulic fractures, leading to a slight increase of cumulative gas production.However, for 10,000 natural fractures, there is a much more complex connected fracture network between natural and hydraulic fractures, resulting in a big increase of contact area between fracture and matrix and a large increase of cumulative gas production.The CPU time for the model with well spacing of 300 m is about 5, 7 min 17 and 30 min, corresponding to the number of natural fractures of 100, 1000, 5000, and 10,000, respectively.Hence, the natural fractures play an important role in long-term multiple shale-gas well performance, which can be easily and efficiently handled by the non-intrusive EDFM technology.

Effect of Natural Fractures on Well Spacing Optimization
The effect of natural fracture density on the incremental gas production after 30 years under different well spacing is presented in Figure 11.Here, the natural fracture density refers to the number of natural fractures per 929 m 2 , which is about 0.06, 0.6, 3, and 6 for the total number of natural fractures of 100, 1000, 5000, and 10,000, respectively.The incremental gas production after 30 years means that the difference of cumulative gas production with and without natural fractures under a given well spacing.It can be clearly seen that the incremental gas production of different well spacing with and without natural fractures almost linearly increases with the increasing natural fracture density.In addition, the difference of incremental gas production between different well spacing increases with an increase in natural fracture density.The difference of incremental gas production between well spacing of 200 m and 400 m under the natural fracture density of 0.06, 0.6, 3 and 6 per 929 m 2 is about 3.2 × 10 6 m 3 , 31.2 × 10 6 m 3 , 125.6 × 10 6 m 3 , and 184.7 × 10 6 m 3 , respectively.Consequently, it can be implied that a larger well spacing is preferred for the shale-gas reservoir with a larger natural fracture density.

Effect of Natural Fractures on Well Spacing Optimization
The effect of natural fracture density on the incremental gas production after 30 years under different well spacing is presented in Figure 11.Here, the natural fracture density refers to the number of natural fractures per 929 m 2 , which is about 0.06, 0.6, 3, and 6 for the total number of natural fractures of 100, 1000, 5000, and 10,000, respectively.The incremental gas production after 30 years means that the difference of cumulative gas production with and without natural fractures under a given well spacing.It can be clearly seen that the incremental gas production of different well spacing with and without natural fractures almost linearly increases with the increasing natural fracture density.In addition, the difference of incremental gas production between different well spacing increases with an increase in natural fracture density.The difference of incremental gas production between well spacing of 200 m and 400 m under the natural fracture density of 0.06, 0.6, 3 and 6 per 929 m 2 is about 3.2 × 10 6 m 3 , 31.2 × 10 6 m 3 , 125.6 × 10 6 m 3 , and 184.7 × 10 6 m 3 , respectively.Consequently, it can be implied that a larger well spacing is preferred for the shale-gas reservoir with a larger natural fracture density.

Effect of Natural Fractures on Well Spacing Optimization
The effect of natural fracture density on the incremental gas production after 30 years under different well spacing is presented in Figure 11.Here, the natural fracture density refers to the number of natural fractures per 929 m 2 , which is about 0.06, 0.6, 3, and 6 for the total number of natural fractures of 100, 1000, 5000, and 10,000, respectively.The incremental gas production after 30 years means that the difference of cumulative gas production with and without natural fractures under a given well spacing.It can be clearly seen that the incremental gas production of different well spacing with and without natural fractures almost linearly increases with the increasing natural fracture density.In addition, the difference of incremental gas production between different well spacing increases with an increase in natural fracture density.The difference of incremental gas production between well spacing of 200 m and 400 m under the natural fracture density of 0.06, 0.6, 3 and 6 per 929 m 2 is about 3.2 × 10 6 m 3 , 31.2 × 10 6 m 3 , 125.6 × 10 6 m 3 , and 184.7 × 10 6 m 3 , respectively.Consequently, it can be implied that a larger well spacing is preferred for the shale-gas reservoir with a larger natural fracture density.

Effect of Natural Fractures on Pressure Distribution
Figures 12 and 13 compare pressure distribution after one and 30 years of production under different well spacing with and without natural fractures.The extreme case with 10,000 natural fractures was plotted for comparison.As shown, different drainage area at early and later time of production can be clearly observed.In addition, a stronger well interference between two wells after 30 years under a smaller well spacing occurs.Especially, the drainage area expands more when considering a bigger number of natural fractures.

Effect of Natural Fractures on Pressure Distribution
Figures 12 and 13 compare pressure distribution after one and 30 years of production under different well spacing with and without natural fractures.The extreme case with 10,000 natural fractures was plotted for comparison.As shown, different drainage area at early and later time of production can be clearly observed.In addition, a stronger well interference between two wells after 30 years under a smaller well spacing occurs.Especially, the drainage area expands more when considering a bigger number of natural fractures.

Conclusions
We applied the non-intrusive EDFM technology as well as a compositional reservoir simulator to perform shale-gas two-phase flow simulations.The impact of natural fractures on multiple shalegas well performance with varying well spacing was investigated based on history matching results of an actual shale-gas well from Marcellus shale formation.The following conclusions can be drawn from this study: (1) The effective matrix and fracture properties were obtained based on good history matching results such as matrix permeability of 0.000525 mD, fracture half-length of 127 m, fracture conductivity of 0.73 mD-m, and fracture water saturation of 37.4%.(2) The EUR and cumulative water production after 30 years of the actual shale-gas well were determined as 413 × 10 6 m 3 and 747 m 3 , respectively.(3) The effect of natural fractures on two shale-gas well performance almost linearly increases with the increasing number of natural fractures.For example, the natural fractures contribute to the incremental EUR after 30 years of 1.4%, 13.9%, 66.7%, and 124.3% for the well spacing of 300 m with fracture number of 100, 1000, 5000, and 10,000, respectively.(4) The CPU time for the model with well spacing of 300 m is about 30 min when dealing with 10,000 natural fractures based on the non-intrusive EDFM technology.(5) A larger number of natural fractures is easier to form a more complex connected fracture network with hydraulic fractures, resulting in a higher well productivity.(6) The difference of well performance between different well spacing increases with the increasing natural fracture density.For example, the difference of well performance between well spacing of 200 m and 400 m is about 3.2 × 10 6 m 3 , 31.2 × 10 6 m 3 , 125.6 × 10 6 m 3 , and 184.7 × 10 6 m 3 corresponding to the natural fracture density of 0.06, 0.6, 3 and 6 per 929 m 2 , respectively.(7) This study provides a better understanding of the impact of complex natural fractures on multiple shale-gas well performance with varying well spacing.A larger well spacing is suggested when the shale-gas reservoir has a larger natural fracture density.

Conclusions
We applied the non-intrusive EDFM technology as well as a compositional reservoir simulator to perform shale-gas two-phase flow simulations.The impact of natural fractures on multiple shale-gas well performance with varying well spacing was investigated based on history matching results of an actual shale-gas well from Marcellus shale formation.The following conclusions can be drawn from this study: (1) The effective matrix and fracture properties were obtained based on good history matching results such as matrix permeability of 0.000525 mD, fracture half-length of 127 m, fracture conductivity of 0.73 mD-m, and fracture water saturation of 37.4%.(2) The EUR and cumulative water production after 30 years of the actual shale-gas well were determined as 413 × 10 6 m 3 and 747 m 3 , respectively.(3) The effect of natural fractures on two shale-gas well performance almost linearly increases with the increasing number of natural fractures.For example, the natural fractures contribute to the incremental EUR after 30 years of 1.4%, 13.9%, 66.7%, and 124.3% for the well spacing of 300 m with fracture number of 100, 1000, 5000, and 10,000, respectively.(4) The CPU time for the model with well spacing of 300 m is about 30 min when dealing with 10,000 natural fractures based on the non-intrusive EDFM technology.(5) A larger number of natural fractures is easier to form a more complex connected fracture network with hydraulic fractures, resulting in a higher well productivity.(6) The difference of well performance between different well spacing increases with the increasing natural fracture density.For example, the difference of well performance between well spacing of 200 m and 400 m is about 3.2 × 10 6 m 3 , 31.2 × 10 6 m 3 , 125.6 × 10 6 m 3 , and 184.7 × 10 6 m 3 corresponding to the natural fracture density of 0.06, 0.6, 3 and 6 per 929 m 2 , respectively.(7) This study provides a better understanding of the impact of complex natural fractures on multiple shale-gas well performance with varying well spacing.A larger well spacing is suggested when the shale-gas reservoir has a larger natural fracture density.

Figure 1 .
Figure 1.The non-intrusive EDFM technology in conjunction with commercial reservoir simulator to efficiently model 3D complex fractures: (a) Physical domain; (b) Computational domain [30].

Figure 1 .
Figure 1.The non-intrusive EDFM technology in conjunction with commercial reservoir simulator to efficiently model 3D complex fractures: (a) Physical domain; (b) Computational domain [30].

Figure 2 .
Figure 2. The field-scale shale-gas reservoir model with two-phase flow including one horizontal well and 97 hydraulic fractures.

Figure 2 .
Figure 2. The field-scale shale-gas reservoir model with two-phase flow including one horizontal well and 97 hydraulic fractures.

Figure 3 .
Figure 3.Comparison of flowing BHP and gas and water flow rates between filed data and simulation results: (a) Flowing BHP; (b) Gas flow rate; (c) Water flow rate.

Figure 3 .
Figure 3.Comparison of flowing BHP and gas and water flow rates between filed data and simulation results: (a) Flowing BHP; (b) Gas flow rate; (c) Water flow rate.

Figure 5 .
Figure 5.Comparison of pressure distribution at different times: (a) After the history-matching period; (b) After 30 years of production.

Figure 6 .
Figure 6.The basic reservoir simulation model including two horizontal wells and 97 hydraulic fractures for each well: (a) Well spacing of 200 m; (b) Well spacing of 300 m; (c) Well spacing of 400 m.

Figure 6 .
Figure 6.The basic reservoir simulation model including two horizontal wells and 97 hydraulic fractures for each well: (a) Well spacing of 200 m; (b) Well spacing of 300 m; (c) Well spacing of 400 m.

Figure 7 .
Figure 7.The basic reservoir simulation model including two horizontal wells with well spacing of 300 m and four different numbers of natural fractures: (a) 100 natural fractures; (b) 1000 natural fractures; (c) 5000 natural fractures; (d) 10,000 natural fractures.

Figure 7 .
Figure 7.The basic reservoir simulation model including two horizontal wells with well spacing of 300 m and four different numbers of natural fractures: (a) 100 natural fractures; (b) 1000 natural fractures; (c) 5000 natural fractures; (d) 10,000 natural fractures.

Figure 8 .
Figure 8.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 200 m (NF in the legend represents natural fractures).

Figure 9 .
Figure 9.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 300 m (NF in the legend represents natural fractures).

Figure 8 .
Figure 8.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 200 m (NF in the legend represents natural fractures).

Figure 8 .
Figure 8.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 200 m (NF in the legend represents natural fractures).

Figure 9 .
Figure 9.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 300 m (NF in the legend represents natural fractures).

Figure 9 .
Figure 9.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 300 m (NF in the legend represents natural fractures).

Figure 10 .
Figure 10.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 400 m (NF in the legend represents natural fractures).

Figure 11 .
Figure 11.Effect of natural fracture density on the incremental gas production after 30 years under different well spacing.

Figure 10 .
Figure 10.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 400 m (NF in the legend represents natural fractures).

Figure 10 .
Figure 10.Effect of different numbers of natural fractures on cumulative gas production for the reservoir model with well spacing of 400 m (NF in the legend represents natural fractures).

Figure 11 .
Figure 11.Effect of natural fracture density on the incremental gas production after 30 years under different well spacing.

Figure 11 .
Figure 11.Effect of natural fracture density on the incremental gas production after 30 years under different well spacing.

Figure 12 .
Figure 12.Comparison of pressure distribution after 1 and 30 years of production under different well spacing without natural fractures: (a) After 1 year of production with well spacing 200 m; (b) After 30 years of production with well spacing 200 m; (c) After 1 year of production with well spacing 300 m; (d) After 30 years of production with well spacing 300 m; (e) After 1 year of production with well spacing 400 m; (f) After 30 years of production with well spacing 400 m.

Figure 12 . 16 5. 3 .Figure 12 .Figure 13 .
Figure 12.Comparison of pressure distribution after 1 and 30 years of production under different well spacing without natural fractures: (a) After 1 year of production with well spacing 200 m; (b) After 30 years of production with well spacing 200 m; (c) After 1 year of production with well spacing 300 m; (d) After 30 years of production with well spacing 300 m; (e) After 1 year of production with well spacing 400 m; (f) After 30 years of production with well spacing 400 m.

Figure 13 .
Figure 13.Comparison of pressure distribution after one and 30 years of production under different well spacing with 10,000 natural fractures: (a) After one year of production with well spacing 200 m; (b) After 30 years of production with well spacing 200 m; (c) After one year of production with well spacing 300 m; (d) After 30 years of production with well spacing 300 m; (e) After one year of production with well spacing 400 m; (f) After 30 years of production with well spacing 400 m.

Table 1 .
Reservoir and fracture properties used in the field-scale shale-gas model.

Table 1 .
Reservoir and fracture properties used in the field-scale shale-gas model.