An integrated model for multiobjective design optimization of hydraulic fracturing

Journal of Petroleum Science and Engineering 31 Ž2001. 41–62 www.elsevier.comrlocaterjpetscieng

An integrated model for multiobjective design optimization of hydraulic fracturing M.M. Rahman, M.K. Rahman ) , S.S. Rahman School of Petroleum Engineering, The UniÕersity of New South Wales, Sydney 2052, Australia Accepted 28 June 2001

Abstract An integrated novel model for hydraulic fracturing design optimization is presented, which recognizes complex interactions between a hydraulically coupled fracture geometry module, a hydrocarbon production module and an investment-return cash flow module. Free design variables are identified and various design constraints are formulated, which must be satisfied so that an optimum design obtained is executable in the field using the specified surface equipment Žpump, tubing, etc.. and that the treatment does not cause any undesirable formation damage by uncontrolled fracture growth andror multiple secondary fracture initiation. The model is formulated within the framework of a multivariate and multiobjective optimization method, which is based on the combined features of Genetic Algorithm and Evolutionary Operation. A 2D fracture model is used to establish relationships between treatment parameters and fracture growth. The potential for hydraulic fracturing design improvement is demonstrated by application to a tight gas reservoir. Results show that the proposed model is instrumental in improving hydraulic fracturing design and achieving a goal-oriented optimum design in a conflicting environment. About 12% compromise with maximum possible production, or net present value ŽNPV., over 10 years can save up to 44% of initial hydraulic fracturing treatment cost. Furthermore, the 88% production, or net present value, as a result of an optimum treatment program is significantly higher than any arbitrary design. The issues of real-time design modification, however, are not included in this study. q 2001 Elsevier Science B.V. All rights reserved. Keywords: Reservoir stimulation; Hydraulic fracturing; Treatment design; Optimization

1. Introduction The petroleum industry has long been applying hydraulic fracturing treatment as a technique to stimulate oil and gas reservoirs. For a given reservoir, one must decide a set of values for various treatment ) Corresponding author. Tel.: q61-2-9385-6746; fax: q61-29385-5182. E-mail address: [email protected] ŽM.K. Rahman..

parameters, such as viscosity of fracturing fluid, injection rate of the fluid, injection time and proppant concentration, etc., such that a favourable hydraulic fracture geometry is created to improve a design objective. The net present value ŽNPV. is usually used as an objective for hydraulic fracturing treatment design ŽMeng and Brown, 1987; Balen et al., 1988; Economides et al., 1994; Aggour and Economides, 1998.; in these, the improvement in NPV is attempted by parametric sensitivity analysis.

0920-4105r01r$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 0 - 4 1 0 5 Ž 0 1 . 0 0 1 4 0 - 1

42

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

The height of the fracture is assumed to be a known value Žoften equal to the pay zone height. and a set of discrete values are taken for each of the treatment parameters. For each combination of these treatment parameters, the NPV is calculated for a given fracture length. The procedure is repeated for a number of fracture lengths and NPV versus length curves are plotted. From these curves, values of stimulation parameters and the resulting fracture length are ascertained as optimum based on a maximum value of NPV. This type of design optimization procedure by parametric sensitivity analysis is tedious. Although the procedure can find a design better than an arbitrary design, it is almost impossible to explore all potential design scenarios and to satisfy various operational limitations Žpump capacity, tubular strength, pressure rating of surface equipment, etc.. and fracture growth control requirements based on formation characteristics. Furthermore, the authors believe that there are other design objectives, such as maximize production, minimize treatment cost, achieve a production target, etc., which should be reflected in the optimum design. Thus, the overall hydraulic fracturing design optimization task may be redefined to find a ‘best possible’ set of values for treatment and other parameters while the revenuerproduction from a well is maximized with minimum treatment cost, and operational limitations and fracture growth control requirements are satisfied. The requirements discussed above are within the framework of various industrial design optimization problems which are usually modeled in terms of a number of design variables, design constraints and design objectives, and are often solved using various optimization methods. In oil and gas engineering, a limited number of such optimization studies have been reported in the area of production system design ŽCarroll and Horne, 1992; Fujii and Horne, 1995. and in development planning and management of petroleum reservoirs ŽMcFarland et al., 1984.. Surprisingly, the use of any formal mathematical optimization method for hydraulic fracture treatment design is rare. Yang et al. Ž1996. made an attempt for procedural optimization of hydraulic fracturing considering only injection rate, injection time and proppant concentration as free variables. Similar to other works, Yang’s work also did not consider any operational limitations and fracture growth control

requirements. The three variables were optimized within specified ranges to maximize NPV systematically varying the fracture length with a fixed height. Thus, the work did not exploit the full benefit of optimization procedure; rather it became another improved version of design improvement by parametric sensitivity analysis. Mohaghegh et al. Ž1996, 1999. recently reported a hybrid Neuro-genetic approach for hydraulic fracture treatment design optimisation processing historical data. Mohaghegh et al. thus avoided the mathematical formulations of issues addressed here, however, correctly acknowledging that their methodology is not a substitution for physicsbased approaches, which leaves room for the kind of work presented in this paper. Mahrer Ž1999. reviewed studies of hydraulic fracture geometry and treatment practice and he concluded the review with Aoptimizing treatment procedures based on in-situ conditions is an area of technology that requires research and developmentB. An optimization model was developed in this work for hydraulic fracturing design integrating various treatment and fracture parameters, operational limitations, fracture growth control requirements, potential design objectives and a suitable optimization algorithm. The model thus allows the hydraulic fracture itself to assume a suitable geometry to improve the specified design objectiveŽs. to its fullest extent. Time and place-dependent reservoir properties and price data can be entered into the model in order to obtain an automated optimum design using selected fracturing fluid and proppant type. The principal objectives of this paper are to present the integrated model and to demonstrate its capability to investigate various design issues for improved decision-making by a series of applications to an example gas reservoir.

2. The optimization algorithm The central problem of mathematical programming is usually stated as the optimization Žin the mathematical sense of maximization or minimization. of a function f Ž x 1 , x 2 , . . . , x N . of several variables x 1 , x 2 , . . . , x N subject to a set of con-

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

straints. There exists a voluminous literature on optimization methods and their applications to linear and nonlinear problems Že.g. Pierre, 1969; Wismer, 1971; Jacoby et al., 1972; Wilson et al., 1981; Bunday, 1984; Polyak, 1987.. To solve linear optimization problems, the simplex algorithm is a powerful tool and its applications are widely documented. Most industrial design problems, however, involve nonlinear objective and constraint functions. Manipulating basic mathematics of optimization methods, various algorithms are developed, including sequential linear programming ŽSLP. ŽMistree et al., 1981., sequential quadratic programming ŽSQP. ŽSpellucci, 1998., sequential unconstraint minimization technique ŽSUMT., normal-boundary intersection algorithm ŽNBI. ŽDas and Dennis, 1998., genetic algorithm ŽGA. ŽGoldberg, 1989. and polytope algorithm ŽPA. ŽNelder and Mead, 1965.. Based on the method of root finding, most of the methods are grouped into two broad categories: Ž1. derivative based methods and Ž2. direct search methods. For problems with differentiable smooth functions, derivative based methods are reliable and computationally efficient. Attempts have also been made to apply these methods to nondifferentiable problems using various numerical differentiation techniques. Direct search methods, such as GA and PA, are generally slow in convergence but are successful in finding reliable optimum solutions to problems with a high degree of variable noise, including discontinuity and nondifferentiability in functions. The objective and constraint functions involved with hydraulic fracturing optimization are highly nonlinear and nondifferentiable. These functions are also subjected to a certain number of discontinuities. These include the discontinuous production estimation function at the transition of transient and pseudo-steady-state flow regimes and numerically erroneous sub-functions, such as the complementary error function, bi-section solution, etc. Initial attempts were made to solve the problem using a number of derivative based algorithms including SLP ŽRahman, 1998., SQP ŽSpellucci, 1998. and MATLAB routines with their in-built numerical differentiation features. With exhaustive efforts to achieve noise reduction, it was finally possible to obtain optimum solutions for maximum NPV. With other design objectives, however, success was very limited

43

and unreliable. Finally, a direct search based algorithm called an ‘INTElligent Moving OBject’ ŽINTEMOB., algorithm was developed. The underlying principle of INTEMOB is developed based on the combined concepts of PA, GA and Evolutionary Operation ŽEVOP., ŽGhani, 1989.. General formulations of INTEMOB are: Minimize f Ž x . where, xT s Ž x 1 , . . . , x N .

Ž 1.

subject to bound constraints l1 . . . lN

F

F

x1 . . . xN

F

F

¶ • ß

u1 . . . uN

Ž 2.

and design constraints Cl1 . . . Cl M

F

F

C1 Ž x . . . . CM Ž x .

F

F

¶

Cu1 . . . Cu M

•

ß

Ž 3.

where, x represents the vector of free design variables Žthe superscript ‘T’ for transpose.; l i ’s and u i ’s are constants or functions of x Žin the latter case the bound constraints constitute moving boundaries. representing the ranges of x i ’s in which their corresponding values are searched for; C li ’s and Cui’s are constants or functions of x representing the ranges of design constraints, Ci Ž x.’s, which must not be violated by the optimum design; N is the total number of free design variables and M is the total number of design constraints. Maximization of a function, f Ž x., can be achieved by minimizing yf Ž x .. The optimization procedure starts by generating ‘K ’ number of random vertices in the N-dimensional space bounded by the ranges of design variables, where K s 2 N for N F 5 and K s N q 1 for N ) 5. The coordinates of a random vertex are generated by: x i s l i q ri Ž u i y l i . ;

i s 1, . . . , N

Ž 4.

where ri is a pseudo-random deviate rectangularly distributed over the interval Ž0, 1. which is controlled by a known value, x in , for the i-th coordinate. Therefore, the procedure starts with a given initial point Žvertex. and Eq. Ž4. itself ensures that any

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

44

randomly generated new point remains in the space bounded by the ranges of variables defined by Eq. Ž2.. Executing Eq. Ž4. K times, K different random points are generated. Any of the generated points may initially violate any of the constraints defined by Eq. Ž3.; therefore, a technique is required to move such points towards satisfying Eq. Ž3.. For simplicity in description, a two-dimensional space Ž N s 2. is considered in which four random vertices, a, bX , cX and dX are generated in Fig. 1. Obviously, vertices bX , cX and dX violate Eq. Ž3.. Vertices are modified by moving successively towards the centroid, c, by: xs

1 2

Ž cqxX .

Ž 5.

until the new point, x , satisfies Eq. Ž3.. The coordinates of the centroid, c, are calculated using vertices that have already satisfied Eq. Ž3. as follows: ci s

1 N

N

Ý

xXi

values of the objective function, f Ža., f Žb., f Žc. and f Žd. are calculated and the point that corresponds to the minimum objective function value is preserved and new random points are generated around this point. The procedure is repeated as long as the objective function can be minimized further, i.e. a pre-specified convergence criterion is not met. In a sense, the same object Žpoint. is moving to a new location in every step, using intelligent information until the minimum objective function value Žthe optimum point. is established; hence, the name intelligent moving object algorithm, INTEMOB. There are other fine-tuning features in the procedure.

3. The integrated design optimization model

Ž 6.

is1

where n is the number of vertices which have already satisfied Eq. Ž3.. The modified points are a, b, c and d which satisfy both Eqs. Ž2. and Ž3.. Using these points, the

Essential modules and inter-modular interactions in the integrated model are shown in Fig. 2. In-situ reservoir properties and operational limitations are derived outside the fracturing design optimization process. The optimization process starts with initial

Fig. 1. Generation and modification of random vertices in a two-dimensional space.

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

45

Fig. 2. Integrated model of hydraulic fracturing design optimization.

values of design variables with specified upper and lower bounds on each variable. The optimization algorithm then performs bi-directional interactions with the current fracture geometry and design variables to improve the design objective and satisfy the design constraints formulated based on fracture growth control and operational limitations. The achievement of the best solution requires iterative interactions between these modules. 3.1. Design Õariables Fracturing parameters, which have been independently varied in this study to find an optimum design, include fracturing fluid viscosity, m ; injection rate of fracturing fluid, qi ; injection time, t i ; proppant concentration, Pc and fracture half-length, x f . Fracture height h f , and width w, are calculated solving coupling relationships developed based on fracture geometry and material balance. The first four design variables are controllable at the surface, and the fracture half-length is considered as the fifth variable to allow the fracture geometry to develop as required for improvement in the objective function.

To develop larger fracture width and cause effective proppant transport, a high viscous fracturing fluid is generally required. An excessively high viscous fluid, however, causes undesired vertical height growth ŽMeng and Brown, 1987.. In addition, a higher viscous fluid incurs greater treatment cost because of a price-viscosity relationship. A higher injection-rate generally yields a larger fracture with greater width but may not yield the most efficient treatment. This is because the higher injection rate results in higher treating pressure and surface pressure. The high treating pressure may exceed the formation critical pressure and thus induce uneconomic fracture growth; whereas, the high surface pressure may damage the surface equipment. Depending on the optimum injection rate, there exists an optimum injection time which is required to develop the optimum fracture size. Any prolonged injection after the optimum time will induce unnecessary fracture growth and thus incur additional treatment cost. Proppant concentration determines the final propped fracture penetration and conductivity. Frac-

46

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

ture permeability is a function of the proppant sizer type, its concentration, and the residual damage to the fracture from the treatment fluids being used. An optimum value of proppant concentration at the end of the job is sought for a given proppant typersize. The proppant concentration, and its scheduling over the treatment period, is related to pad volume and slurry volume using the method proposed by Nolte Ž1986., Appendix A.

3.2. Fracture geometry The fracture geometry is defined by length, height and width of a fracture, and is assumed to develop as a function of three stimulation parameters: Ž1. fracturing fluid viscosity, Ž2. injection rate and Ž3. injection time. Because an optimum value of the fracture half-length is determined by the optimization algorithm as one of the free design variables, the fracture height and width are calculated by solving fracture geometry and material balance relationships. A fracture grows with a complex configuration when it initiates in the nonpreferred direction of a wellbore due to the presence of natural fracturesr flaws, induced perforation orrand highly deviatedr horizontal wells. Usually, a fully fluid flow coupled 3D model is required to predict growth of such fractures in a 3D space. These models are developed based on the fundamental theories of linear elastic fracture mechanics which are coupled with the effects of complex fluid flow patterns inside fractures ŽHossain, 2001.. Fracture growth is simulated sequentially using a mixed-mode fracture propagation criterion Žin terms of the local stress-displacement field around the crack tip. by highly capable finiter boundary element methods. Such a model is computationally intensive and is not suitable for incorporation in a design optimization scheme that involves a large number of iterations. To idealize fracture growth in multi-layered formations, pseudo-three-dimensional ŽP-3D. models are proposed. These P-3D models are called ‘pseudo’, because they do not consider the variation of fracture geometry in a three-dimensional space, rather they modify the 2D models by adding height variation along the fracture length and its effect on the fracture

width. Settari and Cleary Ž1986. introduced the concept of modeling hydraulic fracturing by P-3D model. Due to the lack of simplified closed-form equations, the method is not quite suitable for an automated optimization process, such as the one proposed herein, although the work has been incorporated in many hydraulic fracturing simulators ŽBouteca, 1988; Morales and Abu-Sayed, 1989.. A 2D fracture model, PKN-C, is used herein; Appendix A. The widely used PKN model is shown in Fig. 3. This was developed from the pioneering works of Perkins and Kern Ž1961. and Nordgren Ž1972., and later the fluid leakoff was considered by incorporating Carter Equation II ŽHoward and Fast, 1957. for material balance at a constant injection rate, and hence, the model name PKN-C. The model can be used for both Newtonian and non-Newtonian fracturing fluids. The latter fluid is considered in this study. The PKN-C model was selected for the current work, because its vertical plane strain assumption is physically more acceptable for the proposed height constrained fractures where the fracture length, most likely, becomes considerably greater than the fracture height ŽValko and Economides, 1995.. Rahim and Holditch Ž1995. also reported that for most

Fig. 3. PKN fracture geometry model.

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

problems the PKN-C model predicts fracture lengths closer to those computed by 3D models for correct fracture heights than does the KGD-C model ŽKhristianovitch and Zheltov, 1955; Geertsma and De Klerk, 1969.. 4. Problem formulation for hydraulic fracturing design optimization The overall problem of hydraulic fracturing design optimization is formulated within the framework of the optimization algorithm presented in Section 2 as follows: T Minimize ZŽ x. where, x s Ž m ,qi ,t i , Pc , x f . subject to the following bound and design constraints. 4.1. Design objectiÕe The minimization function, ZŽ x. is defined for multiple design objectives using a combined objective function ŽRahman and Caldwell, 1992; Rahman, 1996. as follows: no < f x y T < iŽ . i minimize Z Ž x . s Ý Pi Ž 7. Di is1 where, no is the total number of objectives to be considered in the design, f i Ž x. is the i-th objective, Ti is a target value for the i-th objective, Pi is the priority factor to achieve the i-th objective and Di is a normalizing factor for achieving the i-th objective. If a certain objective has no target value, Ti is set to zero and the f i Ž x. is minimized. To maximize the i-th objective, yf i Ž x. is minimized. Using the combined objective function, the following four design objective functions are formulated for a hydraulically fractured reservoir: Ž1. Maximize total production, Gp , over 10 years: minimize Z Ž x . s minimize Ž yGp .

Ž 8.

Ž2. Maximize net present value, NPV, over 10 years: minimize Z Ž x . s minimize Ž yNPV. Ž 9. Ž3. Maximize NPV and minimize treatment cost, Ctr : yNPV Ctr minimize Z Ž x . s minimize P1 q P D1 D2 2

Ž 10 .

47

Values of D 1 and D 2 are adjusted such that the value of both terms at the right hand side of Eq. Ž10. approaches to 0.5. P1 and P2 are set to 1 to assign equal priority to maximize NPV and minimize treatment cost. Ž4. Achieve production target, T1 , and minimize treatment cost, Ctr : minimize Z Ž x . s minimize

< Gp y T1 < D1

P1 q

Ctr D2

P2

Ž 11 . Values of D 1 , D 2 , P1 and P2 are set, similar to Eq. Ž10., to reflect equal priorities for achieving the production target and minimizing the treatment cost. The cumulative production, Gp from hydraulically fractured wells in tight-gas reservoirs is estimated combining transient and pseudo-steady-state flow regimes, Appendix B. The NPV is formulated as follows: NY

NPV s

Rn

Ý n y Ctr ns1 Ž 1 q i .

Ž 12 .

where R n is the revenue generated at year n, NY is the total number of years during which revenue to be generated, and i is the discount rate. The revenue, R n in the year n is calculated as the product of the total production, Gp n at the n-th year and an average gas price. The treatment cost, Ctr is estimated as follows: Ctr s Pf l = Vf l q Ppr = Wpr q Ppump = HPav q FC

Ž 13 . where Ctr is the treatment cost Ž$., Pfl is the price of fracturing fluid Ž$rm3 ., Vfl is the volume of fracturing fluid Žm3 ., Ppr is the price of proppant Ž$rkg., Wpr is the weight of proppant Žkg., Ppump is the pumping price Ž$rhp., HPav is the hydraulic power of the pump Žhp., and FC is the fixed and miscellaneous costs. The fracturing fluid price is related to viscosity based on industry data, to reflect cost dependency on the fluid type, as follows: Pfl s 271.81 m q 177.56

Ž 14 .

48

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

4.2. Bound constraints The following bound constraints are formulated with the problem specific data and other considerations: 1. 0.070 F m F 0.800: based on industry data and recommendations for minimum viscosity by Economides et al. Ž1994. and Smith Ž1992. for effective proppant transport. 2. 0.026 F qi F 0.12: from field experience. 3. 600 F t i F 30 000: the upper bound 30 000 s was chosen so that the optimum value of t i is never constrained by injection time. 4. 600 F Pc F 1800: based on industry practice and Meng and Brown Ž1987.. 5. 38 F x f F 792: the upper bound was chosen such that the fracture does not exceed the reservoir boundary.

4.3. Design constraints Design constraints are formulated based on operational limitations and fracture growth control requirements as follows. 4.3.1. Constraints to formulate operational limitations Ž1. 1.0 F C1Ž x. F 10.0: where C1Ž x. s HPav = PeffrHPreqd ; HPav is the horsepower available from the pump to be used and Peff is the pump efficiency factor. The lower bound on this constraint ensures that the horsepower required, HPreqd , to deliver fracture net pressure, Pnet , is within the capacity of the pump and the upper bound is set arbitrarily. Formulations to estimate HPreqd are presented in Appendix C. Ž2. 1.0 F C2 Ž x. F 15: where C2 Ž x. s P burstrŽ Psurf = SF.; P burst is the burst strength of the tube in use, Psurf is the pressure developed at the surface and SF is a safety factor. The lower bound on this constraint ensures that the pressure developed inside the tube at the surface level is below the burst strength of the tube during injection. Psurf is calculated in Appendix C based on the required net fracture pressure ac-

counting for static head and dynamic frictional loss using fluid consistency indexes calculated based on industry data. Ž3. 1.0 F C3 Ž x. F 15: where C3 Ž x. s PseqprPsurf ; Pseqp is the minimum pressure among the rated pressures of various surface equipment in the injection line. The lower bound on this constraint ensures that the pressure developed at the surface does not exceed the pressure capacity of the critical equipment.

4.3.2. Constraints to formulate fracture growth control requirements Fracture growth control requirements are formulated as design constraints from geometric considerations and formation considerations. 4.3.2.1. Geometric considerations. Ž1. 1.0 F C4Ž x. F 10.0: where C4 Ž x. s 1.25hrh f . The lower bound on this constraint ensures that the fracture height, h f , does not migrate to the bounding layers by more than 25% ŽSmith, 1992. of the pay zone thickness, h. The fracture height, h f , is calculated as a function of design variables by solving material balance relationships in Appendix A. Ž2. 1.0 F C5 Ž x. F 50.0: where C5 Ž x. s x frh f . The lower bound on this constraint ensures that the fracture half-length is always greater than the fracture height, which is a basic assumption in the PKN-C fracture geometry model. Ž3. 1.0 F C6 Ž x. F 5.0: where C6 Ž x. s warŽ4 Pd .; wa is the average fracture width calculated according to PKN-C fracture geometry model ŽAppendix A. and Pd is the proppant diameter. The lower bound on this constraint ensures that the average dynamic fracture width is at least four times the proppant diameter for effective proppant transport ŽSchechter, 1992.. Ž4. 0.5 F C 7 Ž x. F 1.0: where C 7 Ž x. s VfrVi ; Vf is the total fracture volume calculated from fracture dimensions and lateral and vertical shape factors, and Vi is the total injected volume ŽAppendix A.. The lower bound on this constraint ensures that the fluid efficiency, VfrVi , is greater than 0.5 and the upper bound ensures that the fracture volume is less than the injected volume. In typical fracturing treatments, the fluid efficiency varies between 0.7 and 0.4 ŽSmith, 1992..

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62 Table 1 Reservoir and wellbore properties Drainage area Average depth Thickness Shape Žsquare. Equivalent drainage radius Porosity Permeability Initial reservoir pressure Reservoir temperature Gas gravity Gas saturation Initial gas compressibility factor Initial gas viscosity Water compressibility Pore compressibility Skin factor Max. horizontal stress Min. horizontal stress Min. horizontal stress Žshale. Wellbore radius Flowing bottomhole pressure Tubing inside diameter Measured depth to mid-perf

Table 3 Proppant selection data 2.59 km2 Ž640 acres. 2286 m Ž7500 ft. 30.5 m Ž100 ft. 1609=1609 m2 Ž5280=5280 ft 2 . 908 m Ž2980 ft. 10% 1.98ey16 m2 Ž0.20 md. 30338 kPa Ž4400 psi. 93.3 8C Ž200 8F. 0.85 0.8 0.89 0.02831 mPa s 4.35ey10 Pay1 Ž3.0ey6 psiy1 . 12.5ey10 Pay1 Ž8.6ey6 psiy1 . 0.0 48265 kPa Ž7000 psi. 41370 kPa Ž6000 psi. 46196 kPa Ž6700 psi. 0.107 m Ž0.35 ft. 11722 kPa Ž1700 psi. 7.6 cm Ž2.992 in.. 2286 m Ž7500 ft.

4.3.2.2. Formation considerations. Ž1. 1.0 F C8 Ž x. F 10.0: where C8 Ž x. s PfcrrPtreat ; Pfcr is the formation critical pressure and Ptreat is the fracture treatment pressure, calculated as a function of design variables and formation in-situ stress ŽAppendix C.. The lower bound on this constraint ensures that the treatment pressure is kept below the formation critical pressure

Table 2 Fracture mechanics and operational limitations data Fracture geometry model Closure stress Young’s modulus Poisson’s ratio Leakoff coefficient Spurt loss coefficient Burst strength of the tube Horse power of the pump Pump efficiency Rated pressure of surface equipment

49

PKN-C 41,370 kPa Ž6000 psi. 34.99e6 kPa Ž5.075e6 psi. 0.20 9.84ey6 mrs 0.5 Ž2.5ey4 ftrmin0.5 . 0.0 89,635 kPa Ž13,000 psi. 10,444 kW Ž14,000 hp. 85% 96,530 kPa Ž14,000 psi.

Proppant type Specific Gravity Diameter Porosity Lab-measured frac-cond. at 9.76 kgrm2 proppant conc. and 41,370 kPa closure stress Conductivity damage factor

20r40 Westprop 3.1 grcc 0.063 cm Ž0.0248 in.. 0.35 2.022ey12 m3 Ž6700 md ft.

0.40

to prevent uncontrolled fracture growth ŽMeng and Brown, 1987; Nolte and Smith, 1981.. Ž2. 1.0 F C8 Ž x. F 10.0: where C8 Ž x. s D s hrPnet ; D s h is the difference between the minimum horizontal stresses in the pay zone and the bounding layers. The lower bound ensures that the fracture net pressure does not induce excessive fracture height growth into the bounding layers ŽNolte and Smith, 1981.. Ž3. 1.0 F C 10 Ž x. F 10.0: where C 10 Ž x. s D srŽ0.7Pnet .; D s is the difference between the maximum and the minimum horizontal in-situ stresses in the pay zone. The lower bound ensures that the fracture net pressure does not cause the initiation of auxiliary Žsecondary. fractures, which result in large fluid loss ŽNolte and Smith, 1981..

5. Examples and results A gas well located in a tight formation is used to illustrate application of the model. A 2.59-km2 Ž640-acre. square reservoir is considered with 30.5-m Ž100-ft. pay zone bounded above and below by shale formations subjected to higher stresses. The reservoir Table 4 Economics data Proppant cost Fracturing fluid cost Pumping cost for 14,000 hhp Fixed cost Gas price Discount rate Number of years

Ž$2.2rkg. $1.0rlb varying with viscosity $26.8rkW Ž$20rhhp. $10,000 $35.31rMSCM Ž$1.0rMSCF. 0.10 10

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

50

Table 5 Three initial designs to initiate optimization process Variable symbol

Initial Design 1

Initial Design 2

Initial Design 3

Max. NPV design

qi Žm3 rs. 0.053 0.053 0.0795 0.0643 t i Žs. 5800 5000 4000 7824 Pc Žkgrm3 . 1318 959 839 1797.4 @EOJ m ŽmPa s. 122 120 210 99.0 x f Žm. 511 610 168 761.2 NPV ŽM$. 14.850 12.040 11.423 16.520 Gp Ž10 6 m3 . 662.958 551.616 519.816 743.832 Ctr ŽM$. 0.6512 0.5525 0.5757 0.9920

and wellbore properties are presented in Table 1. Fracture mechanics data required for analysis is presented in Table 2, which also contains information on various operational limitations. A 20r40 Westprop is selected as the propping agent, which has a prescribed laboratory measured fracture conductivity data at 9.76 kgrm2 Ž2 lbrft 2 . proppant concentration. The laboratory data is modified from the laboratory concentration to the expected in-situ concentration; experience shows that a concentration greater than 4.88 kgrm2 Ž1 lbrft 2 . is difficult to achieve under most conditions ŽSmith, 1992.. A conductivity damage factor Ž0.4–0.5. also has been applied after correcting for in-situ concen-

tration to get reasonable fracture conductivity. Proppant selection data are presented in Table 3. Economic data are presented in Table 4. The price of fracturing fluid is related to viscosity by Eq. Ž14.. Net present value and production were calculated for 10 years. To demonstrate the benefit of using the optimization scheme in hydraulic fracturing design, three different designs are prepared arbitrarily Ži.e. without using an optimization procedure.. These designs are named as Initial Design 1, Initial Design 2 and Initial Design 3, presented in Table 5. Using these three initial designs as starting points, the optimization scheme was run to achieve maximum NPV. The optimization scheme achieved the final optimum design by executing a number of redesign iterations for each initial design. Three optimum designs, starting with the three initial designs, were very slightly different values within the convergence tolerance. One of the optimum designs is presented in the last column of Table 5. Iterative redesigns by the optimization scheme are plotted in Fig. 4 for the three initial designs. First three points on the graph at zero iteration represent NPVs of the three arbitrarily obtained initial designs. It is evident from the graph that these three designs have been improved significantly by the optimization scheme: specifically, Initial Design 1 by 11%, Initial Design 2 by 38% and

Fig. 4. Convergence to optimum design from different initial designs.

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62 Table 6 Optimum designs for four different objectives Variable symbol

Design 1

Design 2

Design 3

Design 4

qi Žm3 rs. 0.0635 0.0634 0.0396 0.0569 t i Žs. 8024 8039 4949 1819 Pc Žkgrm3 . 1797.4 1797.4 1729.1 1404.4 @EOJ m ŽmPa s. 99.66 99.75 89.31 84.0 x f Žm. 762.0 761.5 481.8 350.0 NPV ŽM$. 16.530 16.528 14.749 11.539 Ctr ŽM$. 1.00 0.998 0.551 0.424 Gp Ž10 6 m3 . 744.569 744.483 654.407 522.902

Initial Design 3 by 44%, with corresponding NPV increments by $1.7M, $4.5M and $5.0M, respectively, from initial designs. The way in which the optimum design changes with different design objectives is demonstrated by achieving four optimum designs for the four design objectives discussed in the previous section. The results are presented in Table 6. Design 1, Design 2, Design 3 and Design 4 represent maximum production design, maximum NPV design, maximum NPV and minimum treatment cost design, and target production with minimum treatment cost design, respectively. The maximum NPV design is very close to the maximum production design, whereas other designs are significantly different. A significant per-

51

centage Ž45%. of treatment cost saving is achieved over the maximum NPV design by setting the design objective to maximize NPV and minimize treatment cost ŽEq. Ž10... This saving, however, has resulted in 11% NPV reduction over 10 years. This shows the conflict between the design objectives, NPV and treatment cost; however, it is possible to achieve a compromised design, which may also be biased to a certain degree to a particular objective, by adjusting priority factors to individual objectives. For Design 4 ŽTable 6., 70% of the maximum production Ž744.569 = 10 6 m3 ., was set to achieve Ži.e. T1 s 0.7 = 744.569 = 10 6 m3 .. The optimum design ŽDesign 4. has achieved the full target. The time-dependent production rate and resulting decline in average reservoir pressures, if the reservoir were treated by the four optimum treatments mentioned above, are plotted in Figs. 5 and 6. The figures show that curves for Designs 1 and 2 Ži.e. maximum production and maximum NPV. are identical, because there is not much difference between these designs. The decline of reservoir pressure would be more if the well were produced by higher production rate treatments, such as by maximum production or NPV design. The slight kinks in the production rate curves show the transition between the transient state and the pseudo-steady-state during production. The curves indicate that the transient production period in tight reservoirs is significantly long and if

Fig. 5. Time-dependent production rates for optimum designs in Table 6.

52

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

Fig. 6. Time-dependent decline in reservoir pressure for optimum designs in Table 6.

the pseudo-steady-state production model is used during the entire production period, one will significantly overestimate the production Žtrace back of pseudo-steady-state curves will yield higher production rates compared to transient curves.. Any treatment design using such estimation will yield fractures and other treatment parameters such that when the design is executed in the field, the actual production will be significantly less than estimated. This will happen due to double effects: first, the estimation itself is optimistic, and second, the optimistic

estimation yields a design with a shorter fracture causing further reduction in production. The proppant scheduling designed for the four designs is presented in Fig. 7, which clearly shows that the treatment for higher production requires higher proppant concentration over a longer period but delayed proppant loading. Table 7 presents maximum NPV designs for a range of initial reservoir permeability. With increasing permeability, a short fracture and a high viscous fluid are optimum and consequently the optimum volume of injection Ž qi t i . is low. The optimum proppant concentration, however, remains almost unchanged. Productions and resulting NPVs are significantly high for high permeable reservoirs even with shorter fractures. This implies that larger fractures will be optimum in low permeable reservoirs to enhance production and NPV. Fig. 8 presents a trade-off analysis between treatment cost and production. For the upper curve, the maximum production design was achieved first using the objective function defined by Eq. Ž8. Žthe point at PF s 1.. Values of target production, T1 , were then set to achieve production factors Žthe production factor, PF is defined as the ratio of actual production to the maximum production. of 0.95, 0.9, 0.85, 0.8, 0.75, 0.7 and corresponding optimum designs were

Fig. 7. Designed proppant scheduling for optimum designs in Table 6.

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

53

Table 7 Effect of reservoir permeability on optimum design k Žm2 .

x f Žm.

qi Žm3 rs.

EOJ Pc Žkgrm3 .

mf ŽPa s.

t i Žs.

NPV ŽM$.

Gp Ž10 6 m3 .

9.9e y 17 Ž0.1 md. 1.98e y 16 Ž0.2 md. 3.47e y 16 Ž0.35 md. 4.95e y 16 Ž0.5 md. 7.43e y 16 Ž0.75 md. 9.9e y 16 Ž1.0 md. 1.98e y 15 Ž2.0 md.

761.4459 761.4532 721.2125 576.089 572.3921 547.3904 442.6638

0.066641 0.063348 0.066748 0.057273 0.057613 0.056817 0.052659

1797.396 1797.396 1797.257 1796.685 1797.396 1797.128 1797.396

0.097564 0.099748 0.102894 0.125825 0.126344 0.131340 0.159252

7468.428 8038.790 7102.378 6538.600 6446.192 6222.276 5324.190

15.21219 16.52874 19.00475 20.82530 23.15066 24.72175 28.17964

695.388 744.476 831.947 884.231 948.152 985.141 1041.429

obtained. Note that no effort was made to minimize the treatment cost for these designs, as the objective function was formulated as minimize Ž Gp y T1 ., and full production targets were always possible to achieve. Treatment costs associated with these designs are plotted against production factors Župper curve.. For the lower curve, optimum designs were obtained by setting production target, T1 to achieve production factors of 0.7, 0.75, 0.8, 0.85, 0.9, 0.95 and 1.0 in the objective function defined by Eq. Ž11.. Note that the minimization of treatment cost is attempted in these designs while the production is achieved as close to the target value as possible. The

full production targets were met up to the production factors of 0.75 beyond which a compromise is necessary with the production target in order to minimize the treatment cost. Thus, the actual production factor on the lower curve falls behind compared to the upper curve and it was not possible to achieve more than 0.88 production factor while T1 was set by production factors between 0.9 and 1.0. It is quite sensible that the treatment cost is always lower Žin the lower curve. when an attempt is made to minimize it. It is particularly interesting that only a 12% compromise with the production Ži.e. production factor of 0.88. saves about 44% of the treatment cost.

Fig. 8. Trade-off between treatment cost and production.

54

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

Fig. 9. Trade-off between treatment cost and NPV.

Similar trade-offs between treatment cost and NPV are shown in Fig. 9.

6. Conclusions The integrated optimization model for hydraulic fracturing design has accounted for complex interactions among fracture geometry, treatment parameters, material balance, operational limitations, formation characteristics, production behaviours and various potential design objectives. The study can be concluded as follows: Ž1. The model is capable of finding a true optimum design starting with any initial design, which may or may not satisfy all the design constraints, and the optimum design is significantly better than any arbitrary design. Ž2. An optimum design depends on the specified design objective. A maximum NPV design is almost identical to a maximum production design and both designs involve high treatment cost. By using the multi-objective optimization technique, it is possible to resolve the conflict between productionrNPV maximization and treatment cost minimization.

Ž3. Production from a well treated by maximum productionrNPV design causes production at higher rate and consequently greater reservoir pressure decline with time. Ž4. The transient production period is significantly long in a tight reservoir; therefore, the use of pseudo-steady-state production model only significantly overestimates the production. Because the optimum design is highly dependent on production estimation, the overestimate in production will result in an inadequate treatment design, which will result in significantly lower production in the field than estimated. Ž5. Proppant loading is delayed in maximum productionrNPV design and the proppant scheduling is continued longer with higher concentration. Ž6. A shorter fracture is optimum in a relatively high permeability reservoir. Ž7. Using the technique presented, an optimum treatment can be designed for a target production. Utilizing this capability of the model, trade-off analysis between productionrNPV and associated treatment cost can be performed for a wide range of target productions. A balanced solution in terms of productionrNPV and treatment cost can be identi-

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

fied from such analysis. For the set of reservoir, economic and other data used for such analysis in this study, it is found that 44% of treatment cost-saving is possible with only 12% reduction in productionrNPV over 10 years.

Nomenclature fracture area, m2 Af i-th design constraint Ci leak-off coefficient, mres 0.5 Cl treatment cost, $ Ctr normalizing factor to achieve the i-th objecDi tive E Young’s modulous, Pa plane strain modulous, Pa EX H length of the tubing, m hydraulic power of the pump, kW HPav HPreqd hydraulic power required, kW K consistency index of the fracturing fluid, Pa sn M number of design constraints, dimensionless N number of design variables, dimensionless Reynold’s number, dimensionless Nr burst strength of the tube, Pa P burst proppant concentration, kgrm3 Pc average proppant concentration, kgrm3 Pc proppant diameter, m Pd pump efficiency, dimensionless Peff formation critical pressure, Pa Pfcr price of fracturing fluid, $rm3 Pfl hydrostatic head in the tubing, Pa Phead priority factor to achieve the i-th objective, Pi dimensionless net fracture pressure, Pa Pnet price of proppant, $rkg Ppr Ppump pumping price, $rkW surface pressure generated, Pa Psurf treatment pressure, Pa Ptreat revenue generated at year n, $ Rn SF safety factor, dimensionless specific gravity of fracturing fluid, dimenSgf sionless specific gravity of proppant, dimensionless Sgp spurt loss coefficient, m Sp target value for the i-th objective, m3 Ti fracture volume, m3 Vf volume of fracturing fluid, m3 Vfl

Vi Vl Vpad Vpl Vpr Wpr d f hf i n qi ti t pad v wa wf x xf fp g h m me rp rf rr sh Ds n

55

total Žpad plus slurry. volume injected, m3 fracturing fluid loss, m3 pad Žclean. volume, m3 proppant-laden fluid Žslurry. volume, m3 proppant volume, m3 weight of proppant, kg internal diameter of the tubing, m friction factor, dimensionless fracture height, m discount rate ŽEq. Ž12.., dimensionless power law exponent, dimensionless injection rate, m3rs injection time, s pad time, s average velocity, mrs average width of the fracture, m width of the fracture at the wellbore, m design variables fracture half-length, m proppant porosity, dimensionless fracture shape factor, dimensionless fracturing fluid efficiency, dimensionless fracturing fluid viscosity, mPa s effective viscosity of fracturing fluid, Pa s proppant density, kgrm3 average density of fracturing fluid, kgrm3 relative density, dimensionless minimum horizontal stress, Pa difference of principal stresses, Pa Poisson’s ratio, dimensionless

Nomenclature in Appendix B (in oilfield unit) A drainage area of the reservoir ŽEq. ŽB9.., ft 2 reservoir water compressibility, 1rpsi Cw pore compressibility, 1rpsi Cf D non-Darcy flow co-efficient ŽDrMSCF. E gas expansion factor, scfrrcf initial gas expansion factor, scfrrcf Ein nondimensional fracture conductivity FCD cumulative gas production, MSCF Gp initial reservoir pressure, psi Pin flowing wellbore pressure, psi Pwf gas saturation, dimensionless Sg connate water saturation, dimensionless S wc T reservoir temperature, o R gas compressibility factor at average reserZg voir pressure, dimensionless

56

ae ct ct h h perf k ks k wf p qg re rw r wX s sf sX t t pss xf f gg mg

mg,wf

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

slope of gas expansion factor vs. pressure, dimensionless total system compressibility at average reservoir pressure, 1rpsi total system compressibility of the initial reservoir, 1rpsi thickness of the pay zone, ft perforated thickness, ft reservoir permeability, md near wellbore permeability, md fracture conductivity, md ft average reservoir pressure, psi production rate, MSCFrD drainage radius of the reservoir, ft wellbore radius, ft effective wellbore radius, ft near wellbore skin factor, dimensionless pseudo-skin factor, dimensionless effective skin factor, dimensionless production time, h time after which pseudo-steady state begins, h fracture half-length, ft porosity of the reservoir, dimensionless gas gravity, dimensionless viscosity of gas at average reservoir pressure, cp gas viscosity at flowing bottomhole pressure, cp

Acronyms EVOP Evolutionary operation GA Genetic algorithms INTEMOB Intelligent moving object KGD-C Khristianovich Geertsma deKlerk–Carter fracture model MSCFrD Thousand standard cubic feet per day MSCM Thousand standard cubic meter NBI Normal-boundary intersection algorithm NPV Net present value PA Polytope algorithm PKN-C Perkins Kern Nordgren–Carter fracture model SLP Sequential linear programming SQP Sequential quadratic programming SUMT Sequential unconstrained minimization technique

Acknowledgements Authors acknowledge the assistance of Simon Chipperfield of Santos for his valuable comments on technical and economic analysis. The authors also thank the two anonymous reviewers for their constructive comments and the Managing Editor, Earl C. Donaldson, for his editorial efforts, which combinedly have improved the quality of this paper.

Appendix A. The PKN-C fracture model Referring to Fig. 2, the relationship between fracture width at the wellbore, wf , fracture half-length, x f , treatment variables Žinjection rate, qi , and Newtonian fluid viscosity, m . and rock properties can be expressed as ŽValko and Economides, 1995.: wf s 3.27

m qi x f

ž

EX

1

/

4

Ž A1.

where EX s

E

Ž A2.

Ž1yn 2 .

in which E is Young’s modulus ŽPa., EX is plane strain modulus ŽPa. and n is Poisson’s ratio. In Eq. ŽA1., wf and x f are expressed in meters and qi , and m are in m3rs and Pa s. For fracturing fluids of non-Newtonian behaviour, the width equation is a function of power law parameters as follows: 1

ž w s 9.15

2 nq2

ž =K

2 nq2

f

1

/ 3.98 Ž /

ž

n 2 nq2

n

Ž 2 nq2 .

n

qin h1yn xf f E

1 q 2.14 n

.

X

/

ž

1 2 nq2

/

Ž A3.

where, n is the exponent Ždimensionless. and K is the consistency index ŽPa s n .. Based on a large amount of data provided by an industry, the power law parameters are correlated with viscosity as follows: n s 0.1756 my0 .1233

Ž A4.

K s Ž 0.457m y 0.0159 . 47.880

Ž A5.

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

Considering vertical and lateral width variation, the shape factor, g , is determined to be pr5 for PKN geometry and therefore the average fracture width can be expressed as: wa s Ž pr5 . wf

Ž A6.

During fracture growth, the general material balance relationship is: Vi s Vf q V l

Ž A7.

where Vi is the total fluid volume injected, qi t i Žm3 ., Vf is the fracture volume Žm3 . and V l is the fluid loss volume Žm3 .. Using the Carter II solution ŽHoward and Fast, 1957. for constant injection rate and considering fluid leak-off and spurt loss, the following relationship between fracture geometry and fluid injection can be derived ŽValko and Economides, 1995. from Eq. ŽA7.: xf s

Ž wa q 2 Sp . qi 4C l2p h f

exp Ž b 2 . erfc Ž b . q

2b

'p

y1

Ž A8. where,

bs

(

2C l p t i wa q 2 Sp

Ž A9.

and Sp is spurt loss coefficient Žm. and C l is leak-off coefficient Žmrs 0.5 .. Eqs. ŽA6. – ŽA9. constitute the solutions to the fracture propagation problem. Values of x f and qi are decided by the optimization algorithm, as free variables, in each iteration of the optimization process. The fracture height, h f , is then calculated solving Eq. ŽA8. using an inner iterative procedure. The net fracturing pressure, Pnet ŽPa., is then calculated as: pnet s

EX 2 hf

wf

Ž A10.

The fracture volume, Vf Žm3 ., and fracture area, A f Žm2 ., at any time for the linear propagating PKN-type fracture are defined, respectively, as follows for a two-sided symmetric fracture: p Vf s g h f x f wf Ž A11. 2

57

and Af s 4 x f hf Ž A12. The final propped width after the closure of fracture is: Wpr wp s Ž A13. A f Ž 1 y fp . rp where, Wpr is the weight of proppant Žkg.; fp is proppant porosity Ždimensionless. and rp is proppant density Žkgrm3 .. The fluid efficiency is: Vf hs Ž A14. Vi The relationship between total fluid volume requirements, Vi , and the pad volume, Vpad Žm3 ., based on the fluid efficiency, is given by ŽNolte, 1986; Meng and Brown, 1987.: 1yh Vpad s Vi Ž A15. 1qh The volume of proppant-laden fluid Žslurry., Vpl Žm3 ., which is the summation of proppant volume Ž Vpr . and fracturing fluid volume Ž Vfl ., is: Vpl s Vi y Vpad Ž A16.

ž /

The proppant weight, Wpr Žkg., is: Wpr s

Vpl

 0 1

Ž A17.

1

q

rp

Pc

where Pc Žkgrm3 . is the average proppant concentration. The average proppant concentration, Pc Žkgrm3 ., is: Pc s h Pc Ž A18. where Pc is the end of job ŽEOJ. proppant concentration Žkgrm3 ., which is optimized as a free design variable. The volume of fracturing fluid, Vfl Žm3 ., is: Wpr Vfl s Ž A19. Pc Based on a material balance, the continuous proppant addition Žramped proppant schedule versus time. is expressed by ŽNolte, 1986.: Pc Ž t . s Pc

ž

t y t pad t i y t pad

´

/

Ž A20.

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

58

where Pc Ž t . is the slurry concentration Žkgrm3 . at time t. The Pad time, t pad Žs. is: t pad s

Vpad

Ž A21.

qi

The exponent ´ depends on the fluid efficiency and is given by:

´s

1yh

Ž A22.

1qh

Appendix B. Production model

Ds

All equations are presented in oilfield units and the final production results are converted into SI units using appropriate conversion factors. B.1. Transient production model Production rate in the transient flow period can be estimated as follows ŽDake, 1978; Economides et al., 1994; Lee and Wattenbarger, 1996.: qg s

Ž Pin2 y Pwf2 . kh

ž

.1 h 6 = 10y5gg ky0 s

mg ,wf r w h2perf

k

where gg is the gas gravity, k s is near-wellbore permeability Žmd., h and h perf are net and perforated thickness Žft., respectively, and mg,wf is gas viscosity Žcp. evaluated at the flowing bottomhole pressure. The effective wellbore radius Ž r wX . equivalent to the fracture is ŽVan Everdingen, 1953; Mathews and Russell, 1967; Lee and Wattenbarger, 1996.: r wX s r w eys f

Ž B4.

fmg ct r wX 2

y1 X

y 3.23 q 0.869 s

/ Ž B1.

F s sf q ln

Ž B2.

in which sX is effectivertotal skin factor, s is near wellbore skin from formation damagerstimulation Ždimensionless., D is non-Darcy flow co-efficient ŽDrMSCF. and qg is production rate ŽMSCFrD.. In Eq ŽB1., Zg is gas compressibility factor Ždimensionless., m g is average viscosity of gas Ž cp ., k is reservoir permeability Žmd., h is reservoir thickness Žft., t is production time Žh., f is porosity Ž%., c t is total system compressibility Žpsiy1 ., r wX is effective wellbore radius Žft., Pin is initial reservoir pressure Žpsi., Pwf is bottomhole flowing pressure Žpsi. and T is reservoir temperature Ž8R.. Substitution of Eq.

xf rw

Ž B5.

where

where sX s s q Dq g

Ž B3.

where sf is the pseudo-skin which can be calculated as follows ŽCinco-Ley and Samaniego, 1981; Valko et al., 1997.:

1637TZg mg = log t q log

ŽB2. into Eq. ŽB1. yields a quadratic equation which can be solved for qg if other parameters are known. As mentioned in the main body of the paper, the average reservoir pressure declines with time due to cumulative production and therefore the pressure-dependent parameters Zg , mg , and ct should be evaluated at a reasonable time interval to predict realistic production. Consideration of time interval and adjustment of these parameters and production rate are discussed at the end of this appendix. An empirical relationship for non-Darcy flow coefficient is given by ŽEconomides et al., 1994.:

Fs

1.65 y 0.328u q 0.116 u 2 1 q 0.18u q 0.064 u 2 q 0.005u 3

u s ln Ž FCD .

Ž B6. Ž B7.

and FCD s

k wf kx f

Ž B8.

Here, k wf is the fracture conductivity in md-ft and FCD is the nondimensional fracture conductivity. Eq. ŽB1. is used to estimate the production rate in the transient flow period adjusting the value of t after each time interval up to t pss Žh., at which the

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

pseudo-steady state begins. The value of t pss can be estimated as ŽEarlougher, 1977.: t pss s

fm c t A 0.000264 k

Ž t DA . pss

Ž B9.

where c t is the system compressibility for the initial reservoir, A is the drainage area Žft 2 . and t DA has a characteristic value that depends on the drainage shape and well location. For a regular shape such as a circle or a square, the value of t DA is 0.1 or greater. The initial value of t DA is considered to be 0.1 and then adjusted iteratively by matching the production rate at the end of transient period with that at the start of pseudo-steady-state. B.2. Pseudo-steady state production model At the pseudo-steady-state, the production rate from a post-treatment well is ŽValko et al., 1997.: qg s

2 kh Ž p 2 y pwf .

1

=

1424mg Zg T

ln

ž

0.472 re xf

/ ž

q sf q ln

xf rw

/

i

Ž p . s pin y

i

i

Ž p . s pin q Ž for A / 0 .

yŽ B . "

(Ž

B.

i 2

2Ž A.

i

y 4Ž A. Ž C .

i

i

i

Ž A . s ae

Ž B.

i

Ž for A s 0, B / 0 .

ž

Cw Ž Sw .

i

Ž B . s ae q Ein i

Ž C . s Ein

ž

iy1

y Cf

1 y S wc Cw Ž Sw .

ž

Ž B12.

iy1

/

Ž B13.

y Cf

1 y S wc

/

Ž B14.

iy1

Ž Gp .

/

G

Ž B15.

Here, C w is reservoir water compressibility Žpsiy1 ., Cf is pore compressibility Žpsiy1 ., Ein is initial gas expansion factor Žscfrrcf., S w is water saturation, S wc is connate water saturation, Gp is cumulative gas production ŽMSCF., G is initial gas-in-place ŽMSCF. and ae is the slope of experimental E Žgas expansion factor. vs. P Žpressure. data. Note that ‘i’ superscripts represent time steps, not power. Values of Zg , mg and c t are adjusted at the i-th time step as follows:

žZ /

i

g

i Ž p . q pwf . Ž s 35.37

Ž B16.

i

2Ž E . T

Ž mg .

i

y6

s 4.0 = 10

ž

i Ž p . q pwf

2

/

q 0.0107

Ž B17. i

i Ž ct . s Ž Sg . Ž cg . i

Ž B18.

in which i Ž E . s Ein q ae

Ž Sg .

i

i Ž p . q pwf

½

s 1 y Ž Sw .

2

y p in

5

Ž Sw . s

ž

Ž B19.

i

Ž B20.

S wc 1 q C w Ž p in y Ž p .

i

1y

Ž B11.

i

where

Ž B10. where re is the drainage radius of the reservoir Žft. and p Žpsi. is the average reservoir pressure. The average reservoir pressure and other resulting parameters should be adjusted at a reasonable time interval to use both Eqs. ŽB1. and ŽB10.. One day is used as the time interval, D t. Successive time steps are indexed as i s 0, 1, 2, 3, . . . , n. At i s 0, all parameters correspond to the initial reservoir condition and they are used for production in the first day and then adjusted after each 24-h period. To use Eq. ŽB1., the cumulative value of t has been adjusted daily basis, i.e. at i s 1, t s 24 h, at i s 2, t s 48 h, at i s 3, t s 72 h, etc. The average reservoir pressure, p, can be calculated in the i-th time step Ž i G 1. as follows ŽGuo and Evans, 1993.:

ŽC.

59

Cw Ž Sw .

iy1

1 y Sw

y Cf

/

i

.

p in y Ž p .

i

Ž B21.

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

60 i Ž cg . s

1

Ž p.

Ž B22.

i

The cumulative gas production, Gp and the recovery factor, RF, at the end of i-th time step are:

Ž Gp .

i

s Ž Gp .

i

Ž RF. s

iy1

Ž Gp .

i

q Ž qg . D t

less.. The relative density can be calculated by ŽSmith, 1992.:

rr s

Ž B23.

i

G

Ž B24.

Note that D t in Eq. ŽB23. is expressed in days. The loop for calculations starts at the 0 time step at the beginning of transient flow state and continues in the pseudo-steady-state with 1-day time interval up to the specified number of years for consideration. In order to calculate NPV by Eq. Ž12., the total production at the n-th year, Gp n , is calculated by summing up the productions of 365 days in the n-th year.

Appendix C. Required pump capacity and developed surface pressure To develop the fracture geometry using the treatment parameters as decided by the optimization algorithm, the pump has to deliver the required net fracture pressure, Pnet . The pump capacity required to deliver Pnet can be estimated by considering static head and dynamic friction loss with respect to the surface level. The pressure developed at the surface to deliver Pnet in the fracture can also be computed. Only then can the constraints be formulated to adjust design such that the required pump capacity does not exceed the available pump capacity and the resulting surface pressure does not exceed the tubing strength and the pressure rating of other equipment. The calculation procedure is presented in the following subsections. C.1. Hydrostatic head Hydrostatic head ŽPa. in the tubing while pumping fracturing fluidrslurry is estimated by: Phead s 9.807H r f rr

ž ž

1000 q

1000 q

Pc Sgf Pc Sgp

Ž C2.

where Sgf is specific gravity of fluid Ždimensionless., Sgp is specific gravity of proppant Ždimensionless. and Pc is average proppant concentration Žkgrm3 ., which can be obtained from Eq. ŽA18.. C.2. Frictional pressure loss in the tubing The frictional pressure loss in the tubing is estimated by: 3

y2

D Pfric .loss s f Ž r f 10

Ž 3.281 . H . Õ 25.8 Ž 39.37d . CF = 6895 Ž C3. 2

where D Pfric.loss is frictional pressure drop ŽPa., f is friction factor Ždimensionless., Õ is average velocity in the tubing Žmrs., d is internal diameter of the tubing Žm., and CF is the factor to incorporate increase in frictional pressure loss due to proppant addition. The average velocity, Õ, in the tubing is: qi Õs Ž C4. p d2 For turbulent flow, which is usually the case during injection of fracturing fluid, the frictional factor, f, can be obtained by solving numerically the following equation ŽDodge and Metzner, 1959.:

(

1

4 s

f

Ž n.

0.75

n

log Nr Ž f .

Ž 1y 2 .

0.395 y

Ž n.

Ž C5.

1.2

where Nr is Reynold’s number Ždimensionless. which can be calculated by: 39.37d Nr s 928 Ž r f 10y2 . Ž 3.281Õ . Ž C6. 1000 me The effective viscosity me ŽPa s. in the above equation can be calculated as:

Ž C1.

where H is length of the tubing Žm. from surface to the fracture centre, r f is average density of fracturing fluid Žkgrm3 ., rr is relative density Ždimension-

/ /

me s K Ž 39.37d .

1y n

Ž 3 q ny1 . 96 Ž 3.281Õ .

1yn

n

Ž 0.0416 .

n

Ž C7.

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

where n and K are power law parameters, Eqs. ŽA4. and ŽA5.. The factor, CF, to incorporate increase in frictional pressure loss due to proppant addition ŽSmith, 1992. is: 0 .8 CF s m0.2 r = rr

Ž C8.

where

mr s 1 q 2.5 w q 10.05w 2 q 0.00273e 16 .6 w

Ž C9.

in which Pc

ws

ž 1000S

gp q Pc

/

Ž C10.

C.3. Required pump capacity and surface pressure The treatment pressure required to deliver Pnet is: Ptreat s Pnet q s h

Ž C11.

where s h is the minimum in-situ stress ŽPa. in the zone to be fractured. The pressure required at the surface to develop the above treatment pressure is calculated as follows: Psurf s Ptreat q D Pfric .loss q D Psurf .eq .loss y P head

Ž C12. where D Psurf.eq.loss is the pressure loss in the surface equipment. The pump capacity Žkilowatt. required to develop the above surface pressure is: qi HPreqd s Psurf Ž C13. 1000

References Aggour, T.M., Economides, M.J., 1998. Optimization of the performance of high-permeability fractured wells. SPE 39474, Int. Symp. on Formation Damage Control, Lafayette, LA, Feb., 18–19. Balen, M.R., Meng, H.-Z., Economides, M.J., 1988. Application of the NPV Žnet present value. in the optimization of hydraulic fractures. SPE 18541, Eastern Reg. Meet., Charleston, WV, Nov., 1–4. Bouteca, M.J., 1988. Hydraulic fracturing model based on a three-dimensional closed form: tests and analysis of fracture geometry and containment. SPE Prod. Eng., Nov., 445–454. Bunday, B.D., 1984. Basic Optimisation Methods. Edward Arnold, London, 128 pp.

61

Carroll, J.A., Horne, R.N., 1992. Multivariate optimization of production systems. J. Pet. Technol., Jul., 782–789. Cinco-Ley, H., Samaniego, V.F., 1981. Transient pressure analysis for fractured wells. J. Pet. Technol., Sep., 1749–1766. Dake, L.P., 1978. Fundamental of Reservoir Engineering. Elsevier Scientific Publishing, Amsterdam: 249, 264–286. Das, I., Dennis Jr., J.E., 1998. Normal-boundary intersection: a new method for generating the pareto surface in nonlinear multicriteria optimization problems. SIAM J. Optim. 8, 631– 657. Dodge, D.W., Metzner, A.B., 1959. Turbulent flow of non-Newtonian systems. AIChE J. 5 Ž2., 189–204. Earlougher, R.C., Jr., 1977. Advances in Well Test Analysis. SPE Monograph Series, Vol. 5: 74–89, 192–221. Economides, M.J., Hill, A.D., Ehlig-Economides, C., 1994. Petroleum Production Systems. Prentice Hall PTR, New Jersey, Chaps. 2, 4, 16 and 17. Fujii, H., Horne, R., 1995. Multivariate optimization of networked production systems. SPE Prod. Facil., Aug., 165–171. Geertsma, J., De Klerk, F., 1969. A rapid method of predicting width and extent of hydraulically induced fractures. J. Pet. Technol., Dec., 1571–1581. Ghani, S.N., 1989. A Versatile Procedure for Optimization of a Nonlinear Nondifferentiable Constrained Objective Function. AERE R13714, United Kingdom Atomic Energy Authority, Nuclear Physics and Instrumentation Division, Harwell Laboratory, Oxfordshire, UK. Goldberg, D.E., 1989. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley Publishing, Reading, MA, p. 412. Guo, G., Evans, R.D., 1993. Inflow performance and production forecasting of horizontal wells with multiple hydraulic fractures in low-permeability gas reservoirs. SPE 26169, Gas Tech. Symp., Calgary, Canada, Jun., 28–30. Hossain, M.M., 2001. Reservoir Stimulation by Hydraulic Fracturing: Complexities and Remedies with Reference to Initiation and Propagation of Induced and Natural Fractures. PhD Thesis, The University of New South Wales, Sydney. Howard, G.C., Fast, C.R., 1957. Optimum fluid characteristics for fracture extension. Presented at the Spring Meet. of the MidContinent District, Division of Production, Tulsa, Okla., April. Drilling and Production Practices, API, ŽAppendix by R.D. Carter: Derivation of the general equation for estimating the extent of the fractured area., 261–270. Jacoby, S.L.S., Kowalik, J.S., Pizzo, J.T., 1972. Iterative Methods for Nonlinear Optimization Problems. Prentice-Hall, Englewood Cliffs, NJ, p. 274. Khristianovitch, S.A., Zheltov, Y.P., 1955. Formation of vertical fractures by means of highly viscous fluids. Proc. 4th World Pet. Congr., Rome 2, 579–586. Lee, J., Wattenbarger, R.A., 1996. Gas reservoir engineering. SPE Textb. Ser. 5, Chaps. 5–7. Mahrer, K.D., 1999. A review and perspective on far-field hydraulic fracture geometry studies. J. Pet. Sci. Eng. 24, 13–28. Mathews, C.S., Russell, D.G., 1967. Pressure buildup and flow tests in wells. SPE Monogr. Ser., 1, p. 167. McFarland, J.W., Lasdon, L., Loose, V., 1984. Development

62

M.M. Rahman et al.r Journal of Petroleum Science and Engineering 31 (2001) 41–62

planning and management of petroleum reservoirs using tank models and nonlinear programming. Oper. Res. 32 Ž2., 270– 289. Meng, H.-Z., Brown, K.E., 1987. Coupling of production forecasting, fracture geometry requirements and treatment scheduling in the optimum hydraulic fracture design. SPE 16435, Low Permeability Reservoirs Symp., Denver, CO, May 18–19. Mistree, F., Hughes, O.F., Phouc, H.B., 1981. An optimization method for the design of large, highly constrained complex systems. Eng. Optim. 5, 179–197. Mohaghegh, S., Balan, B., Ameri, S., 1996. A hybrid, neurogenetic approach to hydraulic fracture treatment design and optimization. SPE 36602, Ann. Tech. Conf. and Exh., Denver, CO, October 6–9. Mohaghegh, S., Balan, B., Platon, V., Ameri, S., 1999. Hydraulic fracture design and optimization of gas storage wells. J. Pet. Sci. Eng. 23, 161–171. Morales, R.H., Abu-Sayed, A.S., 1989. Microcomputer analysis of hydraulic fracture behavior with a pseudo-three-dimensional model. SPE Prod. Eng., Feb., 198–205. Nelder, J.A., Mead, R., 1965. A simplex method for function minimization. Comput. J. 7, 308. Nolte, K.G., 1986. Determination of proppant and fluid schedules from fracturing pressure decline. SPE 13278, SPE Prod. Eng., Jul., 255–265. Nolte, K.G., Smith, M.B., 1981. Interpretation of fracturing pressures. J. Pet. Technol., Sep., 1767–1775. Nordgren, R.P., 1972. Propagation of a vertical hydraulic fracture. SPE J., Aug., 306–314. Perkins, T.K., Kern, L.R., 1961. Widths of hydraulic fractures. J. Pet. Technol., Sep., 937–949. Pierre, D.A., 1969. Optimization Theory with Applications. Wiley, New York, p. 612. Polyak, B.T., 1987. Introduction to Optimization. Optimization Software, New York, p. 438.

Rahim, Z., Holditch, S.A., 1995. Using a three-dimensional concept in a two-dimensional model to predict accurate hydraulic fracture dimensions. J. Pet. Sci. Eng. 13, 15–27. Rahman, M.K., 1996. Optimization of panel forms for improvement in ship structures. Struct. Optim. 11, 195–212. Rahman, M.K., 1998. Automated optimization of transverse frame layouts for ships by elastic–plastic finite element analysis. Struct. Optim. 15, 187–200. Rahman, M.K., Caldwell, J.B., 1992. Rule-based optimization of midship structures. Mar. Struct. 5, 467–490. Schechter, R.S., 1992. Oil Well Stimulation. Prentice Hall, Englewood Cliffs, NJ, Chaps. 8–12. Settari, A., Cleary, M.P., 1986. Development and testing of a pseudo-three-dimensional model of hydraulic fracture geometry. SPE Prod. Eng., Nov., Trans. AIME 281, 449–466. Smith, M.B., 1992. Hydraulic Fracturing. 2nd edn. NSI Technologies, Tulsa, OK, Chaps. 3 and 5. Spellucci, P., 1998. An SQP method for general nonlinear programs using only equality constrained subprograms. Math. Program. 82, 413–448. Valko, P., Economides, M.J., 1995. Hydraulic Fracture Mechanics. Wiley, Chichester, England, Chaps. 4, 8 and 9. Valko, P., Oligney, R.E., Economides, M.J., 1997. High permeability fracturing of gas wells. Gas TIPS ŽFall. 3 Ž3., 31–40. Van Everdingen, A.F., 1953. The skin effect and its influence on the productive capacity of a well. Trans. AIME 198, 171–176. Wilson, A.G., Coelho, J.D., Macgill, S.M., Williams, H.C.W.L., 1981. Optimization in Locational and Transport Analysis. Wiley, Chichester, UK, p. 283. Wismer, D.A., 1971. Optimization Methods for Large-Scale Systems. McGraw-Hill, New York, p. 335. Yang, Z., Crosby, D.G., Khurana, A.K., 1996. Multivariate optimization of hydraulic fracture design. Aust. Pet. Prod. Explor. Assoc. J., 516–527.