Damage Model for Reservoirs with Multisets of Natural Fractures and Its Application in the Simulation of Hydraulic Fracturing

Abstract

The presence of natural fractures can significantly affect the quality of hydraulic fracturing operations in tightsand and shale or oil/gas formations. This paper describes the procedure used to model natural fractures with continuum damage tensor and the resulting orthotropic permeability tensor. A damage model that uses damage variable in tensor form is presented. In the procedure presented, a nonnegligible angle is assumed to exist between directions of principal stresses in the formation and in the natural-fractures-related damage tensor, and this difference in orientation is modeled by introducing local directions in the model. A damage-dependent permeability tensor in tabular form is then proposed. A second case scenario when the directions of principal stresses and natural fractures align is also analyzed. Numerical results of fracture distribution are presented for both cases, and differences can be seen from the computed contour of the damage variable. The results indicate that the model can effectively simulate the fracture propagation phenomena during hydraulic fracturing.

1. Introduction

Modeling of porous flow that accurately predicts changes in permeability is important to the hydraulic fracturing operation in oil/gas reservoirs. Natural fractures in reservoir formations can significantly affect stimulation-induced fracture propagation and hence the production of oil/gas with hydraulic fracturing. For more than fifty years natural fractures and their role in the modeling of porous flow in oil/gas reservoirs have been investigated by various researchers [1,2,3,4,5,6,7]. Porosity, permeability and connectivity are the three properties of natural fractures that are more commonly used during modeling of porous flow and are referred to as porous flow properties. For a hydraulic fracturing process, these properties can be derived from logging data or be derived from inverse analysis using well production testing flow and pressure tests [5]. The accuracy of these properties is most relevant to the effectiveness of hydraulic fracturing practice.

Various models are available that measure the porosity and permeability of reservoirs with natural fractures based on the permeability values within the matrix and within the natural fractures. In total three kinds of models can be distinguished, the single-porosity and single permeability model, the dual porosity and single permeability model and the dual porosity and dual permeability model [6]. Because the formation matrix of a reservoir with natural fractures (i.e., fractured reservoir) is cut into small blocks by the fractures, the fractures consequently have high permeability but low porosity, and the matrix instead has low permeability but high porosity. The resulting formation model is usually one of dual porosity and dual permeability. Determining the permeability in such situation includes separately computing the porous flow within the porous matrix and then computing the porous flow within the natural fractures and the process is often complex. The dual-porosity and single-permeability reservoir model is instead a simplified model. As a step to simplify, it considers the matrix porosity as a measure for storage only, not for porous flow, and porous flow exists only within the fractures. Then a single-porosity and single-permeability reservoir model can further simplify fractured reservoir. In this case, porosity and permeability values are determined for the equivalent reservoir model that represents their functions of porosity and permeability in the matrix and fractures as a whole.

Continuum damage mechanics (CDM) uses continuum variables to describe cracks and joints, thus providing a means of modeling natural fractures. CDM-based modeling of fracturing is referred to as a damage model and it is effective for simulating the development of fractures caused by well stimulation [7,8,9,10,11]. The Barcelona model based on CDM was proposed by Lubliner et al. [12] originally and Lee and Fenves [13] enhanced its features and applicability to various engineering problems. It is a continuum damage model based on scalar plasticity. Damage arising from tensile cracking and damage resulting from compressive crushing are both included in the Barcelona model’s damage development process. Two hardening parameters influence the evolution of plastic loading: equivalent plastic strain and tensile load.

This paper leverages the CDM method and presents a damage model for simulating the porous flow during hydraulic fracturing in single-porosity single-permeability reservoir with natural fracture. To the knowledge of the authors of this paper, the proposed methodology has only ever been applied by the authors to simulate the hydraulic fracturing in the Dina 2 gas field in the Tarim Basin in Xinjiang, China, and that this is the first time a continuum-damage model has been simulated to help optimize injection in such geological conditions. In theory the procedure presented in this work can be applied to any tightsand formations with fractures, or shale oil/gas formations with natural fractures. The permeability value and porous flow predicted using the proposed damage model is for the type of reservoir that consists of matrix and fractures as a whole.

In the following, first a damage model is presented and the numerical approach for calibrating the proposed model with measurable data from natural fracture multisets (acquired in previous work) is presented. Then the suggested model is used to demonstrate hydraulic fracturing in reservoirs with several sets of natural fractures. Numerical results are presented to show that the hydromechanical finite-element model is an appropriate tool for predicting stimulated volume in tightsand reservoirs.

2. Materials and Methods

2.1. Expression of Natural Fractures with Continuum-Damage Variable

Geometrical information of a set of natural fractures is usually characterized by their volumetric density, azimuth angle, and inclination angle. Furthermore, fracture properties related to porous flow can be characterized by aperture, permeability, and porosity. A vector form is used to describe natural fractures. A second-order tensor form is used here to describe the continuum-damage variable. The second-order tensor is used because the tensor form damage variable is easy to connect with orthotropic permeability. The principal orthotropic permeability value has a unique holonomic relationship with the principal damage value in its related direction.

For a set of natural fractures, the continuum damage tensor ωcomprises one principal value ω and a principal directional vector n, which points to the same direction as the natural fracture. The principal value, which is the magnitude of the continuum-damage tensor ω is determined by the volumetric density value of the natural fractures ψ. Various versions define natural fracture density (e.g., number per foot or fracture area per unit volume) because of the constraints of the various measuring methods. Therefore, the unit of fracture density ψ is not the same for each version. Consequently, there is no direct explicit relationship between the value of the continuum-damage variable ω and the fracture density ψ. However, there is a functional relationship betweenω and ψ, which is defined as f(ψ) in Equation (1)

$$\mathbf{\omega }=\mathit{f}\left(\psi \right)$$

(1)

The f(ψ) should be determined by a phenomena-match method that uses either reservoir pressure data measured in the field or in specific porous flow testing data measured in the laboratory. The phenomena match method is an empirical method first mentioned by Swoboda et al. [8]. The method relies on the geology experts’ estimation of the initial damage variable of observed cracks. It has not been studied well enough to derive an analytical expression, so in practice initial conditions of numerical computations have to rely on educated guesses based on geological measurements.

The direction of a set of natural fractures can be described by a unit directional vector $\mathit{n}$, which is also used as a principal directional vector of damage tensor. The direction of this vector $\mathit{n}$ has three directional components (l, m, n). Therefore, the damage variable tensor in terms of vector of coordinate basis ($\mathit{i},\mathit{j},\mathit{k}$) is defined in Equation (2) as:

$$\mathbf{\omega }=\mathbf{\omega }\mathit{n}\otimes \mathit{n}, \mathrm{where} \mathit{n}=l\mathit{i}+m\mathit{j}+n\mathit{k}$$

(2)

Equation (2) can be rewritten in component form as:

$${\omega }_{ij}=\omega n_i{n}_j, i=1,3; j=1,3 \mathrm{for} 3\mathrm{D} \mathrm{problem}$$

(3)

In practice, these directional components (l, m, n) can be calculated in terms of the values of azimuth angle, inclination angle, and x-direction of the coordinate system adopted in a calculation. Figure 1 (illustrations made by the authors for clarity) shows the spatial relationship between these angles and the directional components of the vector that represents the natural fractures set. In Figure 1, θ is the azimuth angle, β is the inclination angle, and the x-direction of the coordinate system is taken in the north direction.

These relationships are expressed mathematically in Equation (4):

$$\{\begin{array}{c}l=sin\beta cos\theta \\ m=-sin\beta sin\theta \\ n=cos\beta \end{array}$$

(4)

Every set of fractures has a corresponding damage tensor, because each have their own inclination angle and other model-relevant properties.

A multiset of fractures that contains many fractures may share the same value for certain parameters such as inclination angle. While there may be many natural fractures distributes in every direction, in practice, it is suggested to have at maximum 3 sets of fractures, so that there are no more than 3 major directions of natural fractures.

We define the synthetic damage tensor ${\mathbf{\omega }}_{\mathit{T}}$ to be equal to the summation of the damage tensors for each single natural fracture when there are multiple sets of natural fractures:

$${\mathbf{\omega }}_{\mathit{T}}={\mathbf{\omega }}_{\mathbf{1}}+{\mathbf{\omega }}_{\mathbf{2}}+\dots {\mathbf{\omega }}_{\mathit{N}}$$

(5)

Equation (5) can be rewritten in component form as:

$${\omega }_{ij}^{\left(T\right)}={{\displaystyle \sum }}_{r=1}^N{\omega }_{ij}^{\left(r\right)}$$

(6)

where N is the number of sets of natural fractures, superscript r varies from 1 to N.

2.2. Damage Initiation Condition

The damage initiation condition is the condition when the damage variable will start to evolve. The condition adopted in this work is the condition proposed by Mazars et al. in [9,10]. For convenience of reading, it is rewritten in Equation (7):

$$f_t={\epsilon }_t-Y_t\le 0$$

(7)

Subscript t represents tension. The condition can be translated in plain English as the condition for damage initiation is reached when the equivalent deformations favoring tensile strain becomes less than or equal to the thermodynamic variable that characterizes the extreme loading state in the tensile part. The calculation of the damage initiation condition using Equation (7) is performed on each of the principal directions of the orthotropic damage model. Only the case of tensile strain is accounted for in the calculation of damage initiation. For the case of shear fracture, estimation on each one of the three principal directions of strain tensor will be done first. If damage loading criterion expressed in Equation (7) be met in any one of the three principal directions, damage variable will start to evolve.

2.3. Damage Evolution Law

The damage evolution law adopted in this work is a holonomic damage model proposed by Mazars et al. in [11]:

$$\omega =1-Y_{t1}/Y_t,$$

(8)

$Y_{t1}$ is the maximum value of the damage conjugate force that corresponds to the maximum fracture opening. The primary advantage of the Mazars’ holonomic damage model is that the damage variable value can be calculated for a given state of strain. Iteration at the local level of a material point is not necessary.

2.4. Damage-Dependent Permeability

In this work, a damage-dependent permeability model is used. In general, the value of permeability $K_i$ is assumed to depend on the value of the continuum-damage variable ${\omega }_i$:

$$K_i=K_i\left({\omega }_i\right),$$

(9)

where subscript i represents the i-th component of orthotropic principal direction, and i = 1, 2, 3 for a 3D problem. $K_i$ represents the value of permeability in the i-th direction of an orthotropic permeability model. ${\omega }_i$ is the variable of continuum damage in the same i-th direction as $K_i\prime s$. ${\omega }_i$ is also interpreted as the component of damage tensor ω in its i-th principal direction.

In calculation, the damage-dependent permeability model is introduced through tabular form of a data series, which are calibrated by using the phenomena-match method. These three orthotropic principal directions of permeability are determined by the directions of the natural fracture, accordingly.

2.5. Finite Element Model (FEM): Procedure and Software

The procedure to build an FEM for hydraulic fracturing analysis includes the following major steps:

• Build the mesh of the FEM with geometry of the object to be analyzed;
• Assign material properties to the mesh with values of parameters of material properties obtained from 1D geomechanics analysis based on logging data etc;
• Establish the initial geostress field with aforementioned 1D geomechanics solution as well as measured data of minimum horizontal stress component if any;
• Introduce boundary conditions and loading condition includes gravity. Fracturing fluid injection will be modeled as in-flow to the nodes which represent perforation sections in production wells;
• Initial geostress will be introduced into the FEM as prescribed initial conditions.
• Natural fractures will be modeled as initial damage in a form of orthotropic second-order damage tensor.
• Orthotropic permeability tensor will be used in the model’s porous flow equations.
• In order to build the orientation of the damage tensor and the permeability tensor, it is necessary to build an orientation in advance in the Abaqus model before the material properties being assigned to the model.

The procedure is then carried out with the software Abaqus by Dassault Systemes. The model is implemented into the software through writing customized user subroutines with all user-defined properties, specifications and computational steps.

3. Results

The procedure and models proposed in previous paragraphs are used in this section to determine the initial value of the continuum-damage tensor of a tightsand oil reservoir with three sets of natural fractures. Damage tensor and orthotropic permeability are further applied to simulate hydraulic fracturing of a zipper fracture of two horizontal wells in a tightsand oil formation. Here, the authors present the model of hydraulic fracturing of the formation and numerical results of fracture generation and propagation represented by the damage variable.

3.1. Geological Properties of Natural Fractures

The reservoir is located in the Dina 2 gas field in Tarim Basin, Xinjiang, China. The rose plot shown below is a summary of cracks measured from 20 wells in the area. Three sets of natural fractures were identified from image logging data (acquired during a previous project and unpublished elsewhere) of the given fractured tightsand reservoir formation as shown in Figure 2.

Red colored circular areas represent the location where values of inclination and orientation of three sets of natural fractures were determined by image logging data and field experience. The three locations with measurements of fractures are named Frac Set 1, Frac Set 2 and Frac Set 3.

Table 1. Geological properties of natural fractures.

Location	Azimuth Angle (θ/°)	Inclination Angle (β/°)	Fracture Density (1/m)	Aperture (mm)
Frac Set 1	275	72	1	0.5
Frac Set 2	350	68	0.6	0.5
Frac Set 3	80	70	0.4	0.5

provides the azimuth angle, inclination angle, fracture density and aperture values of the given sets of natural fractures. The values listed in Table 1 are taken from real world engineering application in a case of tightsand gas formation in the Tarim Basin.

3.2. Initial Damage Variable Calculation

Equation (10) defines the initial principal damage variable values for the three sets of natural fractures given in Section 3.1 as:

$${\omega }_1=0.015,{\omega }_2=0.03,{\omega }_3=0.01.$$

(10)

Due to the highly nonlinear relationship between the damage variable and the properties of natural fractures, there is no analytical expression between them. Equation (10) are given on the basis of data listed in Table 1 and experience in numerical simulation. The values in Equation (10) should be calibrated with phenomena match by using observed microseismic data related to hydraulic fracturing operations.

The derivation of the initial damage variable tensor is given in Equations (1)–(6). In particular, with reference to the Equations (2)–(4) and (10), damage tensors for the Set 1 to Set 3 fractures are calculated in Equations (11)–(14):

$${\mathbf{\omega }}_1=0.015\left(\begin{array}{c}0.083\\ 0.947\\ 0.309\end{array}\right)\left(\begin{array}{ccc}0.083& 0.947& 0.309\end{array}\right)=\left[\begin{array}{ccc}0.0001& 0.0012& 0.0004\\ & 0.0135& 0.0044\\ Sym& & 0.0014\end{array}\right]$$

(11)

$${\mathbf{\omega }}_2=\left[\begin{array}{ccc}0.025& 0.0044& 0.0102\\ & 0.0008& 0.0018\\ Sym& & 0.0042\end{array}\right]$$

(12)

$${\mathbf{\omega }}_3=\left[\begin{array}{ccc}0.0003& -0.0015& 0.0006\\ & 0.0086& -0.0032\\ Sym& & 0.0012\end{array}\right]$$

(13)

The synthetic damage tensor for the given three sets of natural fracture is:

$${\mathbf{\omega }}_T=\left[\begin{array}{ccc}0.0255& 0.0041& 0.0112\\ & 0.0228& 0.0031\\ Sym& & 0.0068\end{array}\right]$$

(14)

A set of direct solution of principal values and principal directions are solved using the synthetic damage tensor ${\mathit{\omega }}_{\mathit{T}}$ defined in Equation (14). The principal values are:

$${\omega }_1=0.0332,{\omega }_2=0.0015,{\omega }_3=0.0204$$

(15)

Correspondingly, the principal directions of ${\omega }_T$ are:

$$n_1=\left(\begin{array}{c}0.808\\ \begin{array}{c}0.437\\ 0.3945\end{array}\end{array}\right), n_2=\left(\begin{array}{c}-0.415\\ \begin{array}{c}-0.052\\ 0.908\end{array}\end{array}\right), n_3=\left(\begin{array}{c}-0.4716\\ \begin{array}{c}-0.897\\ -0.139\end{array}\end{array}\right)$$

(16)

In Equations (4) and (16), the values of azimuth angle and inclination angle are defined as:

$${\theta }_1=61.56°, {\beta }_1=66.76°, {\theta }_2=277.89°, {\beta }_2=24.77°, {\theta }_3=335.06°, {\beta }_3=82.11°$$

(17)

The above principal directions are orthotropic in 3D space. Values are derived from the initial damage variable in Equation (10) and the natural fractures information from Table 1. Angles defined in Equation (17) will be used as the principal directions of orthotropic permeability tensor.

3.3. Numerical Simulation of Hydraulic Fracturing of a Formation with Natural Fractures

3.3.1. Simplification of Orthotropic Permeability and Damage Tensor and Directions of Principal Stress

Orthotropic permeability caused by multisets of natural fractures introduces orthotropic properties to the entire numerical model. In general, principal stress directions are not the same as the initial permeability tensor derived from natural fractures. In this case, a global system of numerical model coordinates should be selected from the two, along either the principal stress direction or the principal orthotropic permeability tensor direction.

In the following, a simplified model will be used to show the principle of initial damage’s calculation. In order to keep the model as simple as possible, while retaining its functions and principles, a model is chosen consisting of two sets of natural fractures with directions different from those of the principal stresses.

For the purpose of brevity, the values of the initial damage tensor parameters used in Equations (10) and (17) are modified to the following values with only two sets of fractures:

$${\omega }_1=0.0132, {\theta }_1=73°, {\beta }_1={\beta }_2=90°; {\omega }_2=0.0015, {\theta }_2=163°.$$

(18)

Thereby, the chosen numerical model coordinates can be in the same principal direction as the initial damage tensor and initial orthotropic permeability tensor.

Next, the direction of minimum horizontal stress (Sh) is assumed at an azimuth angle ${\theta }_1=90°$. Therefore, there is a 73° angle between the directions of principal stresses and the orthotropic permeability tensor.

The damage model proposed by Mazars et al. [10,11] is adopted in the following. The damage initiation criterion is set by a threshold value of damage conjugate force $Y_t$ as follows

$$Y_t=0.0005.$$

(19)

The maximum value of damage conjugate force $Y_{t1}$ when the damage variable reaches 1 is given as:

$$Y_{t1}=0.015$$

(20)

Figure 3 illustrates the damage-dependency of hydraulic conductivity. It is calibrated using the observed phenomena of microseismic data during hydraulic fracturing of offset wells. Both the vertical and horizontal axes in Figure 3 are shown in logarithm scale.

Values of the data points in Figure 3 are given in the figure. With the set of parameters values, numerical results of shear strain distribution can match the distribution of measured microseismic data. Damage is a dimensionless value.

3.3.2. FEM Mesh, Boundary Conditions and Initial Conditions

The finite element model shown in Figure 4a is 200 m in both width and length and 25 m in height. A 20 m-thick reservoir is covered with a 5 m nonpermeable overburden layer. Injection loading is a point flow in the center of the model. Figure 4b shows a 73° angle between the principal directions of permeability tensor and x-axis. The direction of the x-axis is taken in the direction of initial minimum horizontal principal stress. The geostress level is set for a value as self-weight at true vertical depth (TVD) = 1600 m, and it is described with effective stress. Equations (21) and (22) define the initial conditions.

The following values of stress components are defined for cap rock:

$$S_v=-36.368\times {10}^6\mathrm{Pa}, S_h=-28.59\times {10}^6\mathrm{Pa}, S_H=-31.66 {10}^6\mathrm{Pa}$$

(21)

Values of effective stress components, pore pressure, and void ratio in the reservoir are:

$$S_v=-15.7\times {10}^6\mathrm{Pa}, S_h=-7.93\times {10}^6\mathrm{Pa}, S_H=-11\times {10}^6\mathrm{Pa},P_p=20.66\times {10}^6\mathrm{Pa}, VR=0.15.$$

(22)

The void ratio value represents the total porosity of the formation that includes both the matrix and natural fractures.

For the purpose of convenience of illustration as well as possible comparison between numerical results of damage contour and the microseismic data of a hydraulic fracturing operation, the concept of damage intensity $\overline{\omega }$ is proposed and is defined as:

$$\overline{\omega }=\sqrt{{\omega }_1^2+{\omega }_2^2+{\omega }_3^2}$$

(23)

In Figure 4a the mesh of the finite element model is given. In Figure 4b, the blue markers (right angle line segments) show the principal directions (orthogonal to each other) of the orthotropic damage tensor. These values are determined by the properties of the natural fractures through Equations (10)–(17). The principal directions of the permeability tensor are the same as that of the damage tensor. Because we are modeling natural fractures in the framework of CDM, natural fractures in this model cannot be shown in a form of a figure, but can only be shown mathematically, such as in the form of Table 1 and Equation (10).

Injection loading is applied to the central point of the model, which is the red dot within Figure 4a. It represents injection that occurs at an 8 m long perforation section of a horizontal well, with the value of injection rate vs. time, as shown in Figure 5. In Figure 5, the vertical axis shows negative values because it indicates “injection flow” instead of “outflow” rate, and the negative sign adopted as the sign convention for fluid flow loading.

3.3.3. Numerical Results When Principal Directions of Natural Fracture Differ from Those of Geostress

Numerical simulation of the hydraulic fracturing operations was performed by finite element (FE) software Abaqus by using the information presented previously. Specifically, the model is implemented into Abaqus by writing user subroutines. As previously described, the principal directions of natural fractures are different from those of geostress. Figure 6, Figure 7 and Figure 8 show the numerical solution of pore pressure distribution, damage propagation, and the evolution of shear strain intensity γ_T for two moments in time t₁ = 17.53 and t₂ = 26.64 (minutes) during injection loading.

Figure 6 shows the contour of pore pressure at two moments in time. Figure 6a shows that the pore pressure migrates along the direction of the principal permeability value at the beginning. However, a low pore pressure area is formed ahead of the region of higher pore pressure. Consequently, as shown in Figure 6b, at time moment t₂, the path of the pore pressure migration varies to another direction. Figure 7 shows the contours of ${\omega }_T$ and the contours of damage intensity $\overline{\omega }$ at these two moments in time. The value of damage intensity $\overline{\omega }$ is the total value of continuum damage caused by the natural fracture and hydraulic fracturing operations. Figure 8 shows the contour of shear strain intensity ${\gamma }_T$, which is defined as ${\gamma }_T=\sqrt{{\gamma }_{12}^2+{\gamma }_{13}^2+{\gamma }_{23}^2}$ and used for shear strain intensity illustration purposes.

Figure 6 shows a boundary effect in this pore contour generated by hydraulic fracturing operations. The pore pressure contour that is near the boundaries is significantly nonuniform. Therefore, only the solution of pore pressure some distance away from the boundary should be used for analysis.

Figure 7b shows that in time the synthetic damage contour forms a complex fracture network in the near field of the injection point. The direction of major damage distribution bands is along the two principal directions of damage tensor, which is actually the directions of natural fractures. This phenomenon is of important significance to the work on optimized fracturing design. After a period of fluid injection of fracturing operation, the principal directions of fractures are indicated by the direction of the localization band of damage variables which is bright band in Figure 7. This change is caused by the difference between principal directions of the stress tensor and that of the damage tensor.

The contour evolution of shear strain intensity γ_T is shown in Figure 8a. Figure 8b indicates the location of expected micro seismic activity, which can be recorded by monitoring the shear sonic data during a hydraulic fracturing operation.

3.3.4. Numerical Results When Principal Directions of Natural Fracture Overlap with Those of Geostress

For the purpose of comparison, an example of hydraulic fracturing is presented for the case that no angle exists between the directions of principal stresses and that of natural fractures. The values of all other data are the same as those used in Section 3.3.3.

It is alternatively assumed that there is no angle between the natural fractures and that of the principal stress directions. Figure 9, Figure 10 and Figure 11 show the numerical results. Figure 9 shows that the pore pressure contour develops in the directions of principal stresses and develops in a symmetrical pattern.

Because the natural fracture is homogenized as continuum damage and embodied into permeability, induced fracture networking is related to natural fractures but does not directly come from natural fractures. There is no network of natural fractures in the model presented here.

Figure 10 shows the fracture network generated by hydraulic injection. This phenomenon of networking of the damage contour actually is generated by damage localization. The mechanism behind this phenomenon could be plural. One of the mechanisms should be the material instability related to strain-softening. The other factor is relevant to the porous behaviors that pore pressure takes the path along which the permeability value is higher than that of other places within a limited area which is so-called stimulated reservoir volume (SRV). For a model with a homogeneous and uniform material, and a permeable boundary condition, ideally, there should no pressure build-up in areas of high permeability. But if this area of high permeability is small and is surrounded by materials of very low permeability, this area will have pressure build-up: fluid can flow in but cannot leak off. A simple case to illustrate this phenomenon is a fracture created by hydraulic fracturing: at the tip of the fracture, which is of high permeability, there is pressure build-up. The interaction of these two mechanisms intensifies the phenomenon of fracture networking generated by injection.

Figure 11 shows contour of shear strain intensity ${\gamma }_T$ at two moments. It can be seen that the value of ${\gamma }_T$ reduces in the center with development of fracture-networking to the outer field. Boundary effects are not significant for these sets of numerical suctions with the current scale of 200 m, but they seem to be significant in the results shown in Figure 6, Figure 7 and Figure 8. For the case in which the principal directions of natural fractures are not the same as that of principal stresses, its scale should be larger than 200 m.

4. Discussion

Comparing Figure 6, Figure 7 and Figure 8 with Figure 9, Figure 10 and Figure 11, the changes in results when there is a difference between the principal directions of natural fractures and the geostresses can be seen:

(1) For the results of model with difference between principal direction of natural fracture and stress components, as shown in Figure 6, Figure 7 and Figure 8, the profile of pore pressure distribution kinks off from its original direction of propagation. At the beginning of the injection process, the zones where higher values of pore pressure were generated in the principal direction of natural fractures. But in the far field to the injection point, direction of zones of higher pore pressure values rotates, as shown in Figure 6. For pore pressure profile shown in Figure 9 which is result of model which natural fractures have the same direction as that of principal stresses, it propagates smoothly to the far field without kinking.

(2) Similarly, the localization band of damage variable shown in Figure 7 also rotates in the far field but that shown in Figure 11 does not rotate. Porous flow develops in the direction of higher permeability, which is the principal direction of natural fractures. Due to the higher value of pore pressure along this path of porous flow development, the magnitude of effective stress in this area is smaller than that in other area. This makes the damage evolution process easier to develop. Usually, hydraulic fracturing-induced-fractures develop in the direction of maximum horizontal stress. This is true for the case when principal directions of natural fractures are the same as that of the geostress tensor. In the case when there is a difference between principal directions of natural fractures and that of the geostress tensor, this statement will not be true. The reason is that the principal direction of the permeability tensor is determined by natural fractures via the damage tensor and it further dominates the behavior of porous flow. Therefore, at the beginning of the hydraulic fracturing process, in the direction of natural fractures, which is also the principal direction of permeability tensor, porous flow develops faster than that in the other directions. Consequently, the induced fractures will not develop in the direction of maximum horizontal stress. Rather, they will develop in the principal direction of natural fractures. This is true until the impact of induced fractures overcomes that of natural fractures.

The following remarks can be made on the accuracy of numerical results of fractures near the wellbore and in the far field. Fracture development in the hydraulic fracturing process can be divided into two types: near-wellbore fractures (within 15 ft (approximately 4.5 m) from the borehole) and far-field fractures (15 ft from the borehole) [14,15]. Because of the stress concentration and other factors, such as drilling, near-wellbore fractures are rather complicated: mode-I, mode-II, and mode-II cracks are all possible, and their path of development is most likely in the form of a curve or 3D surface. A detailed description of the near-wellbore fracture, however, is not as important to the final formation stimulation result as that of the far-field fracture. The production quantity is more affected by resultant far-field fractures, rather that the near-wellbore fractures.

Therefore, the numerical results of fracture development presented in previous Section 3.3.3 and Section 3.3.4 are obtained by using a model without a detailed borehole geometry description. The reason for neglecting the wellbore geometry and near-wellbore behaviors of fracture propagation is that the stimulated reservoir volume (SRV) is mainly determined by far field behavior and stress distribution around wellbore has little influence on the fracturing behaviors in far field.

In laboratory environments, fracture testing can only focus on the near-wellbore crack. This makes the lab testing of hydraulic fracturing less useful to production prediction than that of the numerical results. To further understand the fracture development under hydraulic fracture in the far field, Figure 12 and Figure 13 provide contours of damage intensity $\overline{\omega }$ at time moment t = 36.64 and 46.64 min, respectively. Plane view is used in these two figures because the focus is the fracture development in the horizontal plane.

Because the x-direction is the direction of minimum horizontal stress, the fracture network represented by the dark area in the middle of the model formed by a localized band of damage intensity $\overline{\omega }$ develops more in the x-direction than in the y-direction. Consequently, it forms an area of fractures in the elliptical form, as shown by the dashed circle in Figure 12. Those small dark areas near the border in y-direction are considered as results of boundary effect and thus are regarded as not true. The reason for this is that a real fractured area should develop continuously from center to border.

The contour of damage intensity $\overline{\omega }$ at time moment t = 46.64 min shown in Figure 13 illustrates the continuous development of the fracture network formed by a localized band of damage intensity $\overline{\omega }$ The area of fractures in the elliptical form as the dashed curve area develops on the basis of the one shown in Figure 12. The dark area connected to the border in the y-direction in Figure 13 is considered to be a result of the boundary effect (i.e., the simplified displacement constraints at border of the model causes abnormal stress at these boundary nodes) and thus is regarded as not true.

5. Conclusions

This paper introduces the process for modeling hydraulic fracturing of reservoirs with multiple sets of natural fractures applying the framework of continuum-damage mechanics. It also demonstrates that the coupled finite-element hydromechanical model is an effective tool for estimating the volume stimulated by injection in tightsand reservoirs.

There is an initial damage tensor for each set of natural fractures. The total initial damage tensor can be calculated by a simple summation of those for each single natural fracture set. Because the initial principal damage value is significantly less than 1.0 (e.g., 0.015), this simple summation of the damage tensor calculation is proved practical for modeling initial natural fractures. The example presented shows a 73° angle between the direction of minimum horizontal stress and the principal direction of the damage tensor and the orthotropic permeability tensor. In practice, horizontal well trajectory is usually placed in the direction of minimum horizontal stress in the purpose to have a fracture perpendicular to the well trajectory and thus maximize the stimulated reservoir volume. When there is this difference, the fracture may propagate in a direction different from what is expected. In this case, in order to have the best/optimized fracturing, the best direction for the well trajectory of horizontal sections should be adjusted accordingly.

The concept of damage intensity $\overline{\omega }$ is introduced to illustrate the intensity of fracturing generated by hydraulic injection. It does not have any independent mechanical meaning. The numerical solution of pore pressure contour and damage contour generated by injection are presented. The contour of synthetic damage forms a complex fracture network in the near field of the injection point. The direction of the major damage distribution band turns from the direction of minimum horizontal stress to the direction of natural fractures. This indicates that the directions of principal stresses and natural fractures should be considered and accounted for in related calculations to optimize hydraulic fracturing design.