Introduction

According to BP Statistical Review (2017), the global demand for natural gas increased by 1.5% from 2016. Therefore, the production of unconventional gas reservoir must be optimized to cover the decreasing of conventional gas reservoir production (Bentley 2002; Hughes 2013). The unconventional gas reservoir is commonly defined as reservoir with poor permeability (less than 0.1 mD) (Gordon 2012; Holditch 2006). There are four categories of unconventional gas that are becoming important as future source of energy, they are shale gas, tight gas sand, coal bed methane and hydrates (Abdelfattah et al. 2015). Production of unconventional resources is expected to thrive and grow sixfold from 2011 to 2030 (Ruehl and Giljum 2011).

The shale revolution is referred as the bloom of shale production due to the technological improvement of horizontal drilling and hydraulic fracturing. The improvements include hydraulic fracturing in horizontal well reduced from 200 m apart to 50 m while having over 60 stages per well (Aguilera and Ripple 2012; Ignatyev et al. 2011) the USA is the world’s largest dry natural gas producer because of the shale gas revolution. It contributes to 20% of the world consumption, in which 40% of the supply is from shale (Sieminski 2014). The revolution mainly impacts the USA due to the advancement in technologies in current stage; however, it is expected to spread globally to countries with abundant shale reserves, such as Canada, China and Argentina (Aguilera and Radetzki 2014).

Hydraulic fracturing is a well stimulation method used to improve the permeability of tight gas reservoir for economical production (Britt 2012; King 2012). As the fracture grows, it follows the direction set by in situ stresses in the rock, which can be termed as “fracture direction”. The propagation of fracture is perpendicular to plane of least effective stress (Britt 2012; King 2012). Microseismic monitoring is one of the methods to study the characteristics of hydraulic fracturing (Shapiro et al. 2006). Microseismicity is the occurring microearthquakes events caused by the injection of fluid into the borehole (Shapiro et al. 2006). All hydraulic fracturing operations stimulated different magnitude of microseismicity. The cloud of microseismic events can be translated into Stimulated Reservoir Volume (SRV), which serves as crucial parameter to study the characterization of reservoir and hydraulic fracturing (Shapiro et al. 2006). Microseismic monitoring is a passive measurement of microseismic event and provides crucial information such as magnitude, location and time of the event (Maulianda 2016). It is proven to be useful in showing the fracture geometries such as fracture length, fracture height, fracture width and fracture azimuth with high confidence (Warpinski and Wolhart 2016). Other application of microseismic monitoring in hydraulic fracturing includes observation of activated natural faults and faults, permeability (Shapiro et al. 1997), and fracturing characterization using Moment Tensor Inversion (MTI) (Nolen-Hoeksema and Ruff 2001).

Stimulated reservoir volume is defined as a collection of fluid-induced microseismic events that reflects the volume of fracture network created in the reservoir (Mayerhofer et al. 2010). Stimulated Reservoir Volume (SRV) was introduced to act as a correlation parameter for well performance (Mayerhofer et al. 2010). SRV compromised by a complex fracture network with different geometries ranging from curved, planar, slanted and of different lengths (Vera and Shadravan 2015). In horizontal well, SRV and fracture network size can be improved with longer lateral lengths and increasing stages of simulation. The study of SRV characteristics optimizes the hydraulic fracturing design. Analytical modelling and numerical modelling are used for SRV modelling.

Other uses of hydraulic fracturing

Enhanced geothermal system (EGS)

Geothermal energy is the energy stored within the Earth’s crust and it is one of the promising clean renewable energy resources (MIT 2006; Vik et al. 2018). In 2016, it was estimated that only 6-7% of its total global potential had been extracted (GEA 2016; Vik et al. 2018). The common challenge to exploit geothermal reservoir is the low permeability of the reservoir. To enhance permeability, hydraulic fracturing is performed.

Enhanced Geothermal System (EGS) or Hot Dry Rock (HDR) geothermal system as formerly known was pioneered at Los Alamos National Laboratory, New Mexico, USA in the 1970s (Barbier 2002; Xia et al. 2017). In EGS, a pair of wells are drilled to a depth where rock temperatures approach 300C (White et al. 2017). The wells are then hydraulically fractured, thus creating a connection between the wells. Fluid such a brines, compressed CO2 or liquid mixtures are then circulated from one well to the other in a closed system through the created hydraulic connection. At ambient surface pressure, the heated fluid converts to steam which subsequently drives a turbine to produce electricity. In another method, the heated fluid exchanges heat with a working fluid, which will turn into vapour during the heat exchange process and it will drive the turbine. The cooled fluids are then re-injected into the thermal reservoir or cooled further in a secondary recovery system (White et al. 2017).

To improve heat extraction from geothermal reservoir with low permeability, multiple-induced fractures are created to establish multiple flow paths between the injection well and production well, and to create multiple contact surfaces for heat exchange between hot rock and cold fluid (Bataille et al. 2006; Vik et al. 2018).

Block cave mining

Hydraulic fracturing started to be applied in mining industry only in recent decades. The application initially was for methane extraction in coal mining in the 1970s in the USA and to control hard roof rockburst (He et al. 2015). Following that, it was applied in cave mining.

Block Cave Mining refers to a method of underground mass mining method in which ore extraction relies on gravity action (Adams and Rowe 2013). The method to induce caving involves undercutting of block by the means of blasting it in order to destroy its ability to support the overlying rock. Gravity then acts to fracture the block (Eklind et al. 2007). Pre-conditioning is required in the event there are massive, unfractured ore body to initiate caving and to reduce caving material size. One of the favoured pre-conditioning methods is intensive hydraulic fracturing in boreholes drilled into the ore body (Van and Jeffrey 2000). The hydraulic fracturing pressure can reach up to 10,000 psi, and pumped pure water volume in the range of 4000–5000 litres per fracture, however, the figures can be larger based on pump size and pressure response (Adams and Rowe 2013).

Block cave mining objective is to create horizontal radial hydraulic fractures (HFHFs) that are able to propagate across existing vertical or inclined natural fractures, or to create inclines hydraulic fractures in horizontal or sub-horizontal natural fractures dominated region (He et al. 2015).

Geologic carbon sequestration (GCS)

In geologic carbon sequestration (GCS), large volume of CO2 are injected into deep geological formations to prevent it from releasing to the atmosphere (Fu et al. 2017; Haszeldine 2009; Orr 2009). The targeted storage is typically saline aquifers or depleted hydrocarbon reservoirs which are overlain by caprocks with low permeability. Consideration to be taken in the CO2 storage design is to ensure caprock integrity to prevent CO2 leakage. However, it is known that hydraulic fractures may be initiated and propagated in the caprock when injected fluid pressure exceeds the minimum in situ principal stress of the caprock. Fluid flow along open hydraulic fracture is significantly more efficient as compared to fluid flow in porous medium to deliver the fluid to far-field reservoir; hence, the injection through hydraulic fracturing could improve both storage capacity and CO2 injection. The most desirable scenario in GCS is when the hydraulic fracture can be contained within reservoir rock without fracturing the caprock, at injection pressure level lower than the minimum principal stress of the caprock (Fu et al. 2017).

Analytical modelling

The SRV modelling is of the highest interest in the hydraulic fracturing designs to predict the facture growth in terms of length, height and width (Rahman and Rahman 2010). The study of SRV behaviour will enable better optimization of the hydraulic fracturing procedures. Analytical modelling of SRV applies simplified mathematical equations to model the fracture propagation within the rocks (Table 1 and Fig. 1). There have been many 2D and 3D models developed to study the fracture propagation under different conditions and assumptions (Table 2).

Fig. 1
figure 1

Comparison of analytical models (Nordgren 1972; Peirce et al. 2010; Perkins and Kern 1961; Settari and Cleary 1984; Simonson et al. 1978; Warpinski and Smith 1989; Yew and Weng 2014; Yousefzadeh et al. 2015; Yu and Aguilera 2012; Zheltov 1955)

Full size image
Table 1 Summary of the available analytical models (Rahman and Rahman 2010; Xiang 2011)
Full size table
Table 2 Summary of researches conducted on studying analytical models (Grechka et al. 2010; Shapiro et al. 1997, 2006; Yousefzadeh et al. 2016; Yu and Aguilera 2012)
Full size table

Two-dimensional models

The 2D models were developed in the early 1960s as a simple approach for fracture design of industry requirement (Gidley 1990). The models are Kristinaovic-Geertsma-de Klerk (KGD) model (Geertsma and De Klerk 1969; Zheltov 1955). Perkins-Kem-Nordgren (PKN) model (Nordgren 1972; Perkins and Kern 1961) and Radial model (Abe et al. 1976). The 2D models have been reasonably successful in practical simulation; however, the assumption of fracture shape and fracture height need to be specified to perform the models limits the practicality of 2D models.

KGD model assumes plane strain in horizontal direction. It represents a fracture with horizontal penetration that is lower than the vertical penetration. The model has cusp-shaped fracture tip and the fracture width is uniform in vertical direction. The model is most suitable to be applied for fractures with proportional length-to-height ratio.

In the PKN model, fracture planes are perpendicular to the vertical plane strain. The model has an elliptical-shaped cross section and the fracture toughness is assumed to be not affecting the fracture geometry (Yew and Weng 2014). The pressure is uniform as there is the absence of fluid flow in vertical section. The model is most suitable to be applied for fracture height that is higher than fracture length.

Both 2D constant height models consider fracture width and length as a function of fracture height, treatment parameters and reservoir properties. Vertical fracture propagation is limited by change in material property and minimum horizontal in situ stress (Yousefzadeh et al. 2015). Fracture deformation is linear plastic process in both models. FracCADE is a fracturing design commercial software from Schlumberger for PKN and KGD models.

Three-dimensional models

The pseudo-3D (P3D) model is proposed to idealize the fracture propagation in multi-layered formations. The constant fracture height assumption is removed in this model (Settari and Cleary 1984; Simonson et al. 1978; Warpinski and Smith 1989). The model is “pseudo” because variation of fracture geometry in three-dimensional space is not considered. The model accounts for height variation by considering the in situ stresses contrast, rock toughness and local net fluid pressure. Some of the commercial softwares available in the market are FLAC3D by Itasca, FracPro by CARBO and FracMan by Golder (Hou and Zhou 2011).

The most recent development by Yu and Aguilera (2012) is to use mass balance derived 3D model to account for fracture geometry variation. The diffusivity model accounts for the spatio-temporal distribution of fluid-induced seismicity. In their study, the front of microseismic events represents the front of pore pressure diffusion. Hydraulic diffusivity coefficient can be determined from the slope of microseismic events versus time plot. The model can determine fracture geometry under different stimulations case after the determination of hydraulic diffusivity coefficient from three-dimensional microseismic events. Diffusivity model can be solved using MATLAB Partial Differential Equation Toolbox (Peirce et al. 2010).

Diffusivity model

The 3D model is developed based on linear diffusion equation, for which the approximate solution is shown by Eq. 1.

$$\frac{{\Delta P_{\text{res}} }}{{\Delta P_{\text{inj}} }} = \frac{{P\left( {x,y,t} \right) - P_{i} }}{{P_{\text{inj}} - P_{i} }} = {\text{erfc}}\left( {\frac{x}{{\sqrt {4\eta_{x} t} }}} \right){\text{erfc}}\left( {\frac{y}{{\sqrt {4\eta_{y} t} }}} \right){\text{erfc}}\left( {\frac{z}{{\sqrt {4\eta_{z} t} }}} \right)$$
(1)

Another variation of 3D model based on nonlinear diffusion equation can be summarized as below, where it is derived from continuity equation and Darcy law.

Equation 2 shows the equation derived from continuity equation and Darcy law with the incorporation of real gas equation and formation compressibility equation.

$$\frac{\partial }{\partial x}\left( {\frac{kp}{2\mu z}\frac{\partial }{\partial x}\left( {P^{2} } \right)} \right) = \frac{1}{0.000264}\varphi \mu c_{t} \frac{1}{2\mu z}\frac{\partial }{\partial t}\left( {P^{2} } \right)$$
(2)
$$u = \frac{{p^{2} }}{\mu z}$$
(3)
$$P = \sqrt {u\mu z}$$
(4)
$$\frac{{{\text{d}}p}}{{{\text{d}}x}} = \sqrt {\mu z} \frac{1}{2}u^{{ - \frac{1}{2}}} \frac{{{\text{d}}u}}{{{\text{d}}x}}$$
(5)
$$\frac{{{\text{d}}p}}{{{\text{d}}y}} = \sqrt {\mu z} \frac{1}{2}u^{{ - \frac{1}{2}}} \frac{{{\text{d}}u}}{{{\text{d}}y}}$$
(6)
$$\frac{\partial }{\partial x}\left( {k_{x} \frac{\partial u}{\partial x}} \right) + \frac{\partial }{\partial y}\left( {k_{y} \frac{\partial u}{\partial y}} \right) = \frac{{\varphi \mu c_{t} }}{0.000264}\frac{{{\text{d}}u}}{{{\text{d}}t}}$$
(7)

The porosity as a function of pressure is defined below (Abe et al. 1976).

$$\varphi = \varphi_{i} \left( {1 - c_{t} \left( {P_{0} - P} \right)} \right)$$
(8)

The permeability as a function of porosity and specific surface is based on the Kozeny–Carman equation is defined below (Gates 2011). Rearranging the equation leads to permeability as a function of pressure.

$$k = \frac{{\varphi_{i} \left( {1 - c_{t} \left( {P_{0} - P} \right)} \right)^{3} }}{{5s^{2} \left( {1 - \left( {\varphi_{i} \left( {1 - c_{t} \left( {P_{0} - P} \right)} \right)} \right)} \right)^{2} }}10^{15}$$
(9)

Numerical models

Modelling and analysis of hydraulic fractures propagation and interaction have gained enormous interest in petroleum engineering area. Hydraulic fracturing is one of the key important drivers in the development of unconventional reservoirs. Hydraulic fracturing demand has increased rapidly. Modelling tools have been developed to estimate the hydraulic fracturing performance in the formation.

Commonly used numerical models

Numerical simulation is the incorporation of the analytical solution to identify the behaviour of the events in the presence of more constituents in the system. There are several types of the models available to assess the hydraulic fracture propagation which are divided into two major groups, continuum method and discontinuum method.

Continuum method implies that the matter on which fracture to be simulated is continuous. The pre-existing fractures and stress distribution can be introduced within the continuity and no necessity to introduce the interaction of the boundaries. This method provided the stress distribution at arbitrary location of the fractures and fracture propagation is deduced based on the stress values obtained. However, this method may not show exact behaviour of the fractures. The summary of strength and limitation of the continuum method can be seen on Table 3.

Table 3 Continuum methods (Arndt et al. 2015; Costabel 1987; Hsiao 2006; Lee et al. 2015; Lei et al. 2017; Li et al. 2015; Lucia and Fogg 1990; McClure and Horne 2013; Postnikov et al. 2017)
Full size table

Discontinuum method on the other hand does show the fracture as separate entity which is introduced as boundary or stress computation limit. These methods give the advantage in simulating the cases with high frequency of cohesive elements within the formation. The summary of strength and limitation of the discontinuum method can be seen on Table 4.

Table 4 Discontinuum methods (Costabel 1987; Lei et al. 2017; Li et al. 2015; Marina et al. 2014; McCluer and Home 2013; Sukumar and Prévost 2003)
Full size table

The summary of the methods based on the research conducted by several authors (Arndt et al. 2015; Costable 1987; Hsiao 2006; Lee et al. 2015; Lei et al. 2017; Li et al. 2015; Lucia and Fogg 1990; Mcclure and Horne 2013; Postnikov et al. 2017) and (Costable 1987; Lei et al. 2017; Li et al. 2015; Marina et al. 2014; Mcclure and Horne 2013; Sukumar and Prévost 2003) is provided in the diagram below.

Extended finite element modelling (XFEM)

XFEM is based on the introduction of the enrichment functions on the previously successfully implemented FEM where it is used to simulate the interaction of solid and liquid during injection of fracturing fluid (Maulianda 2016). The enrichment functions allow the fracture propagation simulation to be computed without the necessity to introduce new meshes. The general form of XFEM is as given in Eq. 10 below (Youn 2016).

$$u\left( x \right) = \mathop \sum \limits_{i = 1}^{N} N_{i} \left( x \right)\bar{u}_{i} + \mathop \sum \limits_{j = 1}^{{M_{1} }} N_{j} \left( x \right)f_{1}^{\text{enr}} \bar{a}_{j} + \mathop \sum \limits_{k = 1}^{{M_{2} }} N_{k} \left( x \right)f_{2}^{\text{enr}} \bar{b}_{k} + \mathop \sum \limits_{l = 1}^{{M_{3} }} N_{l} \left( x \right)f_{3}^{\text{enr}} \bar{c}_{l}$$
(10)

where \(M_{1}\), \(M_{2}\) and \(M_{3}\) are different nodes of enrichment, \(f_{1}^{\text{enr}}\), \(f_{2}^{\text{enr}}\) and \(f_{3}^{\text{enr}}\) are enrichment functions, \(N_{i}\), \(N_{j}\), \(N_{k}\) and \(N_{l}\) are finite element shape functions and \(\bar{a}_{j}\), \(\bar{b}_{k}\) and \(\bar{c}_{k}\) are the degrees of freedom.

The enrichment functions introduced in XFEM are Heaviside and Branch enrichment functions. Heaviside function is derived based on the solutions of the level set function (LSF). LSF defines the fracture location and the fracture tip location. Heaviside function is used to determine the fracture behaviour separating two elements or elements being separated by void (Fig. 2) (Youn 2016).

Fig. 2
figure 2

Heaviside enrichment function (Youn 2016)

Full size image

For the cases with the interaction of two fractures or interaction of fracture with void junction function is implemented. Junction function also uses the LSF to define the fracture location while including the factor which accounts for the secondary fracture.

The Branch enrichment function is applied to identify the behaviour of the fracture tip within the element (Fig. 3).

Fig. 3
figure 3

Junction function and branch function (Youn 2016)

Full size image

By incorporating the Heaviside, Branch and Junction enrichment function on the GFEM function Eq. 11 is produced.

$$\begin{aligned} u\left( x \right) = \overbrace {{\mathop \sum \limits_{i = 1}^{N} N_{i} \left( x \right)\bar{u}_{i} }}^{{{\text{FEM}}\,{\text{displacement}}\,{\text{field}}}} + \underbrace {{\mathop \sum \limits_{j = 1}^{{M_{fr} }} N_{j} \left( x \right)\left\{ {H\left( x \right) - H\left( {x_{j} } \right)} \right\}\bar{a}_{j} }}_{{{\text{Heaviside}}\,{\text{enrichment}}}} + \hfill \\ \overbrace {{\mathop \sum \limits_{k = 1}^{{M_{tip} }} N_{k} \left( x \right)\mathop \sum \limits_{\gamma = 1}^{4} \left\{ {\varPhi_{\gamma } \left( x \right) - \varPhi_{\gamma } \left( {x_{k} } \right)} \right\}\bar{b}_{k} }}^{{{\text{Branch}}\,{\text{enrichment}}}} + \overbrace {{\mathop \sum \limits_{l = 1}^{{M_{jct} }} N_{K} \left( x \right)\left\{ {J\left( x \right) - J\left( {x_{l} } \right)} \right\}\bar{c}_{l} }}^{{{\text{Junction}}\,{\text{enrichment}}}} \hfill \\ \end{aligned}$$
(11)

where \(H\left( x \right) - H\left( {x_{j} } \right)\), Heaviside enrichment, \(\sum\nolimits_{y = 1}^{4} {\{ \varPhi_{\gamma } \left( x \right) - \varPhi_{\gamma } \left( {x_{k} } \right)\} }\), branch enrichment and \(J\left( x \right) - J\left( {x_{l} } \right)\), junction enrichment.

As the enrichment functions are user defined the XFEM gives a prospect to simulation of the fracture propagation behaviour including any factor required. The necessity is to define the enrichment function in terms of the FEM model composites. As the enrichment is applied on the nodes surrounding the fracture the definition of the equation should be computed based on the node characteristics (Bordas and Moran 2006; Feng and Gray 2017; Sepehri et al. 2015; Youn and Griffiths 2014). Enrichment functions are added into FEM model within the ABAQUS by incorporating the external simulator or built-in python compiler. The latest inclusion into the XFEM is the stochastic simulation of the fracture behaviour. Addition of the random field theory to the currently available XFEM which is also referred as XRFEM allows the simulation of the cases with the inclusion of the uncertainty.

The XRFEM includes Monte Carlo simulation and the average fracture behaviour is computed. This assists in accounting for the uncertainties present within the formation. Modelling of the XFEM in ABAQUS software is similar to FEM. This is due to the XFEM using the FEM as base function. Therefore, the model generation is same as FEM until interaction assignment.

Level set function LSF function defines the output to be used in XFEM node enrichment based on the fracture geometry and spatial location. Two-level set function, sign distance function and tangential level set functions are used to identify the distance of the fracture and fracture tip location (Youn 2016).

The distance of the node from the fracture in Eq. 12 is defined as a multiple of the normal from the fracture which is visualized in Fig. 4. The location of the fracture tip is defined as a function of fracture location and tangent unit vector of fracture as shown in Eq. 12 (Youn 2016).

$${\varphi }\left( {\mathbf{x}} \right) = {\text{sign}}\left\{ {{\mathbf{n}} \cdot \left( {{\mathbf{x}} - {\bar{\mathbf{x}}}} \right)} \right\}{ \hbox{min} }{\mathbf{x}} - {\bar{\mathbf{x}}}$$
(12)

where x is the nodal point, \(\bar{x}\) is the normal projection of the point x onto the fracture and n is the unit normal from the fracture at \(\bar{x}\).

$$\varPsi_{\text{i}} \left( {\mathbf{x}} \right) = \left( {{\mathbf{x}} - {\mathbf{x}}_{{\mathbf{i}}} } \right) \cdot \hat{\varvec{t}}_{\varvec{i}}$$
(13)

where \(\hat{t}_{i}\) is the unit vector tangent to the fracture at its tip and \(x_{i}\) is the location of the ith fracture tip.

Fig. 4
figure 4

LSF function parameter visualization (Youn 2016)

Full size image

Due to the dynamicity of the fracture propagation behaviour it is also necessary to define the algorithm to locate where LSF update is required.

As the direction of the fracture propagation is based on the stress intensity factor (SIF) as shown in Eq. 14 it is possible to define the LSF update location based on SIF.

$$\theta = 2\arctan \frac{1}{4}\left( {\frac{{K_{\text{I}} }}{{K_{\text{II}} }} \pm \sqrt {\left( {\frac{{K_{\text{I}} }}{{K_{\text{II}} }}} \right)^{2} + 8} } \right), - \pi < \theta < \pi$$
(14)

where \(\theta\) is the angle of the fracture propagation relative to the existing fracture direction while \(K_{\text{I}}\) and \(K_{\text{II}}\) are the mode one and two stress intensity factors, respectively.

The domain near the fracture tip is divided into two sections by one planar plane which perpendicular to the fracture tip. Following discretization, the fracture propagation direction is computed based on the stress intensity factors. Then, another perpendicular plane to the unit vector showing the direction of fracture propagation is generated. The LSF parameters are then updated for the nodes behind the newly generated plane (Fig. 5).

Fig. 5
figure 5

Identifying the nodes for which LSF is to be updated (Youn 2016)

Full size image

Initial stress intensity is defined based on the analytical computation as provided in Eq. 15.

$$K_{I} = \sigma \sin^{2} \left( {90 - \theta_{f} } \right)\sqrt {\frac{\pi L}{2}}$$
(15)
$$K_{II} = \sigma \sin \left( {90 - \theta_{f} } \right)\cos \left( {90 - \theta_{f} } \right)\sqrt {\frac{\pi L}{2}}$$
(16)

where \(\theta_{\text{f}}\) is the fracture angle relative to the coordinate system and \(L\) is the fracture length.

For the dynamic fracture behaviour, the SIFs are derived based on Eq. 17 which is later solved to identify the SIFs for the XFEM model provided in Eqs. 18 and 19.

$$G = \frac{{K_{\text{I}}^{2} + K_{\text{II}}^{2} }}{{E_{\text{ff}} }}$$
(17)

where \(G\) is the energy release rate and \(E_{ff}\) is the effective Young’s modulus.

$$K_{\text{I}}^{\text{XFEM}} = \frac{{I_{I}^{\text{int}} E_{\text{eff}} }}{2}$$
(18)

where \(K_{\text{I}}^{\text{XFEM}}\) is the mode I stress intensity factor for XFEM and \(I_{\text{I}}^{\text{int}}\) is the interaction state for mode I fracture propagation in XFEM.

$$K_{\text{II}}^{\text{XFEM}} = \frac{{I_{\text{II}}^{\text{int}} E_{\text{eff}} }}{2}$$
(19)

where \(K_{\text{II}}^{\text{XFEM}}\) is the mode II stress intensity factor for XFEM and \(I_{\text{II}}^{\text{int}}\) is the interaction state for mode II fracture propagation in XFEM.

The assumption present in the derivation of Eqs. 18 and 19 are that the auxiliary SIFs are exhibiting no effect on SIFs computed.

XFEM modelling in abaqus: general steps

The main steps in modelling the XFEM is shown by Fig. 6. There will be a few assumptions done prior to building the model. The assumptions are:

Fig. 6
figure 6

General steps in XFEM modelling

Full size image
  1. (i)

    Fracturing fluid is incompressible.

  2. (ii)

    Fracturing fluid is assumed to be liquid base and not foam base.

  3. (iii)

    Proppant will not be simulated.

The creation of XFEM model begin with the creation of geometry. The geometry of the formation and existing fractures are created at the part section. The geometry is drawn initially at 2D and then extruded to form the 3D object. The fractured formation and existing fracture are defined at this stage. There will be SRV and Non-SRV region in the model. At the assign properties section the physical properties of the formation and fracture are introduced. The formation is defined by its mechanical properties such as density, Young’s modulus, Poison’s ratio, porosity, permeability and damage criteria. Each part will have different properties and will be named accordingly. Following that, section is created. Section is the set which contains one of the materials later to be used for the assignation of the materials to the respective part. The model is then assembled. Within assembly, all the parts are combined to produce single system where all constituents will be simulated. After assembling the model, three steps, namely initial step, geostatic step and pumping step need to be created. In the initial step, initial and boundary and conditions will be applied. The initial condition refers to the formation properties such as saturation, stresses and pore pressure. The boundary conditions refers to the symmetry, restriction of displacement and rotation over certain coordinates. In the geostatic step, loads and geomechanical properties are introduced onto the formation. In the pumping step, the injection of the fracturing fluid is defined. Following pumping step, interaction between the components of the system is then introduced. The fracture surface and propagation algorithm are defined within this step. This stage also includes the introduction of the functions to define the enrichment nodes within the simulation. The nodes to be enriched are identified based on LSF and stress intensity factor (SIF). Meshing is then applied on the model. Mesh is the gridding which discretizes the whole part into sub-segments which will be computed separately later. There are several mesh types available within the software, hex, hex-dominated, tet and wedge. For the simulation of the FEM tet type of mesh must be selected as it the mesh compatible with remeshing rule. However, for XFEM hex type mesh can be implemented. Hex type mesh produced lesser number of elements for the same geometry than tet. Lastly, job file is created. At this section, simulation is initiated by including required input data which was defined in previous sections. Within this section, user defines the step size and step number to be simulated if required to differ from default values.

In summary, a review has been made on hydraulic fracturing modelling through analytical and numerical method. The basic theory of the methods was presented and comparison are made between the strength and weaknesses based on published studies. Analytical method can be classified into 2D modelling and 3D modelling. Two-dimensional modelling has an advantage of obtaining fast solution due to its simplistic calculation, but it is limited to model constant fracture height. This limitation is removed under 3D modelling; however, it come at the expenses of longer computational time due to higher complexity. Numerical method can be classified as continuum modelling and discontinuum modelling, and under it, XFEM has the greatest advantage in simulating hydraulic fracturing. The significance of XFEM is its ability to simulate arbitrarily propagating fracture, whereas other methods require the fracture trajectory to be pre-defined. This paper has included general steps to create XFEM model for the purpose of simulating hydraulic fracturing.