Hydraulic fracturing analysis in fluid‐saturated porous medium

Abstract This paper addresses fluid‐driven crack propagation in a porous medium. Cohesive interface elements are employed to model the behaviour of the crack. To simulate hydraulic fracturing, a fluid pressure degree of freedom is introduced inside the crack, separate from the fluid degrees of freedom in the bulk. Powell‐Sabin B‐splines, which are based on triangles, are employed to describe the geometry of the domain and to interpolate the field variables: displacements and interstitial fluid pressure. Due to their C1‐continuity, the stress and pressure gradient are smooth throughout the whole domain, enabling a direct assessment of the fracture criterion at the crack tip and ensuring local mass conservation. Due to the use of triangles, crack insertion and remeshing are straightforward and can be done directly in the physical domain. During remeshing a mapping of the state vector (displacement and interstitial fluid pressure) is required. For this, a new methodology is exploited based on a least‐square fit with the energy balance and mass conservation as constraints. The accuracy to model free crack propagation is demonstrated by two numerical examples, including crack propagation in a plate with two notches.

have been proposed, such as finite elements, 4 the extended finite element method, [5][6][7] isogeometric analysis, 8 extended isogeometric analysis, 9 embedded strong discontinuities, 10 the phase-field method, 11 a coupled finite element-peridynamics model, 12 interfaces elements equipped with a cohesive zone model 13 and a combined finite-discrete element method. 14 Due to their simplicity and robust performance interface elements have gained popularity for modelling fracture initiation and propagation in a poroelastic medium. 15,16 The interface elements are placed in the mesh a priori, requiring information of the crack location. This limits the application of the interface elements in a general framework. Remeshing can however remove this limitation and has been used successfully in the simulation of fracturing in a porous medium. 17,18 For the interpolation of the field variables several basis functions have been employed in this context. Lagrange basis functions have been used frequently, 15,16 due to their simplicity and ease of implementation. However, the  0 -continuity nature of the basis function deteriorates the accuracy of the simulation. The stress is generally discontinuous at element boundaries and the crack tip. Normally, several elements are needed to well capture the fracturing behaviour, such as the crack initiation and propagation direction. Furthermore,  0 -continuous Lagrange bases lead to a discontinuous interelement pressure gradient. Accordingly, a local mass balance is not be guaranteed.
Due to their higher-order continuity, the basis functions used in isogemetric analysis -NURBS and T-splines -normally avoid the discontinuous stress field and a loss of local mass balance, and have been used for discrete crack modelling, 19 including in simulations of fracture in a fluid-saturated medium. 8,9 However, in discrete crack modelling the isogeometric approach requires the introduction of  0 -lines in cracked elements to confine the influence of cracks locally, 20 thus eliminating locally the advantage of isogeometric analysis, namely higher-order continuity. Moreover, new crack segments are inserted in the parameter domain, rather than in the physical domain, and a reparameterisation of the domain is required to align the mesh with the crack path in the physical domain. This makes it mandatory that the initial mesh is sufficiently aligned with the final crack path a priori. 20 In this contribution, we employ Powell-Sabin B-splines to simulate hydraulic fracturing in a fluid-saturated porous medium. Powell-Sabin B-splines are defined on triangles, holding  1 -continuity throughout the entire domain, even at crack tips. This avoids the inaccuracy of the stress evaluation when employing Lagrange basis functions. Due to the flexibility of triangular elements crack insertion is carried out directly in the physical domain, thus avoiding the limitation adhering to isogeometric analysis. After the crack insertion the domain is remeshed to avoid elements with unsuitable aspect ratio, resulting in new Powell-Sabin B-splines on a new mesh. Then, a mapping of the state vector (displacements and pressure) is performed from the old onto the new mesh.
We start with an introduction of the governing equations for the hydraulic fracturing analysis. Section 3 derives the weak form of the governing equations. Next, we present the Powell-Sabin finite element discretisation. The basis functions and poromechanical interface elements are introduced here. In Section 5, we discuss the algorithm to insert a new crack segment, including the algorithm for remeshing and state vector mapping. In Section 6, numerical examples are given which demonstrate the versatility and accuracy of the method.

GOVERNING EQUATIONS FOR THE POROUS MEDIUM
Hydraulic fracturing in porous media is a complex physical phenomenon, including fluid flow in fractures, the pore fluid flow in the porous medium, and the deformation of the porous medium. In this contribution, we confine the discussion to (1) a fully saturated porous medium, (2) a Newtonian fluid, (3) infinitesimal strains and linear elastic material behaviour, (4) no mass transfer or chemical interaction between the solid and the fluid, (5) no consideration of gravity, and convective and inertia effects and (6) the fluid fully occupying the fracture, thus not considering possible fluid lag. We refer to Figure 1 for a graphical illustration of the investigated problem. The porous medium is split into two parts by an interface Γ . Biot's theory is used to model the porous media. 21 The cohesive-zone model is used to model the fracturing behaviour. 22,23

Deformations of the porous media
The fully saturated porous medium is modelled as a two-phase system with the solid skeleton fully filled with pore fluid. The deformation of the solid develops fast compared to the pore fluid pressure change. Thus, the deformation of the porous media can be considered as a quasi-static process, governed by the conservation of the hydro-static linear momentum 24 : solid pore fluid fracturing fluid F I G U R E 1 (A) A solid body with an internal discontinuity Γ . Γ is an interface boundary with positive and negative sides, Γ + and Γ − , respectively. Boundary Γ is prescribed with a displacement̄; Γ with a prescribed traction̂; Γ with a prescribed fluid pressurē; Γ with a prescribed infloŵ; and Γ with a prescribed inflow of fluid̂; (B) pressure around the internal discontinuity Γ .
where is the total stress, composed of the solid and pore fluid parts in which represents the apparent fluid pressure, denotes the unit tensor and is the Biot coefficient. is the stress inside the solid material, which is linearly related to the strain by with the elastic stiffness tensor, which is a function of two material constants -the Young's modulus and the Poisson's ratio for isotropy.

Interstitial fluid pressure
The governing equations of the pressure in the interstitial fluid are obtained from the mass balance of the mixture 13 in which denotes the effective permeability coefficient of the porous medium, = ∕ . and are the intrinsic permeability of the porous medium and the viscosity of the fluid, respectively, and is the Biot modulus: with and the solid and fluid moduli, respectively, and the porosity.̇represents the velocity of the solid and□ denotes the time derivative:□ = □̇= (6) As starting point for the derivation of Equation (4), Darcy's law is used: with the fluid flux and the velocity of the fluid.

Fluid flow in the fracture
To characterise the fluid flow inside the crack Γ , we consider a fully open crack filled with a Newtonian fluid, 24 see Figure 1A. The crack opening is assumed to be small compared to its length. In a two-dimensional setting the mass balance for the flow inside the crack reads: with = ⋅ and = ⋅ being the tangential and normal components of the fluid velocity in the crack, respectively. and are the tangential and normal vectors at the crack Γ , see Figure 1A, and ( , ) denotes the local coordinate system along the crack Γ . Integrating Equation (8) over the fracture height ℎ leads to the difference in the fluid velocity normal to crack faces: where + and − are the normal fluid velocities at the interface Γ + and Γ − , respectively. In this contribution, the fracture height ℎ is taken as the normal displacement jump in the fracture. We assume that the tangential fluid velocity inside the crack is much higher than the normal fluid velocity . Since the fracture height ℎ is much smaller than its length, we consider a constant fluid pressure inside the crack, see Figure 1B. The balance of momentum in the tangential direction is then derived as with being the shear stress, which is derived from = . Now, we reformulate Equation (10) as and the tangential fluid velocity is obtained by solving Equation (11) with no-slip boundary conditions at = ±ℎ∕2: Substituting Equation (12) into Equation (9) yields To obtain the normal fluid velocity difference in Equation (13), we assign the fluid pressure + to Γ + , − to Γ − and to Γ , respectively. The existence of an independent pressure inside the crack allows to model pressurising the crack. The inflow of the fluid through the interface Γ + and Γ − could be different due to the insertion of the interface Γ and the independent pressure . In principle, the resistances at Γ + and Γ − can be different, but we assume them to be equal, so that the interface permeability reads . Analogous to Darcy's law, the normal fluid velocity is determined as 24 : Substituting Equation (14) into Equation (13) yields the mass balance equation for the flow within Γ

Initial and boundary conditions
The governing equations of the porous medium are completed by the initial and boundary conditions. The Dirichlet boundary conditions read wherēand̄denote the prescribed displacement and pressure, respectively. The Neumann boundary conditions are given by in whicĥand̂represent the prescribed traction and inflow, respectively. For the flow in the fracture, the boundary conditions read is the pressure imposed on Γ .̂denotes the inflow imposed on Γ , see Figure 1A. Finally, the initial conditions give as where 0 ,̇0 and 0 represent initial displacements, velocities and pressures separately.

Cohesive-zone model
On the crack faces, we consider traction boundary conditions: ) are a non-linear function of the displacement jump across the crack interface Γ , which in the local coordinate system ( , ) reads: being the displacement jump across Γ in the local coordinate system ( , ). The traction vector relates to the traction in the global coordinate system via a standard transformation: with the rotation matrix. 25 In this study, an exponential traction-separation law is used, defining the traction in the normal and shear directions as: with the tensile strength of the material,  the fracture energy, int the initial crack shear stiffness (when = 0), and ℎ = ln( =1.0 ∕ int ). is a history parameter, set through a loading function = ( For the unloading ( < 0), the tractions are derived from a secant relation. To prevent interpenetration, a penalty stiffness is specified in the normal direction.

WEAK FORM OF THE GOVERNING EQUATIONS
The weak form of balance equations (1), (4) and (15) is obtained through multiplication by the test functions , and for the solid skeleton, the interstitial pressure and the fluid pressure within the fracture, respectively. Employing the divergence theorem and considering the internal boundaries Γ + and Γ − as well as the conditions at the external boundaries Γ , Γ , Γ and Γ , yields the weak form: Considering force equilibrium conditions at crack faces, we have with = − = − + . Reformulating Equation (25a) with the aid of Equation (26) leads to The fluid transport across the crack interface Γ + and Γ − is formulated in a Darcy-like manner 24 : Substituting this equation into the weak form of the mass balance, Equation (25b), results in

POWELL-SABIN FINITE ELEMENT DISCRETISATION
In this section, we will succinctly summarise the finite element discretisation. Powell-Sabin B-splines are used to approximate the trial functions in the solution space, and to parameterise the geometry. 26 Poromechanical interface elements are used to consider the porous effect on the crack faces. The linearised set of equations for the Newton-Raphson iterative scheme are given as well.

Powell-Sabin finite elements
Powell-Sabin B-splines are employed to discretise Equations (27) and (29). They are defined on triangles, maintaining  1continuity throughout the whole domain, even at crack tips. They describe the geometry and interpolate the displacement field and the fluid pressure in an isoparametric sense: where represent the coordinates of Powell-Sabin triangle corners , see Figure 2B. and denote the degrees of freedom at . The indices = 1, 2, 3 imply that each vertex has three Powell-Sabin B-splines attached, Figure 2B. denotes the total number of vertices on the triangulation  . A cracked domain Ω is considered, cf. Figure 2A. In Figure 2B, it is discretised by a triangulation  , and can be generated by any package for standard triangular elements, such as Gmsh. 27 On the triangulation  , there are = 1, 2, … , triangles and vertices, represented by thick black lines in Figure 2B. Each triangle is split into six ( = 1, 2, … , 6) mini-triangles, cf. Figure 2B. We perform the Powell-Sabin refinement  * and get Powell-Sabin points, green dots in Figure 2B. Then, we define a Powell-Sabin triangle for each vertex , 25,28 drawn in red in Figure 2B. Finally, we follow the procedure given in 29 to define the Powell-Sabin triangles on the boundary.
The Powell-Sabin B-splines over each mini-triangle in Figure 2C are obtained using the Bézier ordinates , , 25 :   26 We now formulate the Powell-Sabin B-splines in a matrix form and implement them in existing finite element codes by Bézier extraction = where are Powell-Sabin B-splines associated with each mini-triangle in element . are Bernstein polynomials. is a matrix filled by Bézier ordinates.
Powell-Sabin B-splines do not hold the Kronecker-delta property and are non-interpolatory at the vertex. 29 Thus, imposing Dirichlet boundary conditions on Γ and Γ is not as straightforward as for standard finite elements. In this contribution, we will employ Lagrange multipliers to weakly impose essential boundary conditions. 29

Poromechanical interface elements in the Powell-Sabin finite element scheme
The fluid inside the crack Γ induces a pressure on the crack faces. To include this effect we will apply interface elements enhanced with porosity. Powell-Sabin B-splines do not satisfy the Kronecker-delta property and are non-interpolatory at the vertex. The augmentation with pressure degrees of freedomis not as standard as in Lagrange basis functions. For the three pressure degrees of freedom (3PDOF) model, three pressure degrees of freedom are added ( Figure 3): one on each side of the crack face, and one inside the crack, which allows the discontinuity in the pressure field across the crack Γ and fluid injections along the crack Γ . The 3PDOF model assigns + to Γ + , − to Γ − and to Γ , respectively. To discretise Equation (25c), Powell-Sabin B-splines cannot be used because of their definition on triangles, not along a curve. In this study, standard finite elements (FEM) are used to describe the geometry of the interface Γ and to interpolate the pressure inside the interface, see Figure 3C. Here, quadratic Lagrange shape functions are considered.
with nodal degrees of freedom, quadratic Lagrange shape functions, and the coordinates of standard FEM nodes. The FEM nodes are the triangle vertices itself and the middle points between two vertices, see Figure 3C.
The mass conservation Equations (25c) and (29) contain time derivatives, which are discretised using the Backward Euler scheme. Considering the Powell-Sabin approximation Equation (30) and standard FEM approximation Equation (33), the weak form Equations (25), (29) and (25c) yields: where + and − are shape functions of related to interfaces Γ + and Γ − , respectively. and are shape functions and their derivatives related to the pressure degree of freedom inside the interface Γ . Linearisation of Equation (34) leads to equations for the Newton-Raphson iterative scheme: with the tangential stiffness matrices given in Appendix A. ext and int are external and internal force vectors, which could be obtained from Equation (34).

ADAPTIVE ANALYSIS FOR CRACK GROWTH
Due to the  1 -continuity of Powell-Sabin B-splines at the crack tip, point in Figure 4A, we can directly assess the fracture criterion at this point. The Rankine criterion has been used here to check crack initiation, comparing the major principal stress 1 with the tensile strength . If 1 ⩾ , a crack is inserted through the entire element, 0 in Figure 4B, in front the crack tip. There is no information about the curvature of the crack segment within the element 0 . Therefore, a straight line is inserted within 0 , shown in Figure 4C. The normal vector 1 of the new crack segment, in Figure 4 (B, can directly be obtained from the stress tensor at the crack tip. A further improvement of the quality of the crack direction prediction can be obtained by averaging the stress tensor over a finite space around the tip. After crack insertion element 0 is divided into two triangles 1 and 2 , Figure 4B. The triangular element next to the new crack tip, 3 in Figure 4B, will have four vertices. It is impossible to define Powell-Sabin B-splines on this element. Thus, remeshing is needed to remove triangles with four vertices or with a bad aspect ratio. We employ the algorithm proposed in 25 to remesh the domain, see Figure 4C. In the figure, the original mesh is denoted by Ω , while after remeshing the mesh is represented by Ω . To determine the domain Ω we stand at the element with the newly inserted crack segment, shaded grey in Figure 4A. Then, a radial marching is done until three elements have been crossed in all directions, see Figure 4. The elements along one side of the crack interface are excluded, which avoids updating the field variables along the crack interface.
After remeshing new elements and vertices can be added, and old elements may be moved to ensure elements with a suitable aspect ratio, see Figure 4. Consequently, the mesh is modified and Powell-Sabin B-splines are required to be computed on new triangles. Furthermore, due to non-interpolatory property of Powell-Sabin B-splines, the state vector, displacement and interstitial fluid pressure, needs to be mapped from the old mesh Ω onto the new mesh Ω at time step and + Δ . We take the mapping at time step + Δ for illustrative purposes. The mapping of the state vector is completed with a least-square fit subject to certain constraints. We first map the displacement +Δ from Ω to Ω : subject to: int, + coh, − int, − coh, = 0 on Ω and Ω =̂on Γ Γ (36) in which the subscript 'b' represents the matrix or vector associated with the old mesh Ω , while the subscript 'r' related to the new mesh Ω . denotes the matrix with the shape functions for the displacements. and are displacement vectors, and Γ is the boundary with prescribed displacement. We fix the degree of freedom on the red polygonal boundary and along the crack path of Ω , see Figure 4C.
In Equation (36), we consider the energy balance as the constraint equation, which matches the energy linked to Ω and Ω . The internal work i and the work c ℎ related to the cohesive traction on the crack surface are considered due to their direct relation with the displacement , and given as For the interstitial fluid pressure, the mapping is performed in a similar way as for the displacement, the optimisation problem being where is the fluid pressure shape function matrix, is the fluid pressure vector, and Γ is the boundary with the prescribed fluid pressure. In the constraint equation, we force the mass to be conserved between Ω and Ω . Only mass terms directly linked to the fluid pressure are considered. For the fluid pressure inside the crack Γ , after the crack insertion, we assign zero pressure values to the new crack segment. For the old crack segments, the value does not change. In this study, the MATLAB function fmincon is used to find the optimum in Equations (36) and (38). Alternatively, one can use optimisation packages like MOSEK 30 or ALGLIB, 31 which may provide a better efficiency for large-scale problems. In general, the constraint equations from the energy balance and mass conservation reduce the error level of state vector update. A detailed error analysis of the state vector update with a constraint equation for energy conservation under dynamic loading has been carried out in elsewhere. 28 The computation efficiency in the proposed method is somewhat lower than that in standard finite element analysis. In the evaluation of Equations (36) and (38), we need to find the state vector of Gauss points on the refined mesh from the old mesh. 32 For each triangular element, there are six mini-triangles used to carry out the integration. Thus, the number of triangles used in the integration is × 6, and the number of Gauss points on the refined mesh is × × 6, where denotes the total number of triangular elements and is the number of Gauss integration points inside each minitriangle. In standard FEM the number of Gauss points is × , which is smaller than that in the proposed method. Thus, the computation time in the proposed method will be increased in comparing to standard FEM.

NUMERICAL EXAMPLES
Below we will consider two examples. The first example deals with a pre-fractured specimen, assessing the accuracy of the method. The last example features crack propagation under mixed-mode loading conditions, demonstrating the ability of the method to analyse the propagation of curved cracks.

Single-edge notched plate
The problem consists of a square plate with dimensions 250 mm × 250 mm, with a horizontal crack through the centre of the plate, shown in Figure 5. The first 50 mm of the crack is pre-fractured, and an inflow of t = 50 mm 2 ∕s is imposed on the left end of the crack. The pressure is zero at the top, bottom and right boundaries. The displacement in the horizontal direction is constrained at the right boundary, while the displacement in the vertical direction is fixed at the top and bottom boundaries, see Figure 5A.  (23) is used with the tensile strength = 2.7 MPa and fracture energy  = 0.095 N/mm. Only mode I fracture is considered, that is, int = 0 in Equation (23). To avoid interpenetration, a penalty stiffness = 10 10 MPa/mm is specified in the normal direction of the crack. The interface permeability, , is set as the effective permeability coefficient of the porous medium . The plate has been discretised by the triangulation presented in Figure 5B. A constant time step size Δ = 0.01s is used in the simulation.  Figure 6A shows the crack opening along the crack interface Γ for different time steps. Obviously, the crack opening increases as the time evolves. This evolution is also evident from the pressure along the crack interface, see Figure 6B. In the figure, the pressure inside the crack is higher than the pressure on either side of the interface. The results match well with those in ref., 8 validating the proposed method. Indeed, in ref. 8 the set-up of the problem is similar to the KGD problem. 33,34 Figure 7 presents the contour plot of the interstitial fluid pressure, the flux in the -direction and the displacement norm. As shown in the Figure 7A, almost the entire crack interface Γ is pressurised due to the flow injection at the left edge, except the small part around the crack tip. The pressurisation along the crack interface can also be observed from the flux profile in the -direction, that is Figure 7B. Due to the symmetric setup of the problem, the displacement profile shows a symmetric distribution along the middle plane of the plate, see Figure 7B.

Arbitrary propagation: A plate with two propagating cracks
We next consider the specimen of Figure 8 to demonstrate the ability of the proposed method to properly analyse mixedmode crack problems. The specimen has a thickness of 50 mm and has two initial horizontal notches. Figure 8 (A shows the geometry and the boundary conditions. In the analysis, the specimen is subjected to a prescribed horizontal velocitẏ̄a nd a vertical velocitẏ. Fluid is injected at the inlet of the initial notches at a constant rate t . The time increment is set as Δ = 0.005s. The  (23) is employed to describe the fracturing process with a tensile strength = 3.0 MPa and a fracture energy  = 0.11 N/mm. Mode-II crack behaviour is considered: int = 10 N/mm and ℎ = 0 in Equation (23). 35 Plane-stress conditions are assumed and the loading condition is set up as:  Step 2 The specimen is then subjected to a prescribed horizontal and vertical velocitẏ=̇= 2 × 10 −2 mm/s, see Figure 8. A displacement control is employed to apply the velocitẏanḋin the simulation, Δ = Δ =̄× Δ =̇× Δ = 1.0 × 10 −4 mm. The displacement is imposed by Lagrange multiplier method in. 29 At this stage, fluid injects are terminated.
The load-displacement diagram in Figure 9A presents the relation between the vertical resultant force and the vertical displacement on the top edge at loading Step 2. In the process of crack propagation, the fluid pressure inside the crack and on the crack interface will gradually impose a tensile stress (negative fluid pressure) on the crack faces, see Figure 10A-middle and 10A-right. In the course of time nearly the entire crack interface will come in tension, as illustrated in Figure 10A-right, which induces the increase of the force on the top panel. The prediction of the crack path is presented in Figure 9B. The curved crack path shows the refinement ability of the Powell-Sabin B-splines. There is a kink on the crack path, due to the change of loading conditions from Step 1 to Step 2. In Figure 10, the profiles of the interstitial fluid pressure, the flux and the displacement are illustrated. At the loading Step 1, the crack interface is pressured due to the fluid injection t , Figure 10A-left. The pressurisation can also be observed in the flux plot, Figure 10B-left. After terminating fluid injection at loading Step 2, the fluid will gradually flow back into the fracture, Figure 10A-middle and 10A-right. The flux profile also presents this flowing back, see Figure 10B-middle and 10B-right. The displacement in the -direction is increases as time evolves due to the monotonic loading conditions, Figure 10C.

CONCLUSIONS
Powell-Sabin B-splines have been used in the analysis of hydraulic fracturing. Cohesive interface elements are employed to model the interface behaviour of the solid part, while a three-pressure degree of freedom model describes the fluid flow inside the fracture. Powell-Sabin B-splines are  1 -continuous throughout the domain, even at crack tips. Such higher-order continuity remedies the inaccurate stress issue in employing Lagrange basis functions. Due to the definition of Powell-Sabin B-splines on triangles, crack insertions and crack path tracking are directly performed in the physical domain, circumventing the initial mesh alignment issue in the isogeometric analysis. After inserting new crack segments, remeshing is required to avoid elements with unsuitable aspect ratios, which necessitates a mappling of the state vector from the old onto the new mesh. A novel least-square fit methodology is introduced in combination with constraint equations from the energy balance and mass conservation.
Numerical examples show that the refinement ability of the Powell-Sabin B-splines is very suitable for the analysis of hydraulic fracturing. When fluids are injected into the fracture, the crack interface will be pressured, rendering an increase of the crack opening and forcing the crack propagate. Supposing that no fluids are injected into the fracture, under tensile loading conditions, the existence of fractures will induce tensile stresses on the crack interface, preventing crack openings. The fluid then flows from the porous medium into the fracture. The extension of Powell-Sabin B-splines to three-dimensional problems is non-trivial due to certain constraints with neighboring tetrahedrons. 36,37 Alternatively, one can construct prisms as a tensor product of two-dimensional Powell-Sabin B-splines and Non-Uniform Rational Basis splines (NURBS) in the third dimension.

A C K N O W L E D G E M E N T
Financial support from the European Research Council (Advanced Grant 664734 "PoroFrac") is gratefully acknowledged.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing not applicable to this article as no datasets were generated or analysed during the current study.
with the rotation matrix, 25 the tangent stiffness of traction-opening law at the interface Γ . 25