Hydraulic Fracture



A Hydraulic fractures (HF) are a class of brittle fractures that propagate in pre-stressed solid media due to the injection of a viscous fluid. These fractures occur naturally when pressurized magma from deep underground chambers form vertical intrusions driven by buoyancy forces, which can in turn form geological structures such as dykes and sills. In the oil and gas industry, HF are deliberately created in reservoirs to enhance the recovery of hydrocarbons by the creation of permeable pathways . The application of HF in geotechnical engineering is growing. For example, in the mining industry they have recently been used to weaken the rock and enhance the so-called block-caving process. Similarly, the extraction of geothermal energy equires the creation of new fractures to increase the existing fracture networks. Hydraulic fractures have also been used for waste disposal and are likely to play a key role in the gravitational trapping of CO2 in deep, low permeability, ocean sediments. The propagation of hydraulic fractures into undesirable locations can have severe safety consequences in the mining industry, can cause considerable loss of hydrocarbons and environmental damage in the oil industry. Likewise, the perforation of the caprock by HF can reverse the costly capture process involved in CO2 sequestration. It is therefore of considerable importance to have accurate models in order to be able to predict the advance of HF to achieve an effective design of the engineering parameters in the injection process.

Mathematical models of HF involve a degenerate system of hyper-singular integro-partial differential equations defined on a domain with a moving boundary. The integral operator may be classified as a fractional Laplacian, which when combined with a free boundary problem forms is a class of problems that has received considerable attention in the theoretical PDE community (see for example the following link to a BIRS Workshop dedicated to this topic). The strong degeneracy of the lubrication equation yields solutions with singular pressure fields close to the tips or near "pinch points," which can occur within the fracture due to local changes in the geological structure. This class of problem has been shown to exhibit a rich multi-scale structure in which the fracture can propagate in a number of different modes each determined by the dominant physical process active at the tip of the fracture. The resolution of these problems requires a unique blend of numerical and analytic techniques to capture the free boundary of the evolving fracture as well as the complex coupling between the fluid which drives the fracture and the ambient elastic medium in which the fracture is being formed. The very different time scales on which the fluid flow and elastic response occur, poses considerable interesting challenges for numerical simulation. My interest in this research is the use asymptotic analysis to identify multi-scale singular solutions and to develop numerical tools to accommodate and exploit these asymptotic solutions.

 


29. An enhanced pseudo-3D model for hydraulic fracturing accounting for viscous height growth, non-local elasticity

by E.V. Dontsov and A. P. Peirce

Engineering Fracture mechanics, 142, p 116-139, 2015.

Comparison between the results from the classical P3D model, the enhanced P3D model (EP3D), and those of the fully planar simulator ILSA. For the zero toughness case KIc=0 the values are shown in blue and those of the experiment by Jeffrey and Bunger (2007) are depicted by the blue circles. The results shown in red and orange correspond to the toughness KIc=0.94 MPa.m^(1/2).

Abstract: The goal of this paper is to develop a more accurate pseudo-3D model for hydraulic fracturing. The primary weaknesses of the classical pseudo-3D model are: (1) its inability to capture viscous height growth and (2) its failure to include lateral fracture toughness.

These flaws are addressed respectively by: (1) introducing an apparent fracture toughness
in the vertical direction and (2) using an approximate non-local elasticity operator. To evaluate the accuracy and the level of improvement of the model we have developed, the results are compared to the predictions calculated using a recently developed fully planar hydraulic fracturing simulator.

 

28. Proppant transport in hydraulic fracturing: Crack tip screen-out in KGD and P3D models

by E.V. Dontsov and A. P. Peirce

Int. J. Solids & Structures, 63, p 206-218, 2015.

Fracture footprints with the proppant concentration indicated by colour calculated for the reference parameters at t=800 s (top left), t=1200 s (top right), t=1700 s (centre left), and t=3000 s (centre right). The case with t=3000 s corresponds to G_s=0.45. The bottom pictures show the pressure at the inlet, the length of the fracture, and the height of the fracture at the inlet versus time.

Fracture footprints and the proppant concentration indicated by colour calculated for the reference parameters and t=4000 s, except C'=0 (top left), C'=0 and mu'=0.24 Pa s (top right), C'=0 and Delta sigma=1.5x10^6 Pa (bottom left) and C'=0 and H=35 m (bottom right). The gravitational settling parameters are G_s=0.67 for both top pictures (notable settling), G_s=5.15 for the bottom left picture (significant settling), and G_s=0.17 for the bottom right picture (almost no settling).
Abstract: The aim of this study is to develop a model for proppant transport in hydraulic
fractures capable of capturing both gravitational settling and tip screen-out effects, while prohibiting the particles from reaching the crack tips by imposing a width restriction based on the particle size. First, the equations that govern the propagation of hydraulic fractures and the proppant transport inside them are formulated. They are based on the solution for the steady flow of a viscous fluid, mixed with spherical particles, in a channel, which is obtained assuming an empirical constitutive model. This proppant transport model is applied to two fracture geometries -- Khristianovich-Zheltov-Geertsma-De Klerk (KGD) and pseudo-3D (P3D). Numerical simulations show that the proposed method makes it possible to capture proppant plug formation and growth, as well as the gravitational settling for both geometries. A dimensionless parameter, whose magnitude reflects the intensity of the settling, is introduced for the P3D fracture.

 

27. Slurry flow, gravitational settling and a proppant transport model for hydraulic fractures,

by E.V. Dontsov and A. P. Peirce

J. Fluid Mech. V 760, pp 567-590, 2014.

Schematic view of a hydraulic fracture and the slurry flow with gravitational settling inside it

Normalized particle velocity profiles and particle volume fractions vs y/w

Abstract: The goal of this study is to analyse the steady flow of a Newtonian fluid mixed with spherical particles in a channel for the purpose of modelling proppant transport with gravitational settling in hydraulic fractures. The developments are based on a continuum constitutive model for a slurry, which is approximated by an empirical formula. It is shown that the problem under consideration features a two-dimensional flow and a boundary layer, which effectively introduces slip at the boundary and allows us to describe a transition from Poiseuille flow to Darcy’s law for high proppant concentrations. The expressions for both the outer (i.e. outside the boundary layer) and inner (i.e. within the boundary layer) solutions are obtained in terms of the particle concentration, particle velocity and fluid velocity. Unfortunately, these solutions require the numerical solution of an integral equation, and, as a result, the development of a proppant transport model for hydraulic fracturing based on these results is not practicable. To reduce the complexity of the problem, an approximate solution is introduced. To validate the use of this approximation, the error is estimated for different regimes of flow. The approximate solution is then used to calculate the expressions for the slurry flux and the proppant flux, which are the basis for a model that can be used to account for proppant transport with gravitational settling in a fully coupled hydraulic fracturing simulator.

 

26. Enrichment strategies and convergence properties of the XFEM for hydraulic fracture problems,

by E. Gordeliy and A. P. Peirce,

Comp. Meth. in Appl. Mech. and Eng., 283, p474-502, 2015.

 

 

 

 

 

 

 

 

Wronskian of the new enrichment basis functions

 

Abstract: In two recent papers (Gordeliy and Peirce, 2013; Gordeliy and Peirce, 2013) investigating the use of Extended Finite Element Method (XFEM) for modeling hydraulic fractures (HF), two classes of boundary value problem and two distinct enrichment types were identified as being essential components in constructing successful XFEM HF algorithms. In this paper we explore the accuracy and convergence properties of these boundary value formulations and enrichment strategies. In addition, we derive a novel set of crack-tip enrichment functions that enable the XFEM to model HF with the full range of power law r^{\lambda} behavior of the displacement field and the corresponding r^{\lambda-1} singularity in the stress field, for \frac{1}{2}\le\lambda<1. This novel crack-tip enrichment enables the XFEM to achieve the optimal convergence rate, which is not achieved by existing enrichment functions used for this range of power law. The two XFEM boundary value problem classes are as follows: i) a Neumann to Dirichlet map in which the pressure applied to the crack faces is the specified boundary condition and the XFEM is used to solve for the corresponding crack width (P->W); and ii) a mixed hybrid formulation of the XFEM that makes it possible to incorporate the singular behavior of the crack width in the fracture tip and uses a pressure boundary condition away from it (P&W). The two enrichment schemes considered are: i) the XFEM-t scheme with full singular crack-tip enrichment and ii) a simpler, more efficient, XFEM-s scheme in which the singular tip behavior is only imposed in a weak sense. If enrichment is applied to all the nodes of tip-enriched elements, then the resulting XFEM stiffness matrix is singular due to a linear dependence among the set of enrichment shape functions, which is a situation that also holds for the classic set of square-root enrichment functions. For the novel set of enrichment functions we show how to remove this rank deficiency by eliminating those enrichment shape functions associated with the null space of the stiffness matrix. Numerical experiments indicate that the XFEM-t scheme, with the new tip enrichment, achieves the optimal O(h^{2}) convergence rate we expect of the underlying piece-wise linear FEM discretization, which is superior to the enrichment functions currently available in the literature for \frac{1}{2}<\lambda<1. The XFEM-s scheme, with only signum enrichment to represent the crack geometry, achieves an O(h) convergence rate. It is also demonstrated that the standard P->W formulation, based on the variational principle of minimum potential energy, is not suitable for modeling hydraulic fractures in which the fluid and the fracture fronts coalesce, while the mixed hybrid P\&W formulation based on the Hellinger-Reissner variational principle does not have this disadvantage.

 

25. Modeling multi-scale processes in hydraulic fracture propagation using the implicit level set algorithm,

by A. Peirce

Comp. Meth. in Appl. Mech and Eng., 283, p881-908, 2015.

 

ILSA-MK fracture widths (red) for various toughnesses compared to a reference solution (dots) and the M-vertex solution (dashed magenta) and K-vertex solution (dased blue)

ILSA-MK pressures (red) for various toughnesses compared to a reference solution (dots) and the M-vertex solution (dashed magenta) and K-vertex solution (dased blue)

Fracture footprints for a stress jump problem at t=604 s. The M-vertex footprint is shown in red and the footprint with a scaled toughness K' = 0.75 MPa.m^(1/2) is shown in blue. The colour represents the length scale used to locate the free boudary.

Fracture footprints for a stress jump problem at t=604 s. The M-vertex footprint is shown in red and the footprint with a scaled toughness K' = 1.5 MPa.m^(1/2) is shown in blue. The colour represents the length scale used to locate the free boudary.

Horizontal cross-section of the ILSA width field (solid) and the P3D width (dashed).

Vertical cross-section of the ILSA width field (solid) and the P3D width (dashed), and the theoretical width assuming a constant pressure (red)..


Abstract: In this paper we describe an implicit level set algorithm (ILSA) suitable for modeling multi-scale behavior in planar hydraulic fractures propagating in three dimensional elastic media. This multi-scale behavior is typically encountered when multiple physical processes compete to determine the location of the fracture free boundary. Instead of having to match the mesh size to the finest active length scale, or having to re-mesh as the dominant length scales change in space and time, the novel ILSA scheme is able to represent the required multi-scale behavior on a relatively
coarse rectangular mesh. This is achieved by using the local front velocity to construct, for each point of a set of control points, a mapping that adaptively identifies the dominant length scale at which the appropriate multi-scale universal asymptotic solution needs to be sampled. Finer-scale behavior is captured in a weak sense by integrating the universal asymptotic solution for the fracture width over partially filled tip elements and using these integrals to set the average values of the widths in all tip elements. The ILSA solution shows good agreement with a multi-scale reference solution comprising a radial solution that transitions from viscosity to toughness dominated ropagation regimes. The ILSA scheme is also used to model blade-like hydraulic fractures that break through stress barriers located symmetrically with respect to the injection point. For the zero toughness case, the ILSA solutionshows close agreement to experimental results. The multi-scale ILSA
scheme is also used to provide results when the material toughness K_{Ic} is non-zero. In this case different parts of the fracture-free-boundary can be propagating in different regimes. It is hoped that the multi-scale ILSA solutions presented here will form a set of reference results that can be used to benchmark simulators that use a propagation criterion based on only one dissipative process (either toughness or viscosity). The multi-scale ILSA solutions at larger times (for which plane strain conditions develop in vertical cross sections) are compared with and show close
agreement to plane strain exact solutions for height-growth and the fracture width in vertical cross sections. This comparison provides some measure of the accuracy of the multi-scale ILSA scheme. The multi-scale ILSA solutions are also used to identify the regimes of applicability of pseudo 3D (P3D) approximate solutions. These ILSA
solutions can also be used to design improved P3D models.

24. A New Technique for Proppant Schedule Design,

by E. V. Dontsov and A. P. Peirce

Hydraulic Fracturing Journal, V1, No. 3, 2014.

Proppant Schedule and placment

 

 

 

Abstract: This study introduces a novel methodology for the design of the proppant pumping schedule for a hydraulic fracture, in which the final proppant distribution along the crack is prescribed. The method is based on the assumption that the particles have relatively small impact on the fracture propagation, unless they reach the tip region. This makes it possible to relate the proppant velocity to the clear fluid velocity inside the fracture, which is calculated assuming no proppant. Having the history of the clear fluid velocity distribution, the prospective proppant motion can be computed. Then, volume balance is used to relate the final concentration at some point inside the fracture to the corresponding input concentration at a speci c time instant, which helps to avoid solving an inverse problem. One exceptional feature of the approach lies in the fact that it is applicable to multiple fracture geometries and can be implemented using various hydraulic fracturing simulators. To verify the technique, two fracture geometries are considered - Khristianovich-Zheltov-Geertsma-De Klerk (KGD) and pseudo-3D
(P3D). It is shown that the developed approach is capable of properly estimating the pumping schedule for both geometries. In particular, the proppant placement along the fracture at the end of the pumping period, calculated according to the adopted proppant transport model, shows close agreement with the design distribution. The comparison with Nolte's scheduling scheme shows that the latter is not always accurate, and cannot capture the essential differences between the schedules for the fracture geometries considered.

23. On the moving boundary conditions for a Hydraulic Fracture
by E. Detournay and A. Peirce
International Journal of Engineering Science, 84, p 147–155, 2014.

Evolving Fluid and Fracture Fronts

Abstract: This paper re-examines the boundary conditions at the moving front of a hydraulic fracture when the fluid front has coalesced with the crack edge. This practically important particular case is treated as the zero fluid lag limit of the general case when the two fronts are distinct. The limiting process shows what becomes of the two boundary conditions on the fluid front: a pressure condition and a Stefan condition equating the front velocity to the average fluid velocity, when the lag vanishes. On the one hand, the pressure condition disappears as the net pressure (the difference between the fluid pressure and the magnitude of the far-field stress normal to the fracture) becomes singular. On the other hand, the Stefan condition transforms into a zero flux boundary condition at the front, and, as a consequence, the velocity of the coalesced front does not appear explicitly in the boundary conditions. However, the front velocity can still be extracted from the near-tip aperture field by a nonlinear asymptotic analysis. The paper concludes with a description of an algorithm to propagate the combined front, which explicitly uses the known multiscale asymptotics of the fracture aperture.

22..Numerical Simulation of Simultaneous Growth of Multiple Interacting Hydraulic Fractures from Horizontal Wells,


by A.P. Bunger and A. P. Peirce, submitted to the ASCE Shale Energy Conference, Pittsburgh, 21-23 July, 2014.

5 hydraulic fractures propagating in the toughness regime with contours corresponding to opening.

Evolution of the radius of 5 hydraulic fractures propagating in the toughness regime normalized by the spacing between them.

5 hydraulic fractures propagating in the viscous regime with contours corresponding to opening

Evolution of the radius of 5 hydraulic fractures propagating in the viscus regime normalized by the spacing between them.

Height limited hydraulic fractures with spacing larger than the height 7 hydraulic fractures with contours corresponding to opening.

Evolution of the maximal radius of the fractures normalized by the spacing between them.


Abstract: The technique of multistage hydraulic fracturing from horizontal wells is universally credited with enabling the economical production of hydrocarbon resources from shale formations. The method almost always entails the injection of fluid through the wellbore with the potential to create hydraulic fractures from multiple reservoir entry points, typically clusters of wellbore perforations, that are spaced out along the wellbore within a section that is colloquially referred to as a “stage”. Arguably the most basic question about this situation is how many perforation clusters within a given fracturing stage can be expected to produce growing hydraulic fractures. This paper presents a numerical investigation of this issue that employs a newly-developed, fully coupled parallel planar 3D hydraulic
fracturing simulator that features: implicit time stepping, an implicit level set scheme to locate the propagating hydraulic fracture fronts that respond to their regimes of propagation and enables highly accurate simulations using a very coarse mesh, and the capability to dynamically partition the fluid among multiple, simultaneously growing hydraulic fractures in parallel, overlapping planes. Our results demonstrate the dependence of the energetically preferred number of growing hydraulic fractures on the length of the isolated zone, the height of the reservoir, and the relative importance of the fluid viscosity. In particular, we show that reservoirs with effective height containment and injection strategies that ensure substantial viscous dissipation will promote growth of multiple simultaneous hydraulic fractures rather than localization to just one or two dominant fractures.

21. Interference Fracturing: Non-Uniform Distributions of Perforation Clusters that Promote Simultaneous Growth of Multiple Hydraulic Fractures

(link to Posting on UBC Circle Repository)


by A. P. Peirce and A.P. Bunger, Submitted to the SPE Journal on the 17 October, 2013.

Numerical model of a 4-fracture stage of length Z awith a stress jump with a pay zone H.

Geometriy for parallel-planar model

Uniform 3-fracture array with stress shadowing

Uniform 5-fracture array with stress shadowing

Non-uniform 5-fracture array in which interference fractring leads to different dynamics that avoids stress shadowing

Quantities sampled at the well-bore for fracures 1 (red), 2 (green), and 3 (blue) (numbering from the end of the array) (a) largest fracture dimension/H; (b) well-bore fluxes qk; (c) well-bore widths; (d) well-bore pressures.

Abstract: One of the important hurdles in horizontal well stimulation is the generation of hydraulic fractures (HFs) from all perforation clusters within a given stage despite the challenges posed by stress shadowing and reservoir variability. In this paper we use a newly-developed, fully coupled, parallel-planar 3D HF model to investigate the potential to minimize the negative impact of stress shadowing and thereby to promote more uniform fracture growth across an array of HFs by adjusting the location of the perforation clusters. In this model the HFs are assumed to evolve in an array of parallel planes with full 3D stress coupling while the constant fluid influx into the wellbore is dynamically partitioned to each fracture so that the wellbore pressure is the same throughout the array. The model confirms the phenomenon of inner fracture suppression due to stress shadowing when the perforation clusters are uniformly distributed. Indeed, the localization of the fracture growth to the outer fractures is so dominant that the total fractured area generated by uniform arrays is largely independent of the number of perforation clusters. However, numerical experiments indicate that certain non-uniform cluster spacings promote a profound improvement in the even development of fracture growth. Identifying this effect relies on this new model’s ability to capture the full hydro-dynamical coupling between the simultaneously evolving HFs in their transition from radial to PKN-like geometries.

20. Implicit level set schemes for modeling hydraulic fractures using the XFEM


by E. Gordeliy and A. P. Peirce
Comp. Meth. in Appl. Mech. and Eng., 266, p125-143, 2013.

(a) Sign enrichment model of a crack and (b) crack model with sign and tip enrichment

Mesh for model of a crack approaching a bi-material interface.

Accuracy coparison between the XFEM solution (blue circles), high reolution BEM solution (black triangles), and Erdogan semi-analytic solution (black diamonds).

Three width and pressure field pairs sampled at three disticnt times for a hydraulic fracture passing through a bi-material interface.

 

Abstract: We describe two novel XFEM schemes for modeling fluid driven fractures both of which exploit an implicit level set algorithm (ILSA) for locating the singular free boundary that occurs when the fluid and fracture fronts coalesce. Both schemes use the mixed P&W XFEM formulation developed in (Gordeliy and Peirce, 2013) to incorporate the singular asymptotic solution in the fracture tips. The proposed level set strategy also exploits the asymptotic solution to provide a robust procedure to locate the free boundary, which is not restricted to symmetric growth of the fracture geometry or to a particular mode of propagation. The versatility of the ILSA-XFEM scheme is demonstrated by sampling different asymptotic behaviors along the so-called MK edge of parameter space (Detournay, 2004) by making use of a universal asymptote (Garagash, 1998; Garagash and Detournay, 2000). The two ILSA-XFEM schemes differ in the enrichment strategies that they use to represent the fracture tips: a scheme with full tip enrichment and a simpler, more efficient, scheme in which the tip asymptotic behavior is only imposed in a weak sense. Numerical experiments indicate that the XFEM-t scheme, with full tip enrichment, achieves an O(h^{2}) asymptotic convergence rate, while the XFEM-s scheme, with only signum enrichment to represent the crack geometry, achieves an O(h) asymptotic convergence rate.

19. .The impact of the Near-Tip Logic on the Accuracy and Convergence Rate of Hydraulic Fracture Simulators Compared to Reference solutions.


by Lecampion B., Peirce A.P, Detournay E., Zhang X., Chen Z., Bunger A.P., Detournay C., Napier J., Abbas S., Garagash D., Cundall, P.,

In: Effective and Sustainable Hydraulic Fracturing, AP Bunger, J McLennan and R Jeffrey (eds.), ISBN 978-953-51-1137-5, (Intech), Chapter 43, p855-873, 2013.


Evolution of the relative error in the fracture radius as a function of the mesh size for the different simulators - radial fracture in the viscosity dominated regime (i.e. zero toughness / early time, t < 1). All simulations displayed here are for t < 10E-6.

We benchmark a series of simulators against available reference solutions for propagating plane-strain and radial hydraulic fractures. In particular, we focus on the accuracy and convergence of the numerical solutions in the important practical case of viscosity dominated
propagation. The simulators are based on different propagation criteria: linear elastic fracture mechanics (LEFM), cohesive zone models/tensile strength criteria, and algorithms accounting for the multi-scale nature of hydraulic fracture propagation in the near-tip region. All the simulators tested here are able to capture the analytical solutions of the different
configurations tested, but at vastly different computational costs. Algorithms based on the classical LEFM propagation condition require a fine mesh in order to capture viscosity dominated hydraulic fracture evolution. Cohesive zone models, which model the fracture process zone, require even finer meshes to obtain the same accuracy. By contrast, when the algorithms use the appropriate multi-scale hydraulic fracture asymptote in the near-tip region, the exact solution can be matched accurately with a very coarse mesh. The different analytical reference solutions used in this paper provide a crucial series of benchmark tests that any successful hydraulic fracturing simulator should pass.

18. Coupling schemes for modeling hydraulic fracture propagation using the XFEM.


by E. Gordeliy and A. P. Peirce
Comp. Meth. in Appl. Mech. and Eng., 266, p125-143, 2013.


Top: The crack path obtained using a DD cde OribiC, XFEM-t, and XFEM
Bottom: A fragment of the crack path within the FEM mesh sampled within the rectangular box shown in the top figure..

Abstract: We describe coupled algorithms that use the Extended Finite Element Method (XFEM) to solve the elastic crack component of the elasto-hydrodynamic equations that govern the propagation of hydraulic fractures in an elastic medium. With appropriate enrichment, the XFEM resolves the Neumann to Dirichlet (ND) map for crack problems with O(h^{2}) accuracy and the Dirichlet to Neumann (DN) map with O(h) accuracy. For hydraulic fracture problems with a lag separating the fluid front from the fracture front, we demonstrate that the finite pressure field makes it possible to use a scheme based on the O(h^{2}) XFEM solution to the ND map. To treat problems in which there is a coalescence of the fluid and fracture fronts, resulting in singular tip pressures, we developed a novel mixed algorithm that combines the tip width asymptotic solution with the O(h^{2}) XFEM solution of the ND map away from the tips. Enrichment basis functions required for these singular pressure fields correspond to width power law indices \lambda >=\frac{1}{2), which are different from the index \lambda=\frac{1}{2}
of linear elastic fracture mechanics. The solutions obtained from the new coupled XFEM schemes agree extremely well with those of published reference solutions.


17. An Integrated Extended Kalman Filter Implicit Level Set Algorithm for monitoring Planar Hydraulic Fracture.
by A.P. Peirce and F. Rochinha
Inverse Problems, 28,(2012) 015009, (22 pp).

Earlier time synthetic (red) and identified (black) fracture footprints for a hydraulic fracture evolving in a stress gradient. Synthetic tiltmeter measurements were sampled from a near-field planar array 10 units away from the fracture plane with noise levels of 1%-5%.

 

Later time synthetic (red) and identified (black) fracture footprints for a hydraulic fracture evolving in a stress gradient. Synthetic tiltmeter measurements were sampled from a near-field planar array 10 units away from the fracture plane with noise levels of 1%-5%.

 

Movie of the fracture footprints for a fracture evolving in a Stress Gradient field identified from a Near-Field Tiltmeter array with 6% Gaussian Noise. The Tiltmeters lie in a planar grid (indicated on the movie by black dots) pararllel to the fracture plane and 10 units away from the fracture plane. The synthetic footprints are shown in red and the identified footprints are shown in blue.

Movie of the fracture footprints for a fracture evolving in a Stress Gradient field identified from a Far-Field Tiltmeter array with 4% Gaussian Noise. The Tiltmeters lie in a planar grid pararllel to the fracture plane and 50 units away from the fracture plane. The synthetic footprints are shown in red and the identified footprints are shown in blue.

Synthetic (red) and identified (black) fracture footprints for a symmetric stress jump

Vertical cross-section of synthetic (red) and identified (black) fracture widths for a symmetric stress jump
Movie of the fracture footprints for a fracture evolving between two symmetric stress jumps. Inversions are obtained from measurments from a Near-Field planar tiltmeter array pararllel to the fracture plane and 10 units away. No noise was added to the measurments. The synthetic footprints are shown in red and the identified footprints are shown in blue. Movie of the fracture footprints for a fracture evolving between two symmetric stress jumps. Inversions are obtained from measurements from a Near-Field planar tiltmeter array pararllel to the fracture plane and 10 units away. A noise level of 2% was added to the measurments. The synthetic footprints are shown in red and the identified footprints are shown in blue.

Synthetic (red) and identified (black) fracture footprints for an Asymmetric stress jump

Vertical cross-section of synthetic (red) and identified (black) fracture widths for an Asymmetric stress jump

Movie of the fracture footprints for a fracture initiated in layer between a higher confining stress region (above) and a lower confining stress region (below).. Inversions are obtained from measurments from a Far-Field planar tiltmeter array pararllel to the fracture plane and 50 units away. No noise was added to the measurments. The synthetic footprints are shown in red and the identified footprints are shown in blue. Movie of the fracture footprints for a fracture initiated in layer between a higher confining stress region (above) and a lower confining stress region (below).. Inversions are obtained from measurments from a Far-Field planar tiltmeter array pararllel to the fracture plane and 50 units away. A noise level of 4% was added to thes measurments. The synthetic footprints are shown in red and the identified footprints are shown in blue.

 

Abstract: We describe a novel approach to the inversion of elasto-static tiltmeter measurements to monitor a planar hydraulic fracture propagating within a three dimensional elastic medium. The technique involves integrating the Extended Kalman Filter (EKF) for state prediction and updates using tiltmeter measurement time-series combined with a novel Implicit Level Set Algorithm (ILSA) for the
solution of the coupled elasto-hydrodynamic equations in order tolocate the unknown fracture free boundary. By means of a scaling argument a strategy is derived to tune the parameters in the algorithm to enable measurement information to compensate for significant un-modeled dynamics. A sequence of numerical experiments are considered in which significant changes are introduced to the fracture geometry by altering the confining geological stress field. Even though there is no confining stress field in the dynamic model used by the new EKF-ILSA scheme, it is able to use synthetic measurements from tiltmeters in order to arrive at remarkably accurate predictions of the fracture opening and the fracture footprint. These experiments also explore the robustness of the algorithm to noise and to placement of tiltmeter arrays operating in the near-field and far-field regimes. In these experiments the
appropriate parameter choices and strategies to improve the robustness of the algorithm to significant measurement noise are explored.

16. A Hermite Cubic Collocation Scheme for Plane Strain Hydraulic Fractures
by A.P. Peirce
Computer Methods in Applied Mechanics and Engineering , 199, Issues 29-32, p1949-1962, 2010.

Hermite Basis Functions and the stresses they generate

 

Abstract: We describe a novel cubic Hermite collocation scheme for the solution of the coupled integro-partial differential equations governing the propagation of a hydraulic fracture in a state of plane strain. Special blended cubic Hermite-power law basis functions, with arbitrary index 0<a<1, are developed to treat the singular behavior of the solution that typically occurs at the tips of a hydraulic fracture. The implementation of blended infinite elements to model semi-infinite crack problems is also described. Explicit formulae for the integrated kernels associated with the cubic Hermite and blended basis functions are provided. The cubic Hermite collocation algorithm is used to solve a number of different test problems with two distinct propagation regimes and the results are shown to converge to published similarity and asymptotic solutions. The convergence rate of the cubic Hermite scheme is determined by the order of accuracy of the tip asymptotic expansion as well as the O(h^4) error due to the Hermite cubic interpolation. The errors due to these two approximations need to be matched in order to achieve optimal convergence. Backward Euler time-stepping yields a robust algorithm that, along with geometric increments in the time-step, can be used to explore the transition between propagation regimes over many orders of magnitude in time.

15. Monitoring Hydraulic Fractures: State Estimation using an Extended Kalman Filter
by F. Rochinha and A. Peirce
Inverse Problems, 26, 2, 025009, 17p, 2010.

Abstract: There is considerable interest in using remote elastostatic deformations to identify the evolving geometry of underground fractures that are forced to propagate by the injection of a high pressure viscous fluid. These so-called hydraulic fractures are used to increase the permeability in oil and gas reservoirs as well as to pre-fracture ore-bodies in preparation for a mining process known as block-caving. The undesirable intrusion of these hydraulic fractures into environmentally sensitive areas or into regions in mines which might pose safety hazards, has stimulated the search for techniques to enable the evolving hydraulic fracture geometries to be monitored. Previous approaches to this problem have involved the inversion of the elastostatic data at isolated time steps in the time series provided by tiltmeter measurements of the displacement gradient field at selected points in the elastic medium. At each time step, parameters in simple static models of the fracture (e.g., a single displacement discontinuity or low order moments of the fracture opening) are identified. The approach adopted in this paper is not to regard the sequence of sampled elastostatic data as independent, but rather to treat the data as linked by the coupled elastic-lubrication equations that govern the propagation of the evolving hydraulic fracture. We combine the Extended Kalman Filter (EKF) with features of a recently developed implicit numerical scheme to solve the coupled free boundary problem in order to form a novel algorithm to identify the evolving fracture geometry. Numerical experiments demonstrate that, despite excluding significant physical processes in the forward numerical model that substantially alter the geometry of the evolving hydraulic fracture, the EKF-numerical algorithm is able to compensate for the un-modeled processes by using the information fed back from tiltmeter data. Indeed, the proposed algorithm is able to provide reasonably faithful estimates of the fracture geometry, which are shown to converge to the actual hydraulic fracture geometry as the number of tiltmeters is increased. The location of tiltmeters can affect the resolution of the technique, which opens the possibility of using the algorithm to design the deployment of tiltmeters to optimize the resolution in regions of particular interest.

14. An Analysis of Classical Pseudo-3D Model for Hydraulic Fracture with Equilibrium Height Growth across Stress Barriers
by J. I. Adachi, E. Detournay, and A. P. Peirce,
Int. J. of Rock Mech. and Min. Sci., 47, p 625-639, 2010. (MATLAB Code: ColP3D.m and .pdf file from the Circle UBC Digital Repository).

Movies

The MATLAB code that implements the 4 th Order Collocation Scheme is given below.

MATLAB Code

or

colP3D posted on Circle UBC Digital Repository

Abstract: This paper deals with the so-called “pseudo three-dimensional” (P3D) model for a hydraulic fracture with equilibrium height growth across two symmetric stress barriers. The key simplifying assumptions behind the P3D model, which was originally introduced by A. Settari and M. P. Cleary (Development and Testing of a Pseudo-Three-Dimensional Model of Hydraulic Fracture Geometry (P3DH), SPE 10505, Proc. 6th SPE Symposium on Reservoir Simulation, pp. 185-214, 1982) and by I. D. Palmer and H. B. Carrol (Numerical Solution for Height and Elongated Hydraulic Fractures, SPE/DOE 11627, Proc. SPE/DOE Symposium on Low Permeability, pp. 249-878, 1983), are that (i) each cross section perpendicular to the main
propagation direction is in a condition of plane-strain, and (ii) the local fracture height is determined by a balance between the effect of the stress jump across the barriers and that of the rock toughness. Furthermore, in the equilibrium height growth P3D models, the pressure is assumed uniform in each vertical cross-section. We revisit this particular model by first formulating the non-linear differential equations governing the evolution of the length, height, and aperture of the hydraulic fracture, in contrast to the numerical formulations adopted in many previous studies. Scaling of these equations shows that the solution depends, besides the dimensionless space and time coordinates, on only two numbers representing a scaled toughness and a scaled leak-off coefficient. Analysis of the governing equations enables us to determine explicitly the conditions under which breakthrough takes place (i.e., the onset of growth into the bounding layers), as well as the conditions of unstable height growth (i.e., the conditions of “runaway height” when the main assumptions of the equilibrium height model become invalid). The mathematical model is solved numerically using a novel implicit fourth order collocation scheme on a moving mesh, which makes explicit use of the fracture tip asymptotics. A complete listing of the MATLAB code for this algorithm is provided in an appendix and can be copied for experimentation. We then report the results of several numerical simulations conducted for different values of the dimensionless toughness and the dimensionless leak-off coefficients, as well as a comparison with closed-form small and large time similarity solutions that are valid under conditions where the fracture remains contained within the reservoir layer. Finally, we discuss possible extensions of the model, including layers of different elastic moduli and power law fluids.

13. Uniform asymptotic Green's functions for efficient modeling of cracks in elastic layers with relative shear deformation controlled by linear springs
by A.Peirce, H. Gu, and E. Siebrits
Int. J. Num. and Anal. Meth. in Geomechanics, 33, 285-308, 2009.

Abstract: We present a uniform asymptotic solution (UAS) for a displacement discontinuity (DD) that lies within the middle layer of a three-layer elastic medium in which relative shear deformation between parallel interfaces is controlled by linear springs. The DD is assumed to be normal to the two interfaces between the elastic media. Using the Fourier transform (FT) method we construct the leading term in the asymptotic expansion for the spectral coefficient functions for a DD in a three layer-spring medium. Although a closed form solution will require a solution in terms of an infinite series, we demonstrate how this UAS can be used to construct highly efficient and accurate solutions even in the case in which the DD actually touches the interface. We compare the results using the Green's function UAS solution for a crack crossing a soft interface with results obtained using a multi-layer boundary element method. We also present results from an implementation of the UAS Green's function approach in a pseudo 3D hydraulic fracturing simulator to analyze the effect of interface shear deformation on the fracture propagation process. These results are compared with field
measurements.

12. An Eulerian Moving Front Algorithm with Weak-Form Tip Asymptotics for Modeling Hydraulically driven Fractures
by A. Peirce and E. Detournay
Communications in Numerical Methods in Engineering, 25, 185-200, 2009.

Exact and numerical Width and Pressure solutions
Abstract: The coupled equations for a propagating hydraulic fracture exhibit a multi-scale structure for which resolution on all length scales becomes computationally prohibitive. Asymptotic analysis is able to identify the dominant physical process active at the computational length scale. This paper describes a novel algorithm which uses weak-form tip asymptotics on a rectangular Eulerian mesh to solve the problem of a propagating hydraulic fracture. The location of the fracture front is determined within each tip element by matching the volume associated with the known asymptotic solution to the flux of fluid into the given element. Even if the fracture front is curved, the algorithm is able to capture the solution on a relatively coarse rectangular mesh by implementing a weak-form of the tip asymptotic solution, based on averaging the volume over a tip element to determine the fracture widths at tip element centres. The fracture is divided into a ``channel'' region made up of elements that are filled with fluid and a ``tip'' region comprising partially filled elements. An iterative procedure is used to determine the fracture widths and fluid pressures in the channel region and the fracture font locations in the tip regions. The algorithm is tested by analyzing a hydraulic fracture propagating in a viscosity dominated regime within an impermeable homogeneous elastic material for which a similarity solution is available. The numerical results show close agreement to the exact solution even though a relatively coarse mesh is used.

11. An Implicit Level Set Method for Modeling Hydraulically Driven Fractures
by A. Peirce and E. Detournay
Computer Methods in Applied Mechanics and Engineering, 197, 2858-2885, 2008.

Movies

Radial Fracture

 

Comparisons with similarity solutions

Viscosity dominated fracture propagation

Width

Pressure

Toughness dominated fracture propagation

Width

Pressure

Abstract: We describe a novel implicit level set algorithm to locate the free boundary for a propagating hydraulic fracture. A number of characteristics of the governing equations for hydraulic fractures and their coupling present considerable challenges for numerical modeling, namely: the degenerate lubrication equation; the hypersingular elastic integral equation; the indeterminate form of the velocity of the unknown fracture front, which precludes the implementation of established front evolution strategies that require an explicit velocity field; and the computationally prohibitive cost of resolving all the length scales. An implicit algorithm is also necessary for the efficient solution of the stiff evolution equations that involve fully populated matrices associated with the coupled non-local elasticity and degenerate lubrication equations. The implicit level set algorithm that we propose exploits the local tip asymptotic behavior, applicable at the computational length scale, in order to locate the free boundary. Local inversion of this tip asymptotic relation yields the boundary values for the Eikonal equation, whose solution gives the fracture front location as well as the front velocity field. The efficacy of the algorithm is tested by comparing the level set solution to analytic solutions for hydraulic fractures propagating in a number of distinct regimes. The level set algorithm is shown to resolve the free boundary problem with first order accuracy. Further it captures the field variables, such as the fracture width, with the first order accuracy that is consistent with the piecewise constant discretization that is used.

10. Asymptotic Analysis of an Elasticity Equation for a Finger-Like Hydraulic Fracture
by J. Adachi and A. P. Peirce
Journal of Elasticity, 90, 43-69, 2008.


Abstract: We derive a novel integral equation relating the fluid pressure in a fingerlike hydraulic fracture to the fracture width. By means of an asymptotic analysis in the small height to length ratio limit we are able to establish the action of the integral operator for receiving points that lie within three distinct regions: (1) an outer expansion region in which the dimensionless pressure is shown to be equal to the dimensionless width plus a small correction term that involves the second derivative of the width, which accounts for the nonlocal effects of the integral operator. The leading order term in this expansion is the classic local elasticity equation in the PKN model that is widely used in the oil and gas industry; (2) an inner expansion region close to the fracture tip within which the action of the elastic integral operator is shown to be the same as that of a finite Hilbert transform associated with a state of plane strain. This result will enable pressure singularities and stress intensity factors to be incorporated into analytic models of these finger-like fractures in order to model the effect of material toughness; (3) an intermediate region within which the action of the Fredholm integral operator of the first kind is reduced to a second kind operator in which the integral term appears as a small perturbation which is associated with a convergent Neumann series. These results are important for deriving analytic models of finger-like hydraulic fractures that are consistent with linear elastic fracture mechanics.

9. Multipole Moment Decomposition for imaging hydraulic fractures from remote elastic data
by B. Lecampion and A. Peirce
Inverse Problems, 23, 1641-1658, 2007.

Abstract: Hydraulic fracturing involves the propagation of a fracture in brittle rock by the intrusion of a high pressure viscous fluid. There is considerable interest in identifying characteristics of these evolving underground fractures via the passive monitoring of remote elastostatic deformations. In this paper, we present a far-field multipole expansion procedure to identify the harmonic moments of the fracture. The harmonic moments are related to fundamental quantities such as fracture volume and fracture asymmetries. We illustrate the efficacy of the multipole moment expansion technique by inverting synthetic displacement data from a hydraulic fracture simulator in order to identify the harmonic moments up to second order. These results are compared to those obtained by identifying the parameters of a dislocation model with a prescribed geometry—a procedure which is commonly used for such problems. The multipole moment expansion technique has the following features: it provides significantly more accurate fracture volume information; it provides accurate estimates of first-order moments that can be used to identify asymmetric fractures; it is possible to adapt the truncation process to optimize the information content of a given set of measurements; it can, in some cases, provide estimates of the higher order moments which can be used to determine geometric attributes of the fracture. Given this last feature, we explore the possibility of using up to second-order harmonic moments to identify the dimensions of a simple polygonal model of the fracture footprint. This procedure is tested by attempting to identify fracture footprints from synthesized hydraulic fracture data.

8. Computer Simulation of Hydraulic Fractures
by J. Adachi, E. Siebrits, A. Peirce, and J. Desroches
Int. J. of Rock Mech. and Min. Sci. & Geomech. Abstr., 44, 739-757, 2007.

Propppant concentration (left) before and after shut-in and the corresponding fracture widths (right)

Abstract: We provide a brief historical background of the development of hydraulic fracturing models for use in the petroleum and other industries. We discuss scaling laws and the propagation regimes that control the growth of hydraulic fractures from the laboratory to the field scale. We introduce the mathematical equations and boundary conditions that govern the hydraulic fracturing process, and discuss numerical implementation issues including: tracking of the fracture footprint, the control of the growth of the hydraulic fracture as a function of time, coupling of the equations, and time-stepping schemes. We demonstrate the complexity of hydraulic fracturing by means of an application example based on real data. Finally, we highlight some key areas of research that need to be addressed in order to improve current models.

7. An asymptotic framework for finite hydraulic fractures including leak-off
by S. L. Mitchell, R. Kuske, and A.P. Peirce
SIAM J. on Appl. Math., 67, Issue 2, 364-386, 2007.

Abstract: The dynamics of hydraulic fracture, described by a system of nonlinear integrodifferential equations, is studied through the development and application of a multiparameter singular perturbation analysis. We present a new single expansion framework which describes the interaction between several physical processes, namely viscosity, toughness, and leak-off. The problem has nonlocal and nonlinear effects which give a complex solution structure involving transitions on small scales near the tip of the fracture. Detailed solutions obtained in the crack tip region vary with the dominant physical processes. The parameters quantifying these processes can be identified from critical scaling relationships, which are then used to construct a smooth solution for the fracture depending on all three processes. Our work focuses on plane strain hydraulic fractures on long time scales, and this methodology shows promise for related models with additional time scales, fluid lag, or different geometries, such as radial (penny-shaped) fractures and the classical
Perkins–Kern–Nordgren (PKN) model.

6. An asymptotic framework for the analysis of hydraulic fractures: the impermeable case
by S. L. Mitchell, R. Kuske, and A.P. Peirce
ASME Journal of Applied Mechanics., 74, 365-372, 2007.

Abstract: This paper presents a novel asymptotic framework to obtain detailed solutions describing the propagation of hydraulic fractures in an elastic material. The problem consists of a system of nonlinear integro-differential equations and a free boundary problem. This combination of local and nonlocal effects leads to transitions on a small scale near the crack tip, which control the behavior across the whole fracture profile. These transitions depend upon the dominant physical process(es) and are identified by simultaneously scaling the associated parameters with the distance from the tip. A smooth analytic
solution incorporating several physical processes in the crucial tip region can be constructed using this new framework. In order to clarify the exposition of the new methodology, this paper is confined to considering the impermeable case in which only the two physical processes of viscous dissipation and structure energy release compete.

5. Localized Jacobain ILU Preconditioners for Hydraulic Fractures
by A. Peirce,
International Journal for Numerical Methods in Engineering", 65, 12, 1935-1946, 2006.
Spectra of Jacobain preconditioned by (3x3) and (5x5) localized Jacobians compared to diagonal preconditioning
Abstract: We discuss the properties of a class of sparse localized approximations to the Jacobian operator that arises in modelling the evolution of a hydraulically driven fracture in a multi-layered elastic medium. The governing equations involve a highly non-linear coupled system of integro-partial differential equations along with the fracture front free boundary problem. We demonstrate that an incomplete LU factorization of these localized Jacobians yields an efficient preconditioner for the fully populated, stiff, non-symmetric system of algebraic equations that need to be solved multiple times for every growth increment of the fracture. The performance characteristics of this class of preconditioners is explored via spectral analysis and numerical experiment.

4. A dual mesh Multigrid preconditioner for the efficient solution of hydraulically driven fracture problems
by A. P. Peirce and E. Siebrits
International Journal for Numerical Methods in Engineering, 63, 13, 1797-1823, 2005.
Grid coarsening to preserve "pinch points" at which width constraints w>wc are active
Abstract: We present a novel multigrid (MG) procedure for the efficient solution of the large non-symmetric system of algebraic equations used to model the evolution of a hydraulically driven fracture in a multi-layered elastic medium. The governing equations involve a highly non-linear coupled system of integro-partial differential equations along with the fracture front free boundary problem. The conditioning of the algebraic equations typically degrades as O(N^3). A number of characteristics of this problem present significant new challenges for designing an effective MG strategy. Large changes in the coefficients of the PDE are dealt with by taking the appropriate harmonic averages of the discrete coefficients. Coarse level Green’s functions for multiple elastic layers are constructed using a single dual mesh and superposition. Coarse grids that are sub-sets of the finest grid are used to treat mixed variable problems associated with ‘pinch points.’ Localized approximations to the Jacobian at each MG level are used to devise efficient Gauss–Seidel smoothers and preferential line iterations are used to eliminate grid anisotropy caused by large aspect ratio elements. The performance of the MG preconditioner is demonstrated in a number of numerical experiments.

3. An efficient multi-layer planar 3D fracture growth algorithm using a fixed mesh approach
by E. Siebrits and A.P.Peirce
International Journal for Numerical Methods in Engineering, 53, 691-717, 2002.

Abstract: We present a planar three-dimensional (3D) fracture growth simulator, based on a displacement discontinuity (DD) method for multi-layer elasticity problems. The method uses a fixed mesh approach, with rectangular panel elements to represent the planar fracture surface. Special fracture tip logic is included that allows a tip element to be partially fractured in the tip region. The fracture perimeter is modelled in a piece-wise linear manner. The algorithm can model any number of interacting fractures that are restricted to lie on a single planar surface, located orthogonal to any number of parallel layers. The multiple layers are treated using a Fourier transform (FT) approach that provides a numerical Green’s function for the DD scheme. The layers are assumed to be fully bonded together. Any fracture growth rule can be postulated for the algorithm. We demonstrate this approach on a number of test problems to verify its accuracy and efficiency, before showing some more general results.

2. The scaled flexibility matrix method for the efficient solution of boundary value problems in 2D and 3D layered elastic media
by A. P. Peirce and E. Siebrits
Computer Methods in Applied Mechanics and Engineering, 190 (45), 5935-5956, 2001.

Abstract: We present a method that extends the flexibility matrix method for multilayer elasticity problems to include problems with very thin layers. This method is particularly important for solving problems in which one or a number of very thin layers are juxtaposed with very thick layers. The standard flexibility matrix method suffers from round-off errors and poor scaling of the flexibility equations which occur when one of the layers in the multilayered medium becomes much smaller than the others. The method proposed in this paper makes use of power series expansions of the various components of the flexibility matrix in order to arrive at a system of equations that is appropriately scaled. The effectiveness of the scaled flexibility matrix method is demonstrated on a number of test problems.

1. Uniform asymptotic approximations for accurate modeling of cracks in layered elastic media
by A. P. Peirce and E. Siebrits
International Journal of Fracture, 110, 205-239, 2001.

Abstract: We present uniform asymptotic solutions (UAS) for displacement discontinuities (DD) that lie within the middle layer of a three layer elastic medium. The DDs are assumed to be normal to the two parallel interfaces between the leastic media, and solutions will be presented for both 2D and 3D elastic media. Using the Fourier transform (FT) method we construct the leading term in the asymptotic expansion for the spectral coefficient functions for a DD in a three layer medium. Although a closed form solution will require an infinite series solution, we demonstrate how this UAS can be used to construct highly efficient and accurate solutions even in the case in which the DD actually touches the interface. We present an explicit UAS for elements in which the
DD fields are assumed to be piecewise constant throughout a line segment in 2D and a rectangular element in 3D. We demonstrate the usefulness of this UAS by providing a number of examples in which the UAS is used to solve problems in which cracks just touch or cross an interface. The accuracy and efficiency of the algorithm is demonstrated and compared with other numerical methods such as the finite element method and the boudary integral method.