You are here
Home > BLOG > Hydrulic Fracture > Modeling hydraulic fractures in finite difference simulators using amalgam local grid refinement (LGR)

Modeling hydraulic fractures in finite difference simulators using amalgam local grid refinement (LGR)

Fig. 2 Overview of the static model


Hydraulic fracturing allows numerous, otherwise unproductive, low-permeability hydrocarbon formations to be produced. The interactions between the fractures and the heterogeneous reservoir rock, however, are quite complex, which makes it quite difficult to model production from hydraulically fractured systems. Various techniques have been applied in the simulation of hydraulically fractured wells using finite difference simulators; most of these techniques are limited by the grid dimensions and computing time and hardware restrictions.


Reda Abdel Azim & Sherif S. Abdelmoneim

School of Petroleum Engineering, University of New South Wales, Sydney, Australia. ENAP Sipetrol, Cairo, Egypt.
Received: 24 February 2012 / Accepted: 15 September 2012 / Published online: 3 October 2012
© The Author(s) 2012.

Most of the current analytical techniques assume a single rectangular shaped fracture in a single-phase homogeneous reservoir, the fracture is limited to the block size and the fracture properties are adjusted using permeability multiplier. The current work demonstrates how to model these systems with a smaller grid block size which allows you to apply sensitivity to the fracture length and model the fracture with enhanced accuracy. It also allows you to study the effect of reservoir heterogeneity on the fractured well performance.

It is proposed to apply amalgam LGR technique to decrease the grid size to the dimensions of the hydraulic fracture without dramatically increasing the number of grid blocks which would cause a great increase in the computing time and the model size with no added value. This paper explains how the amalgam LGR is designed and compares between standard LGRs and the proposed design and a case study is presented from an anonymous field in Egypt to illustrate how to use this technique to model the hydraulically fractured well. The simulation model is matched to available production data by changing fracture lengths. Then the model is used to predict future response from the wells.

The advantage of this technique is that it allows hydraulically fractured reservoirs to be modeled with less grid size which will lead to more realistic models and more accurate predictions; however, the most useful application of this technique may be in the fracture modeling stage. With this tool, various fracture geometries and scenarios can be tested in the simulator, and the most economic scenarios selected and implemented. This will lead to better fracture placement, and ultimately greater production.


Hydraulic fracturing is the process of pumping a fluid into a wellbore at an injection rate that is too high for the formation to accept in a radial flow pattern. As the resistance to flow in the formation increases, the pressure in the wellbore increases to a value that exceeds the breakdown pressure of the formation that is open to the wellbore. Once the formation “breaks-down”, a crack or fracture is formed, and the injected fluid begins moving down the fracture.

DOE research has developed several alternative fracturing techniques designed to accomplish specific tasks such as:

  • Tailored pulse fracturing
  • Foam fracturing
  • CO2, sand fracturing

In general, hydraulic fracture treatments are used to increase the productivity index of a producing well or the injectivity index of an injection well. The productivity index defines the volumes of oil or gas that can be produced at a given pressure differential between the reservoir and the well bore. The injectivity index refers to how much fluid can be injected into an injection well at a given differential pressure.

One of the major problems facing the reservoir engineers in modeling the hydraulic fractures using the finite difference simulators is the wide gap between the grid size of the reservoir model and the fracture dimensions.
The purposes of this paper are to model the flow of the reservoir fluid in hydraulically fractured reservoir using finite difference simulators in a manner that would allow the simulator to mimic the actual fracture geometry without dramatically increasing the number of grid cells and hence increasing the computing requirements and time.

This is achieved as shown in the paper by using amalgam LGR to decrease the fractures dimensions to the size and dimensions required to achieve this goal and leave the number of grid cells only slightly affected which makes a minor change in the required computing time and capabilities. This design was tried on an anonymous field in the western dessert in Egypt, and the results were compared with the actual production data which was recorded after the fracture to verify that the model was capable of modeling the actual reservoir performance.

Hydraulic fracture mechanics

The theory of hydraulic fracturing depends on an understanding of crack behavior in a rock mass at depth. Because rock is predominantly a brittle material, most efforts to understand the behavior of crack equilibrium and growth in rocks have relied on elastic, brittle fracture theories.

However, certain aspects, such as poroelastic theory, are unique to porous, permeable underground formations. The most important parameters are in situ stress, Poisson’s ration, and Young’s modulus.

In situ stresses

Underground formations are confined and under stress. Figure 1 illustrates the local stress state at depth for an element of formation. The stresses can be divided into three principal stresses. In Fig. 1, σ1 is the vertical stress, σ2 is the maximum horizontal stress, while σ3 is the minimum horizontal stress, where σ1 > σ2 > σ3. These stresses are normally compressive and vary in magnitude throughout the reservoir, particularly in the vertical direction (from layer to layer).

Fig. 1 The local stress state at depth for an element of formation

Fig. 1 The local stress state at depth for an element of formation.

The magnitude and direction of the principal stresses are important because they control the pressure required to create and propagate a fracture, the shape and vertical extent of the fracture, the direction of the fracture, and the stresses trying to crush and/or embed the propping agent during production.

A hydraulic fracture will propagate perpendicular to the minimum principal stress (σ3). If the minimum horizontal stress is σ3 the fracture will be vertical and, we can compute the minimum horizontal stress profile with depth using the following equation.

we can compute the minimum horizontal stress profile with depth using the following equation.

Poisson’s ratio can be estimated from acoustic log data or from correlations based upon lithology. The overburden stress can be computed using density log data. The reservoir pressure must be measured or estimated. Biot’s constant must be less than or equal to 1.0 and typically ranges from 0.5 to 1.0. The first two terms on the right hand side of the equation represent the horizontal stress resulting from the vertical stress and the poroelastic behavior of the formation.

Poroelastic theory can be used to determine the minimum horizontal stress in tectonically relaxed areas (Salz 1977). Poroelastic theories combine the equations of linear elastic stress–strain theory for solids with a term that includes the effects of fluid pressure in the pore space of the reservoir rocks.

The fluid pressure acts equally in all directions as a stress on the formation material. The “effective stress” on the rock grains is computed using linear elastic stress–strain theory. Combining the two sources of stress results in the total stress on the formation, which is the stress that must be exceeded to initiate fracturing. In addition to the in situ or minimum horizontal stress, other rock mechanical properties are important when designing a hydraulic fracture. Poisson’s ratio is defined as “the ratio of lateral expansion to longitudinal contraction for a rock under a uniaxial stress condition (Gidley et al. 1989)”.

The theory used to compute fracture dimensions is based upon linear elasticity. To apply this theory, the modulus of the formation is an important parameter. Young’s modulus is defined as “the ratio of stress to strain for uniaxial stress (Gidley et al. 1989)”.

The modulus of a material is a measure of the stiffness of the material. If the modulus is large, the material is stiff. In hydraulic fracturing, a stiff rock will result in more narrow fractures. If the modulus is low, the fractures will be wider. The modulus of a rock will be a function of the lithology, porosity, fluid type, and other variables.

Field case study

Field characterization

We conducted a fracture reservoir simulation study for well A-1 (Appendix 1, Fig 16) over Abu-Roash G reservoir in the field (A), in the western dessert of Egypt.

This field produced light gravity oil (37° API) from Abu-Roash G sand at an average drilled depth of 5,400 ft TVDSS. Only two wells were drilled in the area and had been reviewed in the field study. The interpretation of the well logs shows hydrocarbon bearing in middle A/R G formation which is subdivided into two sand bodies. The two sand bodies were perforated and tested; they showed production with initial rate 370 BOPD by N2 lifting with traces of water and gas production. Production started from well A-1 (Appendix 1, Fig 16) with ESP yielding a production rate of 350 BOPD and 35 % water cut, then the well production rapidly declined to 130 BOPD and water rate started to decline. The ESP then failed several times and has been replaced with sucker rod pump. The last static fluid level measurement showed average static reservoir pressure of 1,247 psig at a reference depth of −5,300 ft TVDSS.

No PVT samples were taken from this field, so calculations were done using both correlations and PVT samples taken from the nearby producing fields. An estimation of the oil in place was done using both material balance and volumetric, showing reserves ranging from 10 to 15 mm STB. Only compression velocity was available from the sonic log and the shear velocity was calculated using the following correlations for both the sand and shale layers. A static model was built using the available data with approximate cell dimensions of 50 × 50 m and a dynamic model was successively prepared to be used to test the different fracture scenarios (Fig. 2).

Fig. 2 Overview of the static model

Fig. 2 Overview of the static model.

The basic simulation model is a three-phase flow single-well model; however, only the oil and gas phase are mobile. The water is not moving as noticed from the well history data, and the amount of produced water at early production period was coming from the water that had been used as a completion fluid, so there is essentially no water production from the well.

Sensitivity runs are completed to determine if grid block size has any impact on the production rates and consequently the fracture length. Static and dynamic models were constructed for this case study; the dynamic model has grid size of 50 × 50 m, consisting of 72,240 cells (80 × 43 × 21). The model is divided into vertical layers. These layers contain productive zones that are fractured and non-productive shale that separate the productive zones. Porosity and permeability for the producing layer is based on the typical values observed in Table 1.

Porosity and permeability for the producing layer is based on the typical values observed in Table 1.

Table 1 Porosity and permeability of the studied field.

39 LGRs were used to model the hydraulic fracture scenarios in three wells; well A-1 (Appendix 1, Fig 16, current well) and the two proposed wells A-2 (Appendix 1, Fig. 17) and A-3 (Appendix 1, Fig. 18). The LGRs were amalgamated in three groups each at the location of each well and were designed so as to match the actual geometry predicted from the rock model as previously mentioned. The original grids were refined depends on the selected length of the fracture to be used in the simulation model.

Six grids with size ~145 ft were refined in all directions around the well as Fig. 3 shows; four of them have been refined regularly and the other cells were refined only in one direction depending on the fracture length. For example, the fracture length used in the current model is 600 ft, so represent this length in grid cells with very small fracture aperture; the grid cells should be refined according to that length. Four grid cells of 145 ft+ the rest of the fracture length from the fifth cell around the wellbore, it gives the desired fracture length with its aperture.

Fig. 3 LGRs around the wellbore

Fig. 3 LGRs around the wellbore.

By using permeability multiplayer keyword only for the fracture, the hydraulic conductivity value of the created fracture increased as shown in Fig. 3. The red line represents the simulated fracture. Coupling of these LGRs was done by amalgam technique. Simulation of hydraulic fracture using this method resulted in decreasing the convergence and stability issue during the running process; also it gave accurate results.

Fig. 4 LGRs around the wellbore

Fig. 4 LGRs around the wellbore.

39 LGRs have been distributed around the main well (1) and the other proposed wells in four directions around the wellbore. Figure 4 shows how the LGRs distributed on the other wells. Table 2 shows some of the LGRs lengths used in the simulation model for the main well.

Table 2 Length of LGRs used in the model

Table 2 Length of LGRs used in the model.

The PVT data were taken from the nearby operating fields and reservoir properties were populated based on the understanding of the depositional environment to match the existing two wells. The model showed a good match with the original production and pressure data. Table 2 and Fig. 5 show the fluid and SCAL data for the current field (Table 3).

Fig. 5 Relative permeabilities curve

Fig. 5 Relative permeabilities curve.

Table 3 PVT and rock data of the studied field

Table 3 PVT and rock data of the studied field.

Rock mechanics properties of the simulation model

The model presented here is to simulate a poroelastic reservoir which is intercepted by a vertical wellbore. The reservoir is horizontal and the Cartesian coordinates of the reservoir are aligned in the same direction of minimum and maximum horizontal in situ stresses.

Formation properties Formation rock properties such as Young’s modulus, Poisson’s ratio, Biot’s coefficient, rock strength are assumed to stay constant throughout the simulation time.

Field case study: results

Determining the fracture length

The simulated model is matched to available history data. First the model is run with no fracture. Then an observation has been taken to the matching process. If the production from this model is greater than the measured production, then an adjustment for the porosity and permeability is applied on the grid block according to the petrophysical data to achieve the matching process, and there is no need for running the model with hydraulic fracture scenario. If it does not happen, then the model with initial fracture length is run, and finally the fracture length is changed to match the measured production rates.

Because the permeability of the studied field is very low, changing the friction factor does not significantly change the production rates. When the fracture permeability is increased, it does not result in increase of the flow rates. This is because the fracture acts like infinite conductivity fracture. History-matching process of the producing well generates fracture length of ~600 ft. Porosity and permeability values are not changed and stay close to the range of the data measured for the well. The result of history match of the well is shown in Fig. 6a. The line with the lowest production is the run was no fracture includes in the model. Figure 6b shows how the simulated reservoir pressure matched with the observed one.

Fig. 6 a Relative permeabilities curve. b Reservoir pressure matching

Fig. 6 a Relative permeabilities curve. b Reservoir pressure matching

Fig. 6 a Relative permeabilities curve. b Reservoir pressure matching.

Although the observed data do not have any recording for the produced gas rates, the model was able to give a trend of the produced gas during the history-matching process by using the solution gas oil ratio (Rs), as Fig. 7 shows.

Fig. 7 Gas producing rates

Fig. 7 Gas producing rates.

Future production rate

The actual production data and model results are from October 2004 to December 2009. In this section, the simulator is run longer. The simulated production rate curves are extended from December 2009 to December 2019. The simulated values are compared to the actual values between October 2004 and December 2009. This period gives an indication of how good the history-matched model can predict production. The extended simulated rates in the period mentions are found in Fig. 8.

The extended simulated rates in the period mentions are found in Fig. 8.

Fig. 8 Prediction of oil rates.

Comparsion the LGRs method with finite element method

A fully poroelastic model is developed to simulate the rock deformation and fluid flow in a reservoir. The governing equations of poroelasticity are solved simultaneously to give change in stress and pressure in the reservoir and it is related to rock deformation by assuming rock to be linearly elastic. The governing equations that guide fluid flow occurring between the reservoir rock and the wellbore are developed from diffusivity equation. These governing equations, which are derived on the basis of mass continuity equations (for both fluids and solids), Darcy law and equation of equilibrium of solids, are as follows (Charlez 1997: Chen et al. 1995).

The governing equations that guide fluid flow occurring between the reservoir rock and the wellbore are developed from diffusivity equation.

where ϕ is the porosity; ct is the total compressibility; p is the pore pressure; α is the Biot’s coefficient; u is the displacement vector; k is the permeability tensor; μ is the fluid viscosity; cf is the fluid compressibility; ∇ is a vector operator; K is the bulk modulus and G is the shear modulus.

Discretizing the governing equations of poroelasticity by using the finite element method (FEM) (Zienkiewicz and Taylor 2000) results in the following coupled linear systems of equations:



Model geometry includes wellbore, reservoir and two pairs of finite conductive hydraulic fractures (see Fig. 9).The initial maximum and minimum horizontal stresses are along x, y axes, respectively.

Fig. 9 Mesh used in finite element method

Fig. 9 Mesh used in finite element method.

Emanuel Martin
Emanuel Martin is a Petroleum Engineer graduate from the Faculty of Engineering and a musician educate in the Arts Faculty at National University of Cuyo. In an independent way he’s researching about shale gas & tight oil and building this website to spread the scientist knowledge of the shale industry.

Leave a Reply

12 + 12 =