Shale Fracturing Characterization and Optimization by ... - OnePetro
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
Shale Fracturing Characterization and Optimization by Using Anisotropic Acoustic Interpretation, 3D Fracture Modeling and Neural Network Deepak Gokaraju, Ming Gu, Dingding Chen, Mehdi E. Far, and John Quirein, Halliburton Society of Petrophysicists and Well Log Analysts Copyright 2015, held jointly by the Society of Petrophysicists and Well Log Analysts (SPWLA) and the submitting authors. This paper was prepared for presentation at the SPWLA 56th Annual Logging Symposium held in Long Beach, CA, USA, July 18-22, 2015.
ABSTRACT Multiple fractures or an extensive fracture network is critical for commercially viable production from low permeability formations, such as shales. Mechanical anisotropy is inherent in shales because of its platy nature. This inherent anisotropy makes fracture prediction in shales more complex, and traditional methods to predict fracture geometry assuming isotropy frequently prove to be inadequate. Current analytical methods boldly assume a constant fracture height and constant mechanical properties for the entire height. Common 3D fracture modeling software are based on isotropic rock models, and models that take anisotropy into account are computationally expensive and time consuming, especially when numerous simulations must be performed by varying the input parameters for parametric study. This paper proposes a workflow to improve the prediction of fracture geometry in anisotropic formations. The workflow involves generating a neural network by using a limited number of 3D fracture modeling cases. After the neural network is obtained from a pilot or offset well, it can be easily embedded into software for optimizing fracture design, identifying geologic sweet spots, and predicting fracture propagation and correlating the results to other horizontal or vertical wells in the same geological area. This process can be divided into three steps. First, the anisotropic models are used to predict horizontal and vertical Young’s modulus (Ehorz and Evert), Poisson’s ratio (horz and vert), and anisotropic minimum horizontal stress (σhmin_ani) from sonic and density log measurements. Second, the elastic moduli properties and σhmin_ani are entered into a 3D fracture modeling simulator to run different cases by varying the completion input parameters. The outputs of the fracture simulator (i.e., the fracture length, height, width, and effective length) serve as a training database to the neural network. In the final step, a neural network is generated based on the training database. After the reservoir-specific neural network is developed, fracture geometry can be predicted or optimized for numerous combinations of completion input parameters in a timely
and cost effective manner. Because the commonly available commercial fracture modeling software assumes isotropy, a new method is presented in this paper to represent mechanical property anisotropy using equivalent Young’s modulus (Eeq) and Poisson’s ratio (eq). Eeq and eq are derived from Ehorz, Evert, horz, and vert and the isotropic (Sneddon and Berry 1958) and anisotropic (Chertov 2012) width functions. This workflow is demonstrated by generating a neural network for two reservoirs using anisotropic elastic moduli as predicted by the dipole sonic log. The fracture geometry predicted by the neural network is compared with the conventional method, assuming the isotropic shale rock. The results show that by assuming an isotropic model the fracture width is overestimated, and the fracture containment and propped length are underestimated. The anisotropic neural network model is further run in a large parametric study to demonstrate how the effective length varies with perforation position, injection volume, and injection rate. The results helped to optimize perforation depth, injection rate, and pumped volume. INTRODUCTION The combination of hydraulic fracturing and horizontal drilling has made production from unconventional reservoirs, such as shales and tight formations possible. However, because of the laminated and platy nature intrinsic to shales, isotropic models built from logs cannot fully describe them. The isotropic model assumes homogeneity across the entire formation and computes a single Young’s modulus and a single Poisson’s ratio from sonic and density measurements. Schoenberg et al. (1996) proposed a simplified anisotropic model for shales (ANNIE model). The ANNIE model computes the stiffness tensor for a transversely isotropic model by assuming Thomsen’s third parameter () to be 0. Another assumption the ANNIE model makes based on observation is in many shales, C12= C13 to compute different elastic moduli in different directions. Quirein et al. (2014) observed the ANNIE model always horz, and the stress predicted using the predicts vert ANNIE model is the same as the isotropic model in absence of any field data for calibration. Quirein et al. (2014) 1
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
proposed the modified ANNIE model allows for cases where vert horz and also predicts the anisotropic closure stress is higher than the closure stress estimated from the isotropic model.
pump time can be optimized using the effective propped area, which greatly influences both short-term and longterm production rates. ELASTIC MODULI FROM ANISOTROPIC MODELS
To predict fracture geometry in shales, it is important they are characterized accurately. Hence, accounting for anisotropy in elastic moduli is equally important as using closure stress computed from an anisotropic model. Most commercial fracture simulation software used today do not account for anisotropy in elastic moduli. Fully 3D finite element models that predict fracture geometry do take anisotropy in elastic moduli into account but are computationally expensive, time consuming and not a practical solution.
To compute elastic moduli (Young’s modulus and Poisson’s ratio), the stiffness tensor must be calculated first. For a transverse isotropic material, five independent stiffness coefficients are required to build the complete stiffness tensor. For a VTI medium (transverse isotropic with a vertical axis of symmetry), they are C11, C33, C44, C12, and C13. The dipole sonic tool in a vertical well in a VTI medium can measure the following velocities: vertical Pwave (C33), two vertically propagating horizontally polarized S-waves (C44 and C55), and the Stoneley wave which is used to derive the horizontally propagating, horizontally polarized S-wave (C66) (Norris et al. 1993).
Hydraulic fracture optimization is tricky and reservoir dependent as very few parameters are under human control. Stress barriers, fracture orientation, and elastic properties of the formation are significant in hydraulic fracture propagation, which cannot be altered or controlled. The parameters under human control are secondary parameters, such as fluid viscosity, pump rate, pump volume, and perforation depth. These parameters have to be chosen carefully on the basis of reservoir prevailing conditions and the targets set, such as fracturing fluid viscosity and the foam quality in fracturing fluids. Not only do they impact fracture geometry but also impact fracture conductivity and effective propped length when combined with proppant density and quantity used. Low viscosity fluids with small proppant injection proved successful in the Barnett Shale but could not reproduce the same amount of success in other formations. The effect of foam quality on effective propped length and success of a hydraulic fracture design has been studied in detail by various authors (Gu et al. 2014; Sani et al. 2001; Powell et al. 1999). The effect of proppant size and density in combination with fluids is not covered in this paper but their impact is described in detail by Gu and Kulkarni (2014), Gu and Dao (2015), and Phatak et al. (2013). To fully understand the impact on fracturing design and fractured well productivity, one needs to conduct a massive parametric study.
Because these four velocities are not sufficient to build the complete stiffness tensor, additional boundary conditions are imposed using the ANNIE model (Shoenberg et al. 1996; Higgins et al. 2008; Waters et al. 2011). C11, C13, and C33 can be estimated by setting the Thomsen parameter to zero. C13 + 2*C44 – C33 = 0 (which gives one C13)
(1)
and the second constraint is that, in many shales, C13 = C12; and therefore, C13 = C12 = C11 – 2*C66
(2)
and combining Eqs. 1 and 2, C66 – C44 = (C11 – C33)/2.
(3)
The stiffness tensor can also be computed by imposing a modified set of boundary conditions using the modified ANNIE model. To account for cases where vert horz, the boundary conditions in Eqs. 2 and 3 have been modified as follows:
This paper proposes a workflow to compute a single Young’s modulus (Eeq) and a single Poisson’s ratio (eq) that take anisotropy in elastic moduli into account and can be used as inputs into the fracture simulation software. This paper also demonstrates the proposed methodology to optimize hydraulic fracturing parameters in unconventional reservoirs more efficiently by combining log interpretation, fracture modeling, neural networks, and a parametric study on two different formations. Optimized solutions for perforation depth, fracturing fluid pump rate, and fluid
C13 = kC12 = k(C11 – 2*C66)
(4)
C11 = k’ (2(C66 – C44)) + C33
(5)
k and k’ are determined using field and core data. Once the stiffness tensor is completed using an appropriate anisotropic model, horizontal and vertical dynamic elastic properties are computed using Eqs. 6–9. 2
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
(6) .
(7)
.
.
(9)
LINEAR ELASTIC FRACTURE MODELS Sneddon and Berry (1958) estimated the width of an elliptical crack in an isotropic medium as shown in Eq. 10. ,
,
(10)
In Eq. 10, w is the maximum fracture width in the center of the elliptical cross section, h is fracture height, σnet is net pressure (fracture pressure minus closure stress), E and v are the isotropic Young’s modulus and Poisson’s ratio. σnet is the distributed pressure inside a fracture and is a function of the distance from the wellbore to the fracture tip and the pumping time. Eq. 10 estimates the width of the fracture at the corresponding location based on σnet. Chertov (2012) proposed a similar equation to estimate fracture width in a transverse isotropic formation using both horizontal and vertical elastic moduli as shown in Eq. 11. ,
√
√
,
1
0, (11)
1
2
(16)
(18)
,
.
(19)
In the equation system above, f(E,v) is the elastic property or depending on the type term and can either be of formation of interest, q0 is the injection rate, μ is the fracturing fluid viscosity, h is the fracture height, t is the injection time, and cw, cl are unit conversion coefficients. The 2D analytical models are simple and straightforward to understand the impact of elastic properties and completion parameters on fracture geometries. A major drawback for such models is they assume a fixed fracture height and constant elastic properties along the height, which is not generally valid for laminated shale reservoirs.
(13) 1
.
(17)
,
(12)
1
√
,
0,
where 1
√
Multiple fracture models have been developed using Eq. 10, which can be extended for anisotropic rocks using Eqs. 1116. Two of the most widely accepted fracture models are the Perkins-Kern-Nordgren (PKN) model (Perkins and Kern 1961; Nordgren 1972) and the Geertsma-de Klerk (GDK) model (Khristianovitch and Zheltov 1955; Geertsma and de Klerk 1969). The PKN and GDK models are similar except for the differences in their basic assumptions. The PKN model assumes an elliptic cross section with a fixed fracture height and is mainly used when the fracture length is much greater than the fracture height whereas the GDK model assumes a rectangular cross section and is mainly used when the fracture height is comparable to or more than the fracture length. These models are developed by combining Eq. 10 with the mass balance and fluid flow physics. The estimated fracture geometry equations given by PKN without considering leakoff are given below:
(8) .
(15)
(14)
where Eh and Ev are horizontal and vertical Young’s modulus, vh and vv are horizontal and vertical Poisson’s ratio, and Gvh is the shear modulus in the x-z plane. The elastic response of the rock in the above equations can be represented using a single term f for simplicity. f for isotropic and anisotropic rocks is given as follows:
Commercial fracture simulation softwares make use of coupled equations and numerical simulations to simulate fracture propagation with fluid and proppant flow in a 3D space (Gu and Leung 1993; Sousa et al. 1993). Some widely used simulators do not take the anisotropy in elastic properties into account although they do take complex stress layering into account. Elastic anisotropy can be accounted 3
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
for during hydraulic fracturing using fully 3D finite element solvers (Gokaraju and Eckert 2014). However, they are time consuming and do not provide a practical solution especially when multiple cases must be run to pick the most effective treatment parameters.
static from dynamic using core measurements. Closure pressure can then be estimated using the elastic moduli as shown below.
H min
THE NEURAL NETWORK APPROACH
p p
An artificial neural network is a powerful tool in recognizing the relation between different parameters and outputs, which can be used to develop an optimized solution. The basic structure of a neural network consists of an input layer, an output layer, and multiple hidden layers. The neural network is usually trained and validated with multiple runs with known outputs; thereby, the model structure and coefficients can be automatically adjusted to optimize the predictions on existing runs and generalize well on the new data for future applications.
E h v ( v (1 ) p p ) E v (1 h ) Eh
1 h
2
H min
E h h 1 h
2
H max
(20)
where,
Hmin is the minimum horizontal stress (psi), v is the overburden stress (psi), pp is the pore pressure (psi), is the Biot’s elastic constant, Hmin is the minimum horizontal strain, Hmax is the maximum horizontal strain, and is the poroelastic constant. f(E,v) is then calculated from the elastic moduli using Eq. 15 or Eq. 16, depending on the formation.
This paper proposes a methodology in which a reservoirspecific neural network is built and trained to quickly and accurately predict fracture geometry, effective propped area for any type of formation, and any given completion parameters (e.g., perforation position, injection rate, injection time, etc.). This can be extended to run a parametric study which helps in selecting optimal hydraulic fracturing parameters to produce the largest EPA (shortterm production) or stimulated reservoir volume (SRV) (long-term production).
Because most simulation software uses isotropic elastic moduli, an equivalent Young’s modulus (Eeq) and Poisson’s ratio (eq) are calculated to best simulate anisotropic elastic moduli using a single E and Eeq can be calculated by the following equation which is a modified and weighted version of the equation used by Chertov (2012).
2
1
(21)
where vvh is the arithmetic averaging for vv and vh and ai is the weight coefficient from 0 to 1 (ah+av+avh=1). Using the calculated f(E,v) and Eeq, an equivalent veq can be calculated by Eq. 21.
1
,
,
, 2
(22)
In the second step, the rock mechanic properties (Eeq, veq) and the closure stress σh from step 1 are input into a 3D fracture modeling simulator. The completion parameters, such as slurry injection rate (qinj), total slurry volume (Qtol), and the perforation depth (TVDperf) are varied to generate the training database with the output results, such as fracture length (Lf), height (Hf), width (wf), and effective propped length (Leff). The training database is then used to produce the neural network. To improve the computation efficiency, multiple modeling cases are run by varying the fracture parameters. To generate an initial training database, each completion parameter is varied by x values equally distributed within the range of interest. The x value is dependent on the level of accuracy sought. For the two
Fig 1. Complete workflow for optimizing completion parameters using sonic measurements.
Fig. 1 illustrates the workflow of generating and training a neural network from sonic measurements and extended to optimize completion parameters. Once the elastic moduli of a rock are computed, they are calibrated and converted to 4
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
organic shales used in this study, x was chosen to be 3-5 (i.e., each parameter was varied up to 5 times). Therefore, if there are n parameters, the total number of training sample combinations is xn. In this study, there are three input completion parameters: slurry injection rate, total injection volume, and perforation position along the horizontal well. The horizontal position is correlated with depth based on geosteering data. Five output results are fracture length, fracture height, fracture width, upper and lower depth of fracture, and effective propped length.
the sonic velocities, Young’s modulus, and closure pressure. The interval can be divided into five zones based on closure stress. Fig. 3 shows a design of a toe-up horizontal well inside the pay. The toe-up strategy has been used for both cases. Fracture simulations will be performed on the lateral part of the shown hypothetical well. Shale 1 is a transversely isotropic media with vertical axis of symmetry (TIV) with horizontal Poisson’s ratio greater than its vertical counterpart. Hence, the modified ANNIE method was used to compute the stiffness tensor from sonic velocities. f(E,) and closure stress are computed using Eqs. 16 and 20 respectively. A 3D isotropic simulator was used to run the simulation models and Eeq and eq were calculated using Eqs. 21 and 22 respectively.
Once the neural network is generated, a group of testing data is run to check the relative error of the outputs for each input parameter. For any input parameter, if the tolerance relative error is not met between Node_i and Node _i+1, one more data point Node_i+1/2 is added in between. The extra cases regarding the added data point are run in fracture modeling to update the current training database, and further update the neural network. Testing is run on the new neural network. If the tolerance error is met for all input parameters, the neural network is left unaltered. Otherwise, it updates itself using the added database until the criterion is met. Once the neural network is obtained, one can either predict fracture geometry and location based on the input of arbitrary completion values, or determine the optimized fracturing design by conducting a massive parametric study with the neural network. For fracturing optimization, the effective propped length (EPL) is one of the best candidates to be the optimization target among all the predicted outputs. It is the propped length within the payzone occupied by infinite relative conductivity (Fcd>50). EPL dominates the short-term production and affects the longterm production. For better optimization, a critical conductivity can be used instead of infinite relative conductivity to define the effective propped length. The critical conductivity is defined as the minimum conductivity needed for fully stimulating a certain propped length during a certain production time. It is a function of propped length, production time, matrix permeability, natural fracture properties, oil API, and other completion and production parameters. (Gu et al. 2014).
Fig. 2. Uranium, sonic, lithology, elastic moduli, and stress data for Shale 1.
STUDY OF FIELD CASES The proposed methodology has been applied to two organic rich shales to develop an optimized EPL and completion parameters. In case 1, the stress barriers are highly defined with the stress difference between the payzone and the bounding layer being as high as 3,000 psi. Fig. 2 is a plot of 5
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
Simulations were performed for a single fracture in Shale 2 that acted as a training database for the generation of a neural network and an optimized solution. Tables 1 and 2 show the variation used for parameter inputs while generating a training database for Shales 1 and 2 respectively. X is the middle of the payzone.
Fig. 3. Toe-up horizontal well inside lower payzone.
For Shale 2, the stress barrier is not as well defined as case 1. Fig. 4 shows uranium, elastic moduli, lithology, and closure pressures. Because of its complex closure stress profile, the interval has been divided into 20 zones, two of which are hydrocarbon-bearing reservoirs.
Parameter
Min value
Max value
Step size
Perforation depth (ft)
X-50
X+50
20
Fluid injection rate (bbl/min)
30
50
10
Fluid injection Volume (M-gal)
10
30
5
Table 1. Input parameter variation for Shale 1.
Parameter
Min value
Max value
Step size
Perforation depth (ft)
X-53
X+72
20
Fluid injection rate (bbl/min)
15
30
5
Fluid injection Volume (M-gal)
10
40
5
Table 2. Input parameter variation for Shale 2.
Fig. 4. Uranium, lithology, elastic moduli, and stress data for Shale 2.
6
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
Shale 2. Fig. 6 is an example showing the target values of Leff and values predicted by the neural network in Shale 1. The solid dots indicate target values, while the hollow dots are predictions. The red ones are 70 training data points, while the green ones are seven random testing cases.
RESULTS AND DISCUSSIONS Fig. 5 shows a typical fracture profile result. This example was run on Shale 1.
Fig. 5. Typical fracture profile plot using fracture design and analysis
Fig. 6. Comparison of target values and values predicted by the neural
software.
network.
The track on the left shows the fracture width-height cross section, while the track on the right shows the length-height cross section. The color shaded contour represents the conductivity distribution. Only the conductivity distribution within the payzone (boundaries marked using yellow colored lines) governs the fracture productivity. The fracture conductivity decreases with increasing propped length as shown by the solid purple line. In the same plot, the minimum conductivity required to fully stimulate different propped length is indicated by red lines. The red lines represent the 1-year, 5-year, and 10-year production using the minimum conductivity criteria (Gu et al. 2014). According to the red lines, as fracture length increases, the minimum conductivity required to fully stimulate the length increases. For the same propped length, the minimum conductivity increases with decreasing production time. Therefore, the intersection between the purple and red lines represents the effective length for the certain production time. For example, for 1-year shot-term optimization, to meet the infinite conductivity criteria, the purple line should always be above the red-dotted line. Their cross is the maximum fracture length with the infinite conductivity. Through this way, the effective length at different production time scale is determined.
Once a neural network is generated for a specific reservoir, it can be used to predict fracture geometry and fracture geometry-based production performance for any combination of input parameters, which in this case are perforation depth, fluid injection rate, and fluid injection volume. Because the anisotropic elastic parameters and complex anisotropic closure pressure profile were both accounted for, the results can be generated instantaneously with a high degree of accuracy. In Fig. 7, this is demonstrated by randomly selecting the slurry injection rate and the total injection volume in Shale 1 to be 48 bbl/min and 23 Mgal. Fracture length, height, width, TVD of the upper and lower fracture boundaries, and effective propped length are all calculated as a function of horizontal well distance. Similarly, Fig. 8 shows fracture length, height, width, and effective propped length are all calculated as a function of depth in Shale 2 using slurry injection rate and the total injection volume to be 30 bbl/min and 40 Mgal respectively.
A neural network is generated using a training database consisting of multiple cases. The different cases are obtained by varying injection rate by three values, injection volume by five values, and perforation position by four values in Shale 1 and by varying injection rate, injection volume and perforation position each by five values in 7
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
position. Leff is 400-460 ft for the first half of the horizontal well, while 360-400 ft for the second half. Injection rate does not have much effect on the effective length. For positions from 0.2 to 1, lower injection rate produces a little longer length. Within the well tip part (
Shale Fracturing Characterization and Optimization by Using Anisotropic Acoustic Interpretation, 3D Fracture Modeling and Neural Network Deepak Gokaraju, Ming Gu, Dingding Chen, Mehdi E. Far, and John Quirein, Halliburton Society of Petrophysicists and Well Log Analysts Copyright 2015, held jointly by the Society of Petrophysicists and Well Log Analysts (SPWLA) and the submitting authors. This paper was prepared for presentation at the SPWLA 56th Annual Logging Symposium held in Long Beach, CA, USA, July 18-22, 2015.
ABSTRACT Multiple fractures or an extensive fracture network is critical for commercially viable production from low permeability formations, such as shales. Mechanical anisotropy is inherent in shales because of its platy nature. This inherent anisotropy makes fracture prediction in shales more complex, and traditional methods to predict fracture geometry assuming isotropy frequently prove to be inadequate. Current analytical methods boldly assume a constant fracture height and constant mechanical properties for the entire height. Common 3D fracture modeling software are based on isotropic rock models, and models that take anisotropy into account are computationally expensive and time consuming, especially when numerous simulations must be performed by varying the input parameters for parametric study. This paper proposes a workflow to improve the prediction of fracture geometry in anisotropic formations. The workflow involves generating a neural network by using a limited number of 3D fracture modeling cases. After the neural network is obtained from a pilot or offset well, it can be easily embedded into software for optimizing fracture design, identifying geologic sweet spots, and predicting fracture propagation and correlating the results to other horizontal or vertical wells in the same geological area. This process can be divided into three steps. First, the anisotropic models are used to predict horizontal and vertical Young’s modulus (Ehorz and Evert), Poisson’s ratio (horz and vert), and anisotropic minimum horizontal stress (σhmin_ani) from sonic and density log measurements. Second, the elastic moduli properties and σhmin_ani are entered into a 3D fracture modeling simulator to run different cases by varying the completion input parameters. The outputs of the fracture simulator (i.e., the fracture length, height, width, and effective length) serve as a training database to the neural network. In the final step, a neural network is generated based on the training database. After the reservoir-specific neural network is developed, fracture geometry can be predicted or optimized for numerous combinations of completion input parameters in a timely
and cost effective manner. Because the commonly available commercial fracture modeling software assumes isotropy, a new method is presented in this paper to represent mechanical property anisotropy using equivalent Young’s modulus (Eeq) and Poisson’s ratio (eq). Eeq and eq are derived from Ehorz, Evert, horz, and vert and the isotropic (Sneddon and Berry 1958) and anisotropic (Chertov 2012) width functions. This workflow is demonstrated by generating a neural network for two reservoirs using anisotropic elastic moduli as predicted by the dipole sonic log. The fracture geometry predicted by the neural network is compared with the conventional method, assuming the isotropic shale rock. The results show that by assuming an isotropic model the fracture width is overestimated, and the fracture containment and propped length are underestimated. The anisotropic neural network model is further run in a large parametric study to demonstrate how the effective length varies with perforation position, injection volume, and injection rate. The results helped to optimize perforation depth, injection rate, and pumped volume. INTRODUCTION The combination of hydraulic fracturing and horizontal drilling has made production from unconventional reservoirs, such as shales and tight formations possible. However, because of the laminated and platy nature intrinsic to shales, isotropic models built from logs cannot fully describe them. The isotropic model assumes homogeneity across the entire formation and computes a single Young’s modulus and a single Poisson’s ratio from sonic and density measurements. Schoenberg et al. (1996) proposed a simplified anisotropic model for shales (ANNIE model). The ANNIE model computes the stiffness tensor for a transversely isotropic model by assuming Thomsen’s third parameter () to be 0. Another assumption the ANNIE model makes based on observation is in many shales, C12= C13 to compute different elastic moduli in different directions. Quirein et al. (2014) observed the ANNIE model always horz, and the stress predicted using the predicts vert ANNIE model is the same as the isotropic model in absence of any field data for calibration. Quirein et al. (2014) 1
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
proposed the modified ANNIE model allows for cases where vert horz and also predicts the anisotropic closure stress is higher than the closure stress estimated from the isotropic model.
pump time can be optimized using the effective propped area, which greatly influences both short-term and longterm production rates. ELASTIC MODULI FROM ANISOTROPIC MODELS
To predict fracture geometry in shales, it is important they are characterized accurately. Hence, accounting for anisotropy in elastic moduli is equally important as using closure stress computed from an anisotropic model. Most commercial fracture simulation software used today do not account for anisotropy in elastic moduli. Fully 3D finite element models that predict fracture geometry do take anisotropy in elastic moduli into account but are computationally expensive, time consuming and not a practical solution.
To compute elastic moduli (Young’s modulus and Poisson’s ratio), the stiffness tensor must be calculated first. For a transverse isotropic material, five independent stiffness coefficients are required to build the complete stiffness tensor. For a VTI medium (transverse isotropic with a vertical axis of symmetry), they are C11, C33, C44, C12, and C13. The dipole sonic tool in a vertical well in a VTI medium can measure the following velocities: vertical Pwave (C33), two vertically propagating horizontally polarized S-waves (C44 and C55), and the Stoneley wave which is used to derive the horizontally propagating, horizontally polarized S-wave (C66) (Norris et al. 1993).
Hydraulic fracture optimization is tricky and reservoir dependent as very few parameters are under human control. Stress barriers, fracture orientation, and elastic properties of the formation are significant in hydraulic fracture propagation, which cannot be altered or controlled. The parameters under human control are secondary parameters, such as fluid viscosity, pump rate, pump volume, and perforation depth. These parameters have to be chosen carefully on the basis of reservoir prevailing conditions and the targets set, such as fracturing fluid viscosity and the foam quality in fracturing fluids. Not only do they impact fracture geometry but also impact fracture conductivity and effective propped length when combined with proppant density and quantity used. Low viscosity fluids with small proppant injection proved successful in the Barnett Shale but could not reproduce the same amount of success in other formations. The effect of foam quality on effective propped length and success of a hydraulic fracture design has been studied in detail by various authors (Gu et al. 2014; Sani et al. 2001; Powell et al. 1999). The effect of proppant size and density in combination with fluids is not covered in this paper but their impact is described in detail by Gu and Kulkarni (2014), Gu and Dao (2015), and Phatak et al. (2013). To fully understand the impact on fracturing design and fractured well productivity, one needs to conduct a massive parametric study.
Because these four velocities are not sufficient to build the complete stiffness tensor, additional boundary conditions are imposed using the ANNIE model (Shoenberg et al. 1996; Higgins et al. 2008; Waters et al. 2011). C11, C13, and C33 can be estimated by setting the Thomsen parameter to zero. C13 + 2*C44 – C33 = 0 (which gives one C13)
(1)
and the second constraint is that, in many shales, C13 = C12; and therefore, C13 = C12 = C11 – 2*C66
(2)
and combining Eqs. 1 and 2, C66 – C44 = (C11 – C33)/2.
(3)
The stiffness tensor can also be computed by imposing a modified set of boundary conditions using the modified ANNIE model. To account for cases where vert horz, the boundary conditions in Eqs. 2 and 3 have been modified as follows:
This paper proposes a workflow to compute a single Young’s modulus (Eeq) and a single Poisson’s ratio (eq) that take anisotropy in elastic moduli into account and can be used as inputs into the fracture simulation software. This paper also demonstrates the proposed methodology to optimize hydraulic fracturing parameters in unconventional reservoirs more efficiently by combining log interpretation, fracture modeling, neural networks, and a parametric study on two different formations. Optimized solutions for perforation depth, fracturing fluid pump rate, and fluid
C13 = kC12 = k(C11 – 2*C66)
(4)
C11 = k’ (2(C66 – C44)) + C33
(5)
k and k’ are determined using field and core data. Once the stiffness tensor is completed using an appropriate anisotropic model, horizontal and vertical dynamic elastic properties are computed using Eqs. 6–9. 2
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
(6) .
(7)
.
.
(9)
LINEAR ELASTIC FRACTURE MODELS Sneddon and Berry (1958) estimated the width of an elliptical crack in an isotropic medium as shown in Eq. 10. ,
,
(10)
In Eq. 10, w is the maximum fracture width in the center of the elliptical cross section, h is fracture height, σnet is net pressure (fracture pressure minus closure stress), E and v are the isotropic Young’s modulus and Poisson’s ratio. σnet is the distributed pressure inside a fracture and is a function of the distance from the wellbore to the fracture tip and the pumping time. Eq. 10 estimates the width of the fracture at the corresponding location based on σnet. Chertov (2012) proposed a similar equation to estimate fracture width in a transverse isotropic formation using both horizontal and vertical elastic moduli as shown in Eq. 11. ,
√
√
,
1
0, (11)
1
2
(16)
(18)
,
.
(19)
In the equation system above, f(E,v) is the elastic property or depending on the type term and can either be of formation of interest, q0 is the injection rate, μ is the fracturing fluid viscosity, h is the fracture height, t is the injection time, and cw, cl are unit conversion coefficients. The 2D analytical models are simple and straightforward to understand the impact of elastic properties and completion parameters on fracture geometries. A major drawback for such models is they assume a fixed fracture height and constant elastic properties along the height, which is not generally valid for laminated shale reservoirs.
(13) 1
.
(17)
,
(12)
1
√
,
0,
where 1
√
Multiple fracture models have been developed using Eq. 10, which can be extended for anisotropic rocks using Eqs. 1116. Two of the most widely accepted fracture models are the Perkins-Kern-Nordgren (PKN) model (Perkins and Kern 1961; Nordgren 1972) and the Geertsma-de Klerk (GDK) model (Khristianovitch and Zheltov 1955; Geertsma and de Klerk 1969). The PKN and GDK models are similar except for the differences in their basic assumptions. The PKN model assumes an elliptic cross section with a fixed fracture height and is mainly used when the fracture length is much greater than the fracture height whereas the GDK model assumes a rectangular cross section and is mainly used when the fracture height is comparable to or more than the fracture length. These models are developed by combining Eq. 10 with the mass balance and fluid flow physics. The estimated fracture geometry equations given by PKN without considering leakoff are given below:
(8) .
(15)
(14)
where Eh and Ev are horizontal and vertical Young’s modulus, vh and vv are horizontal and vertical Poisson’s ratio, and Gvh is the shear modulus in the x-z plane. The elastic response of the rock in the above equations can be represented using a single term f for simplicity. f for isotropic and anisotropic rocks is given as follows:
Commercial fracture simulation softwares make use of coupled equations and numerical simulations to simulate fracture propagation with fluid and proppant flow in a 3D space (Gu and Leung 1993; Sousa et al. 1993). Some widely used simulators do not take the anisotropy in elastic properties into account although they do take complex stress layering into account. Elastic anisotropy can be accounted 3
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
for during hydraulic fracturing using fully 3D finite element solvers (Gokaraju and Eckert 2014). However, they are time consuming and do not provide a practical solution especially when multiple cases must be run to pick the most effective treatment parameters.
static from dynamic using core measurements. Closure pressure can then be estimated using the elastic moduli as shown below.
H min
THE NEURAL NETWORK APPROACH
p p
An artificial neural network is a powerful tool in recognizing the relation between different parameters and outputs, which can be used to develop an optimized solution. The basic structure of a neural network consists of an input layer, an output layer, and multiple hidden layers. The neural network is usually trained and validated with multiple runs with known outputs; thereby, the model structure and coefficients can be automatically adjusted to optimize the predictions on existing runs and generalize well on the new data for future applications.
E h v ( v (1 ) p p ) E v (1 h ) Eh
1 h
2
H min
E h h 1 h
2
H max
(20)
where,
Hmin is the minimum horizontal stress (psi), v is the overburden stress (psi), pp is the pore pressure (psi), is the Biot’s elastic constant, Hmin is the minimum horizontal strain, Hmax is the maximum horizontal strain, and is the poroelastic constant. f(E,v) is then calculated from the elastic moduli using Eq. 15 or Eq. 16, depending on the formation.
This paper proposes a methodology in which a reservoirspecific neural network is built and trained to quickly and accurately predict fracture geometry, effective propped area for any type of formation, and any given completion parameters (e.g., perforation position, injection rate, injection time, etc.). This can be extended to run a parametric study which helps in selecting optimal hydraulic fracturing parameters to produce the largest EPA (shortterm production) or stimulated reservoir volume (SRV) (long-term production).
Because most simulation software uses isotropic elastic moduli, an equivalent Young’s modulus (Eeq) and Poisson’s ratio (eq) are calculated to best simulate anisotropic elastic moduli using a single E and Eeq can be calculated by the following equation which is a modified and weighted version of the equation used by Chertov (2012).
2
1
(21)
where vvh is the arithmetic averaging for vv and vh and ai is the weight coefficient from 0 to 1 (ah+av+avh=1). Using the calculated f(E,v) and Eeq, an equivalent veq can be calculated by Eq. 21.
1
,
,
, 2
(22)
In the second step, the rock mechanic properties (Eeq, veq) and the closure stress σh from step 1 are input into a 3D fracture modeling simulator. The completion parameters, such as slurry injection rate (qinj), total slurry volume (Qtol), and the perforation depth (TVDperf) are varied to generate the training database with the output results, such as fracture length (Lf), height (Hf), width (wf), and effective propped length (Leff). The training database is then used to produce the neural network. To improve the computation efficiency, multiple modeling cases are run by varying the fracture parameters. To generate an initial training database, each completion parameter is varied by x values equally distributed within the range of interest. The x value is dependent on the level of accuracy sought. For the two
Fig 1. Complete workflow for optimizing completion parameters using sonic measurements.
Fig. 1 illustrates the workflow of generating and training a neural network from sonic measurements and extended to optimize completion parameters. Once the elastic moduli of a rock are computed, they are calibrated and converted to 4
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
organic shales used in this study, x was chosen to be 3-5 (i.e., each parameter was varied up to 5 times). Therefore, if there are n parameters, the total number of training sample combinations is xn. In this study, there are three input completion parameters: slurry injection rate, total injection volume, and perforation position along the horizontal well. The horizontal position is correlated with depth based on geosteering data. Five output results are fracture length, fracture height, fracture width, upper and lower depth of fracture, and effective propped length.
the sonic velocities, Young’s modulus, and closure pressure. The interval can be divided into five zones based on closure stress. Fig. 3 shows a design of a toe-up horizontal well inside the pay. The toe-up strategy has been used for both cases. Fracture simulations will be performed on the lateral part of the shown hypothetical well. Shale 1 is a transversely isotropic media with vertical axis of symmetry (TIV) with horizontal Poisson’s ratio greater than its vertical counterpart. Hence, the modified ANNIE method was used to compute the stiffness tensor from sonic velocities. f(E,) and closure stress are computed using Eqs. 16 and 20 respectively. A 3D isotropic simulator was used to run the simulation models and Eeq and eq were calculated using Eqs. 21 and 22 respectively.
Once the neural network is generated, a group of testing data is run to check the relative error of the outputs for each input parameter. For any input parameter, if the tolerance relative error is not met between Node_i and Node _i+1, one more data point Node_i+1/2 is added in between. The extra cases regarding the added data point are run in fracture modeling to update the current training database, and further update the neural network. Testing is run on the new neural network. If the tolerance error is met for all input parameters, the neural network is left unaltered. Otherwise, it updates itself using the added database until the criterion is met. Once the neural network is obtained, one can either predict fracture geometry and location based on the input of arbitrary completion values, or determine the optimized fracturing design by conducting a massive parametric study with the neural network. For fracturing optimization, the effective propped length (EPL) is one of the best candidates to be the optimization target among all the predicted outputs. It is the propped length within the payzone occupied by infinite relative conductivity (Fcd>50). EPL dominates the short-term production and affects the longterm production. For better optimization, a critical conductivity can be used instead of infinite relative conductivity to define the effective propped length. The critical conductivity is defined as the minimum conductivity needed for fully stimulating a certain propped length during a certain production time. It is a function of propped length, production time, matrix permeability, natural fracture properties, oil API, and other completion and production parameters. (Gu et al. 2014).
Fig. 2. Uranium, sonic, lithology, elastic moduli, and stress data for Shale 1.
STUDY OF FIELD CASES The proposed methodology has been applied to two organic rich shales to develop an optimized EPL and completion parameters. In case 1, the stress barriers are highly defined with the stress difference between the payzone and the bounding layer being as high as 3,000 psi. Fig. 2 is a plot of 5
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
Simulations were performed for a single fracture in Shale 2 that acted as a training database for the generation of a neural network and an optimized solution. Tables 1 and 2 show the variation used for parameter inputs while generating a training database for Shales 1 and 2 respectively. X is the middle of the payzone.
Fig. 3. Toe-up horizontal well inside lower payzone.
For Shale 2, the stress barrier is not as well defined as case 1. Fig. 4 shows uranium, elastic moduli, lithology, and closure pressures. Because of its complex closure stress profile, the interval has been divided into 20 zones, two of which are hydrocarbon-bearing reservoirs.
Parameter
Min value
Max value
Step size
Perforation depth (ft)
X-50
X+50
20
Fluid injection rate (bbl/min)
30
50
10
Fluid injection Volume (M-gal)
10
30
5
Table 1. Input parameter variation for Shale 1.
Parameter
Min value
Max value
Step size
Perforation depth (ft)
X-53
X+72
20
Fluid injection rate (bbl/min)
15
30
5
Fluid injection Volume (M-gal)
10
40
5
Table 2. Input parameter variation for Shale 2.
Fig. 4. Uranium, lithology, elastic moduli, and stress data for Shale 2.
6
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
Shale 2. Fig. 6 is an example showing the target values of Leff and values predicted by the neural network in Shale 1. The solid dots indicate target values, while the hollow dots are predictions. The red ones are 70 training data points, while the green ones are seven random testing cases.
RESULTS AND DISCUSSIONS Fig. 5 shows a typical fracture profile result. This example was run on Shale 1.
Fig. 5. Typical fracture profile plot using fracture design and analysis
Fig. 6. Comparison of target values and values predicted by the neural
software.
network.
The track on the left shows the fracture width-height cross section, while the track on the right shows the length-height cross section. The color shaded contour represents the conductivity distribution. Only the conductivity distribution within the payzone (boundaries marked using yellow colored lines) governs the fracture productivity. The fracture conductivity decreases with increasing propped length as shown by the solid purple line. In the same plot, the minimum conductivity required to fully stimulate different propped length is indicated by red lines. The red lines represent the 1-year, 5-year, and 10-year production using the minimum conductivity criteria (Gu et al. 2014). According to the red lines, as fracture length increases, the minimum conductivity required to fully stimulate the length increases. For the same propped length, the minimum conductivity increases with decreasing production time. Therefore, the intersection between the purple and red lines represents the effective length for the certain production time. For example, for 1-year shot-term optimization, to meet the infinite conductivity criteria, the purple line should always be above the red-dotted line. Their cross is the maximum fracture length with the infinite conductivity. Through this way, the effective length at different production time scale is determined.
Once a neural network is generated for a specific reservoir, it can be used to predict fracture geometry and fracture geometry-based production performance for any combination of input parameters, which in this case are perforation depth, fluid injection rate, and fluid injection volume. Because the anisotropic elastic parameters and complex anisotropic closure pressure profile were both accounted for, the results can be generated instantaneously with a high degree of accuracy. In Fig. 7, this is demonstrated by randomly selecting the slurry injection rate and the total injection volume in Shale 1 to be 48 bbl/min and 23 Mgal. Fracture length, height, width, TVD of the upper and lower fracture boundaries, and effective propped length are all calculated as a function of horizontal well distance. Similarly, Fig. 8 shows fracture length, height, width, and effective propped length are all calculated as a function of depth in Shale 2 using slurry injection rate and the total injection volume to be 30 bbl/min and 40 Mgal respectively.
A neural network is generated using a training database consisting of multiple cases. The different cases are obtained by varying injection rate by three values, injection volume by five values, and perforation position by four values in Shale 1 and by varying injection rate, injection volume and perforation position each by five values in 7
SPWLA 56th Annual Logging Symposium, July 18-22, 2015
position. Leff is 400-460 ft for the first half of the horizontal well, while 360-400 ft for the second half. Injection rate does not have much effect on the effective length. For positions from 0.2 to 1, lower injection rate produces a little longer length. Within the well tip part (
