A Data-Driven Workflow Approach to Optimization of Fracture Spacing in Multi-Fractured Shale Oil Wells

A data-driven workflow approach is presented in this study for optimizing fracture spacing of multifractured horizontal wells (MFHW) in shale oil reservoirs. The workflow employs a simple well productivity model for the initial design of hydraulic fracturing well completions. This provides a transparent approach to the identification of key fracturing parameters affecting well productivity. The workflow uses transient pressure or production data to identify fracture interference. This offers a reliable and cost-effective means for assessment of well production potential in terms of optimization of fracture spacing in the MFHW. Result of a field case study indicated that three wells were drilled in an area with dense natural fractures, and the fracture spacing of MFHW in this area was short enough to effectively drain the stimulated reservoir volume (SRV), while the other three wells were drilled in an area with less natural fractures, and the fracture spacing of MFHW in this area could be shortened to double well productivity.


Introduction
Thanks to multifractured horizontal well (MFHW) technology, development of shale oil fields in the past five years has made United States energy self-sufficient. However, the results of the technology in other regions of the world are mixed. This is partially attributed to inadequate effort in optimizing a completed design of MHFW.
Modern MFHW techniques include simul frac, zipper frac, and modified zipper frac. Simul frac creates simultaneous hydraulic fractures in symmetric stages in two wells. Zipper frac creates alternate hydraulic fractures in symmetric stages in two wells. Modified zipper frac creates hydraulic fractures in nonsymmetric stages in two alternate wells [1]. These fracturing techniques create symmetric, transverse fractures of ideally the same length and spacing. This common structure of multifractured horizontal wells allows for using general mathematical models to predict well productivity in the stage of well completion design.
The available methods for predicting shale well productivity include analytical transient flow models [2][3][4][5][6] and numerical computer models [7][8][9][10][11]. Analytical transient flow models are mostly used in pressure transient test analyses rather than in optimizing well completion design. Numerical computer models are flexible in handling systems with nonsymmetrical fractures and multiphase flows [12]. In addition to model complicity, the uncertainty in locations of natural fractures is a concern regarding the accuracy of computer simulation results. Based on the work of Guo et al. [13], Zhang et al. [14] presented an analytical solution for estimating long-term productivity of multifractured shale gas wells by coupling the linear flow in the reservoir and the linear flow in hydraulic fractures.

Workflow for Optimization
The proposed workflow to optimize fracture spacing is based on industry practice [16], with the addition of data analysis, to identify fracture interference. It is summarized in the following steps: 1. Design a fracturing job for the first well in the area for a desirable production rate using a well productivity model. Fracture spacing is selected based on horizontal wellbore length, volumes of fracturing fluid and proppant, and completion tools. 2. Execute the fracturing job using the designed fracture spacing and other parameters. 3. Run a pressure transient test on the well if possible, and put the well into production. 4. Perform transient pressure or transient production rate data analyses to identify fracture interference. 5. If well completion permits, refracture the well based on the identified fracture interference. 6. Modify the well completion design including the fracturing treatment design and/or the spacing between perforation clusters for the next well on the basis of fracture interference in the previous well. Some of these workflow steps are further outlined in the subsections that follow.

Fracturing Design
The physics of fluid flow in shale oil reservoirs is adequately described in the literature. The assumptions employed in this study were oversimplified, as the desorption of hydrocarbon from rock matrix and non-Darcy flow were not considered. However, it is a fact that pressure transient and production data support the mathematical models derived from Darcy's law. This implies that fluid flow in shale oil reservoirs is dominated by Darcy's law.
Fracturing design can be guided by mathematical models derived from Darcy's law for shale gas and oil well production under pseudosteady flow conditions. For multifractured horizontal oil wells, Li et al.'s model for oil production rate is expressed as [15]:

Workflow for Optimization
The proposed workflow to optimize fracture spacing is based on industry practice [16], with the addition of data analysis, to identify fracture interference. It is summarized in the following steps:

1.
Design a fracturing job for the first well in the area for a desirable production rate using a well productivity model. Fracture spacing is selected based on horizontal wellbore length, volumes of fracturing fluid and proppant, and completion tools.

2.
Execute the fracturing job using the designed fracture spacing and other parameters.

3.
Run a pressure transient test on the well if possible, and put the well into production.

4.
Perform transient pressure or transient production rate data analyses to identify fracture interference.

5.
If well completion permits, refracture the well based on the identified fracture interference. 6.
Modify the well completion design including the fracturing treatment design and/or the spacing between perforation clusters for the next well on the basis of fracture interference in the previous well.
Some of these workflow steps are further outlined in the subsections that follow.

Fracturing Design
The physics of fluid flow in shale oil reservoirs is adequately described in the literature. The assumptions employed in this study were oversimplified, as the desorption of hydrocarbon from rock matrix and non-Darcy flow were not considered. However, it is a fact that pressure transient and production data support the mathematical models derived from Darcy's law. This implies that fluid flow in shale oil reservoirs is dominated by Darcy's law.
Fracturing design can be guided by mathematical models derived from Darcy's law for shale gas and oil well production under pseudosteady flow conditions. For multifractured horizontal oil wells, Li et al.'s model for oil production rate is expressed as [15]: where c is expressed as: where q o is the well production rate in stb/d, n f is the number of perforation clusters (hydraulic fractures), k m is the matrix permeability in md, h f is the fracture height in ft, p is the average formation pressure in psia, p w is the wellbore pressure in psia, µ o is the oil viscosity in cp, B o is the oil formation factor in rb/stb, S f is the fracture spacing in ft, e is an exponential function, x f is the hydraulic fracture half-length in ft, k f is the fracture permeability in md, and w is the average fracture width in inches. Substituting Equation (2) into Equation (1) yields: Although Equation (3) shows that maximum well productivity can be obtained by optimizing the effects of the number of hydraulic fractures, hydraulic fracture spacing, and hydraulic fracture width, this is not true in reality as fracture complexity and additional factors such as frac hits also exist. As pointed out by Potapenko et al. [17], in most of the unconventional reservoirs where the majority of multistage wells are drilled, the fracture geometry is quite complex and interaction between fractures starts very early. In this case, spacing between perforation clusters should not have much impact on well productivity. However, in situations where a significant portion of hydraulic fractures is damaged during well startup, resulting in disconnections of such fractures from the wellbore, having small spacing between perforation clusters still appears to be beneficial because disconnected fractures still may contribute to production through other nondamaged perforation clusters [17]. The issues of fracture complexity, interaction, and disconnection are quite difficult to consider in mathematical modeling because there is a lack of data to explicitly describe their configurations and dynamics [18]. These complication issues are assumed to be negligible in the analysis that follows.
For a given horizontal wellbore length L, well operators tend to increase the number of perforation clusters during the well completion optimization process. It brings some additional considerations about fracture length, the number of frac stages, and the volume of frac treatment. The number of clusters is expressed as: For a given total volume of fracture proppant V f in ft 3 , the following material balance holds: It should be mentioned that w in Equation (4) is in inches, and V f is the bulk volume (not the physical volume) of proppant.
Substituting Equation (4) into Equation (5) and rearranging the latter gives: Substituting Equations (4) and (6) into Equation (3) results in: The average reservoir pressure is [15]: where p e is the reservoir pressure in psia. It can be shown that in the practical ranges of parameters the sum of the two terms in the bracket in the denominator is very close to unity. Therefore, if fracture spacing is reduced, both the number of clusters and the total fracture surface will increase. Under the assumption of reservoir linear flow (RLF) with infinite fracture conductivity, the well production rate is proportional to the total fracture surface area. This supports the concept of massive volume fracturing where many clusters with the shortest possible spacings are used for pumping massive proppant into the created hydraulic fractures.
Although cluster spacing should be theoretically as low as possible for maximizing well productivity, a minimum required cluster spacing (MRCS) has to be considered in practice if temporal change and well lifetime are not considered. MRCS is controlled by (1) well completion design and equipment limitations and (2) hydraulic fracture interference. In well completion design, frac plugs and perforation clusters themselves are spaced apart, and casing couplings have to be avoided in perforation clusters. Perf gun can carry a limited number of charges, which are used for the creation of perf clusters. Adding a higher number of perf clusters may require a higher number of perf trips, which impacts overall cost and the operational efficiency. For the concern of hydraulic fracture interferences, sever connections of hydraulic fractures from different perforation clusters are usually not desirable, although overlapping of the frac networks from different perforation clusters may be quite beneficial to well productivity. This interference can be assessed using transient pressure or transient production data analysis.

Transient Pressure Analysis for Fracture Interference Identification
The theoretical basis of pressure transient data analysis can be found in numerous references [19][20][21]. Modern computer software packages available for pressure transient and production data analyses include F.A.S.T. WellTest [22] and PanSystem [23]. Applications of the pressure transient data analysis theory to shale gas and oil wells are shown by a number of investigators such as Shan et al. [24], Pang et al. [25], and He et al. [26]. Figure 2 presents a diagnostic plot of pressure transient data for flow regime identification where the vertical axis is the pressure derivative defined by: where ∆p is the pressure change, defined as reservoir pressure minus flowing bottom hole pressure for drawdown tests and shut-in bottom pressure minus the last flowing bottom hole pressure before shutting-in for pressure buildup tests. t is the test time defined as the flow time for drawdown tests and shut-in time for pressure buildup tests. Mathematical models describing reservoir transient linear flow to fractures are found in the literature [27,28]. For a well with infinite conductivity fractures, reservoir linear flow (RLF) can be identified by the half-slope of the pressure derivative data versus time data plotted on a log-log scale. For a well with finite conductivity fractures, a bilinear flow may occur in the fracture and in the formation matrix during the initial stages. Pressure derivative data versus time data plotted on a log-log scale should show a straight line with a slope of 0.25 [29]. After the pressure change propagates to the midline between fractures, pressure derivative data versus time data plotted on a log-log scale should show a straight line with the unit slope for boundary dominated flow (BDF) during the "depletion" in the fractured volume. If a slope value greater than 0.5 is observed soon after wellbore storage, interference between fractures is indicated. "depletion" in the fractured volume. If a slope value greater than 0.5 is observed soon after wellbore storage, interference between fractures is indicated.

Transient Rate Analysis for Fracture Interference Identification
Applications of the transient production data analysis theory to shale oil wells are also demonstrated by previous investigators [30][31][32]. Figure 3 shows a diagnostic plot of rate transient data for flow regime identification. For a well with infinite conductivity fractures, the production rate data versus time data plotted on a log-log scale should show a straight line with a slope of -0.5 during the RLF. For a well with finite conductivity fractures, a bilinear flow may occur in the fracture and in the formation matrix during the initial stages. Rate data versus time data plotted on a log-log scale should show a straight line with a slope of -0.25. After the depth of investigation propagates to the midline between fractures, rate data versus time data plotted on a log-log scale should show a straight line with a slope of -1 for BDF. If a slope value between -0.5 and -1 is observed early on, interference between fractures is indicated.

Transient Rate Analysis for Fracture Interference Identification
Applications of the transient production data analysis theory to shale oil wells are also demonstrated by previous investigators [30][31][32]. Figure 3 shows a diagnostic plot of rate transient data for flow regime identification. For a well with infinite conductivity fractures, the production rate data versus time data plotted on a log-log scale should show a straight line with a slope of −0.5 during the RLF. For a well with finite conductivity fractures, a bilinear flow may occur in the fracture and in the formation matrix during the initial stages. Rate data versus time data plotted on a log-log scale should show a straight line with a slope of −0.25. After the depth of investigation propagates to the midline between fractures, rate data versus time data plotted on a log-log scale should show a straight line with a slope of −1 for BDF. If a slope value between −0.5 and −1 is observed early on, interference between fractures is indicated. "depletion" in the fractured volume. If a slope value greater than 0.5 is observed soon after wellbore storage, interference between fractures is indicated.

Transient Rate Analysis for Fracture Interference Identification
Applications of the transient production data analysis theory to shale oil wells are also demonstrated by previous investigators [30][31][32]. Figure 3 shows a diagnostic plot of rate transient data for flow regime identification. For a well with infinite conductivity fractures, the production rate data versus time data plotted on a log-log scale should show a straight line with a slope of -0.5 during the RLF. For a well with finite conductivity fractures, a bilinear flow may occur in the fracture and in the formation matrix during the initial stages. Rate data versus time data plotted on a log-log scale should show a straight line with a slope of -0.25. After the depth of investigation propagates to the midline between fractures, rate data versus time data plotted on a log-log scale should show a straight line with a slope of -1 for BDF. If a slope value between -0.5 and -1 is observed early on, interference between fractures is indicated.

Field Case Studies
The within the Upper Cretaceous Mid-Dip Tuscaloosa trend, which occurs updip of the Lower Cretaceous shelf margin. In the Mississippi, production from the Mid-Dip trend extends 150 mi (240 km) to the east-southeast from the Mississippi River in a belt 30 to 60 mi (50-100 km) wide. Reservoir rocks are fine-to medium-grained sublitharenite in the stratigraphic unit of the Lower Tuscaloosa formation in the Late Cretaceous, Cenomanian, deposited in the fluvial meander belt environment. The productive facies are point bars. Rock total porosity is 21% to 27% with 13% to 15% primary, 3% to 4% secondary dissolution, and 5% to 8% microporosity. Core porosity ranges 10% to 35% with an average value of 24%. The core permeability measured with air ranges from 0.1 to 100 md with an average value of 10 md. The initial water saturation is between 40% and 75%, averaging at 55%. The reservoir is composed of Q and Q2 sandstone members at an average depth of 10,770 ft (3283 m) in an area of 6 × 3 mi 2 (9.7 × 4.8 km 2 ). The productive area is 6200 acres (2510 ha.) The hydrocarbon column height is 100 ft (30. TMS needs to use multifractured horizontal wells to become rejuvenated into a key exploration target for the industry as a leading oil field. Production rate and cumulative production data in May 2018 were gathered from 80 TMS wells in Louisiana and Mississippi. The majority of TMS wells were characterized by transient production behaviors.
Selection of wells for analysis in this study was based on the bubble map of the initial well productivity shown in Figure 4. Wells were selected from the Louisiana side of TMS with analyzable production decline curves. Only wells with initial oil production rates between 300 and 1200 stb/d were selected so wells with severe formation damage could be excluded. Wells with abnormal water cuts and gas-oil ratios were also excluded. The Tuscaloosa Marine Shale (TMS) formation, deposited during the Upper Cretaceous, extends from Louisiana to the southern portion of the Mississippi. Little Creek Field is located in Lincoln and Pike Counties, southwestern Mississippi, on the south rim of the Mississippi Salt Basin. It is within the Upper Cretaceous Mid-Dip Tuscaloosa trend, which occurs updip of the Lower Cretaceous shelf margin. In the Mississippi, production from the Mid-Dip trend extends 150 mi (240 km) to the eastsoutheast from the Mississippi River in a belt 30 to 60 mi (50-100 km) wide. Reservoir rocks are fineto medium-grained sublitharenite in the stratigraphic unit of the Lower Tuscaloosa formation in the Late Cretaceous, Cenomanian, deposited in the fluvial meander belt environment. The productive facies are point bars. Rock total porosity is 21% to 27% with 13% to 15% primary, 3% to 4% secondary dissolution, and 5% to 8% microporosity. Core porosity ranges 10% to 35% with an average value of 24%. The core permeability measured with air ranges from 0.1 to 100 md with an average value of 10 md. The initial water saturation is between 40% and 75%, averaging at 55%. The reservoir is composed of Q and Q2 sandstone members at an average depth of 10,770 ft (3283 m) in an area of 6 × 3 mi 2 (9.7 × 4. TMS needs to use multifractured horizontal wells to become rejuvenated into a key exploration target for the industry as a leading oil field. Production rate and cumulative production data in May 2018 were gathered from 80 TMS wells in Louisiana and Mississippi. The majority of TMS wells were characterized by transient production behaviors. Selection of wells for analysis in this study was based on the bubble map of the initial well productivity shown in Figure 4. Wells were selected from the Louisiana side of TMS with analyzable production decline curves. Only wells with initial oil production rates between 300 and 1200 stb/d were selected so wells with severe formation damage could be excluded. Wells with abnormal water cuts and gas-oil ratios were also excluded.  Figure 5 shows the statistical results of slopes of 44 TMS wells that were used to identify RLF. It should be mentioned that 11 TMS wells were not analyzed owing to their abrupt changes in production rate over time. We saw that the slopes of 15 wells were in the range of -0.52 and -0.48. Based on the log-log plot of data from 44 oil wells in the TMS trend, a slope of -0.5 was observed for some wells but not all wells. The correlation coefficient between the initial production rate and the slope was -0.076, which meant that there was no linear correlation.  Figure 5 shows the statistical results of slopes of 44 TMS wells that were used to identify RLF. It should be mentioned that 11 TMS wells were not analyzed owing to their abrupt changes in production rate over time. We saw that the slopes of 15 wells were in the range of −0.52 and −0.48. Based on the log-log plot of data from 44 oil wells in the TMS trend, a slope of −0.5 was observed for some wells but not all wells. The correlation coefficient between the initial production rate and the slope was −0.076, which meant that there was no linear correlation.  Figure 6 shows the log-log diagnostic plot of production rate data from Well #1. This well was drilled and completed with a horizontal length of 9102 ft and a 29-stage fracturing operation. Four perforation clusters were made in each stage, giving an average cluster spacing of 69 ft. Time data up to five months showed a -0.51 slope (very close to -0.5), indicating RLF. Late time production rate data versus material balance time followed a linear trend line with a slope of -0.78. A plot of the same production data versus actual production time formed a linear trend line with a slope of -1.03, indicating typical behavior for the boundary-dominated flow. Figure 7 shows the log-log diagnostic plot of production rate data from Well #2. This well was 11,300 ft deep with an effective lateral length of 4791 ft. Time data up to five months showed a -0.51 slope (very close to -0.5), indicating RLF. RLF ended at a production rate of about 226 stb/d. Late time production rate data versus material balance time followed a linear trend line with a slope of -0.76. A plot of the same production data versus actual time formed a linear trend line with a slope of -1.03, indicating typical behavior for BDF. Figure 8 demonstrates the log-log diagnostic plot of production rate data for Well #3. This well was drilled and completed with a 6099 ft lateral length and a 24-stage fracturing operation. Time data up to four months showed a -0.50 slope, indicating RLF. Late time production rate data versus material balance time followed a linear trend line with a slope of -0.72. A plot of the same production data versus actual production time formed a linear trend line with a slope of -1.00. The pseudosteady production rate was about 126 stb/d at the beginning of the BDF.     Figure 6 shows the log-log diagnostic plot of production rate data from Well #1. This well was drilled and completed with a horizontal length of 9102 ft and a 29-stage fracturing operation. Four perforation clusters were made in each stage, giving an average cluster spacing of 69 ft. Time data up to five months showed a -0.51 slope (very close to -0.5), indicating RLF. Late time production rate data versus material balance time followed a linear trend line with a slope of -0.78. A plot of the same production data versus actual production time formed a linear trend line with a slope of -1.03, indicating typical behavior for the boundary-dominated flow. Figure 7 shows the log-log diagnostic plot of production rate data from Well #2. This well was 11,300 ft deep with an effective lateral length of 4791 ft. Time data up to five months showed a -0.51 slope (very close to -0.5), indicating RLF. RLF ended at a production rate of about 226 stb/d. Late time production rate data versus material balance time followed a linear trend line with a slope of -0.76. A plot of the same production data versus actual time formed a linear trend line with a slope of -1.03, indicating typical behavior for BDF. Figure 8 demonstrates the log-log diagnostic plot of production rate data for Well #3. This well was drilled and completed with a 6099 ft lateral length and a 24-stage fracturing operation. Time data up to four months showed a -0.50 slope, indicating RLF. Late time production rate data versus material balance time followed a linear trend line with a slope of -0.72. A plot of the same production data versus actual production time formed a linear trend line with a slope of -1.00. The pseudosteady production rate was about 126 stb/d at the beginning of the BDF.   Figure 7 shows the log-log diagnostic plot of production rate data from Well #2. This well was 11,300 ft deep with an effective lateral length of 4791 ft. Time data up to five months showed a −0.51 slope (very close to −0.5), indicating RLF. RLF ended at a production rate of about 226 stb/d. Late time production rate data versus material balance time followed a linear trend line with a slope of −0.76. A plot of the same production data versus actual time formed a linear trend line with a slope of −1.03, indicating typical behavior for BDF. Figure 8 demonstrates the log-log diagnostic plot of production rate data for Well #3. This well was drilled and completed with a 6099 ft lateral length and a 24-stage fracturing operation. Time data up to four months showed a −0.50 slope, indicating RLF. Late time production rate data versus material balance time followed a linear trend line with a slope of −0.72. A plot of the same production data versus actual production time formed a linear trend line with a slope of −1.00. The pseudosteady production rate was about 126 stb/d at the beginning of the BDF.   Table 1 presents a summary of the estimated reservoir and well completion/fracture data for these three wells. It was assumed that fractures were created from all perf clusters for simplicity. The fracture height used in this analysis was discounted from the pay zone thickness to the net pay. It should be mentioned that typical values were taken from the area for all wells because there was a lack of fracture job data. The formation pressure gradient and the bottom hole pressure gradient were assumed to be 0.52 and 0.35 psia/ft, respectively. Matrix permeability was estimated using the following equation [33,34]: where tehs is the time at the end of the half slope in days. The average fracture spacing can be estimated using the following equations [33,34]:   Table 1 presents a summary of the estimated reservoir and well completion/fracture data for these three wells. It was assumed that fractures were created from all perf clusters for simplicity. The fracture height used in this analysis was discounted from the pay zone thickness to the net pay. It should be mentioned that typical values were taken from the area for all wells because there was a lack of fracture job data. The formation pressure gradient and the bottom hole pressure gradient were assumed to be 0.52 and 0.35 psia/ft, respectively. Matrix permeability was estimated using the following equation [33,34]: where tehs is the time at the end of the half slope in days. The average fracture spacing can be estimated using the following equations [33,34]:  Table 1 presents a summary of the estimated reservoir and well completion/fracture data for these three wells. It was assumed that fractures were created from all perf clusters for simplicity. The fracture height used in this analysis was discounted from the pay zone thickness to the net pay. It should be mentioned that typical values were taken from the area for all wells because there was a lack of fracture job data. The formation pressure gradient and the bottom hole pressure gradient were assumed to be 0.52 and 0.35 psia/ft, respectively. Matrix permeability was estimated using the following equation [33,34]: where t ehs is the time at the end of the half slope in days. The average fracture spacing can be estimated using the following equations [33,34]: where η is the diffusivity in md·psia·cp −1 , m is the slope of the inverse of production rate versus √ t in day 1/2 ·stb −1 , and c t is the total compressibility in psi −1 . Figure 9 illustrates the calculated well productivity curve generated using the data as input to Equation (7). It showed that the model-predicted well production rate at the beginning of BDF was 121 stb/d (at S f =69 ft), which was 0.83% higher than the observed value of 120 stb/d. The production rate at the beginning of BDF of Well #1 for S f = 15 ft was 89% higher than that for S f = 69 ft. This was explained by Wattenbarger's solution for the initial production rate [34]. It demonstrated that well production rate in the RLF period was proportional to the total fracture surface area. For a given horizontal wellbore length, the total fracture surface is dependent on the number of hydraulic fractures. Therefore, decreasing fracture spacing can increase production rate if the fracture half-length is held constant. Similarly, this indicated that the model-calculated well production rate at the beginning of BDF of Well #2 was 94 stb/d (at S f = 63 ft), which was 2.08% lower than the observed value of 96 stb/d. When the fracture spacing was reduced to 15 ft, the production rate at the beginning of BDF should be 71% higher than that the original one, assuming no other complex factors were missing in the model development. Of course, the decrease in perf cluster spacing may have less than ideal impacts on well productivity as a result of other complex factors that were not considered in the model development.


where η is the diffusivity in md•psia•cp -1 , m is the slope of the inverse of production rate versus t in day 1/2 •stb -1 , and ct is the total compressibility in psi -1 .  Figure 9 illustrates the calculated well productivity curve generated using the data as input to Equation (7). It showed that the model-predicted well production rate at the beginning of BDF was 121 stb/d (at Sf =69 ft), which was 0.83% higher than the observed value of 120 stb/d. The production rate at the beginning of BDF of Well #1 for Sf = 15 ft was 89% higher than that for Sf = 69 ft. This was explained by Wattenbarger's solution for the initial production rate [34]. It demonstrated that well production rate in the RLF period was proportional to the total fracture surface area. For a given horizontal wellbore length, the total fracture surface is dependent on the number of hydraulic fractures. Therefore, decreasing fracture spacing can increase production rate if the fracture halflength is held constant. Similarly, this indicated that the model-calculated well production rate at the beginning of BDF of Well #2 was 94 stb/d (at Sf = 63 ft), which was 2.08% lower than the observed value of 96 stb/d. When the fracture spacing was reduced to 15 ft, the production rate at the beginning of BDF should be 71% higher than that the original one, assuming no other complex factors were missing in the model development. Of course, the decrease in perf cluster spacing may have less than ideal impacts on well productivity as a result of other complex factors that were not considered in the model development.  Figure 10 demonstrates the diagnostic plot of production rate data for Well #4. Early production data showed a trend with a slope of −0.57, which was between −0.5 and −1, not indicating RLF. This could be due to the existence of natural fractures or the interference of hydraulic fractures. Figures 11 and 12 present the diagnostic plot of production rate data for Well #5 and Well #6, respectively. Early production data showed trends with slope values of −0.58 and −0.62, respectively, not indicating RLF. This could be due to the existence of natural fractures or the interference of hydraulic fractures. The bubble map in Figure 4 was compared with geological maps for natural fracture/shale quality identification. A further examination of locations of Wells #4, #5, and #6 confirmed that these wells were drilled in an area where multiple natural fractures were found. It was generally believed that not all the hydraulic fractures were equal in length, owing to heterogeneity within the reservoir or stress shadow. Ambrose et al. [33] indicated fracture interference had a great effect on the performance of a well. They demonstrated that fracture interference tended to increase the cumulative production. However, the higher the heterogeneity, the less the recovery factor. It was therefore not recommended to reduce fracture spacing for these wells or wells in the areas. Figure 9. Model-calculated productivity for TMS wells. Figure 10 demonstrates the diagnostic plot of production rate data for Well #4. Early production data showed a trend with a slope of -0.57, which was between -0.5 and -1, not indicating RLF. This could be due to the existence of natural fractures or the interference of hydraulic fractures. Figure 11 and Figure 12 present the diagnostic plot of production rate data for Well #5 and Well #6, respectively. Early production data showed trends with slope values of -0.58 and -0.62, respectively, not indicating RLF. This could be due to the existence of natural fractures or the interference of hydraulic fractures. The bubble map in Figure 4 was compared with geological maps for natural fracture/shale quality identification. A further examination of locations of Wells #4, #5, and #6 confirmed that these wells were drilled in an area where multiple natural fractures were found. It was generally believed that not all the hydraulic fractures were equal in length, owing to heterogeneity within the reservoir or stress shadow. Ambrose et al. [33] indicated fracture interference had a great effect on the performance of a well. They demonstrated that fracture interference tended to increase the cumulative production. However, the higher the heterogeneity, the less the recovery factor. It was therefore not recommended to reduce fracture spacing for these wells or wells in the areas.   Figure 10 demonstrates the diagnostic plot of production rate data for Well #4. Early production data showed a trend with a slope of -0.57, which was between -0.5 and -1, not indicating RLF. This could be due to the existence of natural fractures or the interference of hydraulic fractures. Figure 11 and Figure 12 present the diagnostic plot of production rate data for Well #5 and Well #6, respectively. Early production data showed trends with slope values of -0.58 and -0.62, respectively, not indicating RLF. This could be due to the existence of natural fractures or the interference of hydraulic fractures. The bubble map in Figure 4 was compared with geological maps for natural fracture/shale quality identification. A further examination of locations of Wells #4, #5, and #6 confirmed that these wells were drilled in an area where multiple natural fractures were found. It was generally believed that not all the hydraulic fractures were equal in length, owing to heterogeneity within the reservoir or stress shadow. Ambrose et al. [33] indicated fracture interference had a great effect on the performance of a well. They demonstrated that fracture interference tended to increase the cumulative production. However, the higher the heterogeneity, the less the recovery factor. It was therefore not recommended to reduce fracture spacing for these wells or wells in the areas.  Material balance time Actual time Figure 11. Plot of production rate data for Well #5.