## Abstract

Recent analysis of gas outflow histories at wellheads shows that the hydraulic crack spacing must be of the order of 0.1 m (rather than 1 m or 10 m). Consequently, the existing models, limited to one or several cracks, are unrealistic. The reality is 10^{5}–10^{6} almost vertical hydraulic cracks per fracking stage. Here, we study the growth of two intersecting near-orthogonal systems of parallel hydraulic cracks spaced at 0.1 m, preferably following pre-existing rock joints. One key idea is that, to model lateral cracks branching from a primary crack wall, crack pressurization, by viscous Poiseuille-type flow, of compressible (proppant-laden) frac water must be complemented with the pressurization of a sufficient volume of micropores and microcracks by Darcy-type water diffusion into the shale, to generate tension along existing crack walls, overcoming the strength limit of the cohesive-crack or crack-band model. A second key idea is that enforcing the equilibrium of stresses in cracks, pores and water, with the generation of tension in the solid phase, requires a new three-phase medium concept, which is transitional between Biot’s two-phase medium and Terzaghi’s effective stress and introduces the loading of the solid by pressure gradients of diffusing pore water. A computer program, combining finite elements for deformation and fracture with volume elements for water flow, is developed to validate the new model.

This article is part of the themed issue ‘Energy and the subsurface’.

## 1. Introduction

Although fracking technology is now highly developed [1–7], the fracture mechanics of fracking is not. The existing theories, which deal in detail with the growth of one or a few parallel hydraulic cracks, with stress and flow singularities at the crack tip, represent only the micromechanics of the problem, while no scale bridging to the macroscale is in sight.

Based on the equations relating the measured gas flux histories at the wellhead (at Fayetteville shale basin in Oklahoma) to the typical permeability of shale [8–10], recent detailed analysis [11] showed that the crack spacing must be of the order of 0.1 m. If the spacing were 1 m or 10 m (as inferred by some from the typical 10 m spacing of perforation clusters), then the gas extraction would take about 100 times and 10 000 times longer, respectively, than it does. The inescapable conclusion is that currently the fracked portion of a fracking stage (only about 5–15% of the total volume) must involve about 10^{5} vertical branched intersecting cracks, and that about 10^{6} cracks would be required to frac the whole fracking stage (typically 70 m long, 150 m deep and 500 m wide).

How much is the estimate of crack spacing affected by errors in permeability? Not too much, because the spacing depends mainly on the time *τ*_{g} for the gas flux at the wellhead to drop from the peak flux to, say, 20% of the peak (called the ‘halftime’). This halftime, *τ*_{g}, varies roughly as the square of the spacing but depends on the permeability only linearly, i.e. much more weakly.

And how much is the spacing estimate affected by errors in the estimated delay (or transport halftime *τ*) of the gas transport from the hydraulic crack faces to the wellhead on the surface? Hardly at all. Should this transport be considered immediate (i.e. *τ*≈0), the rise of gas flux at the wellhead to the peak will not be captured, but the time, *τ*_{g}, for the flux at the wellhead to drop to 20% of the peak would change by a few per cent only. This time depends virtually only on the crack spacing, and more weakly on permeability (this is confirmed by simulations, which yield a curve that descends monotonically from the early maximum value of 1 and becomes almost identical with the optimum curve at the time of 2 years).

Because the typical depth of a gas- or oil-bearing stratum of shale is about 2.5–3 km [12,13] (figure 1), the overburden pressure generally exceeds the fluid pressure. Consequently, all the cracks and joints must be essentially vertical, and no horizontal cracks can form. One system of parallel distributed cracks is expected to be normal to the minimum principal tectonic stress *T*_{1}. The second system of parallel distributed cracks is expected to be nearly orthogonal to the first (and thus normal to the maximum principal tectonic stress *T*_{2}; figure 2). Why? If the cracks in one system of parallel cracks are open and thus free of shear stresses, the open cracks in a second intersecting system can also be free of shear stresses only if the systems are orthogonal, as required for local equilibrium (a similar question was discussed in the 1970s for distributed cracking in concrete; see appendix A). The same argument applies to the geological genesis of pre-existing systems of rock joints, which usually are almost orthogonal [5,14]. And, of course, the hydraulic cracks are expected to follow the rock joints wherever possible [15–17]. There are also frictional slip faults inclined to the principal stresses, but they must be tightly closed and thus cannot serve as conduits for water.

Often the crack branching has been pictured as V-shaped, occurring at the crack tip (figure 2). But such branching is possible only for dynamic cracks propagating at more than 40% of the Rayleigh wave speed, whereas fracking is generally a static process (except for small crack jumps due to heterogeneities that cause acoustic emissions). Hence, no V-shaped branching at crack tips is here considered.

Given that the spacing (figure 3) is about 0.1 m, the number of vertical near-orthogonal cracks to fill a typical fracking stage must be about a million. Yet, so far, all the previous studies based on fracture mechanics [4,19–24] have dealt with the propagation of one hydraulic crack or, very recently, a few parallel (i.e. non-intersecting) ones spaced by about 10 m [25–29]. These studies have been important for clarifying the micromechanics of the problem, especially the complicated interaction of crack tip singularity with viscous flow and capillarity near the tip. But a multiscale bridging of these results to the macroscopic problem of gas or oil extraction from the bulk of the shale stratum would be a formidable unsolved problem, theoretically interesting but probably unimportant for the global behaviour.

Large numbers of cracks can apparently be simulated by a black-box commercial program based on the discrete-element method, recently marketed in the fracking industry. This program apparently simulates hydraulic cracks by a narrow band of interparticle separations, or cracks (see [30], fig. 7). However, these narrow bands propagate from the casing perforations of a horizontal wellbore along planes that are widely separated, located roughly 10 m apart. Unfortunately, the simulated bands do not branch laterally and do not connect. Besides, the discrete-element method seems unable to capture the fracture mechanics energy release when the imagined discrete elements are scaled up to simulate large enough domains. Nor does it seem capable of reproducing the general triaxial constitutive damage behaviour and orthotropy of shale.

The hydraulic fracture system doubtless includes fractures running along pre-existing rock joints [31–33]. Although the joints are usually filled by calcite, their tensile strength and fracture energy are certain to be smaller, probably much smaller, than those of shale. Hence, the hydraulic cracks must be expected to run along the rock joints wherever possible. That they do is supported by the fact that the typical spacing of rock joints appears to be of the order of 0.1 m, which is roughly equal to the aforementioned spacing of hydraulic cracks deduced from observations of gas flux on the surface.

The ultimate goal is to simulate the whole fracking stage, which is typically 70 m long, 150 m deep and 500 m wide. The present discretization of fracture relies on the crack band model [34,35], in which the crack is smeared as equivalent damage over the finite element of a width equal to the crack band width. A tensorially enriched cohesive softening law is used to simulate the energy dissipation by fracture growth. The finite-element size is taken equal to the crack spacing, about 0.1 m. With this element size, the local stress and flow fields at the crack tip cannot, of course, be captured, but the effective global energy dissipation at the front of a zone of many growing hydraulic cracks can.

Although all the cracks are considered to be vertical, the stress and flow fields are three-dimensional (3D). About a billion finite elements would thus be needed to analyse the whole fracking stage, with its 3D stress field and flow field. Thus, the full problem would eventually require either the largest supercomputer, or a model based on larger-scale continuum smearing (or coarse-graining), which is beyond the scope of this paper. However, to gain understanding of the process and clarify the simulation methodology, it will suffice here to analyse a limited zone, which is the goal of this study.

## 2. Modelling approach to interactive fracture and water flow

Because of the previous conclusion [11] that the spacing of water-filled parallel cracks or joints must be of the order of 0.1 m, and that the cracks must be nearly vertical and mutually orthogonal, the shale mass will be subdivided into six-node cubic (hex) finite elements (with linear shape functions; figure 4). Each element may contain two vertical, mutually orthogonal (or orthorhombic) cracks normal to the *x*_{1}- and *x*_{2}-axes, coinciding with the principal tectonic stresses *T*_{1} and *T*_{2}. The opening widths of the cracks normal to these stresses, denoted as *h*_{1} and *h*_{2}, respectively, govern the components of the vectors (*q*_{1},*q*_{3}) and (*q*_{2},*q*_{3}) of the flux of water along these cracks, *q*_{3} being the vertical flux component. The water pressure, *p*, is considered to be uniform within each element, and since the cracks intersect, the pressure *p* in each element must be the same in both intersecting cracks. Pressure *p*, of course, varies from one element to the next.

The orthogonality of cracks is convenient for using the crack band model. The reason is that the mesh lines can be aligned with the cracks, which is the case in which the crack band model works best.

## 3. Flow of fracking water in shale elements with cracks

Initially, the fracking water was assumed to flow only along the opened cracks. However, with this simplification, it turned out that the computer simulations would not lead to any secondary branched cracks, although they must exist. Often the crack branching was pictured as V-shaped, at the crack tip. But such branching is possible only for dynamic cracks propagating at more than 0.4 of the Rayleigh wave speed, whereas fracking is generally a static process.

Therefore, the secondary hydraulic cracks must initiate orthogonally from the faces of the primary hydraulic crack (figure 2). According to the cohesive crack model or the crack band model [34–36], such cracks should initiate if the horizontal normal stress, *σ*_{11}, parallel to the face of a vertical crack exceeds the tensile strength limit, *σ*_{0}, of the cohesive crack or crack band model. Computations showed, however, that tensile stresses cannot be produced only by the hydraulic pressure in the primary crack. A similar conclusion is suggested by the analytical solution based on Westergaard’s complex stress function [36] for an infinite plane, which indicates that the pressure in a line crack in an infinite elastic plane causes no longitudinal tensile stress *σ*_{11} along the face of the line crack.

On the other hand, in the simple case of a uniformly pressurized porous spherical zone in infinite elastic space, the tangential normal stress at the hole surface is *ϕp*/2 [37], p. 395, where *ϕ*=porosity and *p*=pore water pressure. This stress is quite significant. So, we must consider that a significant volume gets pressurized before the secondary cracks branch out. This can be caused only by diffusion of water into the pores in the shale mass, which is equivalent to a leak-off.

Does the shale have sufficient porosity and permeability? It does, which is evidenced by various observations [38], e.g. (i) the porosity of a typical shale (e.g. Marcellus) is *ϕ*≈9%, even though the pores (and many pre-existing cracks) are of nanometre scale and (ii) only about 15% of the water injected from the wellhead returns to the surface after the fracking ends, which means that there must be a large leak-off into the pores in the shale.

Two kinds of water transport must be considered: (i) Darcy diffusion through the pre-existing pores and flaws in the shale and (ii) Poiseuille (viscous) flow along the created hydraulic cracks. The combined water transport is the essential idea of this paper, as discussed next.

## 4. Flow of water through hydraulic cracks and pores in shale

The flow along the hydraulically created cracks is assumed to follow the Reynolds equations of classical lubrication theory, which are based on the Poiseuille law for viscous flow. For the components *q*_{k} of the flux vectors in the vertical (*x*_{1},*x*_{3}) and (*x*_{2},*x*_{3}) crack planes (*x*_{1}≡*x*, *x*_{2}≡*y*, *x*_{3}≡*z*), this law reads
4.1
4.2
4.3where are the volumetric Poiseuille flow rates per unit cross-sectional area (dimension metre per second) in plane (*x*_{1},*x*_{3}); likewise , in plane (*x*_{2},*x*_{3}); pressure in the cracks; *h*_{1},*h*_{2}=openings of the cracks in (*x*_{2},*x*_{3}) and (*x*_{1},*x*_{3}) planes (i.e. in *x*_{1}- and *x*_{2}-directions); and *μ*′=effective viscosity of water, which differs from viscosity *μ* of water with gellants but without proppants, that permeates the pores (the flow nonlinearity due to the proppant is here neglected, although generalization is possible). Equation (4.3) gives the combined vertical flux in both cracks. The flow through the pores of shale is governed by the Darcy law [39]:
4.4where *k*=1,2,3; th component of the vector of volumetric flow rate through the pores, per unit cross-sectional area (dimension metre per second); pressure in the pores; *μ* = water viscosity [4,40] (taking into account the effect of gellants but not proppants); and *κ*_{h},*κ*_{v}=permeabilities of porous shale along the horizontal bedding planes and in the vertical direction (dimension square metre) [8]. The permeability properties of shale are orthotropic (or transversely isotropic) because, in the bedding planes, the clay mineral platelets have a predominantly horizontal orientation. Therefore, *κ*_{v}≪*κ*_{h} [38,41], and in computations *κ*_{v} is neglected.

Strictly speaking, it would not be difficult to introduce also a local flow rate of water, *q*_{L}, from the pores into the cracks within each finite element, governed by a scalar equation , where *b*_{L} is a local permeability constant and *l*^{−2} is the scaling factor for the element size, *l*. However, this transport of water is doubtless much faster than the global diffusion. Therefore, we will consider it as immediate, for the sake of simplification.

So we may add the fluxes through the pores and through the cracks (figure 4), and set
4.5
4.6
4.7where *Q*_{k} (*k*=1,2,3) are the combined flow vector components and *p* is the water pressure common to the pores and the hydraulic cracks.

## 5. Combined diffusion equation for flow of fracking water through shale

Because the water pressures in the cracks and in the pores are considered to be the same within each finite or volume element, the water mass conservation equation must be written for the combined mass of water in the pores and in the cracks of both orientations. When the combined Darcy and Poiseuille fluxes from equations (4.1) and (4.4) are substituted, the combined diffusion equation ensues as follows:
5.1where *l*=side of cubic element; *ϕ*=*porosity* in shale; *h*_{1}, *h*_{2}=opening widths of cracks normal to axes *x*,*y* (or *x*_{1},*x*_{2}) openings; *κ*_{h},*κ*_{v}=permeabilities of shale in horizontal and vertical directions (the vertical permeability, normal to the bedding layers, is much smaller); and, for the sake of simplicity, *μ*′ is considered equal to *μ*.

Although water was considered as incompressible in previous studies of hydraulic fracture, it is about 10–30 times more compressible than rocks. Therefore, the mass density of the fluid must be considered to depend on pressure *p*. This dependence may be linearized as
5.2where *C*_{f}=compressibility of water with proppant, *ρ*_{0}=mass density at reference pressure, and *p*_{0}=reference pressure, taken as gravity pressure of water at the depth of the shale layer before fracking. Water compressibility is, of course, unimportant for the overall deformation, because the transverse compliance of a crack is anyway negligible compared with the compliance of shale between two parallel cracks. But it does matter for the pressure changes in water due to rock deformation changes, and these affect the flow rate. And it also matters for the flow and local stress field near the crack front.

## 6. Calculation of crack opening using crack band model

In general, three mutually orthogonal (or orthorhombic) cracks could appear in one element. But horizontal cracks are prevented because, at the typical depth of 2.5–3 km, the vertical overburden pressure (about 80 MPa) generally exceeds the pressure of the fracking water (less than or equal to 70 MPa), and so only two vertical orthogonal cracks are allowed. The crack band width is taken equal to the material characteristic length *l*, which is equal to the element size, and is defined as the minimum possible spacing of stable parallel cracks; *l* must be distinguished from Irwin’s material length , which characterizes the length (rather than width) of the fracture process zone and is implied by the cohesive crack model, particularly its tensile strength *σ*_{0} and fracture energy *G*_{f} (and also Young’s modulus *E*). The flow of fluid alters the distribution of stresses in the solid and, vice versa, the fluid pressure affects the stresses and crack widths in the shale, which in turn affect the fluid flow and pressure. The crack openings, *h*_{i} (*i*=1,2), are, in the crack band model, calculated from the continuum damage strain or across the crack band, obtained from the softening constitutive law:
6.1and6.2where subscripts 1 and 2 refer to Cartesian coordinates *x*_{1},*x*_{2} in the horizontal plane; (or damage) strain, and *C*_{ijkl}=orthotropic elastic compliance matrix of the shale. The constitutive law for softening damage of the solid (shale), i.e.
6.3must also be orthotropic (its formulation is another challenge which still awaits optimal resolution and will be addressed in a separate paper using the microplane model); ** σ**,

**are the stress and strain tensors with Cartesian components**

*ϵ**σ*

_{ij},

*ϵ*

_{ij}(

*i*,

*j*=1,2,3), and

**is the tensorial function of strain tensor that defines the constitutive law.**

*f*The distributed microcracking damage is characterized by transversely isotropic damage parameters, *ω*_{1} and *ω*_{2}, considered to be functions of the inelastic parts of strains, and :
6.4and
6.5where *ϵ*_{b1}, *ϵ*_{b2}=specified breaking strain limits (figure 5). For cracks forming in intact shale (figure 5*a*), their values are taken as 10^{−5} and, for cracks running along pre-existing cemented joints (figure 5*b*), as 100 times smaller. The joints are treated as planes (or zero-volume layers) with a cohesive crack (or crack band) model of reduced strength *σ*′_{0} and reduced fracture energy *G*′_{f}. Because the planes have zero volume, the constitutive law for the shale mass, ** σ**=

**(**

*f***), remains unaffected by the presence of joints. The joints merely represent planes of reduced strength and reduced fracture energy. These two reduced properties are embedded in the crack band model [34,35], which, for computational purposes, serves to spread out and smear the cohesive crack into a band of finite width.**

*ϵ*## 7. Three-phase medium for the effect of fluid pressure on the deformation of shale

The shale with cracks and pores represents an amalgamation and generalization of both Biot’s two-phase medium and Terzaghi’s effective stress concept [8]. There are three phases: (i) the solid shale, (ii) the cracks, and (iii) the pores. The shale is assumed to be transversely isotropic in the horizontal planes. Considering no horizontal cracks to exist, one can formulate the equilibrium of the three phases on the crack planes as follows:
7.1
7.2
7.3
7.4where *ϕ*=natural pre-existing porosity of shale (including the pre-existing microcracks or flaws); *ω*_{1},*ω*_{2}=additional porosities for the resultants of fluid pressure in the *x*_{1}- and *x*_{2}-directions (figure 6); and *S*_{ij} are the total stresses in the solid–fluid system, which are used in finite-element analysis to calculate the nodal forces. So, the contributions to the total stress, *S*_{ij}, are three:

(1)

*σ*_{ij}is the stress tensor in the solid, which is calculated from the constitutive law of shale in equation (6.3), including the natural pores (since the stress in the constitutive law for the solid is defined as the force per total unit area, the pores included,*σ*_{ij}is not multiplied by 1−*ϕ*);(2)

*ω*_{1}*p*,*ω*_{2}*p*are the cross-sectional resultants of water pressure in the cohesive (bridged) hydraulic cracks (figure 6); and(3) (1−

*ω*_{1})*ϕp*and (1−*ω*_{2})*ϕp*are the cross-sectional resultants of the pressure of water in the pores.

Note that one may also write *σ*_{11}=(1−*ϕ*)*s*_{11} where *s*_{11}=stress resultant from cross section of solid without the natural pores and pre-existing microcracks, whereas *σ*_{11}=stress that is spread out over the full cross section of solid including the pores. The use of *σ*_{11} instead of *s*_{11} is preferred because *σ*_{11} is the stress measured in material tests and calculated from the material constitutive model.

It is interesting to note that, if there are no cracks with damage, i.e. *ω*_{1}=*ω*_{2}=0, equation (7.2) reduces to
7.5which is the well-known equilibrium relation for the total stress in Biot’s two-phase medium. On the other hand, if the damage is complete, i.e. if *ω*_{1}=*ω*_{2}=1 (no bridges across the crack), equation (7.1) reduces to
7.6which is the well-known Terzaghi effective stress. Thus, the present three-phase medium represents a continuous transition between these two classical concepts of soil mechanics.

## 8. Numerical simulation of fracking

Although the stress analysis of the solid, the shale, is handled by the finite-element method, the fluid flow is better handled by the finite-volume method [42–44]. To combine both methods, the fluid flow elements overlap the finite elements for stress analysis, as shown in figure 7. While the unknowns in the finite-element method are the nodal displacements, the unknowns in the finite-volume method are the water fluxes *q*_{1} and *q*_{2} through the finite-volume sides. The key advantage of the finite-volume method is that it enforces the mass balance condition exactly, which is particularly important at the front of the water infiltration zone.

Strictly for the purpose of verifying the method, a two-dimensional analysis of fracking of a square portion of a horizontal layer of shale, of unit height and size 5×5 m, is conducted first (figure 7). The boundaries are mechanically represented by elastic springs at each boundary node, such that their elastic stiffness would approximately characterize the resistance to a uniform expansion in an infinite horizontal elastic plane. The overburden pressure is 80 MPa, corresponding to the depth of 3 km, and the tectonic stresses are *σ*_{h}=30 MPa and *σ*_{H}=40 MPa. In figure 7, a horizontal wellbore runs at the bottom of the square shown, where detonations of shape charges create three perforation clusters. Water, which is considered to contain enough proppants to prevent any crack closing, is injected at these three clusters at prescribed pressure history *P*_{0}(*t*), starting with a rapid rise to 70 MPa (rapid but static, not emitting waves).

The mechanical behaviour of shale, which is transversely isotropic in the horizontal plane, is described by a microplane model, to be presented in a separate paper, and by the crack band model [45] with post-peak softening, corresponding to material characteristic length *l*=0.1 m (the creep of shale, and the embedment creep of proppant grains, is not yet included in the modelling, but eventually will have to be, for long-time predictions). The porosity of shale is considered as *ϕ*=0.09, and the water viscosity is 8.9×10^{−4} Pa s. The tensile strength of shale and fracture energy are here considered to be *σ*_{0}=4.17 MPa and *G*_{f}=0.7 N mm^{−1}, but 3.62 MPa and 0.5 N mm^{−1} when the cracks run along rock joints (these values are considered much higher than measured without confining pressure on specimens brought to the surface).

Two cases are run: (i) no joints, only potential fractures in intact shale, and (ii) predefined joints, assumed to run through the centre planes of all elements, initially closed, representing potential cohesive cracks. Figures 8 and 9 show the snapshots (at various times, as marked) of the distribution of excess hydraulic pressure. The finite elements with the highest pressure reveal the location of hydraulic cracks. For fluid flow, the boundaries are considered as impervious. Because the elastic boundary representation is only approximate, the simulations are stopped before the growing fracked zone of high pressure approaches the boundary. The difference between fracturing through the intact shale and along the pre-existing joints indicates that the tensile strength and fracture energy of shale are not very important. Much more important is the tectonic stress. Figure 10 shows the evolution of hydraulic crack width. The width (exaggerated) is approximately proportional to line thickness and the greatest line width shown corresponds to the width of 10 mm. The proppant is assumed to prevent any reduction of crack width.

Subsequently, three-dimensional simulations, which are computationally far more demanding, are conducted for a cubic block of a smaller side, 0.9 m (figure 11). Here only one perforation cluster, i.e. one inlet of water from the horizontal wellbore, is considered. All the other parameters are the same. Four subsequent snapshots of the pressure patterns are shown in figure 12. The first shows the pattern seen on the surface of the block, the second on a centric horizontal cut of the block, the third on a centric vertical cut normal to the wellbore, and the fourth on a centric vertical cut parallel to the wellbore. One joint is considered in this case. Finally, figure 13 shows the crack width seen on the block surface and on the horizontal cut.

Studies of the effects of various parameters, made on a large shale volume and run on a larger supercomputer cluster, have to be relegated to a subsequent paper. Of particular interest will be identification of the control parameters that maximize the area of the hydraulic cracks per unit volume of shale, and allow it to frac the largest volume of the fracking stage. The challenge of ‘coarse-graining’ of the element subdivision while keeping the correct fracture energy release and capturing the crack-localization instabilities and bifurcations should be addressed to reduce the heavy computational cost.

## 9. Conclusion

(i) The hydraulic crack system must reasonably be expected to consist of roughly orthogonal vertical cracks (whose spacing should, according to previous analysis [11], be of the order of 0.1 m). This crack system geometry simplifies the analysis.

(ii) Secondary cohesive cracks must initiate nearly orthogonally from the walls of primary hydraulic cracks, triggered by tension along the crack walls that overcomes the tensile strength limit of the cohesive crack or crack band model [34–36]. V-shaped branching at the crack tip is not a realistic hypothesis for static crack growth.

(iii) The

*key idea*to explain and model the tension along crack walls, which gives rise to crack branching, is that a sufficient volume of shale surrounding the crack must get pressurized by Darcy diffusion of water into the micropores and microcracks in the shale, occurring simultaneously with Poiseuille flow of water along the hydraulic cracks.(iv) A

*second key idea*is that the equilibrium of interacting cracks, pores and water must be enforced through a three-phase medium concept, which is transitional between Biot’s two-phase medium and Terzaghi’s effective stress concept, and amalgamates these two classical concepts.(v) Many or most of the vertical hydraulic cracks probably follow pre-existing rock joints, which usually form two nearly orthogonal systems. Even though cemented by calcite, the joints doubtless have a much lower tensile strength and fracture energy than the intact shale.

(vi) The problems of softening strain localization and spurious mesh sensitivity can be effectively avoided by the crack band model, in which each element is imagined to contain at the outset a potential cohesive crack, with fracking water flowing through if the crack is formed.

(vii) Numerical examples document that the present model is capable of simulating a large system of branched and intersecting hydraulic cracks.

## Competing interests

We declare we have no competing interests.

## Funding

Funding from the US Department of Energy through subcontract no. 37008 of Northwestern University with Los Alamos National Laboratory is gratefully acknowledged. Crucial initial funding was provided under Grant 36126 from the Institute for Sustainable Energy (ISEN) of Northwestern University. Ramification of the fracturing analysis from concrete to shale was partially supported also under ARO grant no. W911NF-15-101240 to Northwestern University.

## Acknowledgements

Prof. Marco Salviato of the University of Washington and Prof. Gregory J. Wagner of Northwestern University deserve thanks for valuable discussions.

## Appendix A. The problem of initially non-orthogonal cracks

The present model implies the cracks to be always orthogonal, even during the initial cohesive crack opening stage. This is, of course, a simplification. Before the crack fully opens, i.e. before *ω* becomes 1, there are cohesive bridging stresses and thus shear and frictional stress can be transmitted across the cracks, and so non-orthogonal partially formed cracks are certainly possible.

A generalization to the frac water flow through a 3D system of non-orthogonal cracks is easy. The following equations have been derived:
A 1where
A 2Here ** I**=unit tensor,

**∇**=3D gradient vector,

**=**

*n**unit*vectors of the normals of

*n*

_{c}parallel crack systems

*k*=1,2,…,

*n*

_{c},

*h*

_{k}their opening widths,

**=combined frac water flux vector and**

*q***=permeability tensor, which is generally anisotropic. In addition, to make analysis of the full fracking stage possible, a larger-scale homogenization allowing computations with a coarser discretization will have to be formulated.**

*B*However, initially non-orthogonal cohesive cracks, gradually transiting to orthogonal fully opened cracks, bring about a major complication from the solid mechanics viewpoint. This problem was discussed in depth several decades ago for reinforced concrete, and an artifice called ‘rotating crack’, which implies the initial cracks to close and cracks of a new direction to open, has been used.

Anyway, the analysis of a rotating direction of dominant cracks and its relation to pre-existing rock joints of fixed direction is beyond the scope of this paper, and is probably unimportant for the global picture.

## Appendix B. Environmental benefits and side effects

Since denser and narrower cracks generally make smaller dynamic jumps, emitting weaker waves, the seismicity would also be reduced. The cause of these sound-emitting jumps is that the material is not homogeneous. The fracture energy of shale fluctuates up and down along the crack path. At every drop down, the crack temporarily loses stability and makes a dynamic jump depending on the magnitude of strain energy that can be released. For denser cracks, this magnitude is smaller; hence smaller jumps, smaller microseismic events. While seismicity in shale fracking is generally benign, this point may be of interest for ramifications to deep sequestration of CO_{2} or various toxic fluids, for which a much higher seismicity than in fracking has been experienced.

A realistic model for hydraulic crack growth will probably allow increasing the fracked volume portion and thus the percentage of gas extraction. This should further lead to cracks of denser spacing, with smaller and more uniform opening width, which will have a welcome effect—a decrease in the amount of water per crack surface area, which will reduce the amount of water required for injection, decrease the return of contaminated water per unit amount of gas produced, and thus bring about a beneficial environmental side effect.

## Footnotes

One contribution of 12 to a theme issue ‘Energy and the subsurface’.

- Accepted May 21, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.