## Abstract

Reliable prediction of fracture process in shale-gas rocks remains one of the most significant challenges for establishing sustained economic oil and gas production. This paper presents a modeling framework for simulation of crack propagation in heterogeneous shale rocks. The framework is on the basis of a variational approach, consistent with Griffith’s theory. The modeling framework is used to reproduce the fracture propagation process in shale rock samples under standard Brazilian disk test conditions. Data collected from the experiments are employed to determine the testing specimens’ tensile strength and fracture toughness.

#### Authors

Mohaddeseh Mousavi Nezhad^{1}, Quentin J. Fisher^{2}, Elia Gironacci^{1}, Mohammad Rezania^{1}

^{1}School of Engineering, The University of Warwick, Coventry CV4 7AL, UK. ^{2}School of Earth and Environment, The University of Leeds, Leeds LS2 9JT, UK

Received: 24 February 2017 / Accepted: 30 January 2018 / Published online: 3 March 2018

© The Author(s) 2018.

To incorporate the effects of shale formation heterogeneity in the simulation of crack paths, fracture properties of the specimens are defined as spatially random fields. A computational strategy on the basis of stochastic finite element theory is developed that allows to incorporate the effects of heterogeneity of shale rocks on the fracture evolution. A parametric study has been carried out to better understand how anisotropy and heterogeneity of the mechanical properties affect both direction of cracks and rock strength.

## 1 Introduction

Hydraulic fracturing, commonly known as fracking, is an important process in extracting gas from shale formations. Accurate predictions of the directions and the extent of fractures are required in order to conduct safe and economically viable operations for gas production. Fractures tend to propagate along the direction of the least resistance path. The direction and extent of this path is a complex function of the in situ stress condition, anisotropic mechanical properties of the rock and pore and fracture fluid pressures (Warpinski and Smith 1989). It is possible to evaluate the influence of the in situ stress condition and pore-pressure variations on the evolution of fracture networks. This could be done using data from relevant laboratory tests, deep borehole logs from the field and coupled fluid–solid computational codes. However, the role of spatial variability and anisotropy of mechanical properties on the crack propagation remains poorly understood.

Numerous experimental and numerical studies, using a variety of techniques, have been performed for understanding the behavior of different types of rocks (Wong and Einstein 2009; Morgan et al. 2013). Less effort has been placed on the characterization of the mechanical properties of shale rocks. Study of the effect of bedding direction on crack propagation is a topic of great interest; recently, Morgan and Einstein (2014) conducted an experimental investigation on the effects of bedding laminations on crack propagation in the Opalinus shale samples and concluded that the shear strength and elastic modulus of shale samples as well as fracture advancement strongly depend on bedding laminates.

The experimental observations showed that the direction of stress-induced cracks which are generally perpendicular to the maximum principal stress changes toward the bedding orientations. These tests are crucial to understand the effect of material composition and local heterogeneities on the mechanical response of rocks, which, because of their natural genesis, significantly differ from one to another (Bobko and Ulm 2008; Veytskin et al. 2017). Hou et al. (2016) compared the behavior of shale, coal and sandstone specimens subjected to high-pressure fluid-driven fracturing, and observed how different microscopic composition of these materials leads to important differences in the measured values of mechanical properties, strains and crack patterns. Kim et al. (2012) also performed a comparative study between shale, gneiss and schist rocks samples, investigating the influence of their composition and direction of bedding planes on their thermo-mechanical properties.

The proper determination of mechanical properties of fracturing rock plays a fundamental role in predicting the mechanical response of the rock. Fracture toughness, the resistance of the rock to crack propagation, is considered to be a significant factor in determining the artificially created fracture network patterns. In particular, there is an extreme scarcity of published data and knowledge on the fracture toughness and elastic modulus of sedimentary geomaterials such as shales. It is particularly crucial to gain reasonable understanding of the fundamental processes controlling fracture propagation in shale rocks. This necessitates the undertaking of experimental and numerical studies to correctly assess the strength of these materials and identify the geomechanical characterization of these heterogeneous rocks.

Tests conducted on Brazilian disks specimens which are traditionally used to indirectly calculate rock’s tensile strength prove significant influence of laminations inclination on the mechanical response of the materials (Wang et al. 2016; Mokhtari and Tutuncu 2016; Mousavi Nezhad et al. 2018). Zhong et al. (2015) also studied the effect of damage at microscale on the overall macroscopic response of the shale samples subjected to triaxial compression and highlighted that both bedding inclinations and different confining pressures influence on the compressive strength, elastic modulus, crack morphology and opening modes. Holt et al. (2015) also studied the effect of local heterogeneities on the material response and rock brittleness. Mahanta et al. (2017) investigated the importance of understanding the effect of the applied load on shale specimens and how, for varying strain rates, the fracture toughness, and consequently the energy required for fracturing the material, changes.

Duan and Kwok (2015) and Zhou et al. (2016) modeled the Brazilian disk problem and a fluid-driven fracture problem, respectively, in shale using discrete element method (DEM). However, the main issue with DEM model is the computational cost, which becomes prohibitive if large-scale problems need to be analyzed. Classical finite element method (FEM) (Dokhani et al. 2016; Zeng and Wei 2017) and extended FEM (Sun et al. 2016) have been successfully applied to fracture propagation and fracking simulations on shale rocks. Anisotropy has been included in some of these numerical simulations by decreasing the values of mechanical properties, such as elastic modulus, cohesion coefficients or fracture toughness, in precise locations of the domain in order to create weak layers which affect fracture extents and shapes. However, Chen et al. (2016) showed that the introduction of weak bands is not enough to entirely reproduce the real anisotropic mechanical response of the shale rocks. It should be considered that in addition to the presence of the voids and microfractures, there are many types of inclusions, such as quartz, calcite and dolomite, forming a locally heterogeneous structure for the shales. In order to capture the real behavior of such heterogeneous materials, an in-depth knowledge of their microstructure and comprehensive numerical implementations are required (Chen et al. 2016; Mousavi Nezhad et al. 2016).

In the context of numerical modeling, strategies for characterizing the random heterogeneity in materials can be categorized as either multi-scale homogenization or spatially varying random fields. The former is attractive from the point of view that the different phases in a material such as the matrix, inclusions and the interfaces are modeled through identifying a representative volume element and computing the effective properties of complex microstructures (Belytschko et al. 2008). However, the applicability of the method is limited to the cases that a clear separation of scales exists and a closed-form macroscopic equation can be driven to describe the relevant microstructural details, e.g., morphology or constituent material properties (Geers et al. 2010; Novák et al. 2012). Spatially varying random field-based approaches are able to directly approximate the randomness through generation of spatial realizations of the properties associated with a given correlation structures, and the crack trajectory is modeled using FEM within the context of Monte Carlo simulation (Yang and Xu 2008; Stefanou 2009).

An attempt to use probability theory for modeling layered materials has been proposed by Bossi et al. (2016), in which the authors modeled a simplified slope made up by clay with horizontal gravel layers. A probabilistic technique has been used to assign the same deterministic value of gravel’s cohesion and friction angle to randomly selected elements forming the domain. This approach improved the reliability of numerical results, but mechanical properties have been deterministically defined. Borghi et al. (2015) also used a stochastic approach for introducing fractures with random dimensions and locations in the model. A noteworthy approach has been recently outlined by Li et al. (2016a, b) where material heterogeneity has been considered through a combined stochastic methodology accounting for different types of natural materials on a soil profile. Although this approach appears to be very promising, it does not account for material microstructures, which is well known to considerably influence the macroscopic behavior (Bobko and Ulm 2008; Gironacci et al. 2017).

This paper focuses on geomechanical characterization of shales and investigates the anisotropic mechanical behavior of shale specimens under Brazilian test conditions, as well as the effects of random orientation of shale layers on the failure mechanisms. The paper also reports original laboratory measurements on the mechanical properties of core samples of shale-gas reservoir rocks collected from the Dove’s Nest site on the east coast of North Yorkshire, in England. The data are presented on the general mechanical behavior of these rocks to delineate the basic parameters that control the propagation of cracks in shale-gas rocks. Elastic modulus and fracture toughness of shale specimens with different microstructural formations and layering orientations are calculated using Brazilian tests. A new computational modeling framework for simulation of crack propagation in heterogeneous layered materials is developed. It combines random field and continuum damage theories, and Weibull distribution and its spatial autocorrelation length information are used for generating samples of stochastic fields representing the layering structure as well as local heterogeneity of the shales. The numerical data are also discussed to better understand how anisotropy of the properties affects the morphology and direction of cracks.

## 2 Fracture Advancement Methodology

### 2.1 The Variational Approach

The basis of the variational approach to fracture mechanics relies on the association of a potential energy consisting of stored elastic energy, the work of external forces and the energy released through fracture to any crack and deformation configuration. A reference configuration ⊂R^{N}, N = 2, of a homogeneous elastic body is considered which contains a generic crack as visualized in Fig. 1.

**Fig. 1** Schematic configuration of a body with initial crack and applied boundary conditions

The total energy of the body is defined as

where f is the body deformation, K ⊂Ω¯ is the fractured zone, W: R^{NxN}→R is the stored energy function of a hyperelastic material, F is the deformation gradient, γ is the fracture energy, and H^{N−1} is the Hausdorff measure of K which provides the measure of the length of the crack for sufficiently regular fractured zone. The first and the second terms on the right hand side of Eq. (1) represent bulk and surface energy of the body, respectively.

Equation (1) can be minimized, so crack growth is deduced by successive minimization of energy at fixed time steps. The minimization of Eq. (1) with respect to any kinematically admissible displacement and any set of crack curves introduces a high level of complexity for the variation calculus of free-discontinuity problems, particularly due to the presence of non-smooth values of the K parameter. In fact, unless topological constrains are added, it is not usually possible to deduce compactness properties from the only information that such kind of energies are bounded. Following the methodology proposed by De Giorgi and Ambrosio (1988), we introduce K as a set of discontinuity points S_{f} of the function f and set the problems in a space of discontinuous functions. The weak formulation of the energy, by replacing the term K with a set of discontinuity points S_{f} of deformation in a Sobolev space SBV (Ω;R^{N}), is therefore given by

The presence of the term H^{N−1}(S_{f}) creates challenges in terms of finite element discretization of the functional. To overcome such challenges, Eq. (2) has been approximated, in the sense of Γ-convergence (Bourdin et al. 2000), by a family of numerically more tractable functionals defined over a Generalized Sobolev space GSBV. GSBV consists of all functions whose truncations are in Sobolev space (Ω;R^{N}) and allows extending the definition of Eq. (2) to L^{1} functions, which are not of bounded variation. Γ-convergence makes it possible to achieve a variational convergence. This means if a minimizer ϑ_{ε} for a function F_{ε} exists for every ε> 0 and if there is a sequence h↦ε_{h} such that ε_{h}→ 0 and the corresponding ϑ_{εh} converges to ϑ, then ϑ is a minimizer for F_{ε}. Based on the regularized formulation of the energy function for brittle fracture problems presented by Bourdin et al. (2000), an auxiliary variable s, which is called the damage parameter, is introduced. s is a regularized representation of the fractured zone defining the discontinuity set in Eq. (2). Therefore, a functional space X can be considered which its elements are pairs (f, s). Then, we take the functional X → [0, + ∞] defined by

with E(f) as in the problem, D is the domain of the functions ∈ GSBV, the problem (P_{0}) defined as min{F(f,s):(f,s)∈X} is considered. Damage parameter provides a picture of the damage state of the body; for an undamaged and intact body, s is equal to 1 everywhere, while it goes to zero in proximity of the discontinuity set S_{f}.

The functional formulation for a generic p > 1 is expressed in the form provided by Ambrosio and Tortorelli (1992)

where p is the pth power of the norm of the function defined in the Sobolev space, p’=p/(p-1),

is the normalization constant, κ_{ε} is a positive regularization parameter, and ε is related to the material length scale. Bulk and surface terms are two integrations over two different domains Ω and Ω′, physical domain and logical domain, respectively. Ω′ is defined as open set such that where ∂_{1}Ω and ∂_{2}Ω are the two disjoint parts of the boundary of Ω, and int∂_{1}Ω is the interior of ∂_{1}Ω relative to ∂Ω. The choice of the size of the logical domain is made on the consideration that it has to be big enough to avoid underestimation of the fracture energy when the crack reaches the boundary ∂_{1}Ω. For two-dimensional problem where p = 2, the total energy formulation for the body can be represented as

where Ω_{0} and Ω′_{0} represent initial unfractured and stress-free configuration of the body in the physical and logical domain, respectively.

## 2.2 Stored Energy Formulation

Following Del Piero et al. (2007), we used an isotropic, compressible neo-Hookean type stored energy model that is defined as

where μ is the Lamé’s second parameter, and C is the right Cauchy–Green tensor. In the above equation, the first term represents the classical formulation of an incompressible neo-Hookean material and the second term is a convex function that is defined as (Del Piero et al. 2007)

where μ is the Lamé’s first parameter and j=e^{(λ+μ)}/λ. Equation (7) is directly related to surface deformation, as it is function of the Jacobian of the deformation gradient (J). As J goes to zero the stored energy function goes to infinity, penalizing the extreme compression. To account for the tension–compression asymmetry of damage behavior of material, the methodology proposed in Li et al. (2016a, b) is followed, and the energy function is decomposed into two parts; a positive part which is considered to contribute to damage, and a negative part that resists to damage:

In the above equations, it can be noticed that the damage parameter appears only in the positive part of the energy function, the part associated with the elements that increase in surface (i.e., in the elements with J > 1), the value for damage in the elements is kept as calculated. The elements that decrease in surface (i.e., the elements with J < 1) do not contribute in damage. In this way, different behaviors for tension and compression are explicitly taken into account.

## 2.3 Numerical Solution Strategy

An approximation solution of the minimization of Eq. (5) is achieved using an iterative procedure, shown in algorithm 1, which consists of imposing a stationary condition to one of the deformation and damage variables, while keeping the other variable fixed. For all v ∈ W^{1,d }(Ω,R_{n}), w ∈ W^{1,d} (Ω_{0}), we look for a deformation that satisfies the stationary condition of

and then for the scalar field s stationary condition of

where S(F)=∂/∂F(W(F)) is the first Piola–Kirchhoff stress tensor. By taking the updated Lagrangian formulation of Eq. (11) and its linearization (Del Piero et al. 2007) the following is obtained

where using Eq. (7)Applying the integration by parts, the final weak form of Eq. (12) is derived as

MATLAB Partial Differential Equation Toolbox together with the Newton–Raphson iteration scheme is used to solve the above equations. The iteration stops when two consecutive pairs of solution (f_{n−1},s_{n−1}) and (f_{n},s_{n}) are close enough according to an identified convergence criterion. In order to avoid the healing of the cracks, an approximation method was used to consider irreversibility condition for damage evolution. We followed the methodology proposed by Del Piero et al. (2007). Based on that, irreversibility condition of s_{n}(x)=s_{n−1}(x) if s_{n}(x)>s_{n−1}(x) is set, and the value of damage parameter associated with each point in the body cannot exceed the one calculated at the previous time step. We leave the development of more advanced and rigorous methods for incorporating the irreversibility condition to future works.

An adaptive h refinement strategy (Del Piero et al. 2007) is used to automatically refine the elements with values of s lower than the given thresholds. At new nodes generated through the remeshing strategy, the values of displacement and damage are calculated by linear interpolation from the existing nodes.

### 2.4 Stochastic Approach

In order to incorporate the effects of material heterogeneity in the computational model, a statistical technique is employed. The Weibull distribution function (Weibull 1939) is used to generate random distributions of the properties, as it has a simple structure and its applicability for modeling failure of brittle materials has been verified (Fang and Harrison 2002; Yang and Xu 2008; Gorjan and Ambrožič 2012). In this study, the fracture energy is considered as a random variable.

Cumulative density function (CDF) and probability density function (PDF) for Weibull distribution, plotted in Fig. 2, take the form of

where ξ is the random parameter, m is the shape parameter, and ξ_{0} is the scale parameter.

**Fig. 2** Weibull function with different values of the shape parameter: a probability density function (PDF), and b cumulative distribution function (CDF).

In order to define both shape and scale parameters, maximum likelihood estimator (MLE) is used. Let ξ_{1}, ξ_{2}, … ξ_{n} be a random sample of size n with a PDF fξ_{i}(ξ_{i},m,ξ_{0}), where ξ_{0} and m are unknown parameters. The method considers the joint density function for all observations. For an independent sample with known ξ_{0}and m likelihood, function L(ξ_{0},m,ξ_{1},…ξ_{n}) is the joint density function of the n random variables defined as

On the basis of the MLE, the maximum likelihood of ξ is achieved by maximizing L or, for convenience, its logarithm

In order to apply the MLE to estimate the Weibull parameters, we can substitute Eq. (17) into Eq. (18) and apply Eq. (19); after a straightforward calculation one gets the following pair of equations

Equation (20a) is solved using the Newton–Raphson method to calculate the value of m as

where n is the Newton–Raphson iteration number,

The calculated value for m is substituted into Eq. (20b) in order to obtain the scale parameter ξ_{0} and define the Weibull distribution. Then, the inverse cumulative distribution function is used to generate random realizations over the simulation domain that is derived as

where R is a random number between 0 and 1 representing the probability of an occurrence. The full stochastic procedure is summarized in Algorithm 2.

The scale parameter can be related to the size of the elements forming the material. It has been shown in Guy et al. (2012) that, for a specific range of material length scale and crack length over a domain, the material length scale itself can be used as scale parameter of the Weibull distribution. In our work, the material length scale is accounted in the generation of the hypothetical sample used successively to generate the shape and scale parameters by means of the MLE method. In this way, material length scale has been implicitly taken into account in the model.

## 3 Laboratory Tests Over Shale Rocks

Theoretical equations to assess the strength of rocks, which are mainly based on the assumptions of isotropy and homogeneity, are not applicable for sedimentary rocks. Therefore, it is necessary to perform experimental studies to directly measure the mechanical properties of the shale samples. These experimental studies were conducted to investigate the failure behavior of the shale rocks under compressive pressure, and their mechanical properties were measured. Then, the experimental data have been used to verify the proposed computational method and to identify the influence of factors such as heterogeneities and isotropy on both fracture propagation and strength.

### 3.1 Experimental Setup

The Brazilian tests, used to evaluate the tensile strength, were conducted on cylindrical samples taken from a depth of 246 m. The experimental setup is shown in Fig. 3. The testing equipment consisted of a high-speed camera to take pictures of the fracture propagation, a strain gauge and a loading rig.

**Fig. 3** Experimental setup in the current study

As loading condition significantly influences damage and stress distribution inside the specimen, following the typical methodology of the Brazilian tests, the loads were applied to the specimens in two different ways: (1) with flat platens and (2) with curved jaws (Fig. 4). Failure initiates directly under loading points if flat steel plates are used, and the specimen crack initiates at the center of the disk if curved jaws are employed (Hobbs 1964).

**Fig. 4** Schematics of Brazilian disk test configurations, a with flat plates, and b with curved jaws

### 3.2 Experimental Method and Measurements

The fracture tests were carried out on specimens with different angles between loading direction and bedding plane orientations. The mechanical strength was determined by crushing the cylindrical specimens between steel plates. The variations of load, axial displacements and time were recorded from measurement devices connected to the loading frame (e.g., load cell and linear variable differential transformer or LVDT). Details of different testing specimens and measured data are presented in Tables 1 and 2. For all specimens, fracture toughness is calculated using the method proposed in Guo et al. (1993). This method relies on the value of P_{min}, which is the post-residual load being maintained by each fractured disk. The elastic modulus, E, was determined by attaching an electrical resistance strain gauge at the center of a disk to measure the lateral strains until failure.

**Table 1** Details of testing specimens and measured data using flat steel plates

**Table 2** Details of testing specimens and measured data using curved jaws

All test specimens are shown in Fig. 5; they were all tested in an air-dried state. In this figure it can be observed that all specimens had a well-defined layered structure, some with rather visible inclusions (e.g., in specimens C and D).

**Fig. 5** Tested shale specimens: specimens a and b tested using flat plates, and specimens c, d and e tested using curved jaws.

The indirect tensile strength (σ_{t}) of the cylindrical specimen under compression along the vertical diameter of the specimen was obtained using (Li and Wong 2013).

where P_{max} is the maximum applied load to the specimen before failure, R is the radius of the specimen, and t is its thickness. The tensile stress–strain plots from tests on all specimens are shown in Fig. 6. For the case of experiments conducted using flat plates, specimen B showed less stiffness, but a higher value of ultimate tensile strength and a wider post-elastic behavior. For the experiments conducted with curved jaws, the higher value of resistance was obtained for the specimen with a clear layered structure. To ensure the repeatability of the experiments, tests related to each case were repeated three times on specimens with similar characteristics and similar ranges for stress–strain variations have been observed for the test associated with each case. Therefore, the apparent differences among the plots presented in Fig. 6 can be associated with the variations of loading direction with respect to bedding direction.

*Fig. 6 Stress–strain plots from Brazilian disk tests on a specimens tested using flat plates, and b specimens tested using curved jaws.*