# Back Analysis of Rock Hydraulic Fracturing by Coupling Numerical Model and Computational Intelligent Technology

Authors: Shaojun Li; Hongbo Zhao; Xiating Feng; Quan Jiang; Qiancheng Sun
Journal: Geofluids
 (2017)

---

## 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.


---

## Body

## 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 [3–9]. 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, 11–13]), 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 (Algorithm1). 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 [15–20]. 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].Algorithm 1:The code of MSVM.
Sub MSVM()
Dim N As Integer＇ The number of training samples
Dim Dim_x As Integer    ＇ The dimension of input variables
Dim Dim_y As Integer    ＇ The number of output variables
Dim C As Double＇ Pentalty factor of SVM
Dim epsilon As Double
Dim sigma As Double
DimX_input() As Double        ＇ The input of training samples
DimY_output() As Double      ＇ The output of traning samples
＇ The one sample
Dimxi()    As Double
Dimxj()    As Double
＇The weights of MSVM
DimW() As Double
Dimb() As Double
DimW_k() As Double
Dimb_k() As Double
DimW_s() As Double
Dimb_s() As Double
＇ the error of each sample
Dimu() As Double
Dimai() As Double
Dimui() As Double
Dimu_new() As Double
＇The coefficient of matrix for computing Ws and bs
DimA() As Double
Dim BB() As Double
DimA_last_row() As Double
DimB_last_row() As Double
＇The Descending direction
DimP_W() As Double
DimP_b() As Double
Dimk() As Double
DimD_a() As Double
DimD_a_1() As Double
Dim kf() As Double
DimI As Integer
Dimj As Integer
Diml As Integer
＇The step size eta_k
Dim eta As Double
＇The control parameters of algorithm convergence
Dim delta_u As Double
Dimt As Integer
＇The value of the numbers of parameters of training samples
N = Range(＇＇N＇＇).Cells.Value
Dim_x = Range(＇＇Dim_x＇＇).Cells.Value
Dim_y = Range(＇＇Dim_y＇＇).Cells.Value
＇The parameters of SVM
C = Range(＇＇_C＇＇).Cells.Value
epsilon = Range(＇＇epsilon＇＇).Cells.Value
sigma = Range(＇＇sigma＇＇).Cells.Value
ReDimX_input(1 To N, 1 To Dim_x) As Double
ReDimY_output(1 To N, 1 To Dim_y) As Double
ReDimxi(1 To Dim_x) As Double
ReDimxj(1 To Dim_x) As Double
ReDimW(1 To Dim_y, 1 To N) As Double
ReDimb(1 To Dim_y) As Double
ReDimW_k(1 To Dim_y, 1 To N) As Double
ReDimb_k(1 To Dim_y) As Double
ReDimW_s(1 To Dim_y, 1 To N) As Double
ReDimb_s(1 To Dim_y) As Double
ReDimu(1 To N) As Double
ReDimai(1 To N) As Double
ReDimui(1 To Dim_y) As Double
ReDimu_new(1 To N) As Double
ReDimA(1 To N+1, 1 To N+1) As Double
＇ReDim X(1 To N+1, 1 To Dim_y) As Double
ReDim BB(1 ToN+1, 1 To Dim_y) As Double
ReDimA_last_row(1 To N) As Double
ReDimB_last_row(1 To Dim_y) As Double
ReDimP_W(1 To Dim_y, 1 To N) As Double
ReDimP_b(1 To Dim_y) As Double
ReDimk(1 To N, 1 To N) As Double
ReDimD_a(1 To N, 1 To N) As Double
ReDimD_a_1(1 To N, 1 To N) As Double
ReDimkf(1 To N) As Double
＇Read the input of training samples
ForI=1 To N
For j=1 To Dim_x
X
_input(I,j) = Range(＇＇xi＇＇).Cells(I,j)
Next
Next
＇Read the output of training samples
ForI=1 To N
For j=1 To Dim_y
Y
_output(I,j) = Range(＇＇yi＇＇).Cells(I,j).Value
Next
Next
＇ The initial value of Wk and bk
ForI=1 To Dim_y
b
(
I
)
=
0
For j=1 To N
W
(
I
,
j
)
=
0
Next
Next
delta_u = 1
t=0
＇The iteration process of algorithm
While (delta_u > 0.001 And t<100)
＇ Replace the value of Wk and bk by the New w and b
For I=1 To Dim_y
For j=1 To N
W
_
k
(
I
,
j
)
=
W
(
I
,
j
)
Next
b
_
k
(
I
)
=
b
(
I
)
Next
＇ Compute the value of ui and ai
For I=1 To N
u
(
I
)
=
0
For j=1 To Dim_x
x
i
(
j
) = X_input(I,j)
Next
For ii = 1 To N
For j=1 To Dim_x
x
j
(
j
) = X_input(ii, j)
Next
k
f(ii) = kernel_fun(xi,xj,sigma)
Next
For ii = 1 To Dim_y
u
i(ii) = b_k(ii)
For j=1 To N
u
i(ii) = ui(ii) + W_k(ii,j)∗kf(j)
Next
Next
For j=1 To Dim_y
u
(
I
)
=
u
(
I
)
+
(
Y
_
o
u
t
p
u
t
(
I
,
j
)
-
u
i
(
j
)
)
∧
2
Next
u
(
I
) = Sqr(u(I))
If (u(I) < epsilon) Then ai(I) = 0
If (u(I)>=epsilon) Then ai(I) = 2∗C∗(u(I)-epsilon)/u(I)
＇
a
i
(
i
)
=
2
∗
C
∗
(
u
(
i
)
-
e
p
s
i
l
o
n
)
/
u
(
i
)
Next
＇compute the Matrix K and Da
For I=1 To N
For j=1 To N
For l=1 To Dim_x
x
i
(
l
)
=
X
_input(I,l)
x
j
(
l
)
=
X
_input(j,l)
Next
k
(
I
,
j
) = kernel_fun(xi,xj,sigma)
If (I=j) Then D_a(I,j)=ai(I) Else D_a(I,j) = 0
Next
Next
For I=1 To N
For j=1 To N
If (I=j) Then A(I,j)=k(I,j)+1/D_a(I,j) Else A(I,j)=k(I,j)
Next
Next
＇Compute Transpose(a)∗K
For I=1 To N
A
_last_row(I)=0
For j=1 To N
A
_last_row(I)=A_last_row(I)+ai(j)∗k(j,I)
Next
Next
A
(
N
+
1
,
N
+
1
)
=
0
For I=1 To N
A
(
N
+
1
,
I
)
=
A
_last_row(I)
A
(
N
+
1
,
N
+
1
)
=
A
(
N
+
1
,
N
+
1
)
+
a
i
(
I
)
A
(
I
,
N
+
1
)
=
1
Next
For I=1 To Dim_y
B
_last_row(I)=0
For j=1 To N
B
_last_row(I)=B_last_row(I)+ai(j)∗Y_output(j,I)
Next
Next
For I=1 To N
For j=1 To Dim_y
BB(I,j)=Y_output(I,j)
Next
Next
For j=1 To Dim_y
BB(N+1,j) = B_last_row(j)
Next
＇Compute Ws and bs
With Application.WorksheetFunction
x
1  =  .MMult(.MInverse(A), BB)
End With
For I=1 To Dim_y
For j=1 To N
W
_
s
(
I
,
j
) = x1(j,I)
Next
b
_
s
(
I
) = x1(N+1,I)
Next
＇Compute the descending direction
For I=1 To Dim_y
For j=1 To N
P
_
W
(
I
,
j
) = W_s(I,j)-W_k(I,j)
Next
P
_
b
(
I
)
=
b
_
s
(
I
)
-
b
_
k
(
I
)
Next
eta = 1
Dim delta_Lp As Double
Dim Lp_k_1 As Double
Dim Lp_k As Double
delta_Lp = 1
＇Update the solution of W and b
While (delta_Lp>0.0001)
For I=1 To Dim_y
For j=1 To N
W
(
I
,
j
)
=
(
1
-
e
t
a
)
∗
W
_
k
(
I
,
j
)
+
e
t
a
∗
P
_
W
(
I
,
j
)
Next
b
(
I
)
=
(
1
-
e
t
a
)
∗
b
_
k
(
I
)
+
e
t
a
∗
P
_
b
(
I
)
Next
For I=1 To N
u
_new(I)=0
For j=1 To Dim_x
x
i
(
j
) = X_input(I,j)
Next
For ii=1 To N
For j=1 To Dim_x
x
j
(
j
) = X_input(ii,j)
Next
kf(ii) = kernel_fun(xi,xj,sigma)
Next
For ii=1 To Dim_y
u
i
(
i
i
)
=
b
(
i
i
)
For j=1 To N
u
i
(
i
i
)
=
u
i
(
i
i
)
+
W
(
i
i
,
j
)
∗
k
f
(
j
)
Next
Next
For j=1 To Dim_y
u
_new(I)=u_new(I)+(Y_output(I,j)-ui(j))∧2
Next
u
_new(I) = Sqr(u_new(I))
If (u_new(I) < epsilon) Then ai(I)=0
If (u_new(I) >= epsilon) Then ai(I)=2∗C∗(u_new(I)-epsilon)/u_new(I)
delta_u = delta_u + u_new(I)
Next
L
p
_
k
_
1
=
0
L
p
_
k
=
0
For I=1 To Dim_y
For j=1 To N
L
p
_
k
_1 = Lp_k_1 + WI,j∧2/2
L
p
_
k = Lp_k + W_kI,j∧2/2
Next
Next
For j=1 To N
If u_new(j) >= epsilon Then Lp_k_1 = Lp_k_1 + C∗(u_newj∧2-2∗u_new(j)∗epsilon+epsilon∧2)
If u(j)>=epsilon Then Lp_k = Lp_k+C∗(uj∧2-2∗u(j)∗epsilon+epsilon∧2)
Next
delta_Lp = Lp_k_1 - Lp_k
eta = 0.5 ∗ eta
Wend
delta_u=0
For I=1 To N
delta_u = delta_u + u_new(I)
Next
delta_u = delta_u/N
t
=
t
+
1
Wend
ForI=1 To N
For j=1 To Dim_y
Range(＇＇wi＇＇).Cells(I,j) = W(j,I)
Next
Next
Forj=1 To Dim_y
Range(＇＇bi＇＇).Cells(j)=b(j)
Next
End Sub
＇ Kernel function of RBF
Function kernel_fun(xx, yy, sigma2) As Double
Dim temp As Double
Dim temp1 As Double
Dim Dim_x As Integer
Dim_x = Range(＇＇Dim_x＇＇).Cells.Value
temp = 0
ForI=1 To Dim_x
temp = temp +xxI-yyI∧2
＇temp = temp + xx(i)∗yy(i)
Next I
temp1 = Sqr(temp)/(2∗sigma2∧2)
kernel_fun = Exp(-temp1)
＇Kf=temp+1∧sigma2
End Function
＇ Compoute the performance    function value using the MSVM
Sub Perffunc()
DimI As Integer
Dimii As Integer
Diml As Integer
Dimj As Integer
DimN As Integer
DimNp As Integer
Dim sigma As Double
Dimxi() As Double
Dimxj() As Double
DimW() As Double
Dimb() As Double
DimYy_p() As Double
Dimkf() As Double
N = Range(＇＇N＇＇).Cells.Value
N
p = Range(＇＇Np＇＇).Cells.Value
Dim_x = Range(＇＇Dim_x＇＇).Cells.Value
Dim_y = Range(＇＇Dim_y＇＇).Cells.Value
sigma = Range(＇＇sigma＇＇).Cells.Value
ReDimxi(1 To Dim_x) As Double
ReDimxj(1 To Dim_x) As Double
ReDimW(1 To Dim_y, 1 To N) As Double
ReDimb(1 To Dim_y) As Double
ReDimYy_p(1 To Np, 1 To Dim_y) As Double
ReDim kf(1 ToN) As Double
ForI=1 To Dim_y
For j=1 To N
W
(
I
,
j
) = Range(＇＇wi＇＇).Cells(j,I)
Next
b
(
I
) = Range(＇＇bi＇＇).Cells(I)
Next
ForI=1 To Np
For j=1 To Dim_x
x
i
(
j
) = Range(＇＇x_p_input＇＇).Cells(I,j)
Next
For l=1 To N
For j=1 To Dim_x
x
j
(
j
) = Range(＇＇xi＇＇).Cells(l,j)
Next
kf(l) = kernel_fun(xi,xj,sigma)
Next
For ii=1 To Dim_y
Y
y
_
p
(
I
,
i
i
)
=
b
(
i
i
)
For j=1 To N
Y
y
_
p
(
I
,
i
i
)
=
Y
y
_
p
(
I
,
i
i
)
+
W
(
i
i
,
j
)
∗
k
f
(
j
)
Next
Next
Next
ForI=1 To Np
For j=1 To Dim_y
Range(＇＇y_p_output＇＇).Cells(I,j)=Yy_p(I,j)
Next
Next
End SubThe rest of this paper is organized as follows. The classic breakdown formula is presented in detail in Section2. 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σv has a magnitude equal to the overburden pressure in the vertical direction. The smallest horizontal principal stress σhmin is usually determined directly in the experiment from the shut-in pressure. The greatest horizontal principle stress σHmax 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 σHmax for hydraulic fracturing in nonporous impermeable rocks [23], ignoring the pore pressure term.(1)Pb=3σhmin-σHmax+T,where T 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]:(2)Pb=3σhmin-σHmax+T-Pp(3)Pb=3σhmin-σHmax+T-α1-2υ/1-υPp2-α1-2υ/1-υ,where Pb is the breakdown pressure, Pp 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]:(4)Pb=3σhmin-σHmax+T-βPp(5)Pb=3σhmin-σHmax+T-α1-2υ/1-υPp1+β-α1-2υ/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 (Section2). 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 samplesX1,Y1,X2,Y2,…,XN,YN, Xi∈Rn, Yi∈RQ, the MSVM model can be built by solving the following optimization problem based on an iterative reweighted least-square algorithm [25].(6)Lp′W,b=12∑j=1QWj2+12∑i=1Naiui2+τai=0uit<ε2Cuit-εuituit≥ε,where N is the number of input, Q is the number of output, W is the weight, b and C are constants, τ is the sum of constant terms that do not depend on either W or b, ε is the tolerant error, and t denotes thetth iteration. A brief description and the MSVM algorithm can be found in the literature [22]. The MSVM model can be expressed as(7)YX=∑k=1nWkX,Xk+b.Based on the above MSVM model, the nonlinear relationship between the borehole pressure and geomechanical parameters can be described as(8)MSVMX:Rn⟶RQY=MSVMX,where X=(x1,x2,…,xn) is the n-dimensional vector of the identified parameter, for example, the in situ stress, Young’s modulus, or Poisson’s ratio. Y=(y1,y2,…,yQ) is the Q-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.(9)xi,j=xminj+rand0,1xmaxj-xminj,where x(i,j) is the candidate solution of the problem; i=1,2,…,SN/2 and SN/2 denotes the size of the population; j=1,2,…,D and D is the dimension number of each solution; rand(0,1) is a random number between [0,1]; xmini and xmaxi 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 Section3.3) in the employed bees phase. A candidate solution is produced according to the following equation:(10)vi,j=xi,j+φijxi,j-xk,j,where k is different from i and is a randomly chosen index from {1,2,…,SN/2}, j is also an index randomly chosen from {1,2,…,D}, and φij is a random number in the range [-1,1] that controls the generation of food sources around x(i,j) 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 as(11)pi=fitnessi∑n=1SNfitnessn,where fitnessi 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,(12)fitness=MSVMX-Y2Q,where MSVM(X) is the predicted pressure using the MSVM model,Y is the vector of the monitored pressure, and Q 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 Figure1 and the algorithm is described as follows.Figure 1
Flowchart of the back analysis process to obtain the rock parameters.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 Step2, 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.


## 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 samplesX1,Y1,X2,Y2,…,XN,YN, Xi∈Rn, Yi∈RQ, the MSVM model can be built by solving the following optimization problem based on an iterative reweighted least-square algorithm [25].(6)Lp′W,b=12∑j=1QWj2+12∑i=1Naiui2+τai=0uit<ε2Cuit-εuituit≥ε,where N is the number of input, Q is the number of output, W is the weight, b and C are constants, τ is the sum of constant terms that do not depend on either W or b, ε is the tolerant error, and t denotes thetth iteration. A brief description and the MSVM algorithm can be found in the literature [22]. The MSVM model can be expressed as(7)YX=∑k=1nWkX,Xk+b.Based on the above MSVM model, the nonlinear relationship between the borehole pressure and geomechanical parameters can be described as(8)MSVMX:Rn⟶RQY=MSVMX,where X=(x1,x2,…,xn) is the n-dimensional vector of the identified parameter, for example, the in situ stress, Young’s modulus, or Poisson’s ratio. Y=(y1,y2,…,yQ) is the Q-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.(9)xi,j=xminj+rand0,1xmaxj-xminj,where x(i,j) is the candidate solution of the problem; i=1,2,…,SN/2 and SN/2 denotes the size of the population; j=1,2,…,D and D is the dimension number of each solution; rand(0,1) is a random number between [0,1]; xmini and xmaxi 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 Section3.3) in the employed bees phase. A candidate solution is produced according to the following equation:(10)vi,j=xi,j+φijxi,j-xk,j,where k is different from i and is a randomly chosen index from {1,2,…,SN/2}, j is also an index randomly chosen from {1,2,…,D}, and φij is a random number in the range [-1,1] that controls the generation of food sources around x(i,j) 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 as(11)pi=fitnessi∑n=1SNfitnessn,where fitnessi 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,(12)fitness=MSVMX-Y2Q,where MSVM(X) is the predicted pressure using the MSVM model,Y is the vector of the monitored pressure, and Q 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 Figure1 and the algorithm is described as follows.Figure 1
Flowchart of the back analysis process to obtain the rock parameters.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 Step2, 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σHmax and σhmin, respectively, Young’s modulusE, 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 wi and constant bi 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.Figure 2
Parameters and model of MSVM in the Excel VBA platform.Figure 3
Comparison of pressure data estimated by MSVM and borehole measured data.
(a)
(b)The ABC code is also written in Excel and VBA. The parameters of the ABC algorithm and the calculation results are shown in Figure4. 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 σHmax, σhmin, E, 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. σhmin agrees well with the value estimated by ANN-GA, (1), and (2). The relative error is only 1.07%. σHmax 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 σHmax (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.Table 1
Comparison of in situ stresses and mechanical parameters between MSVM model and other preexisted models.
Model
σ
H (MPa)
σ
h (MPa)
Young’s modulusE (GPa)
Poisson’s ratioυ
MSVM
24.46
14.33
44.02
0.25
ANN-GA#
35.83
14.968
39.92
0.26
Eq. (1)∗
34.549
15.637
Eq. (2)∗
24.749
15.637
#The model proposed by Zhang and Yin [2]. ∗The model proposed by Hubbert and Willis [23].Figure 4
MSVM-ABC-based back analysis, its parameters, and results in the Excel VBA platform.Figure 5
Comparison of recognized in situ stress using different models (σH, σh are the maximum and minimum in situ stress).Moreover, there are four borehole pressures, namely, formation breakdown pressure (FBP)P1, fracture propagation pressure (FPP) P2, instantaneous shut-in pressure (ISIP) P3, and leak–off pressure (LOP) P4. 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.Figure 6
Comparison of borehole pressures obtained by different prediction methods and monitored data.Figure 7
Fitness variation with increased cycles in the ABC analysis.Figure 8
Fitness distribution in the searching space at different cyclic stages.
(a)
(b)
(c)
(d)


## 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.



---

*Source: 5314628-2017-10-18.xml*