An approach to hydraulic fracture in geomaterials through a porous brittle damage material model

A coupled hydromechanical code based on a distributed brittle damage material model has been developed in order to simulate hydraulic fracturing processes. As speciﬁc feature, the code does not assume the development of a discrete number of sharp cracks but considers a diffused damage that takes place under overall conﬁnement of the system. The coded is validated against a laboratory test reproducing the local effects of a hydraulic fracturing procedure. The numerical model of the laboratory test is then used to conduct a parametric analysis where the inﬂuence of the principal hydromechanical parameters on the hydraulic fracturing performance is assessed. Finally, the code is applied to simulate multiple hydraulic fracturing in-terventions in a large scale shale reservoir. In all applications considered, the model shows a good predictive capability in terms of damage extension.


Introduction
Hydraulic fracturing is a process characterized by the nucleation and propagation of multiple diffused fractures in soil or rock masses as a consequence of a localized solicitation driven by hydraulic pressure. Hydraulic fractures are found frequently in normal geophysics and geology situations, e. g., magma-driven propagation cracks [48,42]. In industrial applications hydraulic fractures are created artificially: heat production from geothermal reservoirs and the estimate of rock mass stress state [3,54,18] are just a couple of examples. In the last decades hydraulic fracturing has been employed massively in the field of extraction of hydrocarbons from unconventional reservoirs [47]. In this context, the United States Environmental Protection Agency defines hydraulic fracturing as a process that creates artificial fractures in the rock formation to facilitate the flow of natural fluids, therefore increasing the volume of recovered hydrocarbons. The first attempts to improve the production of reservoirs through fracturing started on the late 1930s making use of artificial explosions, and later of acidizing treatments. The first hydraulic fracturing process specifically for stimulation was performed in 1947, in western Kansas (USA). Hydraulic fracturing is now used extensively in the petroleum industry to stimulate hydrocarbon wells, in order to increase or, in many cases, to activate their production, cf. [23] and [47]. The oil or gas production as a consequence of a hydraulic fracturing treatment is influenced strongly on length, extension and width of the induced fractures [40]. Although hydraulic fracturing is associated commonly with low permeability gas shales, it is has been applied also to other diverse geological situations, including weakly consolidated offshore sediments, coal beds, and naturally fractured reservoirs [1].
To design satisfactory hydraulic fracturing treatments in a reservoir, it is of outmost importance to estimate qualitatively and quantitatively the evolution of the fracture network as a function of the parameters of the process through analytical or numerical models. Regrettably, modelling of hydraulic fracturing processes is a demanding task, due to the numerous coupled processes involved: (i) the mechanical deformation induced by the localized peak of fluid pressure within the porous medium; (ii) the flow of the fluid within the fracture and the fluid exchange with the surrounding porous medium; (iii) the actual propagation of the fracture surfaces that modify the permeability of the medium [1,10]. The hydro-mechanical process is transient in nature, and its evolution is characterized by a strongly non-linear, non-local and history-dependent response: the actual evolution of fractures in a heterogeneous formation is a complex phenomenon and it is very difficult to obtain reliable predictions through simple models.
Looking back at the literature in the topic, early attempts to model hydraulic fracturing were based on analytical solutions, obtained by assuming that a unique welldefined crack propagates within an elastic medium under the action of a pressure, thus disregarding the hydromechanical coupling and nonlinear effects. Such models are not able to account for the fluid flow within the reservoir, the damage processes around the stimulated area, and the fluid-solid interaction. Moreover, such models do not account for the existence of a moving boundary and for the stress singularity at the crack tip. Classical examples in two-dimensional setting are: the Perkins-Kern (PK) model [39], based on plane strain crack solution; the Perkins-Kern-Nordgren model, an extension of the PK model accounting for the effect of fluid loss [36]; the Kristianovic-Zheltov-Geertsma-de Klerk (KGD) plane strain model [26]. The strong simplifying assumptions at the basis of these models imply severe limitations to their applicability and predictive capabilities, as demonstrated by field measurements [52]. The remarkable differences between model predictions and field observations can be justified possibly by the presence of complex multiple fracture systems induced by the hydraulic solicitation, or the presence of a layer of small cracks around the main hydraulic fracture [5,6], or by the infiltration of the pumping fluid in the porous rock [13]. In spite of such limitations, the use of analytical models is still widespread, because of their negligible computational cost. Morevover, analytical or semi-analytical models proved able to give insights on the competing processes that control the propagation of fluid-driven fractures, for both impermeable and permeable materials (see, e.g., [43,25,8,29,9,19,21,33,22]), evidencing the presence of different timescales and length scales in the evolution of a single fracture. Recently, theoretical models have been used to describe the propagation of arrays of closely spaced fractures [9,32]. Despite scaling provides interesting frameworks to classify different regimes of hydraulic fractures as well as to design small-scale laboratory tests that mimic field conditions, this procedure can be hardly applied for complex injection geometries or in situations where the state of stress is not uniform, or the reservoir material is heterogeneous. Finally, the microstructural features of rocks make scaling procedures rather complex (if not impossible), due to the dramatic dependence of rock hydro-mechanical properties on its microstructure: scaled model tests performed on brittle-elastic materials like glass or polymethyl methacrylate is precluded for real rocks [17].
The necessity of developing fully three-dimensional models thus stems from the typical material and stress anisotropy observed in rock formations, which causes the deviation of cracks from an ideal planar path and the non easily predictable propagation across multiple layers [40]. Thus, pseudo-three-dimensional models (P3D) have been proposed [46,44,24,34,35], and further extended to account for complex fracture networks or natural fractures [51]. More versatile approaches include planar-3D (PL3D), based on a 2D description of the coupled flow equations and of the fracture footprint [44,2], ad fully 3D models, based on elasticity equations and used to describe the changes of the width of an existing fracture opening displacement with the pressure of the fluid. The computational cost, higher than the one requested by P3D models, is balanced by the possibility to obtain reliable results also in geological situations characterized by irregular stress and inhomogeneities [1].
With the improvement of the computational power, in recent years plenty 3D numerical approaches have been proposed to tackle hydraulic fracturing. Approaches are based on Finite Difference (FD), Finite Element (FE), Boundary Element (BE), Distinct Elements Methods (DE), or a combination of the different methods. Most approaches rely on the a priori knowledge of a single crack geometry and topology. A 3D numerical scheme based on a surface integral method to find crack opening for any pressure distribution obtained via a 2D FE discretization of the continuity equation has been discussed in [14].
A FE scheme to simulate the 3D evolution of non-planar hydraulic fractures, accounting for hydro-mechanical coupling, linear elastic fracture mechanics, and lubrication theory was discussed in [11]. Hydromechanical coupling has been treated with a mixed approach, where the evolution of an existing crack was analyzed with a geomechanical FE solver and the flow equations was solved with a FD scheme [30]. The literature testifies a recent effort to develop more reliable models of fracture propagation in the context of hydraulic fracturing [28]. Zero-thickness cohesive elements with proper fracture criteria and constitutive laws have been used in [10]. Cohesive fracture propagation under hydraulic pressure, simulating explicitly the fluid flow across of network of discrete fractures, has been discussed in [45]. The Extended Fi-nite Element Method [31] and the Virtual Element Method [7] have also been used with the aim to increase the predictive capabilities in reproducing the evolution of fracture networks. All these approaches rely on the explicit modelling of pre-existing fractures or defects, a very strong limitation because in field applications the location of fractures is unknown. Further shortcomings are related to the difficulties to model the interactions between natural and induced fractures, the complexity of the adaptive discretization procedures, and the overall computational cost, unacceptable for simulations at the reservoir scale.
Alternatively, by adopting a continuum point of view, fracture networks induced by hydraulic fracturing can be described by a complex pattern of micro-cracks in contrast to a set of well-defined and localized cracks. Using this approach, the modification of the microstructure due to the progression of damage can be rendered readily in terms of modification of hydraulic conductivity, and linked to the reservoir production. Using a continuum approach, it is possible to ascertain the dominant role of shear failure with respect to tensile failure in hydraulic fracturing processes [12]. A continuum approach based on damage mechanics has been considered in the 3D fully coupled thermal-porous-mechanical FE model discussed in [53].
As well known, damage models are affected by mesh dependency and their use requires the introduction of a characteristic length, which can be embedded in nonlocal models. By the way of contrast, in the present study we use a continuum model of distributed brittle damage characterized by multiple length scales, that is not affected by mesh dependency. The material model has been introduced in terms of linearized and finite kinematics in [15,16], and embedded into a in-house finite element code for the solution of coupled hydro-mechanical problems. The approach has the built in capability of dealing with the interaction between closely-spaced fractures, both natural and induced by the treatment, as well as with the leak-off of the fluid into a permeable rock. Complex geological situations can be also easily accounted for, including reservoir heterogeneity, non-uniform far-field stress states and material anisotropy, as generally encountered in shale reservoirs. Finally, the constitutive model used at the material level can be easily calibrated on the basis of classical laboratory tests performed on rock samples, like uniaxial and triaxial tests, proving able to reproduce the most relevant peculiarities of rock behaviour like the confining stress dependence and dilation upon shearing. We discuss the validation of the code against experimental results and use the validated code for a parametric analysis and the simulation of hydraulic fracturing processes in a shale reservoir.

Hydromechanical model
We adopt the standard hydro-mechanical framework for porous media in linearized kinematics. In consideration of the highly nonlinear behaviors involved in fracture and frictional contact, in this study we opt for Therzaghi's approach to porous mechanics, disregarding the more general approach due to Biot. The adoption of the simplified Terzaghi model is based on the assumption of negligible compressibility of solid grains: as well recognized, cf. [20], this ansatz is acceptable for soils, yet debatable for sandstones and mudstones. Our assumption has been forced by the need of obtaining a simplified model for wide field applications and motivated by the following reasons. Biot's theory has been developed for linear elastic porous media and, in particular, the Biot-Wills coefficient α, introduced to deliver a modified effective stress able to account for solids and fluid compressibility, was derived under the assumption of linear elastic behavior. In problems where failure is not of concern and linear elasticity is a good approximation of rock behavior, Biot's approach has been proved to be predictive in terms of volume change, pore pressure built up, and expelled fluid volumes [20,50]. By the way of contrast, the interpretation of failure data on sedimentary rocks, i.e., shales or limestones of main interest for hydraulic fracturing, show that the assumption α=1, corresponding to Terzaghi's approach, is more in line with the experiments [27,41]. The key mechanical issues of the present work are related to diffused fracturing and failure and, additionally, available failure criteria are expressed in terms of Terzaghi's effective stress. Finally, the effects of elastic behavior of the porous matrix are negligible in all the considered applications.
Note that the simplifying assumption of negligible fluid and solid compressibility affects the constrained specific storage coefficient of the material (i.e., the inverse of Biot modulus M), which attains a null value. However, the particular applications of the proposed model render this assumption acceptable from the engineering point of view. In fact, for mass balance purposes, the significant increase in the hydraulic conductivity of rocks due to the opening of multiple fractures and the consequent increase of the fluid flow diminish the relevance of the contribution compressibility of the fluid.
Assuming full saturation of the void space, incompressibility of the solid particles, and incompressibility of the fluid, the linear momentum balance and the fluid mass balance read: where σ is the Cauchy stress tensor, b is the body force vector, v is the seepage velocity, n is the porosity (i.e. the volume of the voids over the total volume of the porous medium), and t the time. Porosity changes are related to volume changes through the relation where ε v is the volumetric strain, assumed positive for expansion. According to Terzaghi's effective stress principle, the deformability and the strength in the porous medium are related to the effective stress σ , defined as where p denotes the pore fluid pressure and I the identity tensor. Accordingly, all the constitutive equations will be expressed in terms of effective stresses: where ε denotes the small strain tensor.
For the fluid flow, we assume a linear relation between fluid velocity v and the gradient of the hydraulic head h, according to the Darcy law: where k is the permeability tensor, ρ f and µ f are fluid density and fluid viscosity, respectively, g is the gravitational acceleration, and h is the hydraulic head Note that in h only the elevation head and the pressure head are considered, while the kinetic contribution is neglected due to the low seepage velocity. Static boundary conditions on the Neumann boundary Γ t of outward normal n require wheret are the surface tractions. On the Dirichlet boundary Γ u the displacement assumes known valuesū u =ū , Hydraulic boundary conditions on the Neumann boundary Γ q of normal n require q · n =q n , whereq n is the assigned flux. On the Dirichlet boundary Γ p the pressure assumes known valuesp p =p .

The constitutive model
We recall briefly the features of the linearized constitutive model, described in detail in [15]. The material is initially homogeneous, linear elastic and isotropic, and reacts to the application of loading according to the Hooke's law. Critical conditions of loading lead to the attainment of ultimate tensile or shear stress states, that induce the inelastic behavior with formation of multiple parallel planar fractures, named 'faults' in reason of their repetitive structure in sharp contrast to the formation of a single, localized fracture. The inelastic behavior manifests itself in terms of reduction of stiffness of the material, therefore the material model falls in the class of damage models and it is named brittle damage model. Unlikely standard damage models, based on a phenomenological description of the material response that leads to pathological mesh dependent effects, the model here considered is microstructured, or complex, and therefore it is not plagued by dependency on the discretization. Brittle damage microstructures assume the aspect of equidistant planes, with an orientation N and a spacing L, where the scale of the spacing is different from the scale of the body. Microstructures are embedded in the elastic matrix and are recursive: a second set of parallel faults, of orientation N 2 and spacing L 2 , can form within an existing set of faults, of orientation N 1 and spacing L 1 , that provide a container for the kinematics of the second set of faults. The outermost set of faults is called rank-1, and the rank of the microstructure correspond to the number of nested sets of faults.
We begin to describe the damaging behavior of model in the case of the preexistence of an unique family of faults, disregarding for the moment how the set of faults forms. To simplify the notation we drop the index 1 from the symbol of orientation and spacing. The kinematics of the damaged material accounts additively for the deformation of the matrix and for fault opening, in the form where u is the displacement vector, ε m is the deformation of the matrix and ε f the deformation due to fault activity. The deformation due to the fault is related to the the relative displacement of the two surfaces of the fault (displacement jump) ∆ according to where ⊗ denotes the dyadic product. The behaviour of the matrix is assumed to be linear elastic, defined by the Young modulus E and the Poisson ratio ν. The behaviour of the faults follows the classical cohesive theories: it is assumed that during the early stages of damage the separation of the fault planes is resisted by the cohesive traction T. The relation between the cohesive traction and ∆ is expressed through a cohesive law T = T(∆, q), where q denotes a set of variables accounting for the irreversibility of the fault separation process. Here we assume a scalar cohesive law, expressed in terms of an effective opening displacement ∆, defined as the weighted average of the normal component, The weight β can be seen as the ratio between the shear and the tensile strength of the material [37]. Adopting a thermodynamically consistent approach, we postulate the existence of a cohesive energy per unit surface Φ(∆, q), dependent on ∆ and on an appropriate set of internal variables, q, necessary to describe irreversibility. The cohesive law follows as through an effective cohesive traction T (∆). We assume a simple effective linear decreasing law, see Fig. 1, defined by two parameters: the tensile strength T c (corresponding to the peak) and the critical energy release rate G c (corresponding to the area enclosed by the cohesive law). Cohesive forces vanish once the critical opening displacement ∆ c = 2G c /T c is attained. Irreversibility is enforced by assuming elastic unloading up to the origin, thus the only internal variable needed is the maximum attained effective displacement q = ∆ max . Upon attainment of full decohesion, friction remains the sole dissipation mechanism for the material. Friction is included in the model through a dual dissipation potential per unit area Ψ * [38], which for the classical Coulomb model reads: where µ = tan ϕ is the friction coefficient and ϕ is the friction angle, N · σN the normal component on the traction vector on the fault plane, and |∆| is the norm of the rate of the displacement jump. Upon time discretization the local state of the material at time t n+1 can be obtained variationally. We introduce an incremental work of deformation per unit of volume E n that accounts for all the energies involved in the process. We assume the state of the material at time t n to be known and the total deformation ε n+1 at time t n+1 to be assigned. The incremental work of deformation is where W m (ε m n+1 ) is the elastic strain energy density per unit volume of the matrix (depending on the deformation of the matrix), Φ/L is the cohesive energy density per unit volume, and the third term accounts for the frictional dissipation density in the time interval ∆t. The inelastic variables characterizing the state of the material at time t n+1 = t n + ∆t, i. e., ∆ n+1 , q n+1 , are obtained through a constrained optimization conducted with respect to the inelastic variables ∆ n+1 and q n+1 , that accounts for the irreversibility of decohesion and for impenetrability of the faults upon closure The constant C φ may account for the presence of granular additive, that can be used to maintain the fracture, produced by pressurization, open even after the dissipation of the fluid pressure.
The constrained minimization, formulated as leads to a set of four equations and two inequalities in the unknowns ∆ n+1 , q n+1 . The solution is sought by setting one of the two inequalities as an equality and solving the equations not automatically satisfied. If also the second inequality is satisfied, the solution is accepted. Otherwise, the procedure is repeated by assuming the satisfaction of the second inequality [37].
As said, at the beginning the material is purely elastic, and faults are not present, thus the orientation N is not defined, and it has to be computed during the calculations on the basis of the applied load. Prior to onset of damage, at every load increment the purely elastic solution is computed. If at a particular integration point the effective stress violates the Mohr-Coulomb (under overall compression) or the Rankine (under a tensile state) criteria, faults are inserted according to the attained failure condition. Notably, in [15] has been proved that the orientation of the faults is the natural outcome of energy optimization, if N is considered as an unknown unit vector in Eq. 16, see the Appendix in [15] for details. Note that, for the Morh-Coulomb criterion, the parameter β coincides with the friction coefficient µ.
By including an additional energetic term, that describes the misfit energy necessary for the accommodation of the faults within the outer fault family, the minimization process provides also the spacing L. The misfit energy represents the energy stored in the boundary layers that form where the faults meet a confining boundary. Since the compatibility between the faults and their container is satisfied only on average, material develops boundary layers that penetrate into the faulted region to a certain depth. The boundary layer can be modelled as an array of dislocations of alternating sign, thus the misfit energy density can be derived from standard dislocation models [37] in the form where L 1 is the size of the container, C is a material constant proportional to the elastic modulus of the matrix, and the coefficient L 0 (with the dimension of a length) plays the role of a core cutoff to bound the lowest scales. The misfit boundary energy is minimized by narrow boundary layers, which are related to small values of L, in competing demand with the cohesive energy, which is minimized by large L. By assuming L as unknown and minimizing the incremental work of deformation with respect to L, the optimization procedure provides also the value of L. Remarkably, for the linear decreasing cohesive law, the optimal fault separation L is found to be independent of ∆ and given by In the presence of multiple fault families, the model can be employed in a recursive manner, by replacing the elastic behavior of the matrix with a new brittle damage behavior. The number of nested fault families defines the rank of the microstructures.

Permeability and porosity
The simple fracture pattern that characterizes the model provides an analytical expression of the porosity and the permeability of the damaged material. The change in porosity due to the fault activity n f can be computed as the first invariant of the fault small strain tensor Thus the total porosity of the material is obtained by adding the fault porosity to the matrix porosity n m , which is assumed to be related to the matrix volumetric strain An additive decomposition holds also for the permeability tensor k where k m is the permeability of the matrix, assumed to be isotropic and dependent on matrix porosity via a Kozeny-Carman type relationship, and k f is the permeability due to the fault activity [15] k where I is the identity tensor and C KC is a material constant. The permeability of a rank-Q material can be estimated as the sum of the contributions of the single fault families, as

Numerical solution
The solution of the coupled problem is achieved through finite elements (FE), using a in-house developed code. The finite element discretization is applied to the weak form of the governing differential equations. The linear momentum balance (1) accounting for the Terzaghi principle reads where s is a test function field and ∇s the corresponding gradient. The weak form of the continuity equation (2) reads where η is test function and ∇η the corresponding gradient. We use a different FE discretization for the displacement field and the pressure field. For the displacements, we use tetrahedral elements with quadratic interpolation for displacement u and test function s, and linear interpolation for p and test function η. The matrix form of the discretized problem is where P collects the nodal pressures, U the nodal displacement components,U the nodal velocity components, F ext (t) the external forces, F int (U) the internal forces, Q ext (t) the imposed nodal hydraulic fluxes. Moreover K is permeability matrix and H the coupling matrix.
In the solution algorithm, coupling is enforced weakly, by solving separately the time discretized equations. For the solution of the continuity equation we adopt a direct solver, while for the solution of the linear momentum balance equation we use an explicit nonlinear solver based on the concept of dynamic relaxation. The code is parallelized with a shared memory algorithm.

Numerical Applications
The brittle damage model has been conceived to describe degradation in geomaterials as a consequence of the modification of the microstructure, due to the onset of diffused cracking with a specific geometry. As such, the model is not able to simulate classic problems of linear elastic fracture mechanics, since it does not feature the presence of an isolated sharp crack tip and it is not possible to monitor the progressive separation of two fracture flanks. Moreover, it has been mathematically and numerically proved [37] that the model is particularly efficient in modelling failure in compression. For these reasons, it is not significant to validate the model against the propagation of single sharp fractures. The validation of the model has been done by simulating a small scale hydraulic fracture laboratory test [4]. The experimental test has been performed on a concrete specimen in controlled, properly applied, boundary and initial conditions. As expected in experiments conducted on opaque geo-materials (e.g., rock or concrete), authors in [4] reported the time history of the measured borehole pressure and the post-mortem diffused fracture pattern. The material was characterized for the elastic properties, while the missing strength properties used in the numerical simulations have been taken from the literature.
The results of the numerical analyses are reported in terms of two significant quantities: the fracture volume V f , and the fracture surface S f . The fracture volume V f is defined as the volume of the opened faults, i. e., where N f is the number of elements that experienced failure, N g the number of integration points of the element e, ∆ N e g and L e g are the normal opening displacement and the spacing, respectively, relevant to the g-th integration point of the element e. Clearly, the fracture volume depends on the opening of the faults and it decreases upon dissipation of the pore pressure. The total fracture surface S f is the sum of the fracture surfaces of the N f element that experienced failure:  Table 1 Laboratory test [4]. Material parameters used in the numerical analysis. The star (*) denotes the parameters estimated from other literature sources.
In all the subsequent calculations we blocked the formation of nested microstructures and allowed the material to develop only one family of faults.

Reference small scale laboratory test
As reported in [4], The solid model of the specimen is discretized into 6,280 10-noded tetrahedral elements and 8,945 nodes. To satisfy the regularity inf-sup conditions, in each element the interpolation functions for the displacement field are quadratic, while for the pressure field are linear. For the sake of simplicity, the borehole is not explicitly modelled and thus the effects of stress concentration are neglected, see Fig. 2b. The discretization is finer in proximity to the cavity, where the pressurized fluid is localized and stress gradients are expected to be higher. The material parameters used for the numerical simulation are listed in Table 3.1. The experimental paper [4] did not provide all the material properties needed for the numerical simulation, therefore the missing data have been estimated from the literature.
We began with a preliminary numerical simulation where the experimental injection flux is applied to the nodes of the mesh that fall within the volume of the cavity. The corresponding pressure of the cavity node is recorded and compared with the experimental data, see Fig. 3.
The predictive ability of the model to describe the pressure history is remarkable. In particular, although faults begin forming at times around 400 s from the beginning of the pressurization, a drop in the pressure is observed only when the faults undergo a sufficient wide opening. Fig. 4 compares the distribution of the damage due to the pressurization obtained in the experiments with the one provided by the numerical simulations. As expected, the damaged zone extends widely in the plane normal to the direction of the minimum compressive stress, Fig. 4a. The extension of the damaged zone in the numerical simulation is visualized in terms of a scalar that describes the extension of the the damaged zone on a quarter of the numerical model that allows to see the different extension of damage on the two symmetry planes, Fig. 4b, and is in agreement with the experimental distribution. Fig. 5 shows the time evolution of fracture surface S f and fracture volume V f . Fig. 5a shows a progressively increasing curve. Clearly, the fracture process is characterized by irreversibility and the extension of the fracture surfaces cannot reduce when the fluid pressure starts to reduce. Indeed, the total fracture surface can extend also when the pressure peak reduce, since the fluid can flow to other regions and induce the formation of additional fractures. This particular aspect of hydraulic fracturing has been verified by means of acoustic emission measurements [49]. Contrariwise, the fracture volume is related to the level of the opening of the faults, strongly influenced by the fluid pressure, and therefore it reduces when the fluid pressure goes down.
We conducted a second analysis by applying the fluid pressure recorded during the experiment, and obtained very similar results in shorter times. Therefore, given the pure speculative and not applicative nature of this study, for the sake of computational cost, the subsequent analyses have been conducted by applying the fluid pressure history.

Sensitivity analysis
In order to assess the influence of mechanical parameters on the mechanism of failure and on the hydraulic fracture efficiency, we have conducted a sensitivity analysis on all the parameters of the model. The study showed that the response of the system in terms of evolution of fracture volume and fracture surfaces is affected only by a few parameters: the material strength (in terms of both the tensile strength T c and the friction angle ϕ), and the initial permeability of the isotropic matrix k m = kI.

Tensile strength
Three different values of tensile strength have been considered: T c =3, 4, and 5 MPa. In terms of Mohr-Coulomb failure criterion, the change in the tensile strength causes the translation of the failure surface according to the σ − τ plots in Fig. 6.
A few significant results of the numerical analysis are collected in Table 2.
The time evolution of fracture surface and fracture volume for these values are shown in Fig. 7. The most evident effect of the increase of T c is the increase in the pressure necessary to activate the first fracture p 1 (cf. see Table 2) and the corresponding delay in the time of first failure, as well as the reduction of the initial extension of the fracture surface. The time and the extension of the initial fracture surface are visualized with black diamonds in Fig. 7b. The fact that higher tensile strengths reduce the extension of the initial fracture surface is related to the higher engagement   requested to the material to attain the strength in the neighborhood of the damaged zone.
Interestingly, the formation of faults does not imply an immediate normal opening of the faults and the fracture volume may be initially null. Faults in general begin to open at a higher pressure p N , see Table 2. The time and the extension of the initial fracture surface are visualized with black diamonds in Fig. 7d. Data in Table 2 show that tensile strength of the material affects significantly the maximum value of the fracture surface and fracture volume extensions.
Note that these analyses have been conducted at a fixed value of the critical energy release rate G c . In our material model the tensile strength is strictly related to the critical opening displacement ∆ c , which is also related to the selection of the scaling parameter L. Therefore a change in T c modifies also the kinematic aspects of the model.

Friction angle
A change in the friction angle does not affect the values of the critical opening displacement ∆ c . We study the effects of friction angles ϕ =30, 38, and 45 degrees, which in the Mohr-Coulomb criterion corresponds to three failure surfaces with different slope, see Fig. 8.
The results of the sensitivity analysis for the friction angle are collected in Table 3 and visualized in Fig. 9. From the data in Table 3 we observe that lower (higher) friction angles induce higher (lower) fracture surfaces and lower (higher) fracture volumes. This behavior can be justifed by the fact that for higher level of friction the sliding between the crack flanks is less prone to occur, and the preferred failure mechanism is the normal opening. Moreover, Fig. 9 show that higher values of friction angle cause a time delay in the formation of the first set of faults, due to a higher fluid pressure necessary to bring the material to failure, and a higher extension of the fracture surface at the fracture onset.

Matrix permeability
At the beginning of the test, the permeability of the block is given by the initial permeability of the porous matrix. Upon fault formation, the permeability of the material increases as a function of fault opening and spacing, Eq. (21). A high initial permeability may facilitate the flow of the fluid in the medium and reduce the local effective stress as so as to anticipate the time of the formation of the first fault set. We performed three analysis using the permeability coefficients k = 2 10 −14 , 2 10 −16 , 2 10 −18 m 2 .
The results of the analyses are illustrated in Fig. 10. Interestingly, simulations show that the initial permeability does not affect the timing of first failure, but it does affect the initial extension of the fracture surface and in general the evolution of both fracture surface and fracture volume.
The lack of monotonicity between initial permeability and evolution of the extension of the fracture surface and fracture volume can be explained as follows. Fracture is the consequence of two concurrent physical phenomena. On one side, a low permeability material does not allow a fast fluid pressure diffusion, therefore high pressures concentrate in a limited region in the neighborhood of the injection point, leading to a concentrated damage for the material. On the other side a high permeability material facilitates the diffusion of the fluid pressure on a larger volume and the reduction of the pressure peaks, limiting the extension of the damage.

Field application
We use the code to simulate a large scale problem, i. e., the artificial increase of the permeability of a low-permeability hydrocarbon reservoir by means of a sequence of hydraulic fracturing processes. We consider a very simplified geometry and homogeneity for the material.
We consider a parallelepiped domain in a shale reservoir of size 3.5 km × 2.15 km × 2.15 km, subject to geostatic loads σ x = 40 MPa, σ y = 50 MPa, and σ z = 60 MPa, see Fig. 11a. The initial pore pressure of the repository is assumed to be uniform and equal to p = 24 MPa. The domain is ideally crossed by a horizontal borehole parallel to the x-axis (not reproduced in the model).
Note that in field application, where hydraulic fracturing is used to create networks of artificial fractures in deep rock embankments to facilitate the extraction of gas/oil, the re-closure of cracks is an undesired event. To maintain the level of opening of the cracks, granular additives named proppant are mixed to the fracturing fluid. Proppant acts as a constraint to the reclosure of the crack flanks and it is widely used in petroleum engineering. We simulate the effect of nine sequential fluid injections, located along the borehole at different sites distant 300 m one from the other, by applying the 600 s fluid pressure history shown Fig. 11b. The domain is discretized in 51,066 10-noded tetrahedral elements, with 69,227 displacement nodes, see Fig. 11c.
The physical properties adopted for the large scale reservoir simulation are summarized in Table 4. The results of the numerical analyses are shown in Fig. 12a which visualizes the extension of the damage induced by the hydraulic fracturing in terms of magnitude of the opening displacement. Fig. 12b shows the time evolution of fractured surface and fractured volume. Both quantities progressively increase in time, following the sequence of fluid injection. The presence of proppant in the injection fluid allows the fractured volume to preserve high values also after the dissipation of the fluid pressures.

Conclusions
We have used a coupled porous-mechanical finite element code, including a brittle damage material model [37,15,16], to the simulation of a hydraulic fracturing laboratory test [4] and to the simulation of an extensive hydraulic fracturing procedure in a hypothetical hydrocarbon reservoir.
We begin by validating the code against the experimental results reported in [4]. For the simulation it is necessary to complete the experimental data with additional material parameters that are taken from the literature. Therefore, after the simulation we conduct a sensitivity analysis on the parameters not provided by the experimental paper (tensile strength T c , friction angle ϕ, scaling ratio L 0 /∆ c , and isotropic matrix permeability). The sensitivity analysis allows to understand the influence of the various parameters on the results of a boundary value problem.
We conclude the study by applying the code to the simulation of a large scale problem, i. e., the simulation of an extensive hydraulic fracturing process in a gas/oil reservoir. The simulation is able to provide important information on two parameters that are usually taken in consideration to estimate the production of a reservoir: the extension of the fracture surface and of the fracture volume.
The simulation is extremely simplified and dwarfs the complexity of the real problem. First, the actual geological information is disregarded. Second, the actual technology requires the drilling of the borehole, which causes a disturb on the geostatic stress state, and the pre-fracturing of the hydraulic fracturing sites through small explosions, which further modifies the original stress state and permeability. Third, the process control quantity is the amount of fluid injected, while the fluid pressure is an output parameter registered at the surface end of the borehole and not the actual (unknown) pressure. To determine the pressure at the hydraulic fracturing location, a complex fluid-dynamical problem involving a mixture of additives and water should be solved.
On the other side, the oversimplification of the approach leads to a numerical efficiency that allows to obtain in a reasonable computational time an reliable estimate of the hydraulic fracturing procedure. The code can be considered as a supporting design tool for technical applications in petroleum engineering.

Consent for publication
Authors provide the consent for the publication of the information included in the manuscript.

Availability of data and material
Data (Geometry and mesh used in the simulations) are available upon request.

Competing interests
Authors have no competing interests on this topic.

Funding
None.