Three-Dimensional poroelastic effects during hydraulic fracturing in permeable rocks

A fully coupled three-dimensional ﬁnite-element model for hydraulic fractures in permeable rocks is pre- sented, and used to investigate the ranges of applicability of the classical analytical solutions that are known to be valid in limiting cases. This model simultaneously accounts for ﬂuid ﬂow within the frac- ture and rock matrix, poroelastic deformation, propagation of the fractures, and ﬂuid leakage into the rock formation. The model is validated against available asymptotic analytical solutions for penny-shaped fractures, in the viscosity-dominated, toughness-dominated, storage-dominated, and leakoff-dominated regimes. However, for intermediate regimes, these analytical solutions cannot be used to predict the key hydraulic fracturing variables, i.e. injection pressure, fracture aperture, and length. For leakoff-dominated cases in permeable rocks, the asymptotic solutions fail to accurately predict the lower-bound for fracture radius and apertures, and the upper-bound for fracture pressure. This is due to the poroelastic effects in the dilated rock matrix, as well as due to the multi-dimensional ﬂow within matrix, which in many simulation codes is idealised as being one-dimensional, normal to the fracture plane.


Introduction
Hydraulic fracturing is the process by which one or more fractures are propagated into a rock formation, driven by the internal flow of a pressurised fluid. While fluid-driven fracturing can occur naturally, it is most often studied within the context of the engineering process of injecting fracturing fluid into a reservoir rock, with the aim of increasing well productivity ( Adachi et al., 2007;Bazant et al., 2014 ). Although the hydraulic fracturing process is currently often thought of in the context of shale gas reservoirs, in current industry practice, almost all oil and gas wells are hydraulically fractured ( Economides and Nolte, 20 0 0 ). Hydraulic fracturing is a complex, multi-physics, multi-dimensional problem, which requires robust models that can simultaneously account for matrix and fracture deformation, fluid flow through the matrix and fractures, fluid exchange between fractures and matrix, and fracture propagation and interaction, all in a fully-coupled, threedimensional setting.
Hydraulic fracturing protocols are designed to control the fracture's surface area and aperture distribution, and also aim to con-trol injection pressure, and the dependence of these variables on fracturing fluid rheology, injection rate, and the hydro-mechanical properties of the rock ( Detournay and Peirce, 2014 ). Analytical and semi-analytical solutions have been developed to quantify hydraulic fracturing variables of interest, such as injection pressure, fracture aperture, and fracture length ( cf . Adachi et al., 2007 ). These solutions provide the foundation for hydraulic fracturing design ( e.g. Cleary, 1980;Cleary et al., 1988 ). These solutions are constructed by combining the equations for laminar flow through the fracture, with the equations for elastic deformation of the adjacent rock. Fluid flow through the fracture is commonly modelled using lubrication theory, which is derived from the general Navier-Stokes equation for flow of a fluid between two parallel plates ( Batchelor, 1967;Zimmerman and Bodvarsson, 1996 ), whereas the fracture aperture is calculated using linear elasticity in conjunction with Linear Elastic Fracture Mechanics (LEFM) to compute the mode I stress intensity factor at the fracture tip ( Geertsma and de Klerk, 1969;Spence and Sharp, 1985 ).
Based on the energy-dissipation mechanism, fracture propagation regimes can be classified as viscosity-dominated or toughnessdominated ( Detournay, 2004 ). In the viscosity-dominated regime, energy dissipation is dominated by the flow of the viscous fluid, whereas in the toughness-dominated regime, energy dissipated is dominated by the creation of new fracture surfaces at the fracture tip. Based on the ability of the rock matrix to dissipate fracturing fluid, two other extremes can be defined: storage-dominated , in which the injected fluid remains mainly inside the fracture, and leakoff-dominated , in which most of the injected fluid dissipates into the surrounding medium. The four resulting combined asymptotic regimes are therefore storage-viscosity, storagetoughness, leakoff-viscosity , and leakoff-toughness ( Garagash et al., 2011 ). Asymptotic solutions that are valid at the end-members of the parameter space provide a fundamental understanding of the hydraulic fracturing process, and provide benchmarking cornerstones for numerical models. However, existing analytical solutions are restricted to simplified fracture geometries in homogeneous rock masses, and are typically constrained to a set of fixed boundary conditions. Standard geometries include the PKN fracture ( Perkins and Kern, 1961;Nordgren, 1972;Mathias and van Reeuwijk, 2009 ), the KGD fracture ( Geertsma and de Klerk, 1969;Spence and Sharp, 1985;Adachi and Detournay, 2008 ), and radial (penny-shaped) fractures ( Savitski and Detournay, 2002;Bunger et al., 2005;Kovalyshen, 2010 ). Moreover, analytical solutions do not exist for cases that are not at the corners of this parameter space.
Numerical models that attempt to simulate hydraulic fracturing include the boundary integral method ( Peirce and Siebrits, 2001 ), the boundary element method ( Simpson and Trevelyan, 2011 ), the distinct element method ( Marina et al., 2014 ), the finite element method ( Carrier and Granet, 2012 ), discrete fracture network ( Fu et al., 2013 ), the embedded fracture model ( Norbeck et al., 2015 ), the lattice approach ( Grassl et al., 2015 ) and the extended finite element method ( Dahi-Taleghani, 2009;Mohammadnejad and Khoei, 2013;Salimzadeh and Khalili, 2015a ). However, in the majority of available models, flow through the rock matrix, and fluid exchange between fracture and rock matrix, are either ignored by assuming an impermeable rock formation ( e.g. Dahi-Taleghani, 2009 ), or simplified by using a one-dimensional analytical leakoff model ( e.g. Zhou et al., 2015 ). Substantial field evidence has proven the impermeable matrix assumption to be an unrealistic assumption ( Economides and Nolte, 20 0 0; Adachi et al.,20 07 ). In onedimensional leakoff models ( Carter, 1957 ), fracture-to-matrix flow is represented as a sink term in the mass balance equation for fracture flow. This approach has several shortcomings, such as the assumption of one-dimensional flow, time-dependency of flow instead of pressure-dependency, and more importantly, this approach cannot model matrix dilation. Although flow from the fracture into the rock matrix is by definition locally one-dimensional at the fracture wall, where the flux vector must be normal to the fracture wall, in a global sense it is three-dimensional, unless the permeability in the direction normal to the fracture plane is significantly higher than in other directions ( Hagoort et al., 1980 ). As time elapses, the leakoff rate predicted by Carter's model, at each position along the fracture, decreases proportionally to square-root of time; consequently, a scenario of fracture arrest is not possible ( Mathias and van Reeuwijk, 2009 ). Finally, this model does not account for the fact that seepage of the fracturing fluid into the rock formation increases the fluid pressure in the matrix, causing dilation of the rock matrix. A dilated matrix applies stresses back onto the fracture, referred to as 'back-stresses' in the hydraulic fracturing literature, which tend to close the fracture ( Kovalyshen, 2010 ). These factors also affect the available semi-analytical solutions for leakoff-dominated regimes that use a simplified onedimensional leakoff model in their formulation. These solutions therefore fail to accurately predict hydraulic fracturing parameters in leakoff-dominated regimes, as shown by Carrier and Granet (2012) , and Salimzadeh and Khalili (2015a) for single-phase flow, and by Salimzadeh and Khalili (2015b) for two-phase flow in two dimensions, as well as in the present study for three dimensions.
In addition to poroelastic effects due to the aforementioned back stress phenomenon, there is a further environmental consequence of fluid seepage through the rock matrix, as it may promote the possible migration of injected fluid towards drinking water aquifers ( Birdsell et al., 2015 ). Therefore, robust modelling of matrix flow is essential for both hydraulic fracture engineering and environmental aspects of subsurface fracturing. Only a few attempts have been made to incorporate flow in the rock matrix, coupled to mechanical deformation and flow in fracture. Rethore et al. ( 2008 ), Mohammadnejad and Khoei (2013) and Khoei et al. (2014) , using the extended finite element method, introduced enriched pressures at the fracture to capture the discontinuous flow velocity at the fracture boundary. However, the enriched pressure represents the fluid pressure in the rock matrix near the fracture, and does not represent the pressure inside the fracture. Therefore, when coupled with mechanical deformation, the enriched pressure will be scaled by the Biot coefficient, whereas the fracture pressure actually does not require such scaling. Carrier and Granet (2012) introduced independent flow through the fracture and the rock matrix into their hydraulic fracture model. Their model was a combination of zero-thickness elements for the propagating fracture, and conventional bulk finite elements with a cohesive zone model. The equality of pressure between fracture and matrix at the fracture walls was enforced in the numerical model using Lagrange multipliers. Salimzadeh and Khalili (2015a, b ) proposed an XFEM model that included two independent flow models in the fracture and the rock matrix, with a leakoff mass transfer between fracture and rock matrix to link the two. The leakoff depends on the pressure gradient in the matrix adjacent to the fracture, as well as on the fluid viscosity and matrix permeability. Norbeck et al. (2015) , using an embedded fracture model, also considered two flow domains for matrix and fracture in two dimensions, and linked them through a similar mass transfer term.
A three-dimensional fully coupled finite element model for hydraulic fracturing is presented in the present paper, validated against known analytical solutions, and subsequently applied to study the influence of fluid exchange between fracture and matrix on fracturing. In particular, 3D diffusion and its related poroelastic effects on the propagation of fractures are investigated. The present model accounts for fluid flow within fracture and matrix, the propagation of the fracture, and fluid leakage into the formation rock. Fluid flow through the permeable rock matrix is modelled using Darcy's law, and is coupled with laminar flow within the fracture. Fracture growth and the direction of growth are estimated using an energy-based criterion that is based on the modal stress intensity factors along the fracture tip ( Paluszny and Zimmerman, 2013 ). This model is validated against available asymptotic solutions for penny-shaped hydraulic fractures. Fifteen cases with varying fluid and rock matrix properties are run, to investigate the impact of fluid and rock matrix properties on the leakoff and fracturing. Numerical simulations conducted over a range of parameter values delineate the limits of validity of the various available asymptotic solutions.

Computational model
Fractures are represented discretely using two-dimensional surfaces embedded in a three-dimensional domain. When deriving the governing equations, each fracture is represented by a discontinuity c in the domain with boundary , as shown in Fig. 1 . The fully coupled model is constructed on three separate yet interacting sub-models, including models for mechanical deformation, fracture flow, and matrix flow. The solid matrix is assumed to be linear elastic, homogeneous, and isotropic, with flow modelled using Darcy's law. An independent fracture flow model is developed based on lubrication theory. The mechanical and fracture flow models are coupled through hydraulic loading on the fracture walls, and by ensuring the compatibility of fracture volumetric strains. The mechanical model is coupled to matrix flow through the effective stress concept, and finally, the fracture flow and matrix flow models are linked to each other through the leakoff mass transfer term. Tension and compression are assumed positive for stresses and pressures, respectively.

Mechanical deformation model
The mechanical deformation model is based on the condition of stress equilibrium for a representative elementary volume of the porous medium. For quasi-static conditions, the linear momentum balance equation for this elementary volume may be written as div σ + F = 0 (1) where F is the body force per unit volume, and σ is the total stress.
The effective stress is defined as the function of total stress and matrix pressure, and controls the mechanical deformation. It is defined exclusively within the rock matrix, linking a change in stress to the change in strain. The effective stress for the rock matrix saturated with a single-phase fluid is defined as ( Biot, 1941 ) where σ is the effective stress, α is the Biot coefficient, p m is the matrix fluid pressure and I is the second-order identity tensor. The Biot coefficient is defined as where K and K S are the bulk modulus of the porous rock and of the rock matrix material ( e.g. mineral grains), respectively ( Zimmerman, 20 0 0 ). The stress and strain relationship of the element is expressed as in which D is the drained stiffness matrix, and ε is the strain tensor in the porous medium. Assuming infinitesimal deformations, strain is related to displacement by where u denotes the displacement vector in the porous medium. Hydraulic loading on the fracture walls is applied as boundary traction, as shown in Fig where p f is the fracture pressure, and n C is the outward unit normal to the fracture wall (on both sides of the fracture). Integrating Eq. (1) over the element, and after some manipulation, the differential equation describing the deformation field for a saturated rock matrix is given by

Fracture flow model
An independent fracture flow model is considered for the hydraulic fractures. This model allows direct computation of fracture fluid pressure, and implicit application of hydraulic pressures onto fracture walls. The objective is to realistically represent fracture flow, instead of smearing it with the flow through the neighbouring matrix. Assuming a planar fracture, in which the area of the fracture plane is much larger than the fracture aperture, the average velocity of fluid along the fracture plane may be calculated using the cubic law as ( Witherspoon et al., 1980 ) where a f is the fracture aperture, μ f is the fluid viscosity, and p f is the fracture fluid pressure. The aperture is given by the differential displacement between two sides of the fracture, a f = ( u + − u − ). n c , where u + and u − are the displacements of the two opposing faces of the fracture. The mass balance equation for a slightly compressible fracture fluid may therefore be written as in which ρ f is the fluid density, and L f is the leakoff flow from the fracture to the matrix. This leakoff leads to mass transfer coupling between the fracture flow and rock matrix flow. Assuming that the fracture fluid is Newtonian, the leakoff flow per unit area of the fracture wall can be written, using Darcy's law, as where k m is the intrinsic permeability of the rock matrix. Substituting the fluid velocity into the mass balance equation, and after some manipulation, it is found that div a f where c f is the fluid compressibility. Note that the term . n c / ∂ t provides direct coupling between the displacement field and the fracture flow field, which is symmetric to the fracture pressure loading term, p f n C . For the case of onedimensional incompressible flow with no leakoff, c f = 0 and L f = 0, Eq. (11) reduces to the lubrication equation ( Batchelor, 1967 ),

Matrix flow model
The matrix flow model that represents flow through the porous matrix is constructed by combining Darcy's law with mass conservation for the fluid. Neglecting inertial and viscous effects, Darcy's law for matrix flow may be written as where v r is the relative velocity vector of the fluid in the matrix, k m is the intrinsic permeability of the rock matrix, μ f is the fluid viscosity, g is the vector of gravitational acceleration, and p m is the matrix fluid pressure. The relative velocity of the fluid with respect to the deforming rock matrix is given by where v s is the rock matrix velocity, defined as v s = ∂u ∂t (15) where u is the displacement vector of the rock matrix. The mass balance equation for the fluid in the rock matrix may be written as where ρ f is the fluid density, φ is the rock matrix porosity, and v m is the matrix fluid velocity. Note that the leakoff only occurs on the boundary of the volume element that is connected to a fracture ( c ). Therefore, a Dirac delta function is applied, where x c represents the position of the fracture. Integrating over the element and after some manipulation, the governing equation for the flow model may be expressed as div where c f is the fluid compressibility. The Biot coefficient α appears in Eqs. (7) and ( 17 ), whereas it does not appear in the fracture flow model ( Eq. 11 ), as the fracture itself is not a "porous medium".
Setting α = 0 will decouple the mechanical deformation model and the matrix flow model, in which case mechanical loading will have no direct effect on the matrix pressures, and vice versa . In contrast, fracture pressures will always be coupled to the mechanical deformation model, irrespective of the value of the Biot coefficient.

Finite element approximation
The governing equations are solved numerically using the finite element method. Spatial and temporal discretisation are accomplished using the Galerkin method and finite difference techniques, respectively. Displacements (three displacements for three dimensions) and fluid pressures (fracture and matrix) are defined as the primary variables. Using the standard Galerkin method, the displacements and pressures within an element may be approximated from the nodal values as where N and N c are the standard shape functions vector for volume and surface elements, respectively. ˆ u , ˆ p m andˆ p f are vectors of nodal values of displacement, matrix pressure, and fracture pressure, respectively. Fracture pressures are only defined for the nodes on the fractures.
Using the finite difference technique, the set of discretised equations can be written as where K is the mechanical stiffness matrix, C f and C m are hydromechanical and poroelastic coupling matrices, respectively, H is the conductance matrix, M is the capacitance mass matrix, L is the leakoff mass matrix, F is the applied load vector, Q is the fluid flux, and ˆ u and ˆ p are the vectors of nodal values of displacement and fluid pressure, respectively. [ B 1 ] 6 ×3 n = ∇ N , [ B 2 ] 1 ×3 n = δ T B 1 , and, [ B 3 ] 3 ×n = ∇N are derivatives of the shape function, δ = { 1 1 1 0 0 0 } T , and ∇ is the gradient operator. Superscript t represents the time at the current step, superscript t + dt represents time at the next step, and dt is the time step. The non-diagonal components of the stiffness matrix are populated with the coupling matrices C f for hydro-mechanical coupling, C p for poroelastic coupling, and L for fracture-matrix flow coupling. The operator ∇ for three-dimensional displacement field is defined as The components of the stiffness matrix depend on the primary unknown variables, i.e. permeability of the fracture depends on the fracture aperture; therefore, an iterative procedure is required to reach the correct solution within acceptable tolerance. The discretised coupled equations are implemented as part of the Imperial College Geomechanics toolkit ( Paluszny and Zimmerman, 2011 ), which interacts with an octree volumetric mesher and the Complex Systems Modelling Platform (CSMP ++ , also known as CSP), an object-oriented application programme interface (API), for the simulation of complex geological processes and their interactions (formerly CSP, cf. Matthäi et al., 2001 ). The set of linear algebraic equations are solved with the algebraic multigrid method for systems, SAMG ( Stüben, 2001 ). Two types of discretisation are used: quadratic tetrahedra for volume elements, and quadratic triangles for surface elements (fractures). The triangles on two opposite surfaces of a fracture are matched with each other, but they don't share nodes, and duplicate nodes are defined for two sides of a fracture. The triangles are matched with faces of the tetrahedra connected to the fractures. Therefore, they share the same nodes. However, the model presented in this study can also be applied to non-matching elements. Fracture flow equation ( Eq. 17 ) is solved only on one-side of fracture ( i.e. matrices H f and M f are accumulated over triangle elements on one side of the fracture); however, the coupling matrices ( C f and L ) are accumulated on both sides of the fracture. Mechanical deformation and matrix flow equations are accumulated over the volume elements.

Stress intensity factors and growth model
The mechanical deformation of the rock matrix leads to concentrations of stress around the fracture tips, which can be quantified locally at each tip by the stress intensity factors (SIFs). The SIFs are key parameters in evaluating and predicting fracture growth, and take into account the effects of fluid pressure and rock properties on the growth of the fracture. The state of stress immediately ahead of the fracture front is known to be singular. Therefore, in contrast to conventional linear polynomial interpolation, quadratic elements are used so as to better approximate the stress tip singularity ( Nejati et al., 2015a ). Two methods for the SIF extraction from the FE solution can be employed. Direct approaches, based on the correlation of the displacements over the crack surface are simple, straightforward, and computationally inexpensive, but require very refined meshes around the crack front in order to yield low approximation errors. Alternatively, energy-based methods that integrate stresses over the region ahead of the crack tip are less prone to numerical error, and yield better approximations, with significantly coarser meshes. Three stress intensity factors for the three modes of fracture opening are computed by computing the energy-based interaction integral ( Yau et al., 1980 ), modified in the present paper for poroelastic media by using the effective stresses in place of the usual stresses, over a set of virtual disk domains distributed along the fracture tip ( Nejati et al., 2015b ). The three SIFs are K I for opening due to tensile loading, K II for in-plane shearing due to sliding, and K III for out-of-plane shearing due to tearing. The crack grows once the equivalent stress intensity factor K Ieq , overcomes the material toughness ( k ic ). The equivalent SIF in the direction of propagation ( θ p ), is calculated as ( Schöllmann et al., 2002;Paluszny and Zimmerman, 2013 ) where K cs = K I cos 2 ( θp 2 ) − 3 2 K II sin ( θ p ) , and θ p is the propagation angle. By setting K III = 0 in the above equation, the equivalent stress intensity factor for two-dimensional space are recovered ( Erdogan and Sih, 1963 ) The propagation angle ( θ p ) is determined using a modified maximum circumferential stress method that takes into account modal stress intensity factors. The equivalent stress intensity factors are computed locally at 100 locations along the fracture front, i.e . tips, and are used to determine if the fracture will or will not advance. Fractures are extended by deforming their geometry, and the mesh is regenerated and optimised at every growth step ( Paluszny and Zimmerman, 2011 ).

Simulation results: penny-shaped hydraulic fracture
Fluid is injected through a horizontal well that perforates the centre of a vertical penny-shaped fracture. The size of the well is assumed to be negligible with respect to the size of the fracture, and so the wellbore is modelled as a point source boundary condition in the simulations. Asymptotic solutions available for this geometric case, under storage regimes (viscosity-storage and toughness-storage), are used to validate the presented numerical model. Further simulations, in which leakoff is modelled by considering a permeable matrix, are performed to investigate the limits of validity of asymptotic solutions under leakoff regime, and to highlight the effects of poroelasticity and the three-dimensionality of flow within the matrix.
A single penny-shaped fracture of initial radius 1 m is located in the centre of a 90 × 90 × 60 m model. The model is spatially discretised using 17,441 tetrahedra and triangles. Convergence is achieved using an average of four iterations to reach a tolerance of 1%. Fracturing fluid is injected at a constant rate of Q = 0.01 m 3 /s into the centre of the fracture. A total of fifteen cases are simulated (seven cases in the viscosity regime, and eight cases in the toughness regime), with varying fluid and matrix properties, including four extreme regimes as well as intermediate cases. The properties for these cases (cases 1-15) are defined in Table 1 . The Young's modulus and Poisson's ratio for all cases are set to 17 GPa and 0.2, respectively. The minimum in situ stress acting in the direction normal to the fracture plane for each case is given in Table 1 , and the in situ stresses in the other two directions are set to 1 × 10 7 Pa. Savitski and Detournay (2002) derived solutions for a pennyshaped hydraulic fracture in an impermeable elastic rock. They defined three hydraulic fracturing variables, fracture aperture a f , fracture net pressure p fnet , and fracture radius R f , as functions of the dimensionless parameters , and γ , as: where ε is a small number, L is a length scale, E = E /(1 − ν 2 ) is the plane-strain elastic modulus, E is the Young's modulus, and ν is Poisson's ratio, respectively. Fracture net pressure is defined as the difference between fracture pressure and the in situ stress normal to the fracture plane. The solution was given for viscositydominated and toughness-dominated regimes for storage cases (impermeable rock matrix). It is worth mentioning that in the viscosity-dominated regime, the rock toughness has a negligible effect on hydraulic fracture growth, whereas in the toughnessdominated regime the influence of the fluid viscosity is negligible. In this latter case, the fracture pressure is essentially uniform within the fracture, and the fracture's mode I stress intensity factor is equal to the fracture toughness, K I = K ic . To distinguish between viscosity and toughness regimes, Savitski and Detournay (2002) defined a dimensionless viscosity as Cases in which M 1 are viscosity-dominated, whereas M 1 represents toughness-dominated cases. In this equation, μ = 12 μ f , and K = 4(2/ π ) 1/2 K ic . Note that the dimensionless viscosity is time-dependent, and so the behaviour moves from viscosity-dominated towards toughness-dominated as time increases.

Viscosity-dominated regime
The fluid viscosity and fracture toughness values for viscositydominated cases are set to μ f = 0.1 Pa s, and K ic = 1 × 10 6 Pa m 0.5 , respectively. For storage-dominated cases, the permeability of the matrix is assumed to be negligible ( k m = 0), whereas for leakoff dominated cases the matrix permeability is taken to be 1 × 10 −12 m 2 . The Biot coefficient α is varied between 0 and 1, corresponding to uncoupled and fully coupled poroelasticity scenarios, respectively, although 0 is an unrealistic value, as the Biot coefficient is bounded below by the porosity ( Zimmerman, 20 0 0 ). For the given parameters, the dimensionless viscosity is M = 39, which corresponds to a viscosity-dominated case ( M 1). The simulation results for injection pressure, fracture aperture at the well, and fracture length (radius) for viscosity-dominated cases 1-7 are shown in Fig. 2 . Included in these figures are the asymptotic solutions for reference cases. Cases 1-3 with zero matrix permeability correspond to the storage-dominated regime, and cases 5-7 with matrix permeability 1 × 10 −12 m 2 , correspond to leakoff-dominated regime. Case 4 is an intermediate case. The simulation time for leakoff dominated cases 5-7 is increased to 200 s. Fig. 2 a shows the net injection pressure versus time when injecting a viscous fluid. Good agreement is found between the analytical solution for the storage case and the present model with an impermeable matrix (case 1). The growth increment in case 1 is 0.5 m, and it has been changed to 1.0 and 0.25 m in cases 2 and 3, respectively. Results show a negligible effect due to the growth increment. When leakoff is allowed by assuming a permeable rock matrix in cases 4-7, the injection pressure increases, compared with the no leakoff cases 1-3. The leakoff is increasing by increasing the permeability, resulting in higher injection pressure, lower fracture aperture, and lower fracture radius. In case 5, the permeability in the matrix is assumed to be "one-dimensional", such that the permeability in direction normal to the fracture is 10 0 0 times higher than the other two directions (one-dimensional leakoff). This case is used for validation against the analytical solution. A very good match is observed between present model results and the analytical solution for the fracture radius for case 5. In case 6, the permeability in all three directions is assumed equal to 10 −12 m 2 (three-dimensional leakoff). Leakoff is also increased by allowing the flow within the matrix to be globally three-dimensional, rather than one-dimensional. Increasing leakoff, again, leads to higher injection pressure, lower fracture aperture, and lower fracture radius. In case 7, the Biot coefficient is increased to 1. Injection pressure increases substantially with increasing Biot coefficient. This is due to the poroelastic coupling, which results in a back-stress on the fracture due to the dilated matrix. The frac-  ture aperture at the injection well is shown in Fig. 2 b as a function of time. For the storage cases 1-3, the results of the present model match the analytical solution very well, whereas for the permeable cases 4-7, the fracture aperture has been reduced, due to the increased leakoff from the increase in matrix permeability. Furthermore, the fracture aperture reduces when considering poroelastic coupling, through the Biot coefficient in case 7. Included in this figure is the simulation of case 5 with a coarser mesh. Very good agreement is found between the results for two different meshes, confirming that the leakoff is mesh-independent (38,142 elements for the fine mesh, versus 20,342 elements for the coarse mesh). The two discretisations for the fracture are shown in Fig. 3.
The fracture radii predicted from the present model are compared with available solutions in Fig. 2 c. Again, good agreement is found between analytical solutions and the present model, for storage cases 1-3. The fracture radius decreases with increasing permeability in cases 4-7. A similar trend is observed between the fracture radius and the Biot coefficient. Increasing the permeabil-ity and the Biot coefficient increases leakoff, and reduces fracture volume. The increase in fracture pressure and decrease in fracture aperture and radius are due to leakoff, which causes dilation of the rock matrix. Good agreement is also found between the present model results in case 5 and the analytical solution for fracture radius. In cases 6 and 7, the present model predicts lower fracture growth than the leakoff-viscosity asymptotic solution; this is due to the combined effects of poroelasticity and multi-dimensional flow through the matrix. The analytical solution applies for onedimensional leakoff without poroelastic effects.

Toughness-dominated regime
The fluid viscosity and fracture toughness values for toughnessdominated cases are set to μ f = 0.0 0 01 Pa s, and K ic = 2 × 10 6 Pa m 0.5 , respectively. Again, for storage-dominated cases, the permeability of the matrix is assumed to be negligible ( k m = 0), whereas for leakoff dominated cases the matrix permeability is taken to be 1 × 10 −15 m 2 . The Biot coefficient α is varied between 0 and 1, corresponding to uncoupled and fully coupled poroelasticity scenarios, respectively. The dimensionless viscosity for given parameters is M = 0.003, which corresponds to a toughnessdominated case ( M 1). In toughness-dominated cases, the effects of matrix permeability, minimum in situ stress normal to the fracture plane, and Biot coefficient on the injection pressure, fracture aperture and fracture radius are illustrated in Figs. 4-6 . Fig. 4 a shows the net injection pressure versus injection time for cases 8-11. These cases show the effect of changing the matrix permeability. Good agreement is found between the present model and the analytical solution for storage case 8. In permeable cases 9 −11, the injection pressure increases. Bunger et al. (2005) proposed solutions for a penny-shaped hydraulic fracture in the leakoff-toughness dominated regime. In their solution, leakoff is modelled using Carter's equation ( Carter, 1957 ). The asymptotic solution for the leakoffdominated regime also shows an increase in the injection pressure. In case 10, the flow through the matrix is considered onedimensional along the direction normal to the fracture. Good agreement is found between the present model results and the analytical solution for this case, validating the present model. Leakoff increases by increasing the flow dimensions within the matrix, as well as by considering poroelasticity effects through the Biot coefficient. Therefore, in case 11, the present model predicts higher leakoff than does the analytical solution. This is due to multidimensional flow in the rock matrix, a feature that is not captured in the analytical solution. In case 11, the growth increment has been reduced to 0.25 m. A smaller growth increment corresponds to more growth steps and more computational time.
In Fig. 4 b, the fracture aperture at the well versus injection time is shown for cases 8-11. Again, good agreement is found between current model's results and the analytical solution for both storage case 8 and one-dimensional leakoff case 10. As permeability increases in cases 9-11, leakoff increases, and the fracture aperture reduces. Similar trends are observed for the graph of fracture radius versus injection time, as shown in Fig. 4 c. In cases 12 and 13, the in situ normal stress is set to zero, which substantially reduces the leakoff, as shown in Fig. 5 . Again, the injection pressure increases with increasing permeability, whereas fracture aperture and radius both decrease. However, the effect of changing the permeability in the absence of an in situ stress is not as significant as changing the permeability in the presence of an in situ stress ( Fig. 4 ). Therefore, the in situ stress has a significant impact on leakoff; meaning that the leakoff could still be important in low permeability reservoirs provided that the minimum in situ stress is high enough. As the fracture pressure is not explicitly computed using Carter's model, but is often approximated as being equal to the minimum in situ stress, Carter's model predicts lower leakoff when the in situ stress is reduced. The fracture net pressure in the present simulations varies between 1 to 10 MPa (higher values for the viscosity-dominated regime, and lower values for toughnessdominated regime), which could be considerable compared to the actual minimum in situ stress.
The effect of poroelasticity on leakoff is investigated through cases 14 and 15, in which the Biot coefficient is increased to 0.5 and 1.0, respectively. The results for net injection pressure, fracture aperture at the well, and fracture radius versus injection time, are shown in Fig. 6 . The increase in injection pressure is more dominant when the poroelastic coupling between the matrix fluid and mechanical deformation is considered, as shown in Fig. 6 a for cases 14 and 15. This is due to the back-stress from the dilated rock matrix. Increased fracture pressure increases leakoff, which in turn reduces fracture aperture and radius, as shown in Fig. 6 b and c, respectively. Included in Fig. 6 c is the simulation of case 15 with stress intensity factors computed using displacement correlation (DC) method. Good agreement is found between the two results, which validates the computation of stress intensity factors. Fig. 7 shows the fracture pressure distribution over the fracture surface, together with matrix pressure contours over a cut plane perpendicular to the fracture plane for cases 7 ( Fig. 7 a), 11 ( Fig. 7 b) and 13 ( Fig. 7 c) in the viscosity, toughness, and intermediate regimes, respectively. Higher viscosity of the injected fluid results in a nonuniform fracture pressure distribution. The matrix pressure distribution in all cases illustrates the three-dimensional nature of flow through the matrix. The time-step shown in Fig. 7 a is the immediate step prior to growth-step 7. However, the minimum fracture pressure is less than the in situ stress, which means that the fracture growth is controlled by the flow of the viscous fluid; the fracture toughness has negligible effect. In the viscosity-dominated regime, the fracture propagates as soon as the fluid reaches the fracture tip. The distance between the fluid front and the fracture tip is called the fluid lag. The size of the fluid lag depends on the in situ stress, and it reduces by increasing the in situ stress. The fracture and matrix pressure distribution for toughness-dominated case 11 is shown in Fig. 7 b. Fracture pressure distribution in this case is almost uniform compared to the fracture pressure distribution in the viscosity-dominated case shown in Fig. 7 a. This is due to the negligible energy required to move the low-viscosity fluid (case 11). This result supports the assumption of a uniform pressure distribution in the fracture in the toughnessdominated regimes. The minimum fracture pressure is above the in situ normal stress in the time-step shown in Fig. 7 b; indicating that the flow has reached the fracture tips, but the fracture is not growing, as its equivalent stress intensity factor is less than the fracture toughness of the rock formation. This verifies the fact that the frac- Fig. 7. Fracture pressure distribution (solid surface) and matrix pressure contours on a cut-plane perpendicular to the fracture. (a) Case 7 ( σ = 7 × 10 6 Pa; k m = 1 × 10 −12 m 2 ; α = 1) in leak-off-viscosity. (b) Case 11 ( σ = 7 × 10 6 Pa; k m = 1 × 10 −15 m 2 ; α = 0.1) in leak-off-toughness. (c) Case 13 ( σ = 0 Pa; k m = 1 × 10 −15 m 2 ; α = 1). Fracture pressure distribution represents the viscosity/toughness regime, while matrix pressure contours shows the multi-dimensional nature of the flow within the matrix. ture growth is controlled by fracture toughness, as is expected for the toughness-dominated regime. For case 13, the matrix pressure distribution is totally different than its distribution in the previous cases 5 and 10. Negative matrix pressures are observed in front of the fracture tips, where the tensile stresses develop matrix expansion. The Biot coefficient is set to 1 in this case, so matrix expansion develops negative matrix pressures due to poroelastic coupling. Zero in situ stress significantly reduces leakoff, such that the leakoff flow cannot compensate the negative pressures due to the matrix expansion.
Effective normal stress contours in a cut-plane perpendicular to the fracture surface for cases 6 and 7 are shown in Fig. 8 . Whereas for uncoupled case 6, the effective normal stress contours follow the fracture pressure pattern shown in Fig. 7 , in coupled case 7, the increased matrix pressure alters the effective stress, so that the contours in the vicinity of the fracture are parallel to the fracture. In the uncoupled case, the maximum effective stress is located at the injection point, and is equal to the maximum fracture pressure, whereas in the coupled case, the maximum effective stress is less than the maximum fracture pressure, and is located away from the fracture. The results shown in Figs. 7 and 8 are for one timestep prior to growth occurrence. In this step, the fluid pressure at the fracture tips is almost equal in magnitude to the in situ stress (7 MPa), and so the stress singularity is very weak, and does not show up at the scale at which the stresses are plotted. In summary, the simulation results show that leakoff increases with increasing the permeability of the rock matrix, in situ stress, the Biot coefficient of poroelasticity, and the dimensionality of the flow within matrix.
Available asymptotic solutions for the leakoff regime do not account for the effects of poroelasticity and multidimensional flow within the matrix, and consequently they predict lower leakoff. Lower leakoff is associated with lower fracture pressure, and higher fracture aperture and radius. Actual field cases may not correspond to any of the extreme cases for which asymptotic solutions are available, warranting the need for a robust numerical model that accounts for all the processes that occur during hydraulic fracturing. Convergence tests has been also performed to investigate the robustness of the model. For this test, the case 2 is considered with three types of mesh: coarse (9540 elements), medium (17,441 elements) and fine (29,080 elements). The error between injection pressures calculated from the present simulations and from the analytical solution at a time of 20.85 s is computed and shown for different meshes in Fig. 9 . The error reduces from 5% for the coarse mesh to close to 1% for the fine mesh.

Conclusions
A fully coupled three-dimensional finite element model has been presented for the simulation of hydraulically driven fractures in poroelastic rocks. The model was validated against available analytical solutions for penny-shaped fractures in different regimes of propagation. The simultaneous impact of fluid and rock matrix parameters on the hydraulic fracturing variables such as injection pressure, fracture aperture, and fracture length, has been studied. It was shown that leakoff increases with increasing matrix permeability, increasing in situ stresses, increasing Biot coefficient, and increasing flow dimensions. The analytical asymptotic solutions for hydraulic fracturing in an impermeable rock matrix provide an upper-bound for fracture radii and apertures, and a lower-bound for fracture pressure. For leakoff-dominated cases in a permeable rock matrix, the analytical asymptotic solutions cannot represent a lower bound for fracture radii and apertures, or an upper bound for fracture pressure. This is because the analytical solutions do not capture the effect of the fracture being compressed by the dilated poroelastic matrix, and underestimate the leakoff due to their neglect of global flow in directions parallel to the fracture wall. For intermediate cases, which are most likely to arise in actual reservoirs, there are currently no analytical solutions available. The extent of leakoff is also important when evaluating environmental impact. In summary, this work has delineated the range of applicability of the available asymptotic analytical solutions, identified the trends (in fracture aperture, fracture length, and injection pressure) due to deviations of the input parameters from the values assumed in the asymptotic analytical solutions, and demonstrated a new code that can, for example, be used for more complex geometries, such as multiple fractures and/or multiple wells (work in progress).
for partially funding this work through the TRUST Collaborative Project, 309067.