Optimization of hydraulic fracturing design under spatially variable shale fracability

Journal of Petroleum Science and Engineering 138 (2016) 174–188

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol

Optimization of hydraulic fracturing design under spatially variable shale fracability Atefeh Jahandideh, Behnam Jafarpour n University of Southern California, United States

art ic l e i nf o

a b s t r a c t

Article history: Received 4 March 2015 Received in revised form 14 August 2015 Accepted 29 November 2015 Available online 2 December 2015

Current hydraulic fracturing designs are based on equal fracture intervals (spacing) along horizontal wells. Typical fracture spacing ranges from 50 to 600 feet while the fracture half-lengths normally vary between 100 and 600 feet. Symmetric spacing designs assume homogeneous shale mechanical property distributions, e.g. rock brittleness or fracability, which may not be the case in many shale plays. Fracability of the shale rock is an important property that controls the efficiency of the fracturing process. Therefore, the heterogeneity in shale fracability can have an important impact on the design and optimization of the hydraulic fracturing process. In particular, the optimal number, location, and length of the fractures can depend on rock fracability distribution. In this paper, we develop an optimization approach for hydraulic fracturing design under geospatial variability in shale fracability and investigate several aspects of the proposed fracture design optimization approach. In particular, we optimize the hydraulic fracturing design by implementing a variant of the Simultaneous Perturbation Stochastic Approximation algorithm to maximize the net present value of the shale asset, including the cost of fracturing and the revenue from gas production. The optimization framework provides a systematic design approach for hydraulic fracturing of unconventional reservoirs and can be applied to improve production efficiency and reduce the cost and environmental impact of hydraulic fracturing. We demonstrate the effectiveness and suitability of the proposed method using a series of numerical experiments in shale gas development. & 2015 Elsevier B.V. All rights reserved.

Keywords: Fracture design Optimization Hydraulic fracturing Shale fracability Rock brittleness

1. Introduction Shale gas plays are ultra-low permeability (10–100 nanodarcy) formations with porosity in the range of 2–8%, typical net thickness of 50–600 feet and a total organic content (TOC) of 1–14%. Because of the compact and impermeable nature of the shale rocks gas production without permeability stimulation is not profitable. Two key technologies that have enabled commercial production from shale gas reservoirs are horizontal drilling and multi-stage hydraulic fracturing. A horizontal well increases the reservoir– borehole contact while hydraulic fracturing stimulation enhances the permeability of the formation to hydrocarbon flow toward the well. To hydraulically induce the fractures, a mixture of fracturing fluid and proppants is injected into the formation at a very high pressure to cause rock failure in shear mode along its plane of weakness. The fracturing process is repeated for multiple locations (stages) along the horizontal well to maximize the stimulated n Correspondence to: Petroleum Engineering, Viterbi School of Engineering, University of Southern California, 925 Bloom Walk, HED 313, Los Angeles, CA 90089, USA. E-mail address: [email protected] (B. Jafarpour).

http://dx.doi.org/10.1016/j.petrol.2015.11.032 0920-4105/& 2015 Elsevier B.V. All rights reserved.

reservoir volume (SRV). The primary production mechanism in shale gas reservoirs is attributed to a combination of induced hydraulic fractures and natural fractures (when present). The shale gas exists as free gas stored in the matrix pores and fractures, and as adsorbed gas to organic material. Initially, the free gas is transported toward the production zones through the fractures. As the pressure decreases during depletion, the adsorbed gas is also released and transported to the matrix–fracture interface and is produced through the fractures. The complexity of the transport processes in fractured shale gas reservoirs make the decline curve analysis inapplicable, and calls for more sophisticated predictive modeling techniques. Simplistic models and assumptions have been proposed for describing the behavior of hydraulic fracturing and for predicting the production from hydraulically fractured shale gas plays. In practice, the fracture initiation points (fracturing stages) are chosen by dividing the wellbore into equally spaced zones, assuming homogeneous geomechanical and flow properties in the shale formation. However, field experience suggests that even under identical fracturing scenarios in the same formation, the production performance can be markedly different. This behavior can, in part, be explained by geo-spatial variability in shale

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

geomechanical and flow properties that identify the fracturing target locations (or sweet spots). To minimize the fracturing cost and environmental impacts, it is important to optimize the design of the fracturing process by understanding the parameters that affect the economic development of unconventional reservoirs (King, 2010; Boulis et al., 2013). Without accounting for the heterogeneity in rock geomechanical properties along the lateral extent of the well, similar treatment efforts may result in different production performance (Daniels et al., 2007). Therefore, to improve shale gas development, it is important to understand and characterize the variability in the spatial distribution of rock mechanical properties and its mineralogy; in particular, the spatial variability in the reservoir rock brittleness, which determines its fracability at different locations (i.e., rocks responsiveness to the fracturing process). Fracability is related to material brittleness and ductility. It is widely accepted that rock brittleness is mainly controlled by its mineralogy, mainly the relative abundance of carbonate and quartz compared to clay content (Buller et al., 2010). Several field experiments in the literature have documented the significance of optimizing hydraulic fracturing design based on brittleness. One example is the BP CGU 13-17H horizontal well located in the Carthage field, Panola County, Texas. Based on LIBS chemostratigraphy data, this well showed a wide range of relative brittleness. However, there was no attempt to use the horizontal well data to optimize the fracturing design. The well was stimulated using equally spaced stages (every 110 ft along the horizontal well), which led to undesirable production outcome in fracture stages that were placed in the ductile areas (Buller et al., 2010). Another example is the Longmaxi black shale in South China, in which the effect of brittleness has been demonstrated through fracture tracing and micro-seismic monitoring. In this field, hydraulic fractures in brittle intervals initiate and propagate more efficiently, resulting in larger fracture apertures and higher production (Li et al., 2013). Whereas fracturing in brittle intervals (lower clay volume) results in wider fracture width, better proppant placement, and hence higher gas production, fracturing in ductile intervals (higher clay volume) can lead to less effective (shorter and narrower) fractures that are subject to high pressure drawdown and proppant embedment (Buller et al., 2010). This important principle is not followed in the current practice in hydraulic fracturing, where equally spaced fracture stages with higher intensity are induced to increase fracture–reservoir contact and gas production. In recent years, the importance of brittleness on fracturing effectiveness is recognized and brittleness properties of the rock have been suggested as fracture initiation points (Parker et al., 2009; Li et al., 2013). In this paper, we propose a systematic framework to optimize fracture design parameters (e.g. fracture location, fracture half-length and number of fracture stages) for a more cost-effective hydraulic fracturing operation. From a geomechanical standpoint, rock brittleness (fracability) is related to the Poisson's ratio and Young's (Tensile/Elastic) modulus. The former is defined as the negative ratio of the transverse to axial strain while the latter represents the ratio of stress to strain along an axis, following the Hook's law. These two measures are combined to reflect rock's ability to fail under stress (Poisson's ratio) and to maintain a fracture (Young's modulus) when it undergoes fracturing treatment. Shale brittleness distribution can be quantified by combining these two mechanical properties (Rickman et al., 2008). Brittle shales have higher Young's modulus and lower Poisson's ratio (usually because of high silica or calcite content). They tend to fracture readily and the induced fractures are more likely to remain stable after treatment (King, 2010). On the other hand, ductile rocks (usually with high clay content) are less amenable to fracturing and tend to heal

175

(close) faster after fracturing. Rock ductility increases the energy required to induce a fracture and can lead to short fractures that are undesirable as they provide less contact with the formation while requiring more proppants. . Given these relations, fracability can be estimated using LWD measurements, core analysis, and chemostratigraphy analysis of drill cuttings. The spatial distribution of fracability, on the other hand, must be inferred from limited available data using spatial modeling and geostatistical interpolation. Since horizontal drilling and hydraulic fracturing stimulation are costly and can have potential environmental impacts, optimizing the fracturing design parameters, including the number, spacing, and geometric attributes of the fracture stages, is critical. Recently, some studies have focused on optimizing well spacing, fracture spacing and fracture length in order to maximize production and minimize the fracturing cost. Yu and Sepehrnoori (2013) used response surface methodology, for two horizontal wells, to optimize the design parameters that are related to reservoir properties and fracture attributes, namely permeability, porosity, fracture spacing, fracture half-length, fracture conductivity and well distance. They reported alternative optimized combinations of these parameters for different gas prices. Boulis et al. (2013) applied a stochastic forward modeling workflow to evaluate the optimal number of wells based on a semi-analytical shale gas model for production prediction in a transverse-fractured horizontal well. They used Monte-Carlo-based stochastic modeling to generate alternative plausible models for forecasting production performance. Ma et al. (2013) applied the Discrete Simultaneous Perturbation Stochastic Approximation (DSPSA) method and genetic algorithm (GA) to optimize hydraulic fracture placement under the spatial variability of permeability throughout the shale formation and compared the performance of these methods with gradient-based finite difference approach for shale gas production. Based on the results of their study, they suggested using hybrid methods that combine the advantages of both types of algorithms. Wilson and Durlofsky (2012) applied a derivativefree generalized pattern search optimization technique to a twodimensional surrogate model to determine optimal fracture location, length and number of stages based on the variability in initial reservoir pressure. The surrogate model parameters were calibrated to reproduce the predictions from a full-physics-based model. Meyer et al. (2010) presented an analytical solution for predicting the behavior of symmetric transverse hydraulic fractures for a horizontal well in a homogenous formation and used the solution to optimize the number of fracture stages and fracture length for equally spaced fractures. An important aspect that has not been considered in optimizing the hydraulic fracturing design is geomechanical rock properties, in particular the spatial variability in rock fracability. Characterization of geospatial variability in geomechanical rock properties of shale formations can provide insight into the fracturing response of the rock, which can, in turn, be used to optimize the fracturing design and implementation. In this paper, we consider geomechanical reservoir parameters and their spatial variability (heterogeneity) for fracture placement optimization. The motivation behind our work is to incorporate rock fracability in optimizing the parameters of hydraulic fracturing design. We use the SPSA algorithm to optimize fracture attributes such as location, length and number of stages in shale formations by considering the geospatial variability in shale fracability.

176

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

is typically above 0.3, with lower values being more favorable. Grieser and Bray (2007) propose the following empirical relation for quantifying rock brittleness based on averaged Young's modulus and Poisson's ratio:

2. Methodology 2.1. Rock fracability Rock brittleness plays an important role in its response to hydraulic fracturing. Brittleness has been related to Young's modulus and Poisson's ratio as two important geomechanical properties (Rickman et al., 2008; Sone, 2012), which are typically estimated from geophysical well logs such as sonic logs (Breyer et al., 2010). Shale brittleness can also be characterized based on composition and mineralogy. Carbonate, clay, and quartz are the main mineral content in shale gas reservoirs (Karastathis, 2007). Quartz and silica show brittle behavior with high fracability while clay is highly ductile with poor fracability. An example of composition-based rock brittleness index is defined by Jarvie et al. (2007), where brittleness index (BI) is quantified based on the relative amounts of the clay, calcium and quartz content of the rock, as follows:

BI =

Quartz Quartz + Clay + Carbonate

(1)

Based on this definition, rock brittleness increases with increasing quartz content and decreasing clay and carbonate contents. Wang and Gale (2009) extended the above brittleness index definition by including dolomite content and the TOC as follows

BI =

Quartz + Dolomite Quartz + Dolomite + Limestone + Clay + TOC

(2)

This definition shows that the impact of dolomite and TOC on brittleness is in opposite directions. However, some studies on Barnett shale, e.g., (Perez Altamar and Marfurt, 2014), report a direct relation between high TOC and high brittleness. Geomechanical definition of rock brittleness based on elastic properties is also quite useful for simulation of fracture generation, propagation, and healing with time. Hydraulic fracturing can induce rock failure in shear and tensile modes. During hydraulic fracturing, shear failure is dominant in reservoirs with natural fractures (Palmer et al., 2009). Due to shear forces, fracture walls are displaced out of their original position, causing surface roughness and asperities. In addition, fracturing in shear mode can open pre-existing natural fractures and planes of weakness through shear slippage (Hossain et al., 2002). In brittle material, when fracturing stress is exceeded tensile (shear) fractures develop if the displacement is perpendicular (tangential) to the displaced surface. On the other hand, ductile material can accommodate relatively large displacements (deformation) without fracturing. For brittle rocks, the shear modulus

G =

F /A FL = ∆x/L A∆x

(3)

takes larger values. For isotropic material, the shear modulus is related to Young's modulus and Poisson's ratio through

G=

E 2 (1+ϑ)

(4)

where ϑ is the Poisson's ratio and E is the Young's modulus. From this equation, an increase (decrease) in Young's modulus (Poisson's ratio) results in an increase in the shear modulus or brittleness. It is important to note that shales are generally considered as transverse isotropic rocks. The shear moduli of important minerals in shale reservoirs are [Quartz ¼44 GPa, Clay¼6.85 GPa, Calcite¼30 GPa, Dolomite¼ 43 GPa], showing that quartz and clays have the highest and lowest shear moduli, respectively. In shale gas reservoirs, Young's modulus values typically vary between 2 × 106 psi and 6 × 106 psi, where any value greater than 4 × 106 psi is desirable. The Poisson's ratio for shale gas reservoirs

BIavg =

EBrit =

ϑbrit =

1 (EBrit +ϑbrit ) 2 E − Emin Emax − Emin

ϑ − ϑmax ϑmin − ϑmax

(5)

(6)

(7)

Here, Emin and Emax are the minimum and maximum of observed Young's modulus, corresponding to the minimum and maximum brittleness, respectively, and ϑmin and ϑmax are the minimum and maximum observed Poisson's ratios that correspond to maximum and minimum brittleness, respectively. Grieser and Bray (2007) verified their model by comparing the brittleness measurements with well production data (based on 12month gas production) for 28 shale wells in Wise and Denton Counties in Texas. They reported a positive correlation between log analysis (empirical BIavg ) and well productivity. In another recent publication, Yang et al. (2015) verified the BIavg effect in predicting natural fractures locations and the targets of hydraulic fracturing for Horn River shale. These authors found a strong positive correlation between hardness and brittleness calculated from well logs (using Grieser and Bray correlation) for different members of the Horn River group. They also reported a positive correlation between the brittleness of the core samples with their natural fracture density. Based on their observations, the authors concluded that, in the absence of borehole image well logs, a brittleness or hardness value can be used to predict the location of natural fractures and potential targets for hydraulic fracturing. In this paper, brittleness is defined based on the method proposed by Grieser and Bray (2007), which require measurements of Young's modulus and Poisson's ratio. Sonic log data have been successfully used to derive Young's modulus and Poisson’ ratios in vertical wells. However, in lateral wells, due to the complexity of layered systems and scarcity of core measurements in horizontal wells, the derivation of geomechanical properties is more challenging. Nonetheless, for more complex wells, dipole sonic logs have been applied to determine rock mechanical properties from directional differences in rock mechanics (Smith et al., 2013). To distinguish ductile shales from brittle shale rocks, a lower bound (0.28 in this paper) is defined for rock brittleness. For any value below this threshold, which is field site-specific and depends on rock properties and responsiveness to fracturing, productive fractures cannot be induced. High values of brittleness imply greater complexity in fracture geometry, resulting in higher conductive flow paths for hydrocarbon. For example, Eagleford and Haynesville or Caney shale plays have rock brittleness as low as 0.31. To include the effect of fracability on gas production, we define the following linear model that relates fracture conductivity to brittleness index,

kf wf =

20 − 0. 1 (BI−0.28)+0.1 0. 7 − 0. 28

(8)

where kf wf denotes fracture conductivity in md-ft. While simple, this relation characterizes a positive correlation between fracture conductivity and brittleness index, as reported based on field studies. While the true relationship between conductivity and brittleness may not be linear, field data and core observations suggest a positive correlation between these two parameters. A more representative relation may be obtained by resorting to

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

detailed coupled flow and geomechanic simulation models. For our optimization purpose in this paper, we use Eq. (8) to approximately capture the relation between the BI and fracture conductivity. With this model, cells with fracability less than 0.28 are treated as matrix (with permeability of 0.00015 md). The dimensionless fracture conductivity (FCD ) is defined as the ratio of fracture to reservoir conductivity, that is,

FCD =

kf . wf k m. xf

Table 1 Parameter used for calculating the objective function (  NPV). Hydraulic Treatment (per stage)

Cost (1000 $)

Parameter

250-ft frac half-length

100

Gas price

500-ft frac half-length

125

750-ft frac half-length 1000-ft frac half-length

150 175

Value

3.5 $/MSCF Water disposal cost 3.3 $/bbl Discount rate 0.1

(9)

which is a key parameter that controls productivity enhancements. In this equation, kf . wf is the fracture conductivity, k m represents the permeability of matrix and xf stands for fracture halflength. For reservoirs with high permeability, short and wide fractures with high values of kf wf (fracpacks) are desirable. For low permeability formations (0.01 md) long and narrow fractures are more desirable as they can produce as much as the reservoir can deliver.

2.2. Hydraulic fracturing optimization method 2.2.1. Optimization objective function We formulate the optimization problem of finding the best fracture length and location along the horizontal wells using the following minimization form:

^ L , u^ = argmin fLϵZ1, uϵZ2 (L (u), u)

177

(10)

where f (·) is a user-specified objective function, L (u) is a vector that contains the fracture length for each stage and u represents the locations of each fracturing stage. For a discretized computer model, the location and length of each fracture are defined by the intersected cells. Hence, the feasible set Z2 for the location of fractures is the cell indices along the horizontal well. For fracture length, the feasible set Z1 consists of integer numbers representing the number of cells intersected by the fractures (in the direction ^ perpendicular to the horizontal well). The solutions L and u^ correspond to the fracture length (after multiplying by the cell dimension) and fracture locations, respectively. This formulation only considers fracture locations and their lengths as decision variables. In general, the number of wells, number of fracture stages, and well spacing could also be considered as additional decision variables. The objective function of our optimization problems represents production performance of the asset. A number of economic/productivity measures may be used to define the objective function. The Net Present Value (NPV) is one of the most frequently used measures to evaluate the economic viability of a project. The NPV accounts for the project cash flow and time value of money. In shale development projects, the NPV includes revenues from gas production as well as capital and operating expenditures, all discounted to obtain their present value. In this paper, to maximize NPV, we minimize the NPV value with a negative sign (i.e.,  NPV). In our problem, the NPV consists of the revenues from gas production and the costs associated with water disposal/recycling as well as the capital expense associated with hydraulic fracturing per unit length of treatment. Since we assume a fixed number of wells, without loss of generality, the objective function does not include the capital cost associated with drilling. Hence, the objective function of our minimization problem takes the form:

f (L (u), u) = − NPV ⎛ K ⎡ ⎜ ⎢ =−⎜ ∑⎢ ⎜ k=1 ⎣ ⎝ P

−

n

P

∑

(Q

⎤ − Q wk , j. rw . ∆tk ⎥ ⎥ tk (1 + b) 365 ⎦

k g, j. rg

j=1

)

⎞

∑ ∑ Li, j. Ci, j ⎟⎟ j=1 i=1

⎟ ⎠

(11)

where Ci, j is the cost of ith fracturing stage in the jth well ($/ft ), L i, j is the length of ith fracturing stage in the jth well (ft), rg is the gas price ($/MSCF); rw denotes water disposal cost , and Q g and Q w are the gas production rate (MSCF/D) and water production rate (bbl/day), respectively. The notation ∆tk, b , n, and p represent, respectively, the k th (time) step-size, the discount rate, the number of fracture stages for each well, and the number of production wells. Table 1 summarizes the values assigned to each quantity in calculating the objective function in the numerical experiments of this paper. The cost per unit length (ft) of each fracture stage was calculated based on the work of Schweitzer and Bilgesu (2009). Fig. 1 shows the cost of a fracturing stage as a function of fracture half-length. It is important to note that the brittleness index has a direct impact on the value of NPV (as described in Eq. (11)). In our model, the BI distribution is used to determine fracture conductivity (using Eq. (8)), which in turn influences the predicted (by commercial simulator, e.g. Eclipse E300) gas production ( Q g ) and water production ( Q w ). Table 2 shows the significant impact of BI value on the NPV based on Eq. (8) and the production simulations results. To solve the optimization problem, we combine a numerical optimization algorithm with Eclipse E300 (for predicting the production and estimating the NPV values). For optimization, we use the SPSA algorithm, which was first introduced by Spall (1992) and later extended to several variants (Spall, 2000, 2003). The SPSA algorithm has also been applied to well placement optimization in conventional reservoirs (Bangerth et al., 2006; Li and Jafarpour, 2012; Li et al., 2013). For brevity of our presentation, the description of the SPSA algorithm for solving our optimization problem is provided in Appendix A.

Fig. 1. Fracture stage cost versus fracture half-length.

178

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

Table 2 Effect of BI value on NPV. BI

NPV ($)

0.28 0.35 0.5 0.7

1.64e þ 5 5.34e þ 7 7.06e þ 7 7.65e þ 7

Fig. 2. Planar fractures with very high permeability to approximate the behavior of fractures.

3. Numerical experiments 3.1. Reservoir model In our numerical experiments, the reservoir model is constructed based on available data for Barnett shale in (Al-Ahmadi and Wattenbarger, 2011; Yu and Sepehrnoori, 2013). Table 3 lists the parameters used in setting up the model. Given the complexity of hydraulically fractured reservoirs, we perform numerical simulation to predict gas production from two horizontal wells with multiple fracturing stages (Cipolla et al., 2010). While analytical and semi-analytical solutions for production from hydraulically fractured horizontal wells are available for simple cases, numerical solutions are better able to incorporate the complexity associated with fracture interference with varying fracture length and spacing. Fractured shale reservoirs have been modeled as either single porosity or dual porosity systems. In the single porosity approach, planar fractures with very high permeability values are placed in the reservoir to approximately simulate the behavior of fractures (Fig. 2), resulting in a single flow system with variable properties. In dual porosity models, on the other hand, the reservoir is divided into two separate systems, fractures and matrix, with distinct flow and transport behaviors. The matrix usually contains a large portion of the hydrocarbon with very low permeability and slow mass transfer while the fractures represent the main conduits for transporting the hydrocarbons to the production wells (Li et al., 2011). In this paper, we assume a single porosity model with a very large contrast between the matrix and fracture permeability values. In single porosity approach, each fracture stage is usually modeled using a logarithmic local grid refinement. With grid refinement the simulation time can increase significantly. To keep the computational demand in our optimization approach at a modest level, we do not implement any grid refinement. However, since we use a uniform coarse grid system with a cell dimension of 50 ft to represent a large reservoir, the fracture cells do not provide a realistic fracture dimensionality (aperture). Therefore, to offset the dimensionality effect we reduce the permeability of the Table 3 Reservoir model properties and simulation information. Parameter Phases Simulation time (t ) Matrix porosity (φm ) Matrix permeability (km ) Initial water saturation (swi ) Gridblock dimensions (x × y × z ) Number of gridblocks (nx × ny × nz )

Value Gas–Water 10 years 0.06 0.00015 md 0.3 50 ft  50 ft  50 ft 80  60  6

Number of production Wells (p) Number of fracture stages for each well (n) Fracture height (hf )

2 15 300 ft

Initial pressure (Pi ) Production wells bottomhole pressure constraint (Pwf )

2950 psi 500 psi

Reservoir temperature (T ) Fracture conductivity (kf × wf )

150°F 0.1–20 md-ft

fracture cells in the model to adjust the conductivity of each fracture. For example if the conductivity of a fracture is 5 md-ft, for a cell dimension of 50 ft, the fracture permeability is reduced to 0.1 md. However, we note that uniform coarse gridding with a fracture width of 50 ft does not represent the true duration of possible flow regimes in fractured reservoirs, especially bilinear or linear flow. Furthermore, the substantial decrease of permeability does not capture the rapid pressure drop that occurs at the vicinity of the wellbore where hydraulic fractures are induced. Nonetheless, we adopt the above approach to predict the long-term cumulative gas production by first validating that the resulting long-term gas production predictions are consistent with those optioned from a finely discretized model. Comparison between predicted gas production in a unconventional reservoir from a model with local grid refinement and the coarser model used in this paper reveals insignificant changes (for this comparison 28 fracture stages with 5 md-ft conductivity were considered). It is also important to note that the proposed optimization framework, i.e., the main contribution of this paper, can be applied with a refined model or even an analytical model so long as the production predictions are reliable. The simulations in the next section were performed using Eclipse E300. Fig. 3 shows a reservoir model with two horizontal production wells, each with 15 stages of fracturing. The shale permeability in this reservoir is set to 0.00015 md. In this paper, hydraulic fracturing design parameters are optimized according to the spatial variability in the shale fracability distribution throughout the reservoir. To describe the brittleness map for the reservoir, the method proposed by Grieser and Bray (2007) is used, i.e. Eqs. (5)– (7). Since we did not have access to real sonic log data, hard data for elastic moduli was synthetized along the horizontal well. The data was then used with variogram-based geostatistical modeling (using the Sequential Gaussian Simulation) to generate conditional distributions for elastic moduli maps, i.e. Young's modulus and Poisson's ratio. Fig. 4(a) and (b) depict the distribution of Young's modulus (ranging from 2 × 106 psi to 6 × 106 psi) and Poisson's ratio (varying between 0.13 and 0.45) generated by SGeMS in the reservoir. The resulting Brittleness Index distribution from these maps is shown is Fig. 4(c). The Brittleness Index has a range of 0.0– 0.7 with higher values showing the sweet spots for hydraulic fracturing. We assume that the locations with the Brittleness Index values smaller than 0.28 are extremely ductile and do not result in producible fractures.

Fig. 3. Matrix and planar fracture permeability distributions for the shale gas reservoir.

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

179

framework is general. Assumptions included in the model used for predicting gas production: – Fracture conductivity enhancement due to improved fracturematrix surface area that can result from shear fracturing is not included. – Permeability reduction during production life of reservoir due to increased net-effective-stress or proppant crush is not considered. – Permeability anisotropy of induced hydraulic fractures is not accounted for. – The well is drilled in the direction of maximum stress. – Planar shape induced fractures are used. – For each fracture location, the BI is averaged over the entire vertical section. – A linear relationship is assumed between rock BI and fracture conductivity. 3.2. Optimization of fracture location

Fig. 4. Spatial distribution of (a) Young's modulus, (b) Poisson's ratio, and (c) brittleness index in the reservoir.

The SPSA optimization algorithm is used to update fracture design parameters (e.g. fracture locations and lengths) with the defined brittleness map. At each optimization iteration for updating the fracture length and location, the brittleness distribution is used with Eq. (8) to calculate the new fracture conductivity for each stage. The updated fracture parameters and conductivity are then assigned to the Eclipse simulation software to predict the corresponding production performance. The process is repeated until no further improvements in the objective function are achieved or when the maximum number of iteration is reached. Before proceeding to optimization results, we summarize our modeling assumptions as follows. Note that these modeling assumptions can be altered by users and the optimization

In this section, we consider fracture stage locations (initiation points) as our decision variables and assume that the fracture lengths and the number of fracturing stages are known. The optimization algorithm is initialized with a symmetric distribution of fracturing stages through the lateral length of the well, as shown in Fig. 5(a). A total of 30 fracture stages (15 stages per each well) are considered and a total production time of 10 years is simulated. We note that the optimization time span is defined by the user and can be shorter or longer depending on the objective of the study and/or the specifics of the field. The corresponding pressure distribution for the initial fracture configuration after 10 years of production is shown in Fig. 5(a). The final solution for the fracture locations and the respective pressure distribution after 10 years are shown in Fig. 5(b). It is clear from the solution that the optimization algorithm attempts to avoid ductile regions with lower fracability (blue areas) and concentrates the fracturing stages around the brittle regions with higher fracability (red areas). The pressure distribution after 10 years of production shows propagation of the pressure front toward neighboring fractures, resulting in fracture interference. Fig. 6(a) shows the increase in the NPV function after each accepted iteration. The NPV value for this example increased from an initial value of

Fig. 5. Location of fracture stages in the reservoir: (a) before optimization, and (b) after optimzization. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

180

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

fracture half-lengths. In other words, Eq. (11) is divided into two optimization: one for finding optimal fracture locations ( u^ ) with ^ fixed fracture length, and one for optimizing fracture length L , assuming a fixed fracture location.

u^ = argmin fuϵZ2 (L′, u)

(fracture location optimization)

^ L = argmin fLϵZ1 L u^ , u^

(() )

(fracture length optimiza tion)

Fig. 6. (a) Evolution of the objective function after each accepted iterations; (b) cumulative gas production before (black) and after (red) optimization. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

$3.2594e þ07 to a final value of $3.6870eþ07, representing a 13% increase. Fig. 6(b) compares the cumulative gas production for the initial and optimized solutions, showing an increase from 14,000 MMSCF to 15,751 MMSCF. The increase in the production is mainly due to the effectiveness of the fracture conductivity in more brittle regions. Additionally, the model used in these examples does not account for production decline due to fracture closing with time. It is expected that fracture closing should work in favor of fractures placed in more brittle areas, which can potentially further increase the contribution of the optimization. 3.3. Optimization of fracture half-length We considered the optimized fracture locations in the previous example (Fig. 4) and use the SPSA algorithm to optimize the

(12)

As shown in Li and Jafarpour (2012), these two optimization problems can be coupled by solving them sequential until no further improvement in the objective function can be achieved. Fig. 7 shows the solution with different fracture half-lengths and pressure propagation around the fractures. In this case, the optimization results show shorter fractures in ductile areas and longer fractures in brittle areas. The final solution provides a NPV value of $3.40e þ 07, which is greater than the NPV value of the initial configuration, $3.15e þ 07. To evaluate this outcome, we considered another solution that maximizes the stimulated reservoir volume (SRV) by assigning maximum fracture length at all locations. The resulting NPV value in that case is $3.69e þ 07, which is higher than the optimized solution. The cumulative gas production plots for all three cases are also shown in Fig. 8. The results suggest that having longer fractures leads to more extensive access to the reservoir (regardless of the conductivity) and hence greater production. To verify this hypothesis, we ran four different cases with fracture half-length values of 100 ft, 200 ft, 400 ft and 800 ft. Fig. 9 displays a comparison between the cumulative gas production for each case, clearly showing an increasing production trend as the fracture half-length increases. Based on these results, the solution obtained for fracture halflength optimization in Fig. 6 has to be a local solution. In general, however, several factors can affect the solution of the fracture half-length optimization problem, including the gas price (set at 3.5 $/MSCF in this example), fracturing cost per unit length and as a function of rock geomechanical properties (or reservoir quality map). Unfortunately, the reported cost of generating a fracture stage in the literature does not differentiate between brittle and ductile regions. However, ductility increases the energy required to propagate a fracture and therefore increases the cost of fracturing. Decreasing the gas price and/or increasing the fracture placement cost can result in solutions that favor fracture halflengths that are shorter than the maximum achievable half-length, especially in ductile regions. In other words, maximum stimulation is only an optimal solution if the added revenue from gas production offsets the additional cost of fracturing, which may not be the case in highly ductile regions of a reservoir. Several underlying assumptions in our optimization problem could also have an impact on the obtained solution. Some of these assumptions include the linear cost of hydraulic fracturing as a function of fracture length, adopting independent homogeneous porosity throughout the formation, the linear relation between fracability and fracture conductivity, using a uniform fracture conductivity. In

Fig. 7. Result of the SPSA algorithm for fracture half-length optimization: brittleness map with the optimized fracture half-lengths (left), and pressure distribution after 10 years of production (right).

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

181

(iv) Pseudo radial/elliptical flow, during which the fluid flow towards the whole fracture-well system is radial. (v) Boundary-dominated flow, which happens when the boundaries of the reservoir are reached.

Fig. 8. Cumulative gas production for the initial configuration (black curve), the optimized fracture half-length values (red curve) and the maximum stimulated reservoir volume, that is maximum fracture half-length value, (blue curve) solutions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9. Cumulative gas production for different fracture half-length values.

general, the heterogeneity in fracability can create regions in the reservoir where the maximum fracture length is not the optimal solution. Our ongoing study is focused on developing optimization methods by using more sophisticated geomechanical models that account for fracture dynamics and important geomechanical effects. 3.4. Optimization of number of fracture stages In this section, we consider the number of fracturing stages as an optimization decision variable. The common practice in hydraulic fracturing is to use intensive fracturing to maximize contact with the sweet spots in the shale formation. However, to minimize the fracturing cost, one should investigate fracturing scenarios that can achieve the same productivity with minimum cost. In this paper, the governing flow equations are solved numerically using a commercial simulator. However, there are five known flow regimes in hydraulically-fractured horizontal wells (Luo et al., 2010): (i) Bilinear or linear flow, where fluid flows within the fracture normal to the fracture planes and the type (linear or bilinear) and duration of the flow depends on fracture conductivity and half-length. During this flow regime, each fracture behaves independently. (ii) Early radial/elliptical flow, where fluid flows across the fracture tips. The duration or existence of this regime depends on fracture spacing and half-length. Fractures still behave independently during this flow regime. (iii) Compound formation linear flow (CFL), in which fractures interact and interfere with each other and the system is characterized by linear flow from formation into the fractures.

Fractures that are placed close to each other can interfere and compete with each other, resulting in a reduced productivity. Therefore, in addition to fracture locations and lengths, the optimal number of fracture stages is unknown. Here, we optimize the location and number of fractures with a fixed length. To optimize the number of fractures, we consider a sequential optimization approach where several optimization subproblems are solved. In each subproblem, one or two fracture stages per well are introduced and their optimal locations are identified by solving the fracture placement optimization problem. The number of fracture stages is continuously increased until adding new fracturing stages results in a decrease in the NPV, which usually happens when the cost of fracturing is not offset by the revenue generated from the incremental gas production. In general, the optimization can start by optimally placing one stage of fracturing and sequentially adding a new fracturing stage and optimizing its location. For brevity of our presentation, however, we consider the solution after introducing 10 fracturing stages for each well. Fig. 10(a) shows the symmetric fracture configuration for the first 10 fractures and the pressure profile after 10 years of production. The optimized locations for these fractures are shown in Fig. 10(b). As expected, the fractures are mainly concentrated around highly brittle regions and the pressure has a greater extension. The optimization continues by adding an additional fracture stage and optimizing its location for up to 30 stages. Fig. 11(a) shows the NPV as a function of number of fracture stages. For this example, the maximum NPV is achieved after including 26 stages. The resulting fracture locations for 26 fracturing stages are shown in Fig. 10(c). Fig. 11(b) displays the cumulative gas production for 22, 24, 26, and 28 fracturing stages. From this figure, the added production due to the increase in fracturing stages from 26 to 28 is incremental and does not offset the cost of fracturing. Examination of Fig. 10(c) indicates that any additional fracturing after 26 stages has to be introduced either into the highly ductile areas of the reservoir or in areas with high fracture concentration, resulting in interference with other fractures and smaller additional production. We also performed an optimization in which the locations of the 26 fracture stages were optimized simultaneously. The initial fracture stages were distributed evenly in the reservoir. Fig. 10 (d) shows the 26 fracturing locations after optimization. In this case, the final NPV value is $3.59e7, which is slightly higher than the final NPV ($3.57e7) that was achieved with 26 fracturing stages using the sequential approach. This result is expected as the solution of the sequential approach is constrained by placing the fracturing stages one at a time without any knowledge about the total number of fractures. The simultaneous optimization on the other hand does not include the number of fracturing stages as an optimization variable. Pressure profile of the simultaneous optimization of 26 fractures also confirms this conclusion since pressure disturbance covers a greater area, and results in a higher gas production. 3.5. Factors influencing production from shale plays The concept of shale capacity has recently been proposed to describe production performance as a normalized product of four shale drivers: TOC, natural fracture density, brittleness and porosity (Ouenes, 2014). However, these four drivers are not independent and can have similar impact on production performance of shale development. The brittleness is mostly due to

182

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

Fig. 10. Optimization of fracture locations and number of stages: (a) initial fracture locations for 10 fracturing stages; (b) optimized fracture locatiosn for 10 fracturing stages; (c) sequential optimization of fracture locations for 26 fracturing stages, where the fracture stages are introduced and their loccations are optimized one at a time; (d) simultaneous optimization of fracture locations for 26 fracturing stages.

quartz content of the shale formation, which is also positively correlated with presence of natural fractures. Shale can adsorb significant amount of gas. The Langmuir isotherm can be used to determine the available area for desorption based on the TOC data. Gas desorption from the shale surface, however, requires high pressure drawdown. On the other hand, Barnett shale (modeled in this paper) has a moderate adsorption capacity, SCF VL = 96 ton , PL = 650psi, (Mengal and Wattenbarger, 2011), implying that desorption of gas from matrix for pressures above 1300 psia is minimal. To account for the effect of TOC and brittleness on fracture spacing optimization, a normalized TOC map (varying between 0 and 1) and a map with a similar scale for the brittleness were generated (Fig. 12). It is noteworthy that, in the published literature, the relation between rock brittleness and TOC remains unresolved. Documented evidence exists to show many high-TOC shales that are ductile, immature and non-productive, e.g. Sylvan shale (Chong et al., 2010). Other reports based on surface seismic and well logs in Barnett formation relate high TOC zones to high quartz content,

leading to high brittleness index (Perez Altamar and Marfurt, 2014). Ouenes (2014) also reports an inverse relation between the TOC and brittleness in Marcellus. We present the optimization outcome for various levels of correlations between the brittleness and conductivity. The controlling parameter for this correlation is m. For m¼ 0, a constant conductivity of 0.1 md-ft is assigned to each stage, regardless of the variability in brittleness. As m increases stronger correlations between brittleness index and fracture conductivity are introduced. For m ¼1, stages located in most brittle regions have the largest conductivity (20 md-ft). For each m, two optimization cases are considered, one with constant TOC and one with spatially varying TOC map (Table 4). The results show that as the dependence between fracture conductivity and brittleness increases desorption also increases due to higher drawdown in high conductivity fracture. Porosity in shale formations has two components: (1) interparticle or matrix porosity and (2) open fractures (Bruner and Smosna, 2011). Marcellus shale exhibits large fracture densities with most of its porosity attributed to natural fractures. Therefore for Marcellus shale there is a positive correlation between

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

183

Fig. 13. Spatial distribution of porosity in the reservoir.

Fig. 11. (a) NPV evolution with increasing number of fracture (b) cumulative gas production for different number of fracturing stages.

stages;

Fig. 12. Spatial distribution of TOC in the reservoir.

Table 4 NPV results for different correlations between conductivity and brittleness with/ without TOC.

m ¼0.0 m ¼0.1 m ¼0.2 m ¼0.3 m ¼0.4 m ¼0.5 m ¼0.6 m ¼0.7 m ¼0.8 m ¼0.9 m ¼1.0

NPV (Constant) ×106

NPV (variable TOC) ×106

Difference ×106

 1.6763 2.4727 5.1373 7.2671 9.0511 10.586 11.93 13.13 14.2 15.166 16.045

 1.6495 2.6204 5.4126 7.6844 9.6117 11.291 12.781 14.118 15.3337 16.442 17.455

0.0268 0.1477 0.2753 0.4173 0.5606 0.705 0.851 0.988 1.1337 1.276 1.41

brittleness and porosity. On the other hand for Barnett formation, fracture porosity is a very small percentage of the total shale porosity. The average porosity in productive portions of the formation ranges from 3% to 6%, of which 4.3% is due to thermal decomposition of the TOC (Bruner and Smosna, 2011). Hence,

unlike Marcellus, there is a positive correlation between the TOC and porosity in Barnett shale. To investigate the impact of porosity, a heterogeneous porosity map with the range 1.7–6% was generated (Fig. 13) assuming a positive correlation with TOC variability. Two extreme cases were considered. In the first case, fracture stages were placed in brittle parts of reservoir with low porosity (Fig. 14a) while in the second case fracture stages were placed in high porosity locations with low brittleness (Fig. 14b). Both cases were simulated for 30 years of production. The cumulative gas production and net present value are shown in Fig. 15. As can be seen in this figure, placing the fracture stages in the brittle parts with lower porosity (first case) will result in a higher cumulative production and NPV values after 30 years. For regions with very high porosity, since fracability is not sufficient, no effective fractures are developed. On the other hand, regions with high fracability are usually associated with lower values of porosity (around 1.7%), and hence lower production. In addition to porosity, TOC, fracture density and brittleness, several other factors contribute to the success of shale gas development. One of the first variables considered before drilling the well is the direction of maximum ( σH ) and minimum ( σ h ) horizontal stresses, which are usually determined from wireline logs in pilot hole. In the same geographic area, the maximum and minimum stresses are typically consistent and do not show significant variations. After determining the maximum and minimum stress orientation, the horizontal well is drilled perpendicular to the direction of maximum stress, subjecting the fracture aperture to the minimum stress direction (Beard, 2011). The shape and geometry of the fractures also depends on the anisotropy of the horizontal stress. In this work, it is assumed that the well is already drilled in the direction of maximum stress and fractures are planar. Another important variable in selecting the potential fracture initiation points is the closure stress or the pressure at which a fracture effectively closes. Lower values of closure stress are more desirable since fracture conductivity will decrease with an increase in closure stress (Cipolla et al., 2010). As closure stress increases and Young's Modulus decreases, a dramatic decrease in fracture conductivity is observed (Cipolla et al., 2010). Therefore, areas with lower clay content, and higher Young's Modulus, exhibit lower closure stress, establishing a negative correlation between brittleness and closure stress. A key assumption in this paper is the knowledge of the true brittleness index distribution. In general, from the collected LWD measurements, core analysis, and chemostratigraphy analysis of drill cuttings, the brittleness along the well can be inferred. Fig. 16 shows the hard data extracted from brittleness map in Fig. 4c along the horizontal wells. In general, however, the exact spatial distribution of fracability is unknown. Geostatistical modeling methods can be used to interpolate between the data and develop full property distribution maps. The simplest assumption for brittleness distribution is extension of the same values at the well locations in the two directions perpendicular to the well, as shown in Fig. 17. The proposed optimization framework, however, does not depend on the manner in which the fracability distribution is developed. Fig. 18 shows the optimal hydraulic fracture

184

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

Fig. 14. Investigating the impact of spatial variability in porosity: fracuring (a) low porosity and (b) high porosity regions. In each case, the fracture stages are shown on the left and pressure distribution after 10 years of production are displayed on the right.

Fig. 17. Random field generated from horizontal well observations (Fig. 16).

Fig. 15. Cumulative gas production and NPV for cases I and II.

Fig. 16. Brittleness hard data collected along the horizontal well.

configuration and pressure propagation for brittleness distribution in Fig. 17. As can be seen, the optimal configuration favors inducing the fractures in brittle areas to maximize the revenue from gas

production. Proper interpolation of brittleness beyond the well location requires a model of spatial variability (e.g., a variogram). Since this model is in general unknown, a significant degree of uncertainty is associated with the distribution of brittleness, especially for locations far from the wellbore. Fig. 19 shows brittleness maps generated assuming different spatial continuity models (variograms) based on the same observations along the well (Fig. 16). Along the wellbore, these maps have the same brittleness values; however, as the distance from the wellbore increases so does the dissimilarity between these models. Even for the same variogram model, multiple distinct brittleness maps can be generated that honor the observed data at the well locations (Fig. 20). Hence, the proposed optimization approach can be extended to a robust optimization method to incorporate the uncertainty in the fracability distribution. Finally, the workflow proposed in this paper is shown in Fig. 21 and can be summarized with the following steps: (1) data

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

185

Fig. 18. Optimal fracture configuration (left) and the corresponding pressure distribution (right).

collection, including logging while drilling (LWD), core analysis and chemostratigraphy analysis of drill cuttings; (2) estimating brittleness index along the well using the available data; (3) generating a brittleness distribution map using the spatial correlation of brittleness in the formation (variogram) or simple extension perpendicular to the well direction; 4) implementing the proposed optimization algorithm to optimize hydraulic fracture design parameters (e.g. fracture location, fracture half-length and number of fracture stages). Implementing this optimization workflow is expected to result in shooting less fracture stages without compromising production.

4. Conclusion Brittleness (fracability) index is an important element in commercial shale gas development that determines rock's response to hydraulic fracturing. Despite its established significance, rock brittleness information is not commonly incorporated in hydraulic fracturing design. The current hydraulic fracturing process is implemented by placing several equally spaced fracture stages to contact a greater area of the reservoir without any attempt to optimize the design parameters, including brittleness data collected from logging and chemostratigraphy analysis. In this work, we developed a systematic framework for optimizing the fracture design parameters to improve fracturing efficiency and reduce the associated cost. We applied the SPSA algorithm to optimize hydraulic fracturing design parameters, including the number of fracture stages, fracture stage locations and half-length for spatially variable shale fracability distribution. The proposed approach is similar to using quality maps for well placement optimization in conventional reservoirs. In this case, the distribution of fracability or brittleness index represents the quality map for hydraulic fracturing. Fracability is an important factor in determining desirable fracture locations or sweet spots. Brittle regions of the reservoir are suitable for fracturing as they respond better to hydraulic fracturing

treatments, while ductile regions tend to resist fracturing and heal (close) after hydraulic fracturing. By relating rock brittleness distribution in the reservoir to fracture conductivity and including the cost of fracturing per unit length in the objective function, we formulated and solved an optimization problem to maximize the NPV of the shale asset development. We presented a series of numerical experiments for optimization of fracturing stage locations, fracture half-length, and number of fracturing stages in a shale reservoir with two horizontal wells. The optimization results for fracture placement show that the method could place the fracture stages in the most brittle sections of the reservoir. Given the local nature of the optimization algorithm, however, local solutions could be found to increase the performance of a given initial configuration. We also presented a brief discussion on the factors that can impact the solution to fracture half-length optimization. We applied this optimization in a sequential manner to identify the number of fracturing stages as well as the location of each stage. Here, the optimization was initialized with a small number of fracture stages, followed by adding a new fracturing stage and optimizing its location. The new fracturing stages were added and optimized continuously until addition of a new fracturing stage (with optimized location) resulted in a decrease in NPV, implying that the gas production revenue from adding the new fracture could not offset the fracturing cost. A key assumption in our study is the knowledge of fracability distribution. While the fracability of the formation around the wellbore can be characterized using well logging techniques, its distribution at locations away from the well has to be inferred. Direct and indirect observations of geomechanical reservoir properties may provide additional insight into the fracability of the reservoir. However, even after integrating various sources of information, shale fracability distribution is likely to be known only with a significant level of uncertainty. The conventional approach to handle the uncertainty in the description of rock property distribution, in this case fracability, is an ensemble-based representation. In that case, robust optimization methods can be used to incorporate the expressed uncertainty in rock fracability

Fig. 19. Brittleness index distribution maps generated using different variograms using observed data.

186

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

Fig. 20. Brittleness index maps generated using similar variograms using observed data.

into the hydraulic fracturing design optimization problem, an approach that has been employed in well placement and production optimization of conventional reservoirs. While fracture design and optimization is important for efficient development, an important step is to ensure that the final design from this optimization is indeed implemented. In practice, the planned treatments could be very different from implemented treatment. This difference calls for a mechanism to modify or adapt the design based on the outcome of previous treatments, which would require real-time implementation. Real-time application can be accomplished by combining the proposed method with fracture monitoring and imaging during the hydraulic fracturing process. Currently, fracture monitoring and imaging is performed using microseismic, pressure and tracer data. If advances in imaging technologies allow for real-time data interpretation and fracture characterization, the optimization approach in this paper can be performed with additional information about the previously completed fractures. This information can inform the need for correcting and repeating the treatment process or even adjusting the fracturing process for the next stage given the outcome of previous stages. Several case studies in the literature show how the use of real-time fracture mapping enables for on-

Core Analysis XDR lithology analysis

the-fly changes in fracture design and perforation strategy (Daniels et al., 2007).

Appendix A. Optimization algorithm The discrete optimization problem that was posed in the text must be solved to find optimal values for fracture locations and lengths. To this end, we use the SPSA algorithm (Spall, 1992, 2000, 2003). The SPSA algorithm has also been applied to well placement optimization in conventional reservoirs (Bangerth et al., 2006; Li and Jafarpour 2012; Li et al., 2013). A main advantage of the SPSA algorithm is its stochastic gradient approximation, which is quite efficient for large dimensional problems where analytical gradient calculation is not possible. The method is particularly useful for stochastic optimization and when the objective function cannot be calculated exactly. For the objective function f (x ) ( NPV in our case), where x = [u L ]T is a m-dimensional vector of unknown parameters over the feasible set of fracture locations and lengths, a minimizer x^ is found by satisfying the first-order op∂f (x ) timality condition g (x^ )= ∂x |x = x^ =0, using an appropriate linesearch method. The following standard recursion algorithm can be

Logging While Drilling (LWD) Data: Spectral Gamma ray Neutron activation spectroscopy sonic log

Chemostratigraphy Analysis of drill cuttings

Estimate Brittleness Index along the well

Generate Brittleness Index distribution map (e.g., with interpolation)

Optimize hydraulic fracturing design parameters Fig. 21. Hydraulic fracture design optimization workflow.

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

implemented to search for a solution:

xk +1 = xk −ak g^k ( xk )

(A1)

where g^k (xk ) is a stochastic approximation of the gradient ∇f (xk ) at k th iteration, and ak is the update step size. Gradient approximation can be computed using a finite difference stochastic approximation (FDSA) or the SPSA algorithm. For the finite difference approximation, the i th component of the gradient vector gk (xk ) at the k th iteration can be written as:

f ( xk + ck ei ) − f (xk − ck ei ) g^ki ( xk ) = (i=1, 2, …,m) 2ck

(A2)

in which ei denotes the standard unit vector (with one in its i th component and zeros everywhere else), and ck represents a small positive number. Considering the dimensionality of x , the FDSA approximation of the gradient is implemented one dimension at a time, requiring a total of 2m function evaluations at each iteration, which is computationally inefficient for large problems. When the gradient is approximated using simultaneous perturbations, all the m components of the gradient vector are randomly perturbed at once to obtain two measurements of f (x ), which are then used to compute the gradient approximation as follows:

f ( xk + ck ∆k ) − f (xk − ck ∆k ) g^k ( xk ) = 2ck ∆k

(A3)

The user-generated m-dimensional random perturbation vector ∆k = [∆k1, ∆k2 , … , ∆km ]T contains independent and symmetrically distributed (about 0) members with finite inverse moment E ( ∆ki −1). One particular distribution for ∆k that satisfies the aforementioned conditions is the symmetric Bernoulli distribution with random values ±1. A main advantage of the SPSA over the FDSA is its efficiency in estimating the approximate gradient, especially in large dimensional problems or when the objective function is noisy. To approximate the gradient, the SPSA requires two function calls in each iteration, regardless of the size of m. On the other hand, the FDSA gradient approximation is achieved through 2m cost function evaluations. Thus, the SPSA uses m times fewer function evaluations than the FDSA. However, since the SPSA does not follow the gradient path exactly, its convergence to the solution may require more iterations than is needed with the FDSA algorithm. Nonetheless, it has been shown that because the approximated gradient with the SPSA is an unbiased estimator of the gradient, on average it nearly tracks the true gradient. Hence, if implemented properly, the SPSA algorithm offers a m-fold savings per iteration over the FDSA algorithm. An important point to consider in fracture placement optimization is that the domain of optimization is discrete, implying that the feasible sets Z1 (fracture length in terms of number of cells) and Z2 (fracture locations) is finite and limited to the indices of gridblocks in the model. Although the SPSA algorithm originally was developed for continuous optimization problems, it is generally applicable to both integer and mixed-integer optimization problems. Examples of its application to well placement and general field development optimization can be found in (Bangerth et al., 2006; Li and Jafarpour 2012; Li et al., 2013).

Appending B. Nomenclature

E ϑ BI G F

Young's modulus, psi Poisson's ratio brittleness Index shear modulus, pa shear force, lbf

A ∆x L kf wf km xf FCD u NPV f (x ) C rg rw Qg Qw ∆tk b n p x g^k (xk ) ak ei ck ∆k

187

area, in2 displacement, ft length, ft fracture permeability, md fracture aperture, ft matrix permeability, md fracture half-length, ft dimensionless fracture conductivity fracture location net present value ( $) objective function cost of fracturing ( $/ft ) gas price ($/MSCF) water disposal cost ($/bbl) gas production rate (MSCF/D) water production rate (bbl/day) time-step ( day ) discount rate number of fracture stages number of production wells vector of unknown parameters Stochastic approximation of ∇f (xk ) step size standard unit vector small positive number random perturbation vector

References Al-Ahmadi, H.A., Wattenbarger, R.A., 2011. Triple-porosity Models: one further step towards capturing fractured reservoirs heterogeneity. In: Proceedings of SPE/ DGS Saudi Arabia Section Technical Symposium and Exhibition (Paper 149054MS), 15–18 May 2011, Al-Khobar, Saudi Arabia. Bangerth, W., Klie, H., Wheeler, M.F., Stoffa, P.L., Sen, M.K., 2006. On optimization algorithms for the reservoir oil well placement problem. Comput. Geosci. 10 (3), 303–319. Beard, T., 2011. Fracture Design in Horizontal Shale Wells–Data Gathering to Implementation. Workshop for the Hydraulic Fracturing Study: Well Construction and Operation, 61. Boulis, A., Jayakumar, R., Rai, R., March, 2013. Application of well spacing optimization workflow in various shale gas resources: lessons learned. In: Proceedings of IPTC 2013: International Petroleum Technology Conference (Paper IPTC17150), 26–28 March 2013, Beijing, China. Bruner, K. R., Smosna, R., 2011. A Comparative Study of the Mississippian Barnett Shale, Fort Worth Basin, and Devonian Marcellus Shale, Appalachian Basin. National Energy Technology Laboratory. DOE/NETL-2011/1478. Buller, D., Hughes, S.N., Market, J., Petre, J.E., Spain, D.R., Odumosu, T., January, 2010. Petrophysical evaluation for enhancing hydraulic stimulation in horizontal shale gas. In: Proceedings of SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Chong, K.K., Jaripatke, O.A., Gierser, W.V., Passman, A., January, 2010. A completions roadmap to shale-play development: a review of successful approaches toward shale-play stimulation in the last two decades. In: Proceedings of International Oil and Gas Conference and Exhibition in China. Society of Petroleum Engineers. Cipolla, C., Lolon, E., Erdle, J., Rubin, B., 2010. Reservoir Modeling in Shale-gas Reservoirs. SPE Reserv. Eval. Eng. 13 (4), 638–653. Daniels, J., Waters, G., Le Calvez, J., Bentley, D., Lassek, J., 2007. Contacting more of the barnett shale through an integration of real-time microseismic monitoring, petrophysics, and hydraulic fracture design. In: Proceedings of SPE Annual Technical Conference and Exhibition (Paper 110562-MS), 12–14 October 2007, Anaheim, California. Grieser, W., Bray, J., 2007. Identification of production potential in unconventional reservoirs. In: Proceedings of Production and Operations Symposium (Paper 106623-MS), 31 March–3 April 2007, Oklahoma City, Oklahoma. Hossain, M.M., Rahman, M.K., Rahman, S.S., 2002. A shear dilation stimulation model for production enhancement from naturally fractured reservoirs. SPE J. 7 (02), 183–198. Jarvie, D.M., Hill, R.J., Ruble, T.E., Pollastro, R.M., 2007. Unconventional shale-gas systems: the mississippian barnett shale of north-central texas as one model for thermogenic shale-gas assessment. AAPG Bull. 91 (4), 475–499. Karastathis, A., 2007. Petrophysical Measurement on Tight Gas Shale (Doctoral dissertation). University of Oklahoma, United States. King, G., 2010. Thirty years of gas shale fracturing: what have we learned? In:

188

A. Jahandideh, B. Jafarpour / Journal of Petroleum Science and Engineering 138 (2016) 174–188

Proceedings of SPE Annual Technical Conference and Exhibition (Paper 133456MS), 19–22 September 2010, Florence, Italy. Li, J., Du, C., Zhang, X., 2011. Critical evaluations of shale gas reservoir simulation approaches: single porosity and dual porosity modeling. In: Proceedings of PE Middle East Unconventional Gas Conference and Exhibition (Paper 141756-MS), 31 January–2 February 2011, Muscat, Oman. Li, L., Jafarpour, B., 2012. A variable-control well placement optimization for improved reservoir development. Comput. Geosci. 16 (4), 871–889. Li, L., Jafarpour, B., Mohammad-Khaninezhad, M.R., 2013. A simultaneous perturbation stochastic approximation algorithm for coupled well placement and control optimization under geologic uncertainty. Comput. Geosci. 17 (1), 167–188. Li, Q., Chen, M., Jin, Y., Zhou, Y., Wang, F.P., Zhang, R., March, 2013. Rock mechanical properties of shale gas reservoir and their influences on hydraulic fracture. In: Proceedings of IPTC 2013, International Petroleum Technology Conference. Luo, S., Neal, L., Arulampalam, P., Ciosek, J.M., January, 2010. Flow regime analysis of multi-stage hydraulically-fractured horizontal wells with reciprocal rate derivative function: Bakken case study. In: Proceedings of Canadian Unconventional Resources and International Petroleum Conference. Society of Petroleum Engineers. Ma, X., Plaksina, T., Gildin, E., 2013. Optimization of placement of hydraulic fracture stages in horizontal wells drilled in shale gas reservoirs. In: Proceedings of Unconventional Resources Technology Conference. Society of Petroleum Engineers (Paper 168780-MS), 12–14 August 2013, Denver, Colorado. Mengal, S.A., Wattenbarger, R.A., 1 January 2011. Accounting For Adsorbed Gas in Shale Gas Reservoirs. Society of Petroleum Engineers. doi: 10.2118/141085-MS. Meyer, B.R., Bazan, L.W., Jacot, R.H., Lattibeaudiere, M.G., January, 2010. Optimization of multiple transverse hydraulic fractures in horizontal wellbores. In: Proceedings SPE Unconventional Gas Conference Society of Petroleum Engineers. Ouenes, A., 25 February 2014. Distribution of Well Performances in Shale Reservoirs and Their Predictions Using the Concept of Shale Capacity. Society of Petroleum Engineers. doi: 10.2118/167779-MS. Palmer, I., Cameron, J., Maschovidis, Z., Ponce, J., 2009. Hydraulic fracturing—1: natural fractures influence shear stimulation direction. Oil Gas J. 107, 12. Parker, M., Petre, E., Dreher, D., Buller, D., May, 2009. Haynesville shale: hydraulic fracture stimulation approach. In: Proceedings of International Coalbed & Shale Gas Symposium, Tuscaloosa, Alabama, USA, pp. 18–22. Perez Altamar, R., Marfurt, K., 2014. Mineralogy-based brittleness prediction from surface seismic data: application to the Barnett shale. Interpretation 2 (4), T255–T271. Rickman, R., Mullen, M., Petre, J., Grieser, W., Kundert, D., 2008. A Practical Use of Shale Petrophysics for Stimulation Design Optimization: All Shale Plays are not Clones of the Barnett Shale. In: Proceedings of SPE Annual Technical Conference and Exhibition, 21–24 September 2008, Denver, Colorado. Schweitzer, R., Bilgesu, H.I., 2009. The Role of Economics on Well and Fracture Design Completions of Marcellus Shale Wells. SPE Eastern Regional Meeting,

23–25 September 2009, Charleston, West Virginia. Smith, C.H., Hamilton, L., Bashkirtseva, A., 2013. Three-dimensional sonic and anisotropy data to provide effective fracture treatments. In: Proceedings of SPE European Formation Damage Conference & Exhibition. Society of Petroleum Engineers. Sone, H., 2012. Mechanical Properties of Shale Gas Reservoir Rocks, and Its Relation to the In-situ Stress Variation Observed in Shale Gas Reservoirs (Doctoral dissertation). Stanford University, United States. Spall, J.C., 1992. Multivariate stochastic approximation using a simultaneous perturbation gradient approximation. IEEE Trans. Autom. Control 37, 332–341. Spall, J.C., 2000. Adaptive stochastic approximation by the simultaneous perturbation method. IEEE Trans. Autom. Control 45, 1839–1853. Spall, J.C., 2003. Introduction to Stochastic Search and Optimization: Estimation, Simulation and Control. Wiley, New Jersey. Wang, F.P., Gale, J.F., 2009. Screening criteria for shale-gas systems. Gulf Coast Assoc. Geol. Soc. Trans. 59, 779–793. Wilson, K., Durlofsky, L., 2012. Computational Optimization of Shale Resource Development Using Reduced-physics Surrogate Models. SPE Western Regional Meeting, 19–23 March 2012, Bakersfield, California. Yang, S., Harris, N., Dong, T., Wu, W., Chen, Z., 2015. Mechanical properties and natural fractures in a Horn River shale core from well logs and hardness measurements. In: Proceedings of EUROPEC 2015. Society of Petroleum Engineers. Yu, W., Sepehrnoori, K., 2013. Optimization of multiple hydraulically fractured horizontal wells in unconventional gas reservoirs. In: Proceedings of SPE Production and Operations Symposium, 23–28 March 2013, Oklahoma City, Oklahoma.

Further reading Biswas, D., October, 2011. Shale Gas Predictive Model (SGPM)—an alternate approach to model shale gas production. Paper 148491-MS Presented at SPE Eastern Regional Meeting, 17–19 August 2011, Columbus, Ohio. Enderlin, M., Alsleben, H., Beyer, J.A. April, 2011. Predicting fracability in shale reservoirs. In AAPG Annual Convention and Exhibition, Houston, Texas, USA, pp. 10–13. Holt, S., 2011. Numerical Optimization of Hydraulic Fracture Stage Placement in a Gas Shale Reservoir (MSc Thesis). TU Deflt University, Netherlands. Spall, J.C., 1998. An overview of the simultaneous perturbation method for efficient optimization. Johns Hopkins APL Tech. Dig. 19 (4), 482–492. Wong, G.K., Fors, R.R., Casassa, J.S., Hite, R.H., Shlyapobersky, J., January 1993. Design execution and evaluation of frac and pack (F&P) treatments in unconsolidated sand formations in the Gulf of Mexico. In: Proceedings of SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.