Abstract

Hydraulic fracturing is widely used to determine in situ stress of rock engineering. In this paper we propose a new method for simultaneously determining the in situ stress and elastic parameters of rock. The method utilizing the hydraulic fracturing numerical model and a computational intelligent method is proposed and verified. The hydraulic fracturing numerical model provides the samples which include borehole pressure, in situ stress, and elastic parameters. A computational intelligent method is applied in back analysis. A multioutput support vector machine is used to map the complex, nonlinear relationship between the in situ stress, elastic parameters, and borehole pressure. The artificial bee colony algorithm is applied in back analysis to find the optimal in situ stress and elastic parameters. The in situ stress is determined using the proposed method and the results are compared with those of the classic breakdown formula. The proposed method provides a good estimate of the relationship between the in situ stress and borehole pressure and predicts the maximum horizontal in situ stress with high precision while considering the influence of pore pressure without the need to estimate Biot’s coefficient and other parameters.

1. Introduction

Hydraulic fracturing is widely used in the recovery of oil, gas, geothermal, and mineral resources [1]. In petroleum engineering it is important to determine the in situ stresses and elastic parameters of the rock mass when using hydraulic fracturing in fracturing operations, wellbore stability analysis, and reservoir simulation [2]. While high accuracy is required for the values of the in situ stress and mechanical parameters of the rock mass, determination of these parameters is still one of the most challenging tasks in hydraulic fracturing.

Hydraulic fracturing tests are considered the most effective method for determining the in situ stress and mechanical parameters of rock mass [39]. The Hubbert and Willis hydraulic fracturing criterion and Haimson and Fairhurst’s hydraulic fracturing criterion are the two classic formulae for hydrofracture breakdown ([10]; Hubbert et al., 1953). However, the pore pressure term, which is a significant factor in deep boreholes, is ignored in Hubbert and Willis’s hydraulic fracturing criterion. Modifications of the original equations were proposed to account for the pore pressure (Detournay et al., 1988; [1, 1113]), but they have not been used in practice because it involves the Biot poroelastic parameters and Poisson’s ratio which are difficult to obtain. Schmitt and Zoback built a more useful generalized form of the hydrofracture breakdown equation by considering the poroelastic effects [1]. It can be used to provide upper and lower bound to the maximum horizontal in situ stress because it depends on the specific pore and microcrack structure. However, this method requires the poroelastic coefficients which are difficult to determine in practice.

Owing to the limitations of the classic breakdown formulae and the complexity of hydraulic fracturing tests, laboratory and field tests have been commonly used to determine the in situ stress and mechanical parameters of rock mass (Algorithm 1). However, these tests may not always produce the poroelastic parameters or may provide inaccurate results because of low-quality core samples [14]. Alternatively, back analysis, associated with the “in situ” approach, has been widely used to determine the mechanical parameters of rock mass in rock engineering [1520]. Zhang and Yin proposed a back analysis method which combined a neural network and a genetic algorithm to simultaneously identify the in situ stresses and elastic parameters [2]; however, this method did not consider the poroelastic effect. To overcome this difficulty, in this paper we extend our proposed displacement back analysis method to determine the in situ stress and mechanical parameters of a rock mass based on measured borehole pressure. The borehole pressure can be easily measured in the field with a pressure gauge installed inside the borehole [21]. Back analysis is implemented following an optimization strategy based on the multioutput support vector machine (MSVM) and artificial bee colony algorithm (ABC) model, which is effective in multiple parameter identification [22].

Sub MSVM()
  
Dim N As Integer                The number of training samples
Dim Dim As Integer     The dimension of input variables
Dim Dim As Integer     The number of output variables
  
Dim C As Double               Pentalty factor of SVM
Dim epsilon As Double
Dim sigma As Double
  
Dim input() As Double         The input of training samples
Dim output() As Double       The output of traning samples
  
  
The one sample
Dim     As Double
Dim     As Double
The weights of MSVM
Dim As Double
Dim As Double
Dim As Double
Dim As Double
Dim As Double
Dim As Double
the error of each sample
Dim As Double
Dim As Double
Dim As Double
Dim new() As Double
  
The coefficient of matrix for computing and
Dim As Double
Dim BB() As Double
Dim lastrow() As Double
Dim lastrow() As Double
  
The Descending direction
Dim () As Double
Dim () As Double
  
Dim () As Double
Dim () As Double
Dim 1() As Double
Dim kf() As Double
  
Dim As Integer
Dim As Integer
Dim As Integer
  
The step size eta
Dim eta As Double
  
The control parameters of algorithm convergence                      
Dim delta As Double
Dim As Integer
  
The value of the numbers of parameters of training samples                  
= Range().Cells.Value
Dim = Range(Dim).Cells.Value
Dim = Range(Dim).Cells.Value
  
The parameters of SVM
= Range().Cells.Value
epsilon = Range(epsilon).Cells.Value
sigma = Range(sigma).Cells.Value
   
   
ReDim input(1 To , 1 To Dim) As Double
ReDim output(1 To , 1 To Dim) As Double
ReDim (1 To Dim) As Double
ReDim (1 To Dim) As Double
ReDim (1 To Dim, 1 To ) As Double
ReDim (1 To Dim) As Double
ReDim (1 To Dim, 1 To ) As Double
ReDim (1 To Dim) As Double
ReDim (1 To Dim, 1 To ) As Double
ReDim (1 To Dim) As Double
ReDim (1 To ) As Double
ReDim (1 To ) As Double
ReDim (1 To Dim) As Double
ReDim new(1 To ) As Double
ReDim (1 To , 1 To ) As Double
ReDim (1 To , 1 To Dim) As Double
ReDim BB(1 To , 1 To Dim) As Double
ReDim lastrow(1 To ) As Double
ReDim lastrow(1 To Dim) As Double
ReDim (1 To Dim, 1 To ) As Double
ReDim (1 To Dim) As Double
ReDim (1 To , 1 To ) As Double
ReDim (1 To , 1 To ) As Double
ReDim 1(1 To , 1 To ) As Double
ReDim (1 To ) As Double
   
Read the input of training samples
For To
For To Dim
input = Range().Cells
Next
Next
   
Read the output of training samples
For To
For To Dim
output = Range().Cells.Value
Next
Next
   
The initial value of and
For To Dim
For To
Next
Next
   
   
delta = 1
   
The iteration process of algorithm
While (delta > 0.001 And )
Replace the value of and by the New and
For To Dim
For To
Next
Next
Compute the value of and
For To
For To Dim
= input
Next
   
For ii = 1 To
For To Dim
= input(ii, j)
Next
(ii) = kernelfun
Next
For ii = 1 To Dim
(ii) = (ii)
For To
(ii) = (ii) +
Next
Next
   
For To Dim
Next
= Sqr
If (u(I) < epsilon) Then = 0
If Then ai(I) =
Next
compute the Matrix and Da
For To
For To
For To Dimx
input
input
Next
= kernelfun
If Then Else = 0
Next
Next
For To
For To
If Then Else                
Next
Next
Compute Transpose(a)
For To
lastrow
For To
lastrowlastrow
Next
Next
For To
lastrow
Next
For To Dim
lastrow
For To
lastrowlastrowoutput
Next
Next
For To
For To Dim
BBoutput
Next
Next
For To Dim
BB = lastrow
Next
Compute and
With Application.WorksheetFunction                        
  =  .MMult(.MInverse(), BB)
End With
   
   
For To Dim
For To
=
Next
=
Next
Compute the descending direction
For To Dim
For To
=
Next
Next
eta = 1
Dim deltaLp As Double
Dim Lpk1 As Double
Dim Lpk As Double
deltaLp = 1
Update the solution of and
While (delta)
For To Dim
For To
Next
Next
For To
new
For To Dim
= input
Next
For To
For To Dim
= input
Next
kf(ii) = kernelfun
Next
For To Dim
For To
Next
Next
   
For To Dim
newnewoutput
Next
new =
If (new() < epsilon) Then
If (new >= epsilon) Then newnew()
delta = delta + new()
Next
For To Dim
For To
1 = 1 +
= + 2/2
Next
Next
For To
If new >= epsilon Then 1 = 1 + newnew)
If Then =
Next
delta = 1 -
eta = 0.5 eta
Wend
   
delta
For To
delta = delta + new()
Next
delta = delta
   
   
Wend
For To
For To Dim
Range().Cells =
Next
Next
   
   
For To Dim
Range(bi).Cells
Next
  
End Sub
Kernel function of RBF
Function kernelfun(xx, yy, sigma2) As Double                      
Dim temp As Double
Dim temp1 As Double
Dim Dim As Integer
   
Dim = Range(Dim).Cells.Value
temp = 0
For To Dim
temp = temp +
temp = temp +
Next I
   
temp1 = Sqr(temp)/()
  
kernelfun =
   
End Function
   
   
Compoute the performance    function value using the MSVM
Sub Perffunc()
Dim As Integer
Dim As Integer
Dim As Integer
Dim As Integer
Dim As Integer
Dim As Integer
Dim sigma As Double
   
Dim As Double
Dim As Double
Dim As Double
Dim As Double
Dim As Double
Dim As Double
  
= Range().Cells.Value
= Range().Cells.Value
Dim = Range(Dim).Cells.Value                          
Dim = Range(Dim).Cells.Value
sigma = Range(sigma).Cells.Value
   
ReDim (1 To Dim) As Double
ReDim (1 To Dimx) As Double
ReDim (1 To Dim, 1 To ) As Double
ReDim (1 To Dim) As Double
ReDim (1 To , 1 To Dim) As Double
ReDim kf(1 To ) As Double
   
   
For To Dim
For To
= Range().Cells
Next
= Range().Cells
Next
    
   
For To
   
For To Dim
= Range(input).Cells
Next
For To
For To Dim
= Range().Cells
Next
kf = kernelfun
Next
For To Dim
For To
Next
Next
Next
   
For To
For To Dim
Range(output).Cells
Next
Next
   
End Sub

The rest of this paper is organized as follows. The classic breakdown formula is presented in detail in Section 2. The formulation and procedure of back analysis based on borehole pressure are presented in detail in Section 3. In Section 4, a numerical example is used to verify the proposed method, and our conclusions are presented in Section 5.

2. Hydraulic Breakdown Equations

Hydraulic fracturing is a widely accepted technology used for determining in situ stress magnitude and direction. The principal stress has a magnitude equal to the overburden pressure in the vertical direction. The smallest horizontal principal stress is usually determined directly in the experiment from the shut-in pressure. The greatest horizontal principle stress must be calculated using a breakdown formula derived from an appropriate hydraulic fracturing model. Hubbert and Willis proposed a classic breakdown formula (1) to calculate for hydraulic fracturing in nonporous impermeable rocks [23], ignoring the pore pressure term.where is the rock tensile strength.

Equations (2) and (3) are the breakdown formulae of porous impermeable rocks and porous permeable rocks, respectively, including the pore pressure [24]:where is the breakdown pressure, is the pore pressure, is the Biot poroelastic parameter, and is Poisson’s ratio. Although (3) may best describe the conditions under which hydraulic fracturing is conducted from an open borehole, (2) is used in practice because of the difficulty of determining and . Schmitt and Zoback proposed a more generalized form for the equations of hydrofracture breakdown for porous impermeable rocks and porous permeable rocks [1]:where β is the poroelastic effect parameter.

3. Back Analysis Model Based on Borehole Pressure

The in situ stress can be estimated based on borehole pressure using the hydraulic breakdown equations (Section 2). However, these equations present some limitations in practice. Therefore, we propose a back analysis method that combines a numerical method and an intelligent computational method. A multioutput support vector machine (MSVM) is used to map the complex, nonlinear relationship between the in situ stress, elastic parameters, and borehole pressure. The numerical method provides the training samples for the MSVM. It is important to use an optimization method in back analysis. Here, we use the ABC algorithm to find the best-fit in situ stress and elastic parameters by comparing the measured pressure data and the MSVM predicted pressure.

3.1. Nonlinear Relationships between Pressure and Geomechanical Parameters

The relationship between the borehole pressure and geomechanical parameters can be derived by the MSVM. The basic idea of MSVM is to extend the single-output support vector machine to a multidimensional output case. Given a set of training samples , , , the MSVM model can be built by solving the following optimization problem based on an iterative reweighted least-square algorithm [25].where is the number of input, is the number of output, is the weight, and are constants, is the sum of constant terms that do not depend on either or , is the tolerant error, and denotes the tth iteration. A brief description and the MSVM algorithm can be found in the literature [22]. The MSVM model can be expressed as

Based on the above MSVM model, the nonlinear relationship between the borehole pressure and geomechanical parameters can be described aswhere is the -dimensional vector of the identified parameter, for example, the in situ stress, Young’s modulus, or Poisson’s ratio. is the -dimensional vector of the borehole pressure.

To build the MSVM model, the necessary training or learning samples are constructed and the MSVM parameters are determined. The samples are constructed by numerical analysis which computes the corresponding borehole pressure for a given set of tentative determined parameters. The MSVM parameters have a strong influence on the performance of the MSVM. In this study, we determined these parameters using the formulation presented by Meza et al. [26].

3.2. Optimization Method

The back analysis ABC algorithm, developed by Karaboga [27], was adopted to search for the optimal geomechanical parameters of the rock mass. In the algorithm, the colony of artificial bees consists of three groups: employed bees, onlookers, and scouts. The ABC algorithm involves a cycle of four phases: the initialization phase, employed bees phase, onlooker bee phase, and scout bee phase.

In the initialization phase, the ABC generates a randomly distributed initial population of SN solutions and calculates the fitness of each solution. where is the candidate solution of the problem; and denotes the size of the population; and is the dimension number of each solution; rand is a random number between ; and are the upper and lower bounds of each solution.

Once initialization is completed, the employed bees search for a solution and calculate the fitness value (see Section 3.3) in the employed bees phase. A candidate solution is produced according to the following equation:where is different from and is a randomly chosen index from , is also an index randomly chosen from , and is a random number in the range that controls the generation of food sources around and represents the comparison of two food positions seen by a bee.

In the onlooker bee phase, the onlooker bees choose a solution based on the fitness value, determine which solution will be abandoned, and allocate its employed bees as scout bees. The probability of being selected for each fitness value can be expressed aswhere is the fitness value of the solution.

Finally, in the scout bee phase the scout bees randomly search for a new solution in the determined ranges. A solution that cannot be improved further through a predetermined number of cycles is assumed to be abandoned by the onlookers.

3.3. The Fitness Function

In order to find the optimal solution, it is necessary to build the fitness function for the ABC algorithm; that is,where is the predicted pressure using the MSVM model, Y is the vector of the monitored pressure, and is the number of monitored points.

3.4. Procedure of the MSVM-ABC Based Method

If the MSVM model can establish the nonlinear relation between the borehole pressures and determined parameters, the model can be used to predict the borehole pressures. Then, the ABC algorithm is utilized to find the optimal parameters through error minimization between the pressures predicted by the MSVM model and the measured pressures. The back analysis flowchart is shown in Figure 1 and the algorithm is described as follows.

Step 1. Determine the general information and data such as the unknown (determined by back analysis) and known parameters of the numerical model, the MSVM and ABC algorithm parameters, and the range of parameters to be determined.

Step 2. Generate the combination of determined parameters, calculate the borehole pressure for each combination, and then build the learning samples for MSVM.

Step 3. Based on the learning samples of Step 2, construct the MSVM model using the MSVM algorithm and activate the ABC algorithm.

Step 4. Search for the optimal determined parameters based on the monitored pressure.

4. Validation and Application

To verify the proposed method, a numerical example is adopted to determine the in situ stress and elastic mechanical parameters of the elastic rock. The numerical experiment is conducted based on 2D hydraulic fracturing model of water injection into a hypothetical deep formation. Further details of the physical and numerical model can be found in the literature published by Schmitt and Zoback [1].

The parameters to be determined are the maximum and minimum horizontal in situ stress and , respectively, Young’s modulus E, and Poisson’s ratio . The rock mass in all the zones is considered to be elastic. The mechanical parameters of the joints and the permeability of the rock mass are known; the parameter values can be seen in the literature published by Schmitt and Zoback [1]. Thirty sets of training samples and ten testing samples derived in previous studies [1, 2] were selected. Based on the MSVM algorithm, the MSVM code was written in Excel and VBA. The MSVM parameters and some of the weight and constant values and samples are shown in Figure 2. Good agreement between the measured data and the pressures estimated by the MSVM is shown in Figure 3, indicating the good performance of the MSVM model. Thus, the proposed model can accurately estimate the borehole pressures, replacing the existing numerical analysis method for calculating borehole pressures. The results also confirm that the MSVM model provides an accurate representation of the nonlinear relationship between the pressures and the determined parameters.

The ABC code is also written in Excel and VBA. The parameters of the ABC algorithm and the calculation results are shown in Figure 4. Based on the proposed method for determining the in situ stress and mechanical parameters of rock mass, the results are shown in Table 1. we obtained the values of , , , and as 24.46 MPa, 14.33 MPa, 44.02 GPa, and 0.25, respectively. The elastic mechanical parameters of the rock agree with the results calculated by the Genetic Algorithm-Neural Network (ANN-GA) [2]. A comparison of in situ stress values calculated using four different formulations is shown in Figure 5. agrees well with the value estimated by ANN-GA, (1), and (2). The relative error is only 1.07%. agrees well with the value estimated by (2), which considers the pore pressure, but differs considerably from the values calculated by (1) and ANN-GA which do not consider the poroelastic effects. The relative error is up to 31.8%. Using (4), we obtain the upper and lower limits of the maximum horizontal in situ stress (14.95–34.55 MPa). The value 24.46 MPa is within this range. Thus, the proposed method can be used in back analysis as an alternative numerical analysis method, which considers the poroelastic effects and provides rational, high-precision results. Note that the proposed method can determine the maximum horizontal in situ stress without estimating the poroelastic coefficient, which is a difficult parameter to obtain.

Moreover, there are four borehole pressures, namely, formation breakdown pressure (FBP) , fracture propagation pressure (FPP) , instantaneous shut-in pressure (ISIP) , and leak–off pressure (LOP) . A comparison of back analysis on borehole pressures obtained by three different methods is presented in Figure 6. The borehole pressure calculated by the proposed MSVM method is very close to the measured pressure. The relative error is less than 3%. On the other hand, the convergence processes of the algorithm and fitness variations are shown in Figures 7 and 8. Initially, the data was distributed randomly in the searching space (Figure 8) and then converged to the solution of the problem in the 500th generation. This indicates that the proposed method can determine both the in situ stress and elastic mechanical parameters of the elastic rock with excellent converging performance.

5. Conclusions

In this paper, a new borehole pressure-based back analysis approach to determine the stress and mechanical parameters of rock mass is proposed. The method combines a coupling numerical model of hydraulic fracturing and a computational intelligent method. The method is applied to a numerical example to successfully determine both the in situ stress and mechanical parameters of a rock mass. In this approach, the MSVM is adopted to represent the nonlinear relationship between the borehole pressure and mechanical parameters of the rock mass, proving more efficient than existing numerical models. The ABC algorithm is used to search for the optimal parameters in the search space. The proposed approach is implemented in Excel with VBA.

In the classic breakdown formula, it is difficult in practice to determine the maximum horizontal in situ stress while considering the poroelastic coefficient. The proposed back analysis method can predict the maximum horizontal in situ stress based on the borehole pressures without the need to obtain the poroelastic coefficient. Thus, it is a more practical method for determining the in situ stress from hydraulic fracturing. The proposed method is practical and accurate and can be conveniently applied to simultaneously determine the in situ stress and mechanical parameters of rock from hydraulic fracturing.

Symbols

:Principal stress in direction
:The greatest horizontal principle stress
:The smallest horizontal principal stress
:Poisson’s ratio
:Poroelastic effect parameter
:The weight vector
MSVM():The predicted pressure using the MSVM model
:Hyper parameter that determines trade-off between the regularization and the error reduction term
:Tolerant error
:Number of output
:A random number in the range
:A constant for classification threshold
:Breakdown pressure
:Pore pressure
:Young’s modulus
:Biot poroelastic parameter
:Rock tensile strength
:Number of input
:Sum of constant terms that do not depend on either or
rand:A random number between
fitness:Fitness value of the solution.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors gratefully acknowledge the National Natural Science Foundation of China (Grant no. 41372315), the State Key Research Development Program of China (Grant no. 2016YFC0600702), and Innovative Research Team (in Science and Technology) in University of Henan Province (no. 15IRTSTHN029).