# Numerical investigation of the fluid lag during hydraulic fracturing

## Abstract

### Purpose

The purpose of this paper is to systematically investigate the fluid lag phenomena and its influence in the hydraulic fracturing process, including all stages of fluid-lag evolution, the transition between different stages and their coupling with dynamic fracture propagation under common conditions.

### Design/methodology/approach

A plane 2D model is developed to simulate the complex evolution of fluid lag during the propagation of a hydraulic fracture driven by an impressible Newtonian fluid. Based on the finite element method, a fully implicit solution scheme is proposed to solve the strongly coupled rock deformation, fluid flow and fracture propagation. Using the proposed model, comprehensive parametric studies are performed to examine the evolution of fluid lag in various geological and operational conditions.

### Findings

The numerical simulations predict that the lag ratio is around 5% or even lower at the beginning stage of hydraulic fracture under practical geological conditions. With the fracture propagation, the lag ratio keeps decreasing and can be ignored in the late stage of hydraulic fracturing for typical parameter combinations. On the numerical aspect, whether the fluid lag can be ignored depends not only on the lag ratio but also on the minimum mesh size used for fluid flow. In addition, an overall mixed-mode fracture propagation factor is proposed to describe the relationship between diverse parameters and fracture curvature.

### Research limitations/implications

In this study, relatively simple physical models such as linear elasticity for solid, Newtonian model for fluid and linear elasticity fracture mechanics for fracture are used. The current model does not account for such effects like leak off, poroelasticity and softening of rock formations, which may also visibly affect the fluid lag depending on specific reservoir conditions.

### Originality/value

This study helps to understand the effect of fluid lag during hydraulic fracturing processes and provides numerical experience in dealing with the fluid lag with finite element simulation.

## Keywords

#### Citation

Chen, B., Cen, S., Barron, A., Owen, D. and Li, C. (2018), "Numerical investigation of the fluid lag during hydraulic fracturing", *Engineering Computations*, Vol. 35 No. 5, pp. 2050-2077. https://doi.org/10.1108/EC-02-2018-0087

### Publisher

:Emerald Publishing Limited

Copyright © 2018, Emerald Publishing Limited

## 1. Introduction

Since the first field test performed in 1947 on a gas well in the Hugoton field (Carter *et al.*, 2000), hydraulic fracturing has become a vital step to achieve commercial or higher production from low permeability rock formations, in particular for tight oil and unconventional gas. To better understand the underlying mechanism and to support the operational practice, a number of analytical and computational techniques have been proposed to simulate the hydraulic fracturing process. Some of these hydraulic fracturing models have been successfully applied in industry, but in many situations, the modeling capacity remains far behind engineering practice, among which is the evolution of fluid lag during the fracturing process of rock formations. Hence, the aim of this work is to develop a versatile and accurate computational model and systematically investigate the fluid lag and its influence on the fracturing process. The development of hydraulic fracturing models is briefly recapped below, followed by a closer examination of the fluid-lag problem.

The earliest hydraulic fracturing models were developed in the 1950s (Zheltov, 1955), following a complex analysis approach with a plane-strain assumption. Two classical models were proposed during the 1960s and the early 1970s: the Perkins–Kern–Nordgren (PKN) model (Nordgren, 1972; Perkins and Kern, 1961) and the Khristianovic–Geertsma–de Klerk (KGD) model (Zheltov, 1955; Geertsma and De Klerk, 1969). The PKN model assumes an elliptical plane-strain fracture in each vertical section and a constant fracture height along the propagation direction. The maximum fracture width in a vertical section is determined by the local fluid pressure and confining stress, and 1D fluid flow along the fracture is simulated. This approach was later extended to the pseudo 3D (P3D) model (Cleary, 1980; Settari and Cleary, 1986) and the planar 3D (PL3D) model (Clifton and Abou-Sayed, 1979; Clifton and Abou-Sayed, 1981). The constant-height assumption is removed in P3D model, while 2D elements were introduced to compute the rock deformation and fluid flow in PL3D model. As for the KGD model, the plane-strain assumption is taken for the horizontal section, the Poiseuille’s law is adopted to describe the one-dimensional fluid flow within the fracture and the elasticity theory is used to describe the relationship between fluid pressure and fracture width. With parameter scaling analysis, Garagash (2006) and Lecampion and Detournay (2007) presented the toughness-dominated and viscosity-dominated propagation regimes of hydraulic fracturing. The influences from rock toughness, formation permeability, fracturing fluid viscosity and other relevant factors were investigated both analytically and numerically.

With further understanding of the field and experimental data, it is proved that hydraulic fracture may propagate in a more complicated manner than the prediction from the aforementioned models. More advanced physical models and diverse numerical methods have been introduced into the hydraulic fracturing simulation to meet the requirements. Physical models used for both solid and fluid phases as well as the fracture propagation criterion have all been improved over the past decade. Biot’s poroelastic theory is more popularly used to replace linear elasticity model to model the mechanical behavior of the rock in the reservoir, which can be regarded as a kind of porous media. For the fluid part, more complex rheology models such as non-Newtonian model have been introduced. Cohesive fracture models, which are more suitable for ductile fracture propagation than the linear elastic fracture mechanics is adopted in discrete fracture approach. Other fracture propagation criteria including strain threshold and Mohr-coulomb criterion are also more commonly used in the smeared fracture approach. On the numerical part, continuum-based methods, non-continuum-based methods and hybrid methods have all been reported. In continuum-based models, the fracture can be simulated with the discrete fracture approach or the smeared fracture approach. For the former one, the fracture surfaces need to be presented by the mesh. Typical numerical methods include finite element method (Secchi and Schrefler, 2012), partition-of-unity methods (e.g. generalized finite element method and extended finite element method) (Pereira *et al.*, 2010; Gupta and Duarte, 2014) and boundary element method (Xu, 1999; Yamamoto *et al.*, 2004; Rungamornrat *et al.*, 2005). With the smeared fracture approach, a hydraulic fracture can be simulated using a fixed mesh by using finite element method (Li *et al.*, 2012; Wangen, 2013) or finite difference method. A significant advantage of smeared fracture approach over discrete fracture approach is that the complex treatment for mesh when fractures intersect with each other is bypassed. Other methods including phase-field models (Lee *et al.*, 2016) have also been adopted to simulate 3D hydraulic fracturing. The non-continuum-based methods typically represent the material with discrete particles and bonds between them (Huang *et al.*, 2013). The fracture propagation is simulated by the breakage of internal bonds. Finite/discrete element method has been proved to be an effective numerical method suitable for rock mechanics and has also been used to simulate hydraulic fracturing (Profit *et al.*, 2016).

The existence of fluid lag has long been recognized. It was considered in Zheltov (1955) that the fracturing fluid might not fill the whole fracture length during a fracturing process, resulting a fluid lag between the fluid front and the fracture tip, and they predicted a strong dependence between the fracture width, the pressure profile and the fluid lag status, which was also confirmed in Geertsma and De Klerk (1969). A much simplified static model was proposed in Advani *et al.* (1997) to relate the fluid lag to the excess pressure encountered in practice, but the dynamic interaction between fluid flow and fracture propagation was ignored. Over the past two decades, the fluid lag effect has been studied in more detail by several research groups. For an ideal state with zero confining stress (Garagash, 2006), an analytical solution was obtained based on the KGD model, showing the lag-fracture ratio remains a constant; later, numerical solutions were obtained under different dimensionless toughness values (Lecampion and Detournay, 2007) showing a monotonic decrease of the lag-fracture ratio. These KGD-based predictions were recently verified by independent finite element method (FEM) simulations (Hunsweck *et al.*, 2013). Using the boundary element method approach, the fluid lag has been considered in conjunction with natural fractures (Zhang and Jeffrey, 2006; Zhang *et al.*, 2009) and multiple hydraulic fractures (Zhang *et al.*, 2011). More recently, based on the KGD model and using the FEM approach, the evolution of fluid lag during fracture initiation was investigated (Shen, 2014; Bao *et al.*, 2016).

Although the existence of fluid lag has long been recognized, the majority of hydraulic fracturing simulations do not model the fluid lag explicitly. In this case, negative fluid pressure may arise near the fracture tip, which is interpreted as the fluid lag (Secchi and Schrefler, 2012; Devloo *et al.*, 2006; Samimi and Pak, 2016; Simoni *et al.*, 2008; Secchi *et al.*, 2007; Cao *et al.*, 2017; Milanese *et al.*, 2016). Several research groups have included the fluid lag in their models but mostly limited to special conditions or relatively simple cases. None of the aforementioned studies have involved the full spectrum of fluid-lag evolution, i.e. initiation, growth, decrease and vanishment coupled with dynamic fracture propagation. It remains unclear to the academic and industrial communities when and where the fluid lag becomes important and how much it affects fracturing operations. The systematic investigation based on comprehensive parameters in the context of the lab experiment and field practice are still missing. Hence, the aim of this work is to systematically investigate the fluid lag phenomena and its influence in the hydraulic fracturing process, including all stages of fluid-lag evolution, the transition between different stages and their coupling with dynamic fracture propagation under common conditions. To achieve this, an accurate and versatile 2D finite element model is developed to simulate the fluid-lag evolution and its effect on hydraulic fracture propagation in different geological and operational conditions. The rest of the paper is organized as follows. In Section 2, the formulation of the strongly coupled solution scheme is derived following the mathematical definition of the problem. More details of the numerical treatment and implementation are given in Section 3. In Section 4, the proposed model is verified by comparing with semi-analytical solutions and numerical benchmarks available in the literature, after which a series of carefully designed numerical studies are presented in Section 5, aiming to provide a clear and complete picture for the fracture and fluid lag evolution under different practical conditions. Finally, concluding remarks are made in Section 6.

## 2. Governing equations and FEM solution

As shown in Figure 1, an extended KGD model is adopted, where the lag evolution is studied on the plane orthogonal to the vertical wellbore. The rock formation is assumed to be linear elastic, under plane-strain conditions and impermeable. The fracturing fluid is modeled as an incompressible Newtonian fluid with a laminar-flow assumption. An initial fracture with a length *L*_{0} and an orientation angle *α* with respect to the maximum horizontal stress is assumed. The FEM simulation domain is set as 50*L*_{0} by 50*L*_{0} to approximate the infinite rock formation. The stress boundary conditions are set according to the minimum and maximum horizontal stresses, while the influence of gravity is ignored, as it is orthogonal to the simulation plane. A constant injection flow rate is imposed on the injection point at the center of the model.

### 2.1 Governing equations

In the context of hydraulic fracturing analysis, linear elasticity and poroelasticity are the two most widely used constitutive models for rock formations. This study is based on the linear elastic model with small deformation assumption. The equilibrium and constitutive equations for quasi-static rock deformation are:

in which ** σ**,

**and**

*ε***denote the stress, strain and elastic tensor, respectively. For linear elastic rock formations,**

*C***is determined by the Poisson’s ratio**

*C**v*and the Young’s modulus

*E*of the rock.

As the width of a hydraulic fracture is much smaller than the other two dimensions, a lubrication theory, known as the Poiseuille’s law (or the cubic law) is commonly adopted to describe the momentum conservation of fracturing fluid:

*q*is the flow rate,

*w*the fracture width,

*μ*the viscosity of the fracturing fluid,

*p*the fluid pressure and

*s*is the local coordinate aligned with the tangential direction to the fracture path. The mass conservation for the one-dimensional fluid flow is expressed as:

*dt*denotes the time derivative. Substituting equation (3) into equation (4) yields:

The equilibrium equation (1) and the fluid equation (5) together form the governing equations for the hydraulic fracturing model in Figure 1. A fracture propagation criterion is required to determine:

when the fracture will propagate; and

where the fracture will propagate to.

Among the most popular fracture criteria are the maximum tensile-stress criterion and the minimum strain-energy density criterion. Based on the former one (Hossain and Rahman, 2008), a fracture will propagate when the following condition is reached:

*K*

_{I},

*K*

_{II}and

*K*

_{Ic}are the stress intensity factors for the Mode I fracture, the Mode II fracture and the fracture toughness, respectively, and the propagation direction

*θ*is determined by Rungamornrat

*et al.*(2005):

The interaction energy integral method is adopted to compute the stress intensity factors in this work (Walters *et al.*, 2005). Specifically, for a virtual unit fracture advance as shown in Figure 2, the energy release along the propagation direction is (Walters *et al.*, 2005):

*χ*is a scalar field with unit value at the fracture tip and a zero value at external boundaries,

*σ*and

*u*are the real stress and displacement fields respectively,

*σ*and

^{a}*u*are the axillary stress and displacement fields, respectively,

^{a}*δ*is Kronecker delta,

_{ij}*τ*is the fluid traction along the fracture surface,

_{i}*A*is the annulus near the fracture tip and

*C*

^{+}and

*C*

^{–}are the fracture surfaces. The stress intensity factors

*K*and

_{I}*K*are related to the energy release rate

_{II}*J*as:

*E*and

*ν*are the elastic modulus and Poisson’s ratio, respectively;

_{I}, and set

_{II}. In our numerical implementation, the second ring of the triangle mesh around the fracture tip is chosen as area A to calculate the first term in equation (8), and three Gaussian integration points are adopted to calculate the integral in each triangle while five-point Gaussian integration is used to compute the second term on related edges.

### 2.2 A strongly coupled FEM solution

In the above hydraulic fracturing model, the rock is assumed to be linear elastic and quasi-static undergoing small deformation that obeys the equilibrium equation (1), the motion of the fracturing fluid is modeled by the fluid equation (5) based on lubrication theory and rock fractures are described by a brittle fracture criterion equation (6). The rock displacement ** u** and the fluid pressure

*p*are coupled in equations (1) and (5), while equation (6) is decoupled from them. Thus, the hydraulic fracturing model is essentially described by

**and**

*u**p*, while other quantities such as the flow rate, strength intensity factors and fracture propagation direction can all be readily computed from the rock displacement and the fluid pressure.

We propose a finite element approach to solve this hydraulic fracturing model: the solid equation (1) is solved with an unstructured triangular mesh that is adaptively updated with the fracture propagation, and the fluid equation (5) is solved on a 1D linear mesh generated consistent with the 2D triangular mesh of the solid. More details about the FE mesh discretization are given in Section 3.1.

To solve the quasi-static equilibrium equation (1) on a FE mesh, we assume the nodal displacement **u*** _{m}* and the nodal pressure

**p**

*at time*

_{m}*t*are both known and have been represented on the current mesh. Then, at time

_{m}*t*

_{m}_{+1}, equation (1) can be discretized as:

**K**is the stiffness matrix,

**u**

_{m}_{+1}and

**p**

_{m}_{+1}are the nodal displacement and nodal pressure at

*t*

_{m}_{+1}respectively,

**F**

*is the nodal force at the domain boundary due to the confining stress of the rock formations, and*

_{external}**F**(

**p**

_{m}_{+1}) is the nodal force on the fracture due to the fluid pressure. Let

**u**

_{m}_{+1}=

**u**

*+Δ*

_{m}**u**, with Δ

**u**denoting the displacement increment from time

*t*to time

_{m}*t*

_{m}_{+1}, equation (10) can be reformulated as:

The stress boundary condition on the fracture can be expressed as:

**T**is a coefficient matrix computed from the 1D fluid mesh. For a linear mesh,

**T**is formed as:

*l*is the element length of the fluid mesh.

^{e}To accurately track the fluid front in relation to the fracture tip, we take an FE approach to solve the fluid equation (5) as well. Specifically, the weak form of equation (5) is:

*p*(0) and

*q*(0) are the fluid pressure and the flow rate at the injection point (i.e. fracture center), and

*p*(

*s*) and

_{f}*q*(

*s*) are the fluid pressure and the flow rate at the fluid front. A constant injection flow rate

_{f}*q*(0) is assumed at the center of the initial fracture. Depending on the lag condition, Dirichlet or Neumann boundary conditions need to be considered at the fluid front. If a fluid lag is present, the pressure at the fluid lag is set to zero, i.e.

*p*(

*s*) = 0; and due to the zero-width at the fracture tip, if the fluid front coincides with the fracture tip, the flow rate at the fracture tip is set to zero, i.e.

_{f}*q*(

*s*) = 0 where

_{t}*s*denotes the fracture tip. To accurately capture the evolution of fluid lag, these two types of boundary conditions must be adaptively set according to the current fluid lag condition.

_{t}Discretizing equation (14) with 1D linear finite elements yields:

*and w*

_{m+1}*are the nodal fracture widths at time t*

_{m}*and t*

_{m+1}*, respectively, L is a coefficient matrix, h(w*

_{m}*, p*

_{m+1}*) is a nonlinear vector function with respect to the nodal fracture width w*

_{m+1}*and the nodal fluid pressure p*

_{m+1}*and q*

_{m+1}*is the nodal flow rate at time t*

_{m+1}*. The fracture path is discretized into*

_{m+1}*N*linear elements of lengths {

*l*

_{1},

*l*

_{2},…l

_{N}} indexed from the fracture center to the fracture tip, and the first N

_{f}elements are occupied by the fracturing fluid. The matrix quantities in equation (15) are defined as follows:

*N*<

_{f}*N*, there exists a fluid lag, and thus the pressure at the fluid front is set to zero, i.e.

*N*=

_{f}*N*, there is no fluid lag, and thus the fluid front overlaps with the fracture tip, i.e.

*q*(

*N*

_{f}_{+1}) = 0.

The solid equation (1) and the fluid equation (5) are discretized into finite element equations (11) and (15), respectively, where the nodal fracture width **w** can be directly represented by the nodal displacement **u**. As equations (11) and (15) are nonlinear and strongly coupled, the Newton–Raphson solution scheme is adopted:

*n*denotes the iteration step,

**M**

*is the Jacobi matrix and*

_{Jocobi}**R**is the residual vector. The Jacobi matrix and the residual vector are given below:

The fluid lag condition affects the coupling between the solid and fluid equations, and hence the formulation of matrix **J**. Specifically, if *N _{f}* =

*N*,

**J**=

**T**; if

*N*<

_{f}*N*,

**J**=

**T**(:, 1:

*N*) and other matrix quantities are accordingly modified by removing the entries associated with the flow rate at the fluid front.

_{f}## 3. Computational issues

Following the mathematical formulation explained in Section 2, this section provides details on several computational issues that affect the implementation and performance of the proposed approach.

### 3.1 Adaptive mesh update

Adaptive remeshing is necessary for fracture simulation based on the standard finite element method, while the specific meshing strategies can differ greatly dependent on the geometry and element type (Figure 3). In this study, to support systematic investigation of fluid lag, it is essential to allow fractures to propagate along arbitrary directions and in arbitrary geometries. In addition, to obtain reliable conclusions, it is critical to ensure the highest analysis accuracy for the rock-fluid interaction. To meet these requirements, we adaptively update the FE meshes with the fracture propagation, during which the nodal-level consistency is maintained between solid and fluid meshes. A convergence test is performed to ensure the mesh density is sufficient for both fluid and solid analysis. The meshing strategies are listed below:

Fracture path determination: When fracture propagation is detected, the center line between the solid nodes on the current fracture surface is calculated, to represent the fracture path as a cubic spline.

Fluid mesh generation: The fluid mesh is generated along the newly calculated fracture path with linear elements, and the node distribution is based on the curvature of the fracture path, the distance to the fluid front and the associated computational cost.

Fracture surface determination: To ensure the consistency between the fluid and solid meshes, the solid nodes along the fracture surface are recalculated from the newly generated fluid mesh, by projecting from each fluid node along the perpendicular direction of the fluid path. For a higher accuracy, the fracture surface is also represented by cubic splines, consistent with the fracture path.

Fracture tip treatment: The fracture path is projected onto the intact rock formation by a predicted length and orientation, and bonded node pairs are placed along the new fracture surface to simulate fracture propagation.

Triangular mesh generation: Aligned with the newly generated solid nodes on the fracture surface, the unstructured triangular mesh is generated for the rock formation, taking into account the stress singularity at the fracture tip, the local geometries and the associated computational cost.

Mapping and re-calculation: After remeshing, the solution information needs to be transferred from the old mesh to the new mesh. First, on the fracture surface, the rock formation is subject to the fluid pressure represented as the nodal forces

**Ku**, which are transferred via a mapping globally smoothed with the same FE interpolation. Then, in the rock formation, the displacement field is recalculated following the newly mapped nodal forces on the fracture surface.

### 3.2 Advancement of the fluid front and fracture tip

One of the key challenges in this study is to accurately capture the positions of the fluid front and the fracture tip. The evolutions of the fluid front and the fracture tip are strongly coupled, due to the interaction between the rock deformation, fluid flow and fracture propagation processes. In addition, to enable systematic investigation of fluid lag, the computational strategy must allow initiation, development and vanishment of fluid lag as well as free transitions between different states. To meet these requirements, a robust and highly accurate advancement algorithm is proposed for the evolution of fluid front and fracture tip. The algorithm has two coupled steps: the fluid front advancement and the fracture tip advancement.

Given the time step Δ*t* and the fluid velocity
*i*th fracture, the fluid front is advanced by a distance of
*t* is determined as
*d* is the advancement limit specified by the user. In this case, the fluid-front advancement for different fracture branches is smaller than the pre-specified limit, and it keeps the simulation stable.

After the fluid front is predicted, we solve the coupled system equation (11) and (15) to obtain the present rock deformation and fluid condition. These state variables are then used to compute stress intensity factors. For any fracture branch, the fracture will propagate when the stress intensity factor K is greater than the given toughness *K _{Ic}*, where the position of the fracture tip (i.e. the propagation length and direction) is determined iteratively depending on the equilibrium of the rock-fluid interaction and the fracture criterion. Multiple advancements of fracture tip may occur in a single time step. During this process, the fluid front will advance together with the fracture tip if and only if the fluid pressure at the fracture tip (

*p*) is positive. An illustration of the proposed fracture propagation strategy is given in Figure 4.

_{tip}### 3.3 Algorithm flow

For clarity, the overall algorithm flow of the proposed simulation strategy is summarized below:

**Initial condition**: Given the initial fracture path and fracture toughness, an iterative approach is used to set the initial condition.

**Repeat**

**Fluid front update**: The fluid front is updated according to the strategy described in Section 3.2.

**Repeat**

**Strongly-coupled FE solution**: Corresponding to the current fluid front and fracture tip, the displacement field and fluid pressure are solved from equations (11) and (15), following the Newton–Raphson [equation (17)].

**Fracture tip update**: The fracture propagation direction and the fracture tip position are iteratively updated as described in Section 3.2.

**Adaptive meshing**: Both the fluid and the solid meshes are updated as necessary, taking into account the requirement of fracture tip advancement and the change of fracture profile. See Section 3.1 for details.

**Until** the SIF is less than the toughness for all fracture branches

**Until** the total simulation time reaches the preset time limit

## 4. Model verification

To check the correctness and accuracy of the proposed hydraulic fracturing model with the presence of fluid lag, we compare our numerical solutions with both analytical and computational results available in the literature and also present a convergence study as part of model verification. Two types of cases are examined here. First, the proposed method is verified by the semi-analytical solutions for early-time stage (Garagash, 2006) and late-time stages (Spence and Sharp, 1985), to prove its effectiveness in capturing fluid lag. Second, the capacity of mixed mode fracture simulation is verified, and a sensitivity analysis of mesh size is also conducted. Other aspects of the model including remeshing and projection strategy are also verified in the meantime.

### 4.1 Khristianovic–Geertsma–de Klerk verification

Within the framework of KGD model, the hydraulic fracture is assumed to propagate along a straight line perpendicular to the minimum principal stress, as shown in Figure 5(a). The solutions for specific propagation regimes can be interpreted in the so-called OMK parameter space, shown in Figure 5(b).

#### 4.1.1 Early-time solution.

In this example, the proposed solution is compared with the early-time solution from (Garagash, 2006). To keep the hydraulic fracturing stay in the early-time stage, *σ*_{1} and *σ*_{2} are both set as zero. In this case, the dimensionless confining stress (or dimensionless time) remains zero. Comparisons of the fracture profile, the pressure distribution and the solution state are shown in Figure 6(a), (b) and (c), respectively, where the solid lines (solid markers) denote our numerical solutions, and the dashed lines (hollow markers) represent the analytical solutions from (Garagash, 2006). Good agreement is observed in all cases. The relation between the fracture width and the fluid pressure and their dimensionless counterparts are given below:

*ξ*denotes the dimensionless coordinate along the fracture path,

*σ*

_{1}is the

*in situ*stress perpendicular to the initial fracture path,

_{m}=

*K*′(

*q*

_{0}

*E*′

^{3}

*μ*′)

^{−1/4},

*K*′ = 4(2/

*π*)

^{1/2},

*q*

_{0}=

*q*(0) is the injection flow rate,

*E*′ =

*E*/(1 –

*ν*

^{2}), and

*μ*′ = 12

*μ*.

#### 4.1.2 Transition stage solution.

In this example, the transient stages from the early-time solution to the late-time solution are checked. Specifically, the dimensionless toughness *K*_{m} is set as 0.5. To obtain the transient solutions corresponding to the lag-fracture ratios at 0.6, 0.7, 0.8, 0.9 and 1 – *η* (0 < *η* ≪ 1) in the same length scale, different confining stresses are adopted. The results are compared with the late-time solution from Hunsweck *et al.* (2013) and Spence and Sharp (1985) and the transition solution from Lecampion and Detournay (2007). Figure 7(a) and (b) compare the fracture profile and the pressure distribution, respectively, where the solid lines represent the transition solutions obtained with our numerical model, and the dashed line indicates the late-time solution from Hunsweck *et al.* (2013) and Spence and Sharp (1985). It can be seen that as the dimensionless time *T* elapses, our numerical solutions converge to the late-time solution in the literature. Figure 7(c) compares our numerical solutions at different time moments to the transition solution from Lecampion and Detournay (2007), and again good agreement is observed in all cases.

### 4.2 Mixed-mode fracture propagation

A mixed-mode fracture propagation driven by a uniform pressure distribution was simulated as a benchmark case in several literatures (Gupta and Duarte, 2014; Dong and de Pater, 2001), where a uniform pressure is assumed as an approximation of the late-time condition and the pressure remains constant during the fracture propagation. The initial fracture length is 0.04 m, and the initial fracture direction is 89^{°} with respect to the maximum horizontal stress. The rock properties are Young’s modulus *E* = 20 GPa, Poisson’s ratio *ν* = 0.2 and fracture toughness *K _{IC}* = 0.6 MPa·m

^{1/2}. The two principal horizontal stresses

*σ*

_{1}= 9.7 MPa and

*σ*

_{2}= 19.4 MPa. Similar to Dong and de Pater (2001), the initial fracture is discretized into eight uniform segments, and three different internal pressures (24.3, 29.1 and 38.8 MPa) are considered. The fracture trajectories are compared in Figure 8, where the lines represent the literature results, and the markers represent our solutions. The comparison shows our solutions match well with the literature results under the same mesh resolution. To check the sensitivity of the simulation with respect to the mesh resolution, a convergence study is also performed as shown in Figure 9, where the fracture trajectories obtained at different mesh resolutions are plotted. As shown in Figure 9, the mesh size of 1/8 initial length is not fine enough to give an accurate fracture path. A finer mesh is recommended to ensure the accuracy.

## 5. Numerical cases

Two sets of numerical studies are presented in this section. The first set of numerical experiments is designed to examine the flexibility and robustness of the proposed hydraulic fracturing model in capturing mixed-mode fracture with coupled fracture and fluid lag propagation. The second set of numerical experiments forms a comprehensive parametric study for mixed-mode hydraulic fracturing under various geological and operational conditions, and it is aimed to provide a clear and complete picture for fracture propagation taking into account of fluid lag evolution.

### 5.1 Straight fracture propagating into different rock zones

To examine the capacity of the proposed model in modeling different lag evolution scenarios, a test case with a straight fracture propagating into different rock zones is designed. As shown in Figure 10, the formation has two rock zones [−*L*_{0}, *L*_{0}] and (−*∞*, −*L*_{0}) ∪ (*L*_{0}, +*∞*), aligned in the direction of the maximum horizontal stress, and an initial fracture with a half-length of *L*_{0}/2 = 0.1 m is placed at the center of the formation. The Young’s modulus and Poisson’s ratio of the formation are set at the same value across both zones, at *E* = 25 GPa and *v* = 0.3, respectively, while the rock toughness varies in different zones. As illustrated in the OMK parameter space in Figure 11, four different lag evolution scenarios are examined, including lag decrease, increase, vanishment and initiation. The physical parameters including rock toughness, confining stress, injection rate and fracturing fluid viscosity are listed in Table I.

#### 5.1.1 Lag decrease and vanishment.

When a hydraulic fracture propagates into a tougher rock zone, the lag may decrease or even vanish and the normalized pressure and width profile can also change. Two cases are designed to examine the decreasing and vanishing scenarios. In the decreasing case, the confining stresses are set to zero, restricting the fracture propagation to the early-time stage. The toughness in the zone [−*L*_{0}, *L*_{0}] and zone (−*∞*, −*L*_{0}) ∪ (*L*_{0}, +*∞*) are set as 0.8879 and 1.2860 MPa·m^{1/2}, respectively, which corresponds to the dimensionless toughness 0.600 and 0.869, respectively. The initial relative position of the fluid front is determined to be 0.5. The width and pressure profiles are shown in Figure 12(a) and (b), while the evolutions of fracture tip and fluid front are shown in Figure 13. It can be observed that the fracture tip stops when it just arrives at the interface and then moves forward again later. During the transition process, the normalized width and pressure profile changes from low-toughness mode to high-toughness mode, as shown in Figure 12(c) and (d). After the fracture tip propagating into the high-toughness zone, the relative position of fluid front changes from 0.6 to 0.8, as shown in Figure 13.

In the lag vanishing case, the two far-field stresses are set as 19.4 and 9.7 MPa, same as the case in Section 4.2. The initial relative position of the fluid front is determined to be 0.8. The evolutions of width and pressure profiles are shown in Figure 14, while the evolutions of the fracture tip and fluid front are shown in Figure 15. It can be observed that the fracture tip starts to move when the relative fluid front reaches about 0.90 and the solution gets into a stable mode several steps later. After that, the lag ratio reduces gradually until the fracture tip approaches to the interface. The increase of rock toughness prevents the fracture tip from moving forward temporally, which leads to a significant decrease of lag ratio. After crossing the interface, the fluid front keeps coincide with the fracture tip and a finite pressure around 5 MPa is observed at the fracture tip.

#### 5.1.2 Lag increase and initiation.

In contrast to the lag decreasing and vanishing cases discussed above, the lag may increase or initiate when a hydraulic fracture propagates into a softer zone. Similarly, two cases are designed to examine the increasing and initiating scenarios separately. For the increasing case, the same toughness parameters as in decreasing case are adopted, but the soft and tough zones are swopped. The initial relative position of the fluid front is determined as 0.7, and the results of fracture and lag evolution are shown in Figures 16 and 17. It can be observed that the fracture tip moves forward a large distance after it arrives at the interface, and during the transition process, the normalized width and pressure profiles change from the high-toughness mode to the low-toughness mode in one step. After the fracture tip propagating into the low-toughness zone, the relative position of the fluid front decreases from 0.8 to 0.6.

For the lag initiating case, the same parameters as in the lag vanishing case are adopted, but the soft and tough zones are swopped. The initial lag ratio is set as 0.2, and the evolutionary history of fracture and fluid lag is shown in Figures 18 and 19. It can be observed that a process of lag vanishment is detected first, after which the fluid front moves forward together with the fracture tip until it reaches the interface. As a result of toughness drop, the fluid lag reappears after the fracture propagating into the soft zone and after some fluctuations in the following three steps, the lag ratio starts to reduce gradually (Figure 19).

### 5.2 Parametric study on mixed-mode fracture propagation

To understand how the fluid lag initiates, evolves and disappears under different geological and operational conditions, a comprehensive parametric study is presented in this section. The initial fracture length is set as *L*_{0} = 0.2 m, and the orientation angle *α* with respect to the maximum horizontal stress is taken to be 89° to consider the influence of anisotropic *in situ* stress. The simulation parameters include:

*geological conditions*:*in situ*stress*σ*_{1}and*σ*_{2}, rock toughness*K*, Young’s modulus_{Ic}*E*and Poisson’s ratio*ν*; and*Operation parameters*: fracturing fluid viscosity*μ*and injection flow rate*q*_{0}.

The influence of *in situ* stress, rock toughness and fracturing fluid viscosity are investigated first, by fixing *E* = 25 GPa, *ν* = 0.3 and *q*_{0} = 0.002 m^{2}/*s*. The fluid lag fractions when the hydraulic fracture propagates from 0.2 to 0.6 m are listed in Table II. The viscosity of fracturing fluid is set as 2, 20 and 500 cp, corresponding to slick water (2-3 cp), linear gel (10-30 cp) and cross-linked gel (100-1000 cp), respectively. The *in situ* stresses are set as (9.7 MPa, 19.4 MPa), (24.7 MPa, 38.0 MPa) and (43.5 MPa, 55.0 MPa), based on typical values reported for lab experiments and field measurement. The toughness of shale is around 0.5-3 MPa·m^{1/2}. As expected, higher fracturing fluid viscosity, lower *in situ* stress and lower rock toughness contribute to a more significant fluid lag, and vice versa. As shown in Table II, there is almost no fluid lag when using slick water and linear gel to fracture rock with a toughness over 1.5 MPa·m^{1/2} under realistic in-situ stress conditions. But when using the cross-linked fracturing fluid, a non-negligible fluid lag is observed in the current length scale even for a very high in-situ stress and rock toughness.

To investigate the influence of injection flow rate, the fracturing fluid viscosity, *in situ* stress and rock toughness are fixed at 20 cp, (24.7 MPa, 38 MPa) and 3 MPa·m^{1/2} (denoted by U20-S38-K3), while the injection flow rate is doubled and halved in two separate simulations. In both cases, there is no fluid lag detected. A lower fluid pressure at fracture tip is observed for a higher injection flow rate, which implies that the fluid lag is prone to initiate at a higher injection flow rate. To investigate the influence of Young’s modulus, the fracturing fluid viscosity, *in situ* stress and rock toughness are fixed at U500-S55-K1.5, while the Young’s modulus is set to 25 GPa, 30 GPa and 35 GPa in three separate simulations. The corresponding fluid lag fraction increase from 1.2-4.5 per cent to 2.2-6.2 per cent and then to 2.2-9.0 per cent.

It is noted that the fluid lag dimension is not only related to the aforementioned five parameters but also related to the current length of the hydraulic fracture. To investigate the influence of fracture length, a series of test cases are chosen from Table II for extended simulations. To reduce the associated computational cost, the fracture propagation is redistricted to a straight line along the direction of maximum horizontal stress. The evolution of fluid lag with respect to the increase of fracture length is shown in Table III. During our numerical experiments, it is also found that ignoring a fluid lag which should be presented in a numerical simulation would lead to extremely negative pressure at the fracture tip, which is nonphysical and significantly affect the computational accuracy of stress intensity factors and worse, affecting the simulation convergence. On the numerical aspect, the appearance of fluid lag is also sensitive to the minimum mesh size. A very small fluid lag described by a fine mesh may disappear when using a coarse mesh. In the parametric study, a minimum mesh size L_{0}/60 is adopted first; if the mesh is too coarse to capture the state of fluid lag, the minimum mesh size is halved adaptively.

Except for the evolution of fluid lag, the aforementioned geological and operational parameters also have a significant influence on the fracture path. Some typical fracture paths are plotted in Figure 20. It can be observed in Figure 20(a) that increasing the fracturing fluid viscosity prevents the significantly curved fracture path caused by anisotropic *in situ* stresses. As shown in Figure 20(b), the most significant bending occurs at the *in situ* stress (38 MPa, 24.7 MPa). Corresponding to the three *in situ* stress cases, the differential stress is 9.7 MPa, 13.3 MPa and 11.5 MPa, respectively, and the stress ratio is 1:2.00, 1:1.53 and 1:1.26, respectively. The trend of the fracture path matches with the differential stress, instead of the *in situ* stress or the stress ratio. As shown in Figure 20(c), (d) and (e), the fracture path is straighter at higher rock toughness, higher flow rate and higher Young’s modulus.

The hydraulic fracture path is a direct result of the complex interaction between *in situ* stress and fluid pressure, which is determined by the rock toughness, fracturing fluid viscosity, injection flow rate, etc. The computational simulation can quantitatively predict the fracture path, but for the practical operation, it is desirable to have a simpler measure for quick evaluation. To achieve this, a much simplified case with the inviscid fracturing fluid and uniform pressure distribution is first considered. In this case, the fracture path is determined mainly by the rock toughness, *in situ* stress and initial fracture, and a classic analytical solution from linear elasticity can be adopted to approximately evaluate the fracture path with the following dimensionless factor (Zhang *et al.*, 2011):

Intuitively, the influence from the non-uniform fluid pressure distribution is related to the fluid viscosity, injection flow rate and Young’s modulus. The fourth power of the dimensionless toughness is adopted to account for the fluid effect:

By combining the solid factor *β _{s}* and the fluid factor

*β*in various ways and checking against the parametric studies, we propose an overall mixed-mode fracture propagation factor as follows:

_{f}To examine the effectiveness of the above fracture propagation measure, four numerical cases with similar *β* values (around 0.022), but very different geological/operational parameters are plotted in Figure 21. It can be seen that despite the significant difference in physical parameters, the fracture paths in all four cases are very similar.

## 6. Conclusions

We proposed a plane 2D hydraulic fracturing model that is capable to accurately capture the fluid lag evolution during the fracturing process. The rock formation is considered as linear elastic, the fracturing fluid is treated as an impressible Newtonian fluid and the fracture propagation is depicted with the linear elastic fracture mechanics. A fully implicit numerical solution scheme is formulated, where the solid and the fluid are solved together in a monolithic finite element scheme. The proposed approach is highly accurate and versatile, and it can capture the full spectrum of lag evolution statuses at various geological and operational conditions. The key technical contributions include:

a novel coupling strategy to accurately track the fracture tip and the fluid front during fracture propagation;

a robust fracture propagation model to enable efficient and highly accurate simulation of mixed-mode hydraulic fractures in arbitrary paths with strong coupling between the rock formation and the fracturing fluid; and

systematic numerical investigation of fluid-lag evolution and its interaction with fracture propagation in various geological and operational conditions.

The numerical simulations predict that the lag ratio is around 5 per cent or even lower at the early stage of the hydraulic fracture in high *in situ* stress conditions (around 50 MPa) and high fracturing fluid viscosity (500 cp). This value will be lower when using a fracturing fluid with lower viscosity but higher (around 20-30 per cent) when lower confining stress is applied in the context of the lab experiments. In addition, with the increase of the fracture length, the lag ratio keeps decreasing and can be ignored at around 10 m for typical geological conditions. In the late stage of hydraulic fracturing, the effect of the fluid lag would only restrict at the area close to the fracture tip. On the numerical aspect, whether the fluid lag can be ignored depends not only on the lag ratio but also the minimum mesh size used for fluid flow. A too small mesh size and high fluid lag ratio tend to result in an extremely negative fluid pressure at fracture tip and convergence problem when the fracture is assumed to be fully filled with fluid and is simulated by the linear elastic fracture mechanics. Based on the numerical results, an overall mixed-mode fracture propagation factor is proposed to describe the relationship between diverse parameters and fracture curvature.

In this study, relatively simple physical models such as linear elasticity for solid, Newtonian model for fluid and linear elasticity fracture mechanics for fracture are adopted for simplicity; therefore, the results are not without limitation. However, considering the form of fluid lag more relies on the coupled system between solid deformation, fluid flow inside fracture and fracture propagation, the results in this study are still reliable. More advanced physical models such as poroelasticity, non-Newtonian fluid model and cohesive fracture model will be introduced into the current model in the future work.

## Acknowledgement

The authors would like to thank the support from the European Community’s Seventh Framework Programme (Marie Curie International Research Staff Exchange Scheme, Grant No. 612607), the Sêr Cymru National Research Network in Advanced Engineering and Materials, the China Scholarship Council, the Welsh Government Sêr Cymru Programme, the Robert A. Welch Foundation (C-0002) and the Royal Academy of Engineering.

## Figures

#### Figure 6.

A comparison between our numerical solution (solid lines and solid markers) and the early-time solution in Garagash (2006)

#### Figure 7.

A comparison between our numerical solutions and the late-time solution (Hunsweck *et al.*, 2013; Spence and Sharp, 1985) and the transition solution (Lecampion and Detournay, 2007)

Simulation parameters for straight fracture propagation in different rock zones

Case index | K_{a} ( MPa·m^{1/2}) |
K_{b} (MPa·m^{1/2}) |
σ_{1} (MPa) |
σ_{2} (MPa) |
q_{0} (m^{2}/s) |
μ (cP) |
---|---|---|---|---|---|---|

① Decrease | 0.888 | 1.286 | 0 | 0 | 0.002 | 1 |

② Increase | 1.286 | 0.888 | 0 | 0 | 0.002 | 1 |

③ Vanishment | 1.5 | 3.0 | 9.7 | 19.4 | 0.002 | 15 |

④ Initiation | 3.0 | 1.5 | 9.7 | 19.4 | 0.002 | 15 |

The fraction of fluid lag under different conditions

Fracturing fluid viscosity μ (cp) | In-situ stress σ_{1} and σ_{2} (MPa) |
Rock toughness K (MPa·m_{Ic}^{1/2}) |
||
---|---|---|---|---|

0.5 (%) | 1.5 (%) | 3.0 (%) | ||

2 | (9.7, 19.4) | 2.5∼4.0 | Zero | Zero |

(24.7, 38.0) | 0∼1.7 | Zero | Zero | |

(43.5, 55.0) | 0∼1.7 | Zero | Zero | |

20 | (9.7, 19.4) | 6.4∼9.1 | 2.3∼4.1 | Zero |

(24.7, 38.0) | 1.7∼3.2 | Zero | Zero | |

(43.5, 55.0) | 0∼1.7 | Zero | Zero | |

500 | (9.7, 19.4) | 23.3∼33.3 | 10.2∼18.3 | 5.7∼10.7 |

(24.7, 38.0) | 11.7∼17.7 | 3.7∼7.0 | 2.2∼6.1 | |

(43.5, 55.0) | 2.2∼6.0 | 1.2∼4.5 | 0.7∼1.7 |

The fluid lag fractions at different fracture lengths

Numerical parameters | 0.1 m (%) | 1m (%) | 10 m (%) | 100 m (%) |
---|---|---|---|---|

U20-S38-K0.5-E25-Q0.002 | 4.7 | Zero | Zero | Zero |

U500-S38-K0.5-E25-Q0.002 | 18.2 | 2.3 | Zero | Zero |

U500-S38-K1.5-E25-Q0.002 | 13.4 | 2.0 | Zero | Zero |

U500-S55-K0.5-E25-Q0.002 | 8.0 | 1.6 | Zero | Zero |

## References

Advani, S.H., Lee, T.S., Dean, R.H., Pak, C.K. and Avasthi, J.M. (1997), “Consequences of fluid lag in three-dimensional hydraulic fractures”, International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 21 No. 4, pp. 229-240.

Bao, J.Q., Fathi, E. and Ameri, S. (2016), “A unified finite element method for the simulation of hydraulic fracturing with and without fluid lag”, Engineering Fracture Mechanics, Vol. 162, pp. 164-178.

Cao, T.D., Milanese, E., Remij, E.W., Rizzato, P., Remmers, J.J.C., Simoni, L., Huyghe, J.M., Hussain, F. and Schrefler, B.A. (2017), “Interaction between crack tip advancement and fluid flow in fracturing saturated porous media”, Mechanics Research Communications, Vol. 80, pp. 24-37.

Carter, B.J., Desroches, J., Ingraffea, A.R. and Wawrzynek, P.A. (2000), “Simulating fully 3D hydraulic fracturing”, Modeling in Geomechanics, Wiley Publishers, New York, NY, pp. 525-557.

Cleary, M.P. (1980), “Comprehensive design formulae for hydraulic fracturing”, *SPE Annual Technical Conference and Exhibition*, Dallas, TX.

Clifton, R.J. and Abou-Sayed, A.S. (1979), “On the computation of the three-dimensional geometry of hydraulic fractures”, *Symposium on Low Permeability Gas Reservoirs*, Colorado.

Clifton, R.J. and Abou-Sayed, A.S. (1981), “A variational approach to the prediction of the three-dimensional geometry of hydraulic fractures”, *SPE/DOE Low Permeability Gas Reservoirs Symposium*, Colorado.

Devloo, P.R.B., Fernandes, P.D., Gomes, S.M., Bravo, C.M.A.A. and Damas, R.G. (2006), “A finite element model for three dimensional hydraulic fracturing”, Mathematics and Computers in Simulation, Vol. 73 Nos 1/4, pp. 142-155.

Dong, C.Y. and de Pater, C.J. (2001), “Numerical implementation of displacement discontinuity method and its application in hydraulic fracturing”, Computer Methods in Applied Mechanics and Engineering, Vol. 191 Nos 8/10, pp. 745-760.

Garagash, D.I. (2006), “Propagation of a plane-strain hydraulic fracture with a fluid lag: early-time solution”, International Journal of Solids and Structures, Vol. 43 Nos 18/19, pp. 5811-5835.

Geertsma, J. and De Klerk, F. (1969), “A rapid method of predicting width and extent of hydraulically induced fractures”, Journal of Petroleum Technology, Vol. 21 No. 12, pp. 1571-1581.

Gupta, P. and Duarte, C.A. (2014), “Simulation of non-planar three-dimensional hydraulic fracture propagation”, International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 38 No. 13, pp. 1397-1430.

Hossain, M.M. and Rahman, M.K. (2008), “Numerical simulation of complex fracture growth during tight reservoir stimulation by hydraulic fracturing”, Journal of Petroleum Science and Engineering, Vol. 60 No. 2, pp. 86-104.

Huang, K., Zhang, Z. and Ghassemi, A. (2013), “Modeling three-dimensional hydraulic fracture propagation using virtual multidimensional internal bonds”, International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 37 No. 13, pp. 2021-2038.

Hunsweck, M.J., Shen, Y. and Lew, A.J. (2013), “A finite element approach to the simulation of hydraulic fractures with lag”, International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 37 No. 9, pp. 993-1015.

Lecampion, B. and Detournay, E. (2007), “An implicit algorithm for the propagation of a hydraulic fracture with a fluid lag”, Computer Methods in Applied Mechanics and Engineering, Vol. 196 Nos 49/52, pp. 4863-4880.

Lee, S., Wheeler, M.F. and Wick, T. (2016), “Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model”, Computer Methods in Applied Mechanics and Engineering, Vol. 305, pp. 111-132.

Li, L.C., Tang, C.A., Li, G., Wang, S.Y., Liang, Z.Z. and Zhang, Y.B. (2012), “Numerical simulation of 3D hydraulic fracturing based on an improved flow-stress-damage model and a parallel FEM technique”, Rock Mechanics and Rock Engineering, Vol. 45 No. 5.

Milanese, E., Rizzato, P., Pesavento, F., Secchi, S. and Schrefler, B.A. (2016), “An explanation for the intermittent crack tip advancement and pressure fluctuations in hydraulic fracturing”, Hydraulic fracturing Journal, Vol. 3 No. 2, pp. 30-43.

Nordgren, R.P. (1972), “Propagation of a vertical hydraulic fracture”, Society of Petroleum Engineers Journal, Vol. 12 No. 4, pp. 306-314.

Pereira, J.P., Duarte, C.A. and Jiao, X. (2010), “Three-dimensional crack growth with HP-generalized finite element and face offsetting methods”, Computational Mechanics, Vol. 46 No. 3, pp. 431-453.

Perkins, T.K. and Kern, L.R. (1961), “Widths of hydraulic fractures”, Journal of Petroleum Technology, Vol. 13 No. 9, pp. 937-949.

Profit, M.L., Dutko, M., Yu, J., Armstrong, J. and Parfitt, D. (2016), “Application of state of the art hydraulic fracture modelling techniques for safe-optimized design and for enhanced production”, *50th US Rock Mechanics/Geomechanics Symposium*, Houston.

Rungamornrat, J., Wheeler, M.F. and Mear, M.E. (2005), “A numerical technique for simulating nonplanar evolution of hydraulic fractures”, *SPE Annual Technical Conference and Exhibition*, TX.

Samimi, S. and Pak, A. (2016), “A fully coupled element-free Galerkin model for hydro-mechanical analysis of advancement of fluid-driven fractures in porous media”, International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 40 No. 16, pp. 2178-2206.

Secchi, S. and Schrefler, B.A. (2012), “A method for 3-D hydraulic fracturing simulation”, International Journal of Fracture, Vol. 178 Nos 1/2, pp. 245-258.

Secchi, S., Simoni, L. and A. Schrefler, B. (2007), “Mesh adaptation and transfer schemes for discrete fracture propagation in porous materials”, International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 31 No. 2, pp. 331-345.

Settari, A. and Cleary, M.P. (1986), “Development and testing of a pseudo-three-dimensional model of hydraulic fracture geometry”, SPE Production Engineering, Vol. 1 No. 6.

Shen, Y. (2014), “A variational inequality formulation to incorporate the fluid lag in fluid-driven fracture propagation”, Computer Methods in Applied Mechanics and Engineering, Vol. 272, pp. 17-33.

Simoni, L., Secchi, S. and Schrefler, B.A. (2008), “Numerical difficulties and computational procedures for thermo-hydro-mechanical coupled problems of saturated porous media”, Computational Mechanics, Vol. 43 No. 1, pp. 179-189.

Spence, D.A. and Sharp, P. (1985), “Self-similar solutions for elastohydrodynamic cavity flow”, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, Vol. 400 No. 1819, pp. 289-313.

Walters, M.C., Paulino, G.H. and Dodds, R.H. (2005), “Interaction integral procedures for 3-D curved cracks including surface tractions”, Engineering Fracture Mechanics, Vol. 72 No. 11, pp. 1635-1663.

Wangen, M. (2013), “Finite element modeling of hydraulic fracturing in 3D”, Computational Geosciences, Vol. 17 No. 4, pp. 647-659.

Xu, G. (1999), “A variational boundary integral method for the analysis of three-dimensional cracks of arbitrary geometry in anisotropic elastic solids”, Journal of Applied Mechanics, Vol. 67 No. 2, pp. 403-408.

Yamamoto, K., Shimamoto, T. and Sukemura, S. (2004), “Multiple fracture propagation model for a three-dimensional hydraulic fracturing simulator”, International Journal of Geomechanics, Vol. 4 No. 1, pp. 46-57.

Zhang, X., Jeffrey, R.G., Bunger, A.P. and Thiercelin, M. (2011), “Initiation and growth of a hydraulic fracture from a circular wellbore”, International Journal of Rock Mechanics and Mining Sciences, Vol. 48 No. 6, pp. 984-995.

Zhang, X., Jeffrey, R.G. and Thiercelin, M. (2009), “Mechanics of fluid-driven fracture growth in naturally fractured reservoirs with simple network geometries”, Journal of Geophysical Research, Vol. 114 No. B12.

Zhang, X. and Jeffrey, R.G. (2006), “The role of friction and secondary flaws on deflection and re-initiation of hydraulic fractures at orthogonal pre-existing fractures”, Geophysical Journal International, Vol. 166 No. 3, pp. 1454-1465.

Zheltov, A.K. (1955), “Formation of vertical fractures by means of highly viscous liquid”, *4th World Petroleum Congress*, Rome.