Integrated optimum design of hydraulic fracturing for tight hydrocarbon-bearing reservoirs

Although hydraulic fracturing is not a new technology, it has not yet been implemented in the United Arab Emirates. Abu Dhabi National Oil Company (ADNOC), the regional producer in United Arab Emirates, has set out to initiate the utilization of this treatment during year 2019. In this work, a systematic design procedure for hydraulic fracturing in tight petroleum-bearing reservoirs is proposed. The design caters for surface and subsurface flow parameters, and it is hoped to provide basic guidelines for ADNOC in this respect. The proposed design process incorporates both unified fracture design (UFD) methodology and the fracture geometry (PKN) model. Excel spreadsheets were developed and utilized to run sensitivity analysis for optimal performance and predict long-term production profiles before and after fracturing. The excel spreadsheets made are flexible in use, in the sense that they resolve issues with infinite/finite fractures, high/low surface injection rate as well as investigate for non-Darcy flow effects. Reliable published data were used to perform the necessary calculations. The results of the performance calculations have shown that it is possible to access commercial quantities of hydrocarbons from a tight reservoir. In addition, improved productivity by 15-folds and increased gas recovery of 1.02 MMMscf over the first 8 years of production can be achieved by proper hydraulic fracturing design and implementation in tight gas reservoirs. The results of calculations of non-Darcy effect revealed a threshold velocity of approximately 0.2 fps above which these effects could become significant in predicting the overall flow efficiency inside the fracture. To the authors’ knowledge, the literature has not fully addressed the hydraulic fracturing design analytically, and the methodology proposed in this work provides a complete design package which incorporates the UFD concept, the PKN model, the non-Darcy model, and long-term prediction of post-fracturing production performance, and applying the proposed approach in a case study.


Introduction
Tight carbonate reservoirs, although less well understood and believed to require higher development costs and risks than conventional reservoirs, have become an important hydrocarbon resource. Historically, tight reservoirs have been unpopular and unfavorable for geologists and reservoir engineers mainly due to difficulty associated with their development with limited commercial productive value. Well stimulation is commonly applied to boost oil or gas production and, consequently, revenues therefrom. Hydraulic fracturing is a well-stimulation technique used commonly in low-permeability rocks like tight sandstone, shale, and some coal beds to increase oil and gas flow to a well from petroleum-bearing rock formations. Hydraulic fracturing has become a critical component in the successful development of unconventional reservoirs; including tight gas, oil and gas-producing shales, and coal bed methane, such resources rely on hydraulic fracturing for commercial viability (Al-Attar and Barkhad 2018). By 2012, hydraulic fracturing had merged as $20 billion annual activity, second only to drilling budgets (Economides et al. 2015). The execution of a hydraulic fracture involves the injection of fluids at a pressure sufficiently high to cause tensile failure of the rock. At the fracture initiation pressure, often known as the "breakdown pressure," the rock opens. As additional fluids are injected, the opening is extended, and the fracture propagates to the design fracture half length.
In recent years, new fracturing designs and techniques have been developed to enhance production of trapped hydrocarbons (Lin and Zhu 2010;Nolte and Economides 1991;Song et al. 2011). Garavand and Podgornov 2018 proposed a modified pseudo-3D model, suggesting a combination of unified fracture design (UFD), 2D fracture propagation models (PKN or KGD) and linear elastic fracture mechanic (LEFM) principle to achieve optimized fracture geometry. Based on two field studies, these authors concluded that their model is applicable to broad spectrum of multilayered oil and gas reservoirs with more accurate estimation of final fracture height and treating pressure.
A primary goal in unconventional reservoirs is to contact as much rock as possible with a fracture or a fracture network of appropriate conductivity. In designing a hydraulic fracture treatment, fracture length, width, and proppant pack permeability must be suited to the formation properties, especially formation permeability. For tight gas reservoirs specifically, which are the target formations of this work, long and narrow fractures are the most appropriate (Economides and Martin 2007). Justification of the previous statement in simple terms is as follows; in such reservoirs, the well would not produce at an economic rate without the hydraulic fracture, and thus it is not a matter of relatively improving the permeability of the matrix system, but rather creating a channel long enough to free as much trapped hydrocarbon as possible, that would otherwise have no means to flow to the well. Tight carbonate reservoirs, although less well understood and believed to require higher development costs and risks than conventional reservoirs, have become an important hydrocarbon resource. The majority of oil and gas fields in the UAE are carbonates of low matrix permeability, and hydraulic fracturing seems to be the most effective technique to develop such reservoirs.
The design of hydraulic fracturing is strongly dependent on reservoir permeability. In high-permeability reservoirs, the fracture accelerates production without impacting the well reserves, and therefore, fracturing relies on net present value economics. The characteristic curves of such reservoirs are shown in Fig. 1. In reservoirs of low or extremely low permeability (less than 1 md for oil and 0.01 md for gas), the hydraulic fracture significantly contributes to well productivity and to the well reserves, as illustrated in Fig. 2. In such reservoirs, the well would not produce an economic Fig. 1 The concept of production and recovery acceleration, high-permeability reservoirs (Holditch 2006) 1 3 rate without induced fracturing. Comparing the two cases yields the following conclusion: Whether our reservoir is of moderate permeability or low permeability, there will be an increase in production rate following the treatment; as can be observed from Figs. 1 and 2 of production versus time. The difference between the two scenarios is the following; in reservoirs of moderate permeability (up to 50 md for oil and 1 md for gas), the amount of recoverable hydrocarbons remains unaffected. However, in low-permeability reservoirs, ultimate recovery will be significantly increased post-fracturing (Economides et al. 2015). Depending on the original reservoir permeability, the objective of our design would vary.
In the past, hydraulic fracturing was limited to lowpermeability reservoirs only, but nowadays fracturing has become an integral part of well and reservoir management rather than a decision. Every hydraulic fracture is characterized by its fracture half conductive length (x f ), dimensionless fracture conductivity (C fD ), and equivalent skin factor (S f ). Agarawal et al. (1979) and Cinco-Ley and Samaniego (1981) introduced the fracture dimensionless conductivity as A graphical presentation of the relationship between C fD and S f was presented by Cinco-Ley and Samaniego (1981) that can be used to estimate the folds of increase in productivity post-fracturing treatment. The concept of equivalent skin effect, however, only applies after pseudoradial flow is established. When pseudoradial flow emerges after bilinear and/or linear flow, the skin determined from a pressure buildup test will be the equivalent skin effect indicated in the graphical presentation of Cinco- Ley and Samaniego (1981). The same authors plotted [S f + ln (x f /r w ) + 0.5 ln (C fD )] versus C fD on a semi-log graph paper and their plot shows a minimum that appears at exactly C fD = 1.6. They concluded that there exists one optimum conductivity that maximizes (1) C fD = k f w∕kx f the productivity index, and this conductivity corresponds to a particular fracture width and length. It should be emphasized at this point that such hypothesis applies whenever x f ˂˂ r e that pseudoradial flow can occur. A methodology provided by Economides, Oligney and Valkó (2002), which assumes that for a given proppant volume, well drainage area and shape there is a fracture half length, width and conductivity that maximizes well productivity. They provided the concept of unified fracture design (UFD) which expands the above optimum productivity approach to include fracture and well drainage dimensions that will not reach pseudoradial flow before the onset of pseudosteady state. These authors developed correlations for optimum fracture dimensions in square drainage areas which will yield maximum dimensionless pseudosteady-state productivity index and as presented in Eqs. 2 through 7. (2) for Np ≥ 100 (3)  Fig. 2 The concept of production and reserves enhancement from hydraulic fracture, low-permeability reservoirs (Holditch 2006) where Np is dimensionless proppant number defined by Eq. 4. For non-square drainage areas, productivity optimization was presented by Daal and Economides (2006). They presented a set of equations similar to Eqs. 2 through 7 but modified to accommodate the drainage shapes other than the square drainage area. Their approach is appropriate for multiple hydraulic fractures intersecting a horizontal well transversely. Perkins and Kern (1961) developed a model for predicting fracture width for a fixed height vertical fracture. This model was later adjusted for laekoff and storage within fracture by Nordgren (1972) and is now known as the PKN model. The average width of the PKN fracture is expressed in Eq. 8.
The PKN model is applicable for cases in which the fracture length is much larger than the fracture height (x f ≫ h f ) and fracturing fluid is Newtonian. This is ideal for the target formation which has an extremely low permeability, thus requiring long and narrow fractures that would create sufficient conductivity to maximize performance and justify cost. In Eq. (8), the term (πγ/4) is approximately equal to 0.59 and γ is a factor equal to 0.75. For non-Newtonian fracturing fluids, Eq. 9 can be used to estimate maximum fracture width.
The quantities n′ and K′ are the power-law rheological properties of the fracturing fluid. The average fracture width can be calculated by multiplying w max by 0.628 (Economides et al. 2015).
Other important design parameters include total injection time (t i ) and volume of treatment (V i ). The approximate relationship between the volume of treatment, and the pad volume, V pad , was provided by Nolte (1986) and Meng and Brown (1987), as expressed in Eq. 10.
where η is the fluid efficiency and is determined by iterative scheme of calculations as follows. A material balance between total fluid injected, created fracture volume, V f , and fluid leakoff, V L , can be expressed as where q i is the injection rate in bpm, t i is the total injection time in minutes, A f is the fracture area in sq.ft, C L is the leakoff coefficient in ft/min 0.5 , and r p is the ratio of the net to fracture height (h/h f ). The variable K L is known as opening time distribution factor and defined by Nolte (1986).
Thus, the value of K L is between 1.33 and 1.57. For a given x f , the average hydraulic fracture width, w bar , can be calculated using the PKN model. Equation (12), which is a quadratic equation, can be solved for t i by trial and error. The product q i t i is equal to the volume of treatment (pad volume plus slurry volume), V i , and The leakoff coefficient, C L , can be obtained from the laboratory or from a fracture calibration treatment as described by Nolte and Economides (1989). During the slurry stage of the fracturing treatment, propping agent is continuously added to the fracturing fluid at time-dependent concentrations until the total injection time is achieved. Therefore, the fracturing design should solve for proppant schedule which depends on the fluid efficiency. Nolte (1986) has shown that where c p (t) is the slurry concentration for any time greater than t pad , in pounds per gallon in ppg, c f is the final slurry concentration in ppg, and the variable ε depends on the fluid efficiency defined in Eq. (18). At this point, it should be emphasized that the average fracture width predicted by the PKN model is the width created during the propagation of the fracture, and is known as hydraulic width. At the end of the fracturing job and when the surface pumps are switched off at the surface, the fracture tends to close due to the natural stresses in earth. The placement of propping agent, however, inside the fracture will only allow partial closing of the fracture which results in a final fracture width, known as propped fracture width denoted by w p , that is a little narrower than the hydraulic fracture width. Assuming that a mass of proppant, M p , has been injected into a fracture (12) of half length x f and height h f and the proppant is uniformly distributed, then where the product 2x f h f w p (1 − ϕ p ) represents the volume of the proppant pack and ρ p is the density of the propping agent in lbm/cu.ft. The proppant concentration in the fracture, C p in lb/sq.ft., may be defined as In the present work, the proposed design targets lowpermeability reservoirs since they pose a challenge to the UAE oil and gas industry. The design combines the concept of UFD, total injection time, volume of treatment, proppant schedule, propped fracture width, non-Darcy flow, and longterm production behavior for a hydraulically fractured in one package.

Design criteria and specifications
The following conditions represent the foundation of the proposed design: Computation of surface injection pressure depends on the value of injection rate. (whether it is greater than 9 bpm) 10. Maximum allowed injection time for efficiency purposes: 24 h. 11. Maximum allowed surface injection pressure for safety purposes: 5000 psi. 12. Maximum possible skin effect due to stimulation: − 8.5 (Flow efficiency goes to infinity at this value). 13. The fracturing fluid efficiency should be greater than 0.5. 14. The following relates formation permeability and fracture width: • Average fracture width is inversely proportional to the shear modulus of elasticity. • Shear modulus of elasticity is directly proportional to Young's modulus of elasticity. • Young's modulus of elasticity is dependent on depth.
• As the depth increases, the rock is more consolidated thus matrix permeability decreases.

Assumptions
The following assumptions were made to simplify the construction of the design and cater for the conditions of the target formation (carbonate tight gas reservoirs): 1. Assumed fracture height equals reservoir height. Height containment is due to the natural stress contrast resulting from the difference in Poisson ratio between the investigated tight carbonate formation and surrounding shale formations. 2. The following range is assumed for Young's modulus of elasticity: (a) 2E05 ˂ E ˂ E07 psi, (b) E ≈ E07 psi in deep reservoirs (narrow fracture), and (c) E ≈ 2E05 psi in shallow reservoirs (wide fracture). 3. For infinite fracture conductivity, pressure drop in the fracture is assumed negligible. 4. Assumed very low value of liquid leak-off coefficient to improve design efficiency. [C L = 0.0001 ft/(min) 1/2 ] 5. Formula of skin effect used is based on radial flow assumption. This parameter, however, is best obtained from post-fracture test analysis. 6. Skin effect before stimulation assumed zero. Hydraulic fracture bypasses the near-wellbore damage zone. Thus, the pretreatment skin effect has little or no impact on the post-fracture equivalent skin effect value. 7. Constant-percentage decline rate. 8. Vertical fractures are induced at depths greater than 3000 ft.
9. The following flow pattern is assumed for each of the mentioned cases (Figs. 3, 4): where ETR: early time region (reflecting wellbore storage and skin effects), MTR: middle time region (infinite-acting system behavior reflecting reservoir properties), LTR: late time region (reflecting boundary effects). Both figures are top views showing fracture extension into the formation. In Fig. 3, the fracture is of infinite conductivity and half-fracture length is equal to the drainage radius of the formation and thus, the pressure transient moves from ETR into LTR without going through MTR. The flow geometry is basically linear from the formation into the fracture. In Fig. 4, the fracture is of finite conductivity and half-fracture length is less than the drainage radius of the formation and thus, the pressure response exhibits ETR then MTR and finally LTR. The flow geometry is bilinear between the wellbore and fracture tip and pseudoradial between fracture tip and formation outer boundary.

Detailed proposed design procedures
ln r e ∕r w −0.75−ln 0.28k f w∕k r w 1 + N Re,nD 15. Based on the results, choose the best M p value. 16. Repeat process, vary k f for chosen M p and fixed k. 17. Conclude the most productive and cost-effective M p and k f values for the design. 18. Predict long-term production performance as presented in the following section.
A flowchart of the proposed design methodology is presented in Fig. 5. In this chart, Chen equation provides a direct solution for the friction factor in terms of Reynolds number and pipe relative roughness. Equation 25 is used to estimate pressure drop due to friction across pipe length. It has been found that this equation loses its accuracy when the fluid injection rate exceeds 9 bpm (Economides et al. 2015).
Applying Eqs. (27 through 32), a series of calculations were performed to investigate the effect of non-Darcy flow in hydraulically induced fractures. The results of these calculations are presented in Figs. 6 and 7, in which the flow velocity in fps is plotted versus pressure drop in psi with correlating parameters x f and k f , respectively. It can be observed that in both plots and regardless of the value of correlating parameter, the relationship is nonlinear and that the dramatic increase in flow velocity takes place at pressure drop approximately equals to 3700 psi. The inflection point in both plots designates threshold conditions in terms of flow velocity and corresponding pressure drop. The threshold velocity is estimated by drawing an asymptote through the higher velocity data points and it begins at approximately 0.2 fps. At velocities higher than 0.2 fps, the non-Darcy flow effects become significant and ignoring these would result in significant errors in the calculations of fracture conductivity.

Predicting long-term production behavior
For an infinitely conductive fracture, with the fractured well in the middle of a rectangular bounded reservoir, Wattenbarger et al. (1998) provided a general formula for the dimensionless production rate: where y e is the distance from the fracture to the drainage boundary in the direction perpendicular to the fracture, and the dimensionless time is Equation (35)

Detailed design calculations
The design methodology described earlier and illustrated in Fig. 5, was implemented in a vertical gas well completed in a carbonate reservoir with low permeability. The data pertinent to the well and the reservoir are listed in Table 1. A complete sample of calculations is presented in "Appendix A" including schedule of proppant as illustrated in Fig. 8. Three scenarios were considered and as listed in Table 2. The last five columns in Table 2 highlights the optimum values of design parameters for each scenario which meet all design criteria. The relationship between optimum design parameters and formation permeability can be deduced from Table 2. A sample of calculations for predicting long-term production profile before and after fracturing is shown in "Appendix B". Since the fracture half length (x f ) and half of the length (0.5 x e ) of square drainage area are approximately equal, therefore the flow geometry is basically linear (for C fD > 50). The time this flow geometry lasts is calculated by Eq. (41), and the production profile is predicted by Eq. (42) (Kelkar 2008).
Q gsc,L = m p i −m p wf (h∕T) x f opt k c t 0.5 (t) −0.5 ∕64.35 The results of calculations of production profile over this period when linear flow geometry dominates are listed in Table 3. For the same reason, there will be no middle time region and pressure transient enters directly into pseudosteady state without going through infiniteacting system behavior. Predicting long-term production profile for times greater than t lf end is achieved by Eq. (43) (Kelkar 2008) and the results of calculations are illustrated in Table 4.

Results and discussion
The proposed design methodology was implemented in three scenarios representing three reservoir permeabilities of 0.01, 0.1, and 1 md, respectively. The data used in all calculations are listed in Table 1. For each scenario, the proppant mass (M p ) was varied between 50,000 lbm and 500,000 lbm, and for each proppant mass, 11 key design parameters were calculated. These parameters include: maximum dimensionless pseudosteady-state productivity index, optimum fracture half length, optimum propped fracture width, optimum fracture conductivity, equivalent skin factor, fracturing fluid injection rate, total injection time, maximum surface injection pressure, fluid efficiency, gas production rate, and dimensionless productivity ratio with and without non-Darcy flow effects. The results of one scenario are presented in Table 5 and it can be observed that the best choice would be for M p = 200,000 lbm because corresponding surface injection pressure, total injection time, fracturing fluid efficiency, and fluid injection rate meet the design criteria set out earlier in this paper. Qualitatively, this value of M p is cost effective in comparison with the higher values in Table 5 (figures highlighted in red color do not meet design criteria). Similarly, the best choices of M p for the other two scenarios were concluded and the results are listed in Table 2. The last column in Table 2 illustrates the ratio J Dmax pss post-frac /J D pss pre-frac for the three selected scenarios and it can be observed that this ratio is highly dependent on the reservoir permeability. A maximum value of 15 for this ratio has been obtained for k = 0.01 md, which is indicative of the effectiveness of hydraulic fracturing treatment in tight reservoirs. The results of the three selected scenarios also show that hydraulic fracturing is less effective in higher permeability reservoirs. The effect of non-Darcy flow was included in the proposed design methodology using Eqs. (26 through 33). To confirm whether non-Darcy effect is important, sensitivity analysis was performed in terms of flow velocity versus pressure drop inside the fracture with fracture half length and fracture permeability as correlating parameters as shown in Figs. 6 and 7, respectively. The resulting relationships shown in both figures exhibit a deflection from linear to nonlinear behavior  which takes place at flow velocity of nearly 0.2 fps. In this paper this velocity is called threshold velocity above which neglecting non-Darcy effects in the calculations of gas flow rate can result in serious errors. In addition, non-Darcy effects can be evaluated by either Eq. (32) or its equivalent Eq. (33), through the calculation of the ratio J nD /J o . A ratio of one indicates negligible non-Darcy effects, while a ratio less than one is indicative of some degree of non-Darcy flow effects. These two criteria, namely, threshold velocity and J nD /J o , have been considered in the proposed design methodology. The last part of the discussion addresses the prediction of long-term production profiles. The results of calculations of this part are presented in Tables 3 and  4 for early time linear flow geometry-transient state and late-time radial flow geometry pseudosteady state, respectively. Detailed calculations can be found in "Appendix B". Low-permeability reservoirs are commonly produced with p wf held constant, so that the well boundary condition for the well production models illustrated in Eqs. (34 through 42) apply (Economides et al. 2015). Consequently, in this work, the bottom-hole flowing pressure was held constant at 2500 psi. The pre-fracturing predictions of production profile over a period of 8 years have shown a significant drop in gas production rate from 92.57 to 8.45 Mscf/d (Fig. 9) and cumulative gas recovery of 0.23 MMMscf at the end of the eighth year. Since there has been no production history available for this case, a constant nominal daily decline rate D d = 0.000417 day −1 was assumed to perform these predictions. The results of post-fracturing predictions for the same period of 8 years indicated a significant drop of gas flow rate from 1500 to 200 Mscf/d (Fig. 10) and a cumulative gas recovery of 1.25 MMMscf. These results are in line with the qualitative field observations shown in Figs. 1 and 2 (Holditch 2006), and that hydraulic fracturing increased the productivity of this particular case 15-folds and increased the cumulative gas recovery by 1.02 MMMscf.

Summary and conclusions
The design methodology proposed in this work encompasses the UFD concept and PKN geometry model. Spread sheets have been prepared to accommodate the necessary design calculations using published data. Three scenarios of reservoir permeability have been selected to investigate the effectiveness of hydraulic fracturing in tight formations. Based on the results obtained in this work, the following conclusions may be drawn.
1. A simple design methodology is presented which makes it possible to compute key parameters in hydraulic fracturing treatment with the limited tools and data. 2. Application of the proposed design package to three scenarios of tight gas reservoirs have shown that hydraulic fracturing is most effective when k = 0.01 md and corresponding optimum M p and maximum ratio J Dmax pss post-frac /J Dpss pre-frac have been found equal to 200,000 lbm and 15, respectively. Certainty in enhancement of gas production justifies initial investment considering the folds of increase in production. 3. Sensitivity analysis of non-Darcy flow effects revealed a threshold flow velocity of 0.2 fps. For any flow velocity below this limit, the non-Darcy flow effects can be neglected from the calculations of gas flow rate and J nD /J o = 1.0. 4. Predictions of production profiles have shown significant improvement in gas production rate and gas recovery over the 8 year time period at 15 folds and 1.02 MMMscf, respectively.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Appendix A
See Table 5

Appendix B
Predicting long-term pre-fracturing production profile See Fig. 9. For p r = 4500 psi, p wf = 2500 psi, and S = 0.