Comparison of XFEM code results with Warpinski and Teufel’s  experiments
Hydraulic fracture growth through naturally fractured reservoirs presents theoretical, design, and application challenges since hydraulic and natural fracture interaction can significantly affect hydraulic fracturing propagation. Although hydraulic fracturing has been used for decades for the stimulation of oil and gas reservoirs, a thorough understanding of the interaction between induced hydraulic fractures and natural fractures is still lacking. This is a key challenge especially in unconventional reservoirs, because without natural fractures, it is not possible to recover hydrocarbons from these reservoirs. Meanwhile, natural fracture systems are important and should be considered for optimal stimulation. For naturally fractured formations under reservoir conditions, natural fractures are narrow apertures which are around 10-5 to 10-3 m wide and have high length/width ratios (>1000:1) .Typically natural fractures are partially or completely sealed but this does not mean that they can be ignored while designing well completion processes since they act as planes of weakness reactivated during hydraulic fracturing treatments that improves the efficiency of stimulation . The problem of hydraulic and natural fracture interaction has been widely investigated both experimentally [3, 4, 5, 6, 7, 8] and numerically [9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. Many field experiments also demonstrated that a propagating hydraulic fracture encountering natural fractures may lead to arrest of fracture propagation, fluid flow into natural fracture, creation of multiple fractures and fracture offsets [19, 20, 21, 22] which will result in a reduced fracture width. This reduction in hydraulic fracture width may cause proppant bridging and consequent premature blocking of proppant transport (so-called screenout) [23, 24] and finally treatment failure. Although various authors have provided fracture interaction criteria [4, 5, 25] determining the induced fracture growth path due to interaction with pre-existing fracture and getting a viewpoint about variable or variables which have a decisive impact on hydraulic fracturing propagation in naturally fractured reservoirs is still unclear and highly controversial. However, experimental studies have suggested that horizontal differential stress, angle of approach and treatment pressure are the parameters affecting hydraulic and natural fracture interaction [4, 5, 6] but a comprehensive analysis of how different parameters influence the fracture behavior has not been fully investigated to date. In this way, in order to assess the outcome of hydraulic fracture stimulation in naturally fractured reservoirs the following questions should be answered:
What is the direction of hydraulic fracture propagation?
How will the propagating hydraulic fracture interact with the natural fracture?
Will the advancing hydraulic fracture cross the natural fracture or will it turn into it?
For the purpose of this study, a 2D eXtended finite element method (XFEM) has been compared with a feed-forward with back-propagation artificial neural network approach to account for hydraulic and natural fracture interaction.
2. Interaction between hydraulic and natural fractures
The interaction between pre-existing natural fractures and the advancing hydraulic fracture is a key issue leading to complex fracture patterns. Large populations of natural fractures are sealed by precipitated cements (Figure 1) which are weakly bonded with mineralization that even if there is no porosity in the sealed fractures, they may still serve as weak paths for the growing hydraulic fractures .
In this way, experimental studies [4, 5, 6] suggested several possibilities that may occur during hydraulic and natural fractures interaction. Blanton  conducted some experiments on naturally fractured Devonian shale as well as blocks of hydrostone in which the angle of approach and horizontal differential stress were varied to analyze hydraulic and natural fracture interaction in various angles of approach and horizontal differential stresses. He concluded that any change in angle of approach and horizontal differential stress can affect hydraulic fracture propagation behavior when it encounters a natural fracture which will be referred to as opening, arresting and crossing. Warpinski and Teufel  investigated the effect of geologic discontinuities on hydraulic fracture propagation by conducting mineback experiments and laboratory studies on Coconino sandstone having pre-existing joints. They observed three modes of induced fracture propagation which were crossing, arrest by opening the joint and arrest by shear slippage of the joint with no dilation and fluid flow along the joint. In 2008  some laboratory experiments were performed to investigate the interaction between hydraulic and natural fractures. They also observed three types of interactions between hydraulic and pre-existing fractures which were the same as Warpinski and Teufel’s observations. The above referenced experimental studies have investigated the initial interaction between the induced fracture and the natural fracture, however, in reality may be the hydraulic fracture is arrested by natural fracture temporarily but with continued pumping of the fluid, the hydraulic fracture may cross (Figure 2) or turn into the natural fracture (Figure 3).
Alternatively, in some cases the hydraulic fracture may get arrested if the natural fracture is long enough and favorably oriented to accept and divert the fluid.
3. Numerical modeling: Extended Finite Element Method (XFEM)
For fracture propagation through numerical modeling an energy based criterion has been considered which is energy release rate, G. The energy release rate, G, is related to the stress intensity factors through Eq. 1 :
where E’ = E for plane stress (E is Young's modulus) and E’ = E/(1- ν2) for plane strain (where ν is the Poisson's ratio). Energy release rate has been calculated by the J integral using the domain integral approach  whereas J integral is equivalent to the definition of the fracture energy release rate, G, for linear elastic medium. If the G is greater than a critical value, Gc, the fracture will propagate critically. The direction of hydraulic fracture propagation will be calculated by Eq. 2 :
During hydraulic and natural fracture interaction at the intersection point the hydraulic fracture has more than one path to follow which are crossing and turning into natural fracture. The most likely path is the one that has the maximum G. So, at the intersection point energy release rate is calculated for both crossing (Gcross) and turning into natural fracture (Gturn), and if (Gturn /Gcross)>1 hydraulic fracture turns into natural fracture while if (Gturn /Gcross)< 1 crossing takes place and hydraulic fracture crosses the natural fracture. To examine the proposed mechanism, eXtended Finite Element method (XFEM) was applied which was first introduced by Belytschko and Black  in order to avoid explicit modeling of discrete cracks by enhancing the basic finite element solution. In comparison to the classical finite element method, XFEM provides significant benefits in the numerical modeling of fracture propagation and it overcomes the difficulties of the conventional finite element method for fracture analysis, such as restriction in remeshing after fracture growth and being able to consider arbitrary varying geometry of fractures . XFEM enhances the basic finite element solution through the use of enrichment functions which are the Heaviside function for elements that are completely cut by the crack and Westergaard-type asymptotic functions for elements containing crack-tips . The displacement field for a point “x” inside the domain can be approximated based on the XFEM formulation as below :
Where NI is the finite element shape function, uI is the nodal displacement vector associated with the continuous part of the finite element solution, H(x)is the Heaviside enrichment function where it takes the value +1 above the crack and –1 below the crack, aI is the nodal enriched degree of freedom vector associated with the Heaviside (discontinuous) function, Fα(x) is the near-tip enrichment function, is the nodal enriched degree of freedom vector associated with the asymptotic crack-tip function, Nu is the set of all nodes in the domain, Nα is the subset of nodes enriched with the Heaviside function and Nb is the subset of nodes enriched with the near tip functions. At the intersection point, instead of Heaviside enrichment function, Junction function will be applied as shown in Figure 4 . By all means, XFEM is well-suited for modeling hydraulic fracture propagation and diversion in the presence of natural fracture.
4. Artificial intelligence: Artificial Neural Network (ANN)
Artificial Neural Network (ANN) is considered as a different paradigm for computing and is being successfully applied across an extraordinary range of problem domains, in all areas of engineering. ANN is a non-linear mapping structure based on the function of the human brain that can solve complicated problems related to non-linear relations in various applications which makes it superior to conventional regression techniques [33, 34]. ANNs are capable of distinguishing complex patterns quickly with high accuracy without any assumptions about the nature and distribution of the data and they are not biased in their analysis. The most important aspect of ANNs is their capacity to realize the patterns in obscure and unknown data that are not perceptible to standard statistical methods. Statistical methods use ordinary models that need to add some terms to become flexible enough to satisfy experimental data, but ANNs are self-adaptable. The structure of the neural network is defined by the interconnection architecture between the neurons which are grouped into layers. A typical ANN mainly consists of an input layer, an output layer, and hidden layer(s) (Figure 5). As shown in Figure 5, each neuron of a layer is connected to each neuron of the next layer. Signals are passed between neurons over the connecting links. Each connecting link has an associated weight which multiplies by the related input.
Also, to diversify the various processing elements, a bias is added to the sum of weighted inputs called net input shown in Eq. (4)
Where n is the net input, w is the weight, p is the input, b is the bias, S is the number of neurons in the current layer and R is the number of neurons in the previous layer.
Each neuron applies an activation function to its input to determine its output signal . Neurons may use any differentiable activation function to generate their output based on problem requirement. The most useful activation functions are as follows:
where a is the neuron layer output. Purelin is a linear activation function (Figure 6A) defined in Eq. (5). Tansig is hyperbolic tangent sigmoid activation function (Figure 6B) mathematically shown in Eq. (6). Radbas is Gaussian activation function (Figure 6C) shown in Eq. (7). In Figure. 7, a one layer network with R inputs and S neurons is shown . The optimum number of hidden layers and the number of neurons in these layers are determined by trial and error during the training/learning process. The hidden layers in the network are used to develop the relationship between the variables. In general, multilayer networks are more powerful than single-layer networks .
4.1. Feed-forward network with back-propagation
The feed-forward network with back-propagation (FFBP) is one of the most eminent and widespread ANNs in engineering applications . In addition, it is easy to implement and solves many types of problems correctly . Usually, FFBP uses tansig and purelin as activation functions in the hidden and output layers, respectively and the net input is calculated the same as Eq. 4. FFBP operates in two steps. First, the phase in which the input information at the input nodes is propagated forward to compute the output information signal at the output layer. In other words, in this step the input data are presented to the input layer and the activation functions process the information through the layers until the network’s response is generated at the output layer. Second, the phase in which adjustments to the connection strengths are made based on the differences between the computed and observed information signals at the output. In this step, the network’s response is compared to the desired output and if it does not agree, an error is generated. The error signals are then transmitted back from the output layer to each node in the hidden layer(s) . Then, based on the error signals received, connection weights between layer neurons and biases are updated. In this way, the network learns to reproduce outputs by learning patterns contained within the data. One iteration of this algorithm can be written as Eq. 8:
where Xk is a vector of current weights and biases, gk is the current gradient, and is the learning rate. Once the network is trained, it can then make predictions from a new set of inputs that was not used to train the network.
4.2. ANN performance criteria
There are several quantitative measures to assess ANN performance that the most usual one in a binary classification test is accuracy  (Fawcett 2006). To understand the meaning of accuracy, some definitions like true positive, false positive, true negative and false negative should be explained. Imagine a scenario where the occurrence of an event is considered. The test outcome can be positive (occurrence of the event) or negative (the event doesn’t occur). According to this scenario:
True Positive (TP): The event occurs and it is correctly diagnosed as it occurs;
False Positive (FP): The event doesn’t occur but it is incorrectly diagnosed as it occurs;
True Negative (TN): The event doesn’t occur and it is correctly diagnosed as it doesn’t occur;
False Negative (FN): The event occurs but it is incorrectly diagnosed as it doesn’t occur.
According to above definitions:
In general, the accuracy of a system is a degree of closeness of the measured values to the actual (true) values .
5. Results and discussions
Physically, modeling hydraulic fracturing is a complicated phenomenon due to the heterogeneity of the earth structure, in-situ stresses, rock behavior and the physical complexities of the problem, hence if natural fractures are added up to the problem it gets much more complex in both field operation and numerical aspects. For simplicity, it is assumed that rock is a homogeneous isotropic material and the fractures are propagating in an elastic medium under plane strain and quasi-static conditions driven by a constant and uniform net pressure throughout the hydraulic fracture system. Fracturing fluid pressure is included in the model by putting force tractions on the necessary degrees of freedom along the fracture. A schematic illustration for the problem has been presented in Figure 8 which shows that hydraulic fracture propagates toward the natural fracture and intersects with it at a specific angle of approach, θ, and in-situ horizontal differential stress, (σ1 -σ3).
So, a 2D XFEM code has been developed to model hydraulic fracture propagation in naturally fractured reservoirs and interaction with natural fractures. For this purpose, firstly Warpinski and Teufel’s  experiments have been modeled to see how much the results of the developed XFEM model for hydraulic and natural fracture interaction, are compatible with them. Table 1, presents the results of XFEM code which can be compared with Warpinski and Teufel’s  experiments. As shown in Table 1, the results of XFEM code indicate that at high to medium angles of approach, crossing and turning into natural fracture both are observed depending on the differential stress while at low angles of approach with low to high differential stress, the predominant case during hydraulic and natural fracture interaction is hydraulic fracture diversion along natural fracture which are in good agreement with Warpinski and Teufel’s  experiments.
|Angle of approach(θo)||Max. horizontal stress(psi)||min. horizontal stress(psi)||Horizontal differential stress(psi)||Experimental results ||Gturn/Gcross||XFEM results|
|30||1000||500||500||Turn into||3.46||Turn into|
|30||1500||500||1000||Turn into||2.05||Turn into|
|30||2000||500||1500||Turn into||1.29||Turn into|
|60||1000||500||500||Turn into||1.948||Turn into|
|60||1500||500||1000||Turn into||1.201||Turn into|
|90||1000||500||500||Turn into||1.013||Turn into|
Meanwhile debonding of natural fracture prior to hydraulic and natural fracture intersection could also be modeled which is a complicated and very interesting phenomena that has been rarely investigated. Figure 9, presents pre-existing fracture debonding before intersection with hydraulic fracture at approaching angles of 30o, 60o, 90o in Warpinski and Teufel’s  experiments. As it is clearly observed in stress maps in Figure 9, a tensile stress is exerted ahead of hydraulic fracture tip for all of the approaching angles which makes the natural fracture debonded. In addition, the length and the position of the debonded zone vary depending on natural fracture orientation and horizontal differential stress.
Figure 10, shows the debonded zone at the intersecting point of hydraulic and natural fracture and Figure 11 presents the result of hydraulic and natural fracture interaction for approaching angles of 30o, 60o, 90o.
In the second step, a FFBP neural network has been applied for predicting growing hydraulic fracturing path due to interaction with natural fracture in such a way that horizontal differential stress, angle of approach, interfacial coefficient of friction, young’s modulus of the rock and flow rate of fracturing fluid are the inputs and hydraulic fracturing path (crossing or turning into natural fracture) is the output whereas tansig is an activation function. The data set used in this model consists of around 100 data based on experimental studies [4, 5, 6]. Table 2, represents the range of the parameters used in the developed ANN.
|Horizontal differential stress (psi)||290||2175|
|Angle of approach (deg)||30||90|
|Interfacial coefficient of friction||0.38||1.21|
|Young’s modulus of the rock (psi)||1.218*106||1.45*106|
|Flow rate of fracturing fluid (m3/s)||4.2*10-9||8.2*10-7|
The data in the database were randomly divided into two subsets. The training subset used 70% and the remaining 30%, was used for the testing subset. For standardizing the range of the input data and improves the training process the data used in network development were pre-processed by normalizing. Normalizing the data enhances the fairness of training by preventing an input with large values from swamping out another input that is equally important but with smaller values . The optimal number of the neurons of a single hidden layer network for the developed ANN using trial and error method based on accuracy is 19 which is shown in Table 3. The developed FFBP neural network represented a high accuracy of 96.66% which was so promising. Also, according to the dataset, around 30 data were assigned for testing subset. Figure 12, show the results of the developed FFBP neural network predictions with actual measurements for testing subset.
|Number of hidden neurons||Accuracy (%)|
As shown in Figure 12, FFBP predictions are in prominent agreement with actual measurements which shows the high efficiency of the developed FFBP neural network approach for predicting hydraulic fracturing path due to interaction with natural fracture based on horizontal differential stress, angle of approach, interfacial coefficient of friction, young’s modulus of the rock and flow rate of fracturing fluid. Finally, both XFEM and ANN approaches have been examined by a set of experimental study data [4, 6] and the results have been compared. The results of a comparison are presented in Table 4. As shown in Table 4 both of the proposed approaches yield quite promising results and in just one case ANN approach result doesn’t agree with the actual case.
|Horizontal Differential Stress(psi)||Angle of approach (deg)||Coefficient of friction||E*106 (psi)||Flow rate of fracturing fluid(m3/s)||Actual Case||XFEM Result||FFBP Prediction|
Two new numerical modeling and artificial intelligence methodologies were introduced and compared to account for hydraulic and natural fracture interaction. First a new approach has been proposed through XFEM model and an energy criterion has been applied to predict hydraulic fracture path due to interaction with natural fracture. To validate and show the efficiency of the developed XFEM code, firstly the results obtained from XFEM model have been compared with experimental studies which shows good agreement. It’s been concluded that natural fracture most probably will divert hydraulic fracture at low angles of approach while at high horizontal differential stress and angles of approach of 60 or greater, the hydraulic fracture crosses the natural fracture. Meanwhile, the growing hydraulic fracture exerts large tensile stress ahead of its tip which leads to debonding of sealed natural fracture before intersecting with hydraulic fracture that is a key point to demonstrate hydraulic and natural fracture behaviors before and after intersection. Then, a FFBP neural network was developed based on horizontal differential stress, angle of approach, interfacial coefficient of friction, young’s modulus of the rock and flow rate of fracturing fluid and the ability and efficiency of the developed ANN approach to predict hydraulic fracturing path due to interaction with natural fracture was represented. The results indicate that the developed ANN is not only feasible but also yields quite accurate outcome. Finally, both of the approaches have been compared and both of them yield promising results. Numerical modeling yields more detailed results which can be used for further investigations and it can explain different observed behaviors of hydraulic fracturing in naturally fractured reservoirs as well as activation of natural fractures. Also, the potential conditions that may lead to hydraulic fracturing operation failure can be investigated through numerical modeling but it is computationally more expensive and time-consuming than artificial neural network approach. In another hand, since artificial neural network approach is mainly data-driven it can be of great use in real-time experimental studies and field hydraulic fracturing in naturally fractured reservoirs. So, as one may conclude easily, numerical modeling and artificial intelligence both have some positive and negative points; hence simultaneous use of these methods will lead to both technical and economical advantages in hydraulic fracturing operation especially in the presence of natural fractures.