Phase-field modelling of fluid driven fracture propagation in poroelastic materials considering the impact of inertial flow within the fractures

This paper presents a computational framework for modelling of fluid pressurised fracture propagation in saturated porous media. The framework rests on the principle of the variational phase-field theory to predict the fracture propagation pathway. The paper sets out the variational formulations and associated weak forms of the partial differential equations describing the pressure-deformation interplays of the fracturing domain, which are solved in the context of the Updated Lagrangian Finite Element method. The proposed formulation reflects the impact of the temporal evolution of the porous media attributes such as porosity, compressibility, permeability, and mechanical stiffness, on the nonlinear hydro-mechanical behaviour of the porous media during the fracture propagation. The inertial effect of the nonlinear flow inside the fracture is resolved using Forchheimer equation. Robustness of the modelling framework is examined by simulating benchmark examples. The effects of poroe-lastic characteristics of porous media such as the compressibility of solid skeleton and drained bulk modulus on the hydro-mechanical and cracking behaviour of porous rocks and on the total energy of the system are addressed. The nonlinearity of the fluid flow is found to be influential on the length of the leak-off and flow-back regions across the fractured zones, and on the amount of the fluid to be exchanged between the fractures and the porous zone, which is important in the prediction of the productivity of the fracking process in engineering applications.


Introduction
Robust modelling of the interplay between solid and fluid bulk phases in porous media is a key for sustainable design of hydraulic fracturing operations 1 with diverse applications ranging from enhanced oil and gas recovery, 2 geothermal energy production, geological carbon dioxide storage, and block cave mining.4][5] With the increasing interest to the hydraulic fracturing technique, more advanced mathematical models were proposed to depict the interactions between the fluid flow and solid deformation in the fracturing zone and crack-tip. 6aboratory experiments also demonstrated the critical role of the fluid characteristics (e.g., fluid density, viscosity, and compressibility) and the flow features (e.g., the leak-off of the pressurised fluid into the porous rock and pore flow velocity) on the formation and evolution of fractures. 7arious numerical methods have been proposed to simulate the hydraulic fracture propagation in porous media.In the context of finite element method, the extended finite element method (XFEM), discrete fracture modelling combined with adaptive remeshing algorithms, and continuum damage modelling have been employed to trace fracture paths. 8The XFEM provides a robust basis upon which enrichment functions can be added to finite elements for modelling displacement jumps and the interaction between fractures. 9However, coupling hydraulic processes to mechanical deformation and cracking is challenging due to the correct choice of enrichment functions, and completely different governing equations might be required for modelling flow inside and outside the fracture, 10 which can cause inconsistency in predicting effective stresses in the solid skeleton.Discrete-based methods lack sufficiency in capturing the fluid exchange between the hydraulic fracture and the surrounding poroelastic medium; therefore, they mostly consider hydraulic fractures in impermeable elastic media. 8However, they show a high capability in modelling large deformations due to cracking and capturing the fluid flow within the discontinuities, also providing high computational efficiency especially in meshfree approaches. 11Continuum damage theory proposes an advantageous approach for solving sophisticated coupled partial differential equations (PDEs) for Multiphysics-multidomain processes, since it offers a mathematical setting to be solved on a unified finite element mesh through borrowing the variational theory.Hybrid methods are developed to benefit from different aspects of numerical methods.For instance, modelling the fluid pressure field and flow throughout the porous rock and within the fracture channels in the framework of discrete-based methods for modelling fractures as displacement discontinuities have been done over a finite element background mesh. 8mong continuum-damage-based approaches, the phase-field method has been successfully applied for modelling fluid-driven fracture propagation in porous media within the context of the linear elasticity [12][13][14][15][16][17] and hyperelasticity. 18,19By modelling hydro-mechanical processes and fracture propagation over a continuous finite element mesh, the phase-field approach provides a reliable framework for capturing fluid flow in the porous solid and inside the damaged zone (fracture) as well as the fluid exchange between these two regions.The phase-field method benefits from the versatility of the variational approach in simulation of the propagation, joining, and branching of the fluid driven fractures with geometrical complexities and potential material heterogeneities. 20,21The inherent characteristic of this approach for modelling fractures as a result of energy minimisation has opened the way for the use of Neural Networks in predicting fracture propagation in structures, which in turn can contribute to enhancing the uncertainty quantification of fracture propagation in heterogeneous materials and the design optimisation due to the development of fractures in structures. 22,23The input parameters in the phase-field method for modelling hydraulic fractures are either material properties such as Lamé parameters, permeability, and porosity, or user-defined parameters such as the size of finite elements and the so-called length-scale parameter which can significantly affect the response of a fracturing body. 24This inherent source of uncertainty has roots in the assumption of modelling fractures as diffusive regions of damage instead of sharp discontinuities, so conducting probabilistic sensitivity analysis to unveil the obscure effects of the user-defined input parameters is vital to ensure the results of phase-field fracture model especially in industrial applications. 25ourdin et al. 12 applied the variational theory to minimise the total energy of a saturated porous system to simulate the hydraulic fracturing.Mikelić et al. 13 adopted the theory of Biot in poroelasticity and defined a staggered method where a fixed-stress state of the damaged solid is considered to solve the associated PDEs for the pressure inside the reservoir and within the cracks; then, a fixed pressure field is used to calculate the mechanical response of the porous system.These works rested on the concept of the linear elasticity and did not account for the changes in material strength due to the mechanical deformations.Building a numerical model based upon the kinematics of large deformations provides the foundation for using complex constitutive laws for modelling the materials which show different brittle/ductile failure behaviour depending on the rate of deformation. 26The phase-field method was employed by Miehe and Mauthe to model the cracking behaviour of the solid skeleton based on hyperelasticity within a Total Lagrangian formulation 18 and based on infinitesimal strain theory in a poroelastic medium. 14They described the fluid flow throughout the porous medium and inside the fracture by Darcy's law and used Poiseuille's law to calculate the evolution of the permeability within the fractures.However, the inertial effect of laminar high-velocity flow inside the fracture was neglected.In their work, there is no clear distinction between compressible or nearly incompressible solid skeletons in the derivation of the bulk and fluid energy functionals.Thus, due to volumetric strain locking effect, the energy functional overestimates the solid effective energy of the materials with low bulk moduli. 27While the previous studies were based on the macroscopic approach in dealing with the porous media, a microscopic approach was adopted in Ref. 15  to derive the weak forms of the balance of linear momentum and mass conservation to be implemented in a set of small strain Finite Element formulations.The phase-field theory has been used by Wilson and Landis 19 who considered the fundamental rules of thermodynamics and the balance of micro-forces on the internal faces of the fractures caused by the fluid pressure in their modelling method.These analyses were restricted to the linear laminar flow regime, using Darcy's law.The cracking behaviour of saturated porous rocks is influenced by poroelasticity coefficients such as bulk modulus of solid constituent and other experimentally measured coefficients like drained and undrained bulk moduli. 28Despite this fact, the previous studies on the phase-field modelling of hydro-fractures have not considered these effects in their simulations, neither clearly stated in the formulations.
It has been attested from the literature that the inertial effect of laminar flow within the fractures and porous media becomes significant when the Reynold's number (Re) is in the range of 1-10; however, the nonlinear laminar regime can persist up to Re = 150. 29Forchheimer equation has been widely used to model this type of non-Darcy flow in fractured porous media containing an inertia resistance factor, so-called β-factor, which can be related to the microscale characteristics of the porous media. 30The inertial effect of the fluid flow causes a reduction in the hydraulic conductivity of the fracture 31 and influences the fluid pressure field in the vicinity of pressurised fractures. 32Hence, a hydro-fracture modelling framework that captures the nonlinear behaviour of inertial flow within the fracture zone is needed to reliably predict the hydro-mechanical interplay and the flow regime in the fracturing zone.This is what has been tackled in this study.Since Darcy's law and Forchheimer equation have the same theoretical base (Navier-Stokes equations), 33 we develop a coupled Darcy-Forchheimer model that can apply the Forchheimer correction of fluid conductivity in the fracture when a significant increase in fluid velocities is detected.
We base our work on a macroscopic view to formulate the potential energy functional of the system 34 and employ the variational principle to propose the finite element weak forms to model the propagation of the fluid driven fractures.The nonlinear response of the system due to complexity of the inertial flow, progressive pore collapses in the solid skeleton, and the evolution of the hydro-fractures, are modelled by describing the microscale hydro-mechanical interplay of the fluid and solid phases. 35The derivation of the energy functionals which mobilise the driving force for the hydraulic fracturing have been done for both compressible and nearly incompressible hyperelastic materials, which have not been formulated thoroughly in previous studies.An updated Lagrangian finite element formulation in the realm of finite strain theory is developed to solve for mechanical deformations and fluid pressure fields monolithically and for the damage in a staggered manner.The generalised concept of the effective stress in a macroscopic sense is used to apply the stiffness degradation criterion of Miehe 36 on the tensile part of the effective strain energy of the solid skeleton.
The paper is structured as described in the following.Mathematical formulations describing the evolution of the damage and hydromechanical mechanisms are presented in section 2, and the finite element weak forms in an updated Lagrangian framework are developed in section 3. Finite element discretisation and computational algorithms are presented in Section 4. Section 5 exposes the verification and key results of the study, and concluding remarks are provided in Section 6.

Methodology
The fluid-driven fracture is modelled by employing the variational principle to minimise the total potential energy of a fracturing solidfluid porous mixture.The mathematical procedure to derive the related governing equations are described in this section.Notice that boldface letters are used to distinguish matrices and vectors from scalar quantities.

Geometrical setting of fluid-solid porous mixture
We consider a physical porous saturated body in Euclidean space R δ , in which the fractures are being developed as subjected to external loading, fluid injection, and deformations.Since the fracturing process N. Sarmadi and M. Mousavi Nezhad occurs in the solid matrix and the motion outside of the body is not within the scope of this study, the motion of the solid skeleton is chosen as the principal, and the motion of fluid is formulated with respect to the motion of the solid.Consider the time interval I = [0, T], the original undeformed state of the body at time t = 0 is recognised by Ω 0 ⊂ R d , d ∈ [2, 3] in the reference coordinates X i , and current displaced positions of the solid particles in space are specified in the current coordinates x i by Ω t at time t = T. Finally, the motion function of the solid matrix φ, shown in Fig. 1a, is prescribed as a mapping function that transfers reference coordinates X i to the associated coordinates in the current configuration x i .
The evolution of the solid deformation in the body can be traced and formulated in subsequent steps using the second-order tensor deformation gradient F defined as Considering a volume element dV within the saturated undeformed porous body Ω 0 , shown in Fig. 1a, that is consisted of a solid phase with mass density ρ s and a fluid phase with mass density ρ f , the following volume fractions are defined in the reference state of the body Ω 0 as where, dV f is the volume of the fluid (equivalent to the volume of the pores in the saturated condition) in the initial undeformed state Ω 0 .It is known that the relationship φ f + φ s = 1 always holds in the current configuration after deformation.The initial mass of a unit volume dV of the saturated porous body is recognised as m 0 = ρ s φ s 0 + ρ f φ f 0 , which is subjected to a deformation φ(X, t).During the process of deformation from the reference to the current state of the body in time t, the rate of the fluid mass m flowing into/out of the unit volume of the undeformed reference body Ω 0 is defined as where, D s ()/Dt operator is the Lagrangian time derivative of the fluid mass with respect to the solid motion, and M represents the Lagrangian fluid mass rate in the space.Using the divergence theorem, the quantity D s (m)/Dt can be defined as the surface integral of the vector field M on the boundaries of Ω 0 as ∯ ∂Ω0 M⋅NdA.Accordingly, the vector m r can be reached using the relationship where, nda = JF − T NdA, J = det(F) is the Jacobian of the deformation gradient, and m r is the mass of the moving fluid phase relative to the solid skeleton movement in the space.The vector m r can be formulated as where, v r = φ f (v f − v s ) is the relative superficial velocity for the fluid phase, and v f and v s are the intrinsic velocities of the fluid and solid phases, respectively, in the current configuration Ω t .

Phase-field modelling of crack propagation in porous media at finite strains
Griffith introduced the required energy to create new crack surfaces in a body as the Critical Energy Release Rate (G c ), which is a material property and introduced as 37 where, D = 2γA is the amount of surface energy, γ is crack surface density, A is half of the newly formed crack surface and Π is the internal potential energy of the body.By forming new surfaces throughout the volume, the corresponding amount of the energy that is released from the system is calculated as ∫ Γ (G c )dA.Since the crack propagates through the body, the boundary Γ is mathematically unknown, so the surface energy can be approximated by using Mumford-Shah-type functional 38 in the manner introduced by Ambrosio and Tortorelli 39 as the regularization functional for brittle fractures.As a result, the surface energy functional ∫ Γ (G c )dA is defined on the reference body Ω 0 by introducing the crack surface density γ(s, ∇s) over the volume as a function of damage parameter s.The parameter s provides a picture of the damaged state of the body, such that s = 1 and s = 0 represent the undamaged intact zone and fully damaged zone, respectively.A Dirac-type function s(x ' ) = 1 − e − |x ' |/2l0 is used to represent the crack over the reference configuration, as depicted in Fig. 1b.A second order approximating functional γ(s,∇s) = l 0 |∇s| 2 + (1 − s) 2 /4l 0 , introduced by Francfort and Marigo, 40 is utilised in this study to formulate the total new surface energy over the reference volume Ω 0 as Ψ sup (s) = G c ∫ Ω0 γ(s, ∇s)dV, where l 0 is the length-scale parameter.Having adopted the variational phase-field approach for modelling the crack propagation in the porous media, the total internal potential energy functional for the system is formed as ) dV (9)   where, the bulk energy of the saturated porous medium Ψ bulk (φ, m) is dependent on the deformation and the fluid mass change in the system. 34

Governing equations in the finite strain hydro-mechanical process
In a closed thermodynamic system, the variation of the internal energy of the system, Equation (10), with respect to an arbitrary change in the fluid mass δm is set to zero.

DΠ(φ,m,s)[δm]=
where , Ψ vol (φ, m) is the volumetric part of the bulk energy and the deviatoric part of the energy Ψ dev (φ) is solely function of the deformation.The pressure function ϑ = ∂Ψ vol /∂m is defined as the work done on the reference porous material by the extraction of a unit mass of fluid causing the pore fluid pressure to change from p 0 to p 34 .
Accordingly, increase of the pore volume per unit initial volume of the bulk material in this process is defined by notation of dη = dm/ρ f , where η is called the fluid content. 34By substituting dη into Equation (11), Equation (10) is reformulated as The above equation results in the derivation of the fluid pore pressure in the saturated porous medium as p = ∂Ψ vol /∂η.We will be able to obtain the volumetric part of the bulk energy by adopting the following relationship for the incremental changes of the fluid pressure 41 where, (J − 1) represents the volumetric strain of the solid matrix, and M b is Biot's modulus which is defined in Ref. 41 as The parameters K and K u are drained and undrained bulk moduli of the mixture, and α is the Biot coefficient, which will be defined later.Using Equation ( 5) and knowing that the rate of the fluid content is related to the rate of the fluid mass via the relationship of ṁ = (ρ f ) η, the final form of the fluid mass conservation in the reference configuration Ω 0 is written as Using Equations ( 6) and (7), the above equation in current configuration is obtained as where, v r = φ f (v f − v s ) and θ = Jp is the Kirchhoff pore pressure.This equation is still recognised in the reference configuration Ω 0 because of the use of the spatial quantity θ and considering dv = JdV.The force equilibrium equation is derived by finding a stationary position of the internal potential energy functional Π(φ, m, s) with respect to arbitrary displacements ζ.

DΠ(φ, m, s)[ζ]
The term ∂Ψ bulk /∂ Ḟ represents the force components acting on the unit area faces of unit volume dV in the reference configuration known as the first Piola-Kirchhoff stress tensor P, that is a work conjugate to the rate of deformation gradient Ḟ27 .Using integration by parts, the balance of linear momentum is resulted as where, DIV( ) is the divergence operator with respect to the reference coordinates X i .The relationship P = Jσ⋅F − T can be used to relate the stresses in the reference configuration Ω 0 to the Cauchy stresses defined in the current configuration Ω t .At this point, two energetically conjugate quantities defined in the reference configuration are introduced as 27 where, C = F T ⋅F is the second-order right Cauchy-Green deformation tensor, S = F − 1 ⋅P is the symmetric second Piola-Kirchhoff stress tensor, and E is the nonlinear Green-Lagrange strain tensor.The concept of the effective stress is adopted in studying the poroelastic behaviour of the porous media, where the Cauchy stress tensor σ affecting the element in the current configuration as where, the tensor σ ′ is the effective Cauchy stress tensor applied on the wetted solid in the current configuration, p is the fluid pore pressure, and I is the identity tensor.A negative value for the effective stress σ ′ implies a compressive state, whereas the compressive fluid pressure is interpreted as a positive value.The Biot's coefficient α is defined as 42 where, K is the drained bulk modulus of the matrix and K s is the bulk modulus of the solid constituents (i.e., grains).Several laboratorial methods can be found in the literature to measure Biot's coefficient α 43 .
In this study, isotropic hyperelastic model is used to model the behaviour of the solid skeleton in a geometrically nonlinear setting.The volumetric response of the solid skeleton subjected to the effective stresses determines the choice of a compressible or nearly incompressible form of the model.In saturated porous rock materials, the volumetric strain of the solid constituent cannot be neglected, 42 so the bulk modulus of solid K s cannot be considered as infinity.According to Ref. 28, the range of variation of α for most of the reservoir rocks are between 0.5 and 0.8, and a compressible hyperelastic material is the appropriate choice to describe their physical behaviour.However, in nearly incompressible materials (α = 1), the volumetric strain is negligible but the volumetric stress on the solid needs to be included in the model to avoid strain locking. 27To provide a comprehensive formulation describing the energy functionals, we present the derivations for both compressible and incompressible hyperelastic material cases.Thus, Equation ( 20) is reformulated as where, σ ′ dev is the deviatoric part of the effective stress, σ ′ vol is the effective volumetric stress affecting on the solid, and p tot = αp + σ ′ vol is interpreted as the total hydrostatic Cauchy stresses affecting on the mixture.The Cauchy stress tensor in the form written in Equation ( 22) is pulled back into the reference configuration as where S ' dev = F − 1 ⋅Jσ ' dev ⋅F − T .Substituting Equation ( 23) into Equation ( 16) and including the external virtual work δW ext , that is assumed independent from the motion, the variation of the total potential energy functional with respect to arbitrary motion ζ is rewritten as The above relationship is indeed the principle of virtual work.From Equations ( 23) and ( 24), the following can be concluded.
To investigate the constituents of the volumetric energy Ψ vol (C, η) caused by the effective pressure (σ ′ vol ), the internal effective potential energy in a drained condition is constructed using the penalty method and treating σ ′ vol as a Lagrange multiplier, 27 notice that the effective stress affecting on the solid skeleton is The variation of the above effective energy functional for the solid skeleton with respect to σ ′ vol reads Enforcing the above equation for all variations of δσ ′ vol results in the effective pressure formulation based on the volumetric strain as 27 Similar to Equation ( 13), the incremental change of the total pressure in terms of the volumetric strain and the fluid content is introduced as 41 Considering the derivatives of the volumetric energy with respect to the fluid content and the volumetric strain in Equations ( 13) and ( 29), the volumetric energy functional is formed below.
By substituting Equation ( 28) into Equation ( 29), the effective volumetric energy due to the solid effective pressure.The fluid energy is obtained by subtracting the effective volumetric energy K(J − 1) 2 /2 from the above equation.
Depending on the inherent properties of the solid skeleton, a compressible or an incompressible hyperelastic material can be chosen.If K s ≫K, an incompressible hyperelastic Neo-Hookean material is used to form the bulk energy as where, Biot modulus M b = K f /φ f and Biot's coefficient α = 1, and Equation ( 30) is used for the volumetric part Ψ vol (C, η).However, if the drained bulk modulus K of the mixture is relatively large compared to the bulk modulus of the solid phase K s , a compressible Neo-Hookean material is chosen to construct the bulk potential energy as where, Equation ( 31) is used for the fluid energy Ψ fluid (C,η).In this case (α < 1), Biot Modulus M b can be defined as 35 1 where, K s and K f are measurable constant parameters of the solid and fluid phases.The evolution of the porosity φ f associated to the change in the deformation causes the bulk response and the fluid pressure evolution to follow a nonlinear state.Using the concept of effective stresses, the principle of virtual work in the reference configuration is written as 27 where, θ = Jp is the Kirchhoff fluid pressure.

The weak forms and constitutive relationships
The energy functional prescribed in Equation ( 9), is modified by affecting a degradation function g(s) = (s 2 +κ ε ) on the tensile part of the solid bulk energy and introducing the fluid energy to the system.Eventually, the internal potential energy of the reference fracturing fluid infiltrated porous body Ω 0 is developed as where, κ ε ≪1 is the regularization parameter, and the effective energy in the solid is split into tensile Ψ + eff (C) and compressive Ψ − eff (C) parts based on the spectral decomposition of the deformation gradient.Due to the complex behaviour and the possible degradations of the overall system stiffness, the energy functional Π(φ, θ, s) of the fracturing saturated porous medium is not always convex with respect to all the variables. 44ence, the Newton-Raphson iterative method is used to solve for (φ, θ) monolithically; then, the crack evolution model is solved for (s) in a staggered manner.

Phase-field crack evolution model in porous media
Following a staggered approach in solving the problem, the energy minimisation (δΠ = 0) with respect to the damage parameter (s) is performed having a fixed state of the deformation and fluid pressure (φ, θ).The evolution of the damage state inside the body is reached by finding a stationary position of the potential energy functional Π(φ, θ, s) with respect to the arbitrary variations of the damage ω ≡ δs as where, Δ is the Laplacian operator with respect to the reference coordinates, and Ĥ (X, t) is a history field to ensure the irreversibility condition for the fractured zone such that where, Ψ + eff (C) is the tensile part of the effective bulk energy related to the fixed state of the deformed solid skeleton, which provides the driving force to breakdown the solid phase.The tensile and compressive parts of the effective bulk energy are separated based on the spectral decomposition of the right Cauchy-Green deformation tensor C, as 27 where N i are the principal directions of the corresponding eigenvalues s 2 i , and the tensile and compressive parts of the stretches are s ± i = 〈s i − 1〉 ± + 1, where 〈x〉 ± = 1 /2(x ±|x|).The tensile part of the effective bulk energy for compressible and nearly incompressible materials are defined in the following, respectively. 27 where J ± = s ± 1 s ± 2 s ± 3 , and ŝ + i = J − 1/3 s + i are distortional stretches.Accordingly, the tensile and compressive parts of the effective Cauchy stress tensor for compressible and nearly incompressible materials are formulated respectively as 27 where n i = (F⋅N i )/λ i and σ ′ i are the principal stresses, defined below for compressible and nearly incompressible solid materials respectively.
Applying Gauss law (Δs = ∇⋅(∇s)) over the reference volume Ω 0 and using the Divergence theorem on Equation (39), the variational form of the phase-field fracture evolution in a total Lagrangian formulation is written as where, l 0 is the length scale parameter, and the zero Neumann type condition GRAD(s)⋅N = 0 is assumed on all boundaries.

Linearised forms of the degraded coupled hydro-mechanical problem
The principle of virtual work presented in Equation ( 37) and the fluid mass conservation law in Equation ( 15) build the governing equations describing the hydro-mechanical behaviour of a fractured saturated porous system.Schematic view of the boundary conditions of a fractured porous system subjected to external loads is depicted in Fig. 2. The functional space of admissible test functions ζ is introduced as where H 1 (Ω 0 ) is the Sobolev functional space of square integrable functions with square integrable weak derivatives.The Dirichlet and Neumann boundary conditions are which represent the external contributions of the total energy.In the following equation, the virtual external work δW ext is replaced by implementing the natural boundary conditions and volume forces.
where, ρ 0 = Jρ is the pull-back mass density, G ≡ g is the acceleration due to gravity in the reference coordinates.The updated Lagrangian form of the principle of virtual work is obtained by pushing forward the constituents into the current configuration using τ ' = F⋅S ' ⋅F T and knowing that grad() = F − 1 ⋅GRAD() as where, τ ′ is the effective Kirchhoff stress tensor, θ is the Kirchhoff pore pressure, and T * = P⋅N represents the natural boundary condition.Utilising the variational principle, the weak form of Equation ( 15) is formed by using the scalar test function ψ belonged to the space of admissible test functions P.
The boundary conditions are defined as θ = θ r on ∂Ω t θ and q⋅n on ∂Ω 0 q (54)   where, n is the unit normal on the natural boundary ∂Ω t q of the current configuration.Substituting the above Dirichlet and Neumann boundary conditions, the weak form of the mass conservation law takes the form Applying the Gauss theorem and integration by parts on the term, Equation ( 55) is rewritten as The pulled-back form of the natural boundary condition V r = JF − 1 v r is substituted in the above equation, and the final variational form of the fluid mass conservation law is formed over the reference configuration in a spatial picture in the following.
The vector Q = − V r ⋅N is the prescribed volumetric rate of flow per unit undeformed area across the boundary ∂Ω 0 .The linearisation of Equations ( 52) and ( 57) are provided in Appendix A. For the undamaged region of the domain, where the Reynold's number is too small, Darcy's law can be used to relate the superficial velocity v r to the fluid pressure as 45 where, ρ f and μ f are the density and dynamic viscosity of the fluid, k is the intrinsic permeability, k = (k /μ f )I, and G is the gravity acceleration vector in the reference configuration.The Forchheimer equation is hired to capture the nonlinear flow within the fully damaged region when the Reynold's number exceeds the limits of linear flow in porous media, the details of which are mentioned in the next section.

Evolution of the fluid transmissivity
Permeability tensor is considered as an evolving quantity conforming with the changes in the volume of the voids in the porous system, caused by formation of the fractures and the deformations in the medium.The flow within the fractured region can be physically assumed as a flow regime between parallel plates, so the well-known Poiseuille's law can be hired mathematically for the modelling.Since the Hagen-Poiseuille equation is theoretically a specific case of Darcy's law by assuming a symmetric permeability tensor, 46 all we need to do is calculating the equivalent symmetric permeability tensor, depending on the openness of the crack, to be used in Darcy's formulation.The following relationship is used as a summation of the symmetric permeability tensor of the solid skeleton k s and the permeability tensor of the cracked zone k cr over a continuous computational domain.
The term (1 − s) ε is used to restrict the increase of permeability on the highly damaged region, and the parameter ε≫1 must be considered significantly large (e.g., ε = 50 has been tested successfully in our analyses) for rock materials to ensure that k cr is only affecting a small portion of the diffusive crack, where the real crack-width is formed.The parameter ε has been introduced as an additional material parameter in Ref. 18, although no experimental method is mentioned for the measurement.It is evident that the permeability remains k = k s for the undamaged region, where (s ≈ 1).The permeability of the solid skeleton may change due to large deformations in highly pressurised regions around the fracture and porosity expansion, especially in weak rock materials where micro-cracks are expected to form. 47This effect can be considered by writing the Kozeny-Carman equation 48 for the permeability value k (k = kI), assuming an isotropic tensor k s = (k)I and a linear Darcy flow regime as where, φ f 0 is the initial porosity, and k 0 = k/μ f is constant depending on the initial undeformed state of the intrinsic permeability k and the fluid dynamic viscosity μ f .For the fractured region, the symmetric permeability tensor of the cracked region k cr is formed as a function of fluid viscosity and the crack openness as where, w is the crack openness and is calculated over the domain using the stretch tensor along the direction perpendicular to fracture. 18The direction perpendicular to the fracture is obtained via N 0 = − GRAD(s)/ |GRAD(s)| in the reference volume, and the spatial stretch vector Λ is calculated from the right Cauchy-Green deformation tensor C. The magnitude of Λ in N direction is recognised by Λ ⊥ and obtained via

GRAD(s)⋅GRAD(s) GRAD(s)⋅C
In the context of Finite Element method, the correct value for fracture opening w is calculated by localizing the objective deformation measure Λ based on the size of the damaged elements.Consequently, the fracture opening is obtained as where, H e ≤ l 0 is the size of the elements for which the damage parameter is lower than the level-set parameter (c≪1).
After the formation of the cracks and significant changes in the fluid velocity within the cracked region, the inertial effect of the laminar fluid flow within the fracture 29 is considered by reformulating Darcy's law.In this study, the nonlinear Forchheimer equation is used to govern the fluid flow inside the crack, when the Reynold's number of the flow inside the cracked region exceeds unity, 49 Re cr > 1.
The Forchheimer equation is formulated as 50 where, β having the unit (m − 1 ) is the non-Darcy flow coefficient.It is evident that the second term adds an extra nonlinearity to the equations, where an iterative approach is required to solve it.Multiplying both sides of the above equation into k revives Darcy's law and the extra the term (F v )v r .The scalar F v , which is called the Forchheimer number, is defined below.
N. Sarmadi and M. Mousavi Nezhad Eventually, we obtain Darcy-Forchheimer formulation as where, the permeability tensor k is defined in Equation ( 59), and the Forchheimer number F v demands subsequent updates in an iterative manner due to the inherent nonlinearity in Equation (66).Notice that Darcy's law in Equation ( 16) is employed for laminar flow in the intact skeleton, whereas Equation ( 66) is activated for the fluid flow in the fractured region when Re cr > 1, representing a laminar flow regime with nonlinear inertial effect.

The linearised weak forms of the coupled poromechanical problem
Following Newton-Raphson method, the time-integrated form of the mass conservation in Equation ( 57) is written as where, Υ is the time-integration parameter, and notice that the terms with no subscript belong to the step t n+1 = t n + Δt.The incremental linearised weak form of the mass conservation law in the direction of δu and δθ is written as The final linearised form of the principle of virtual work in the direction of δu and δθ is where, the fourth-order tensor C (F n , s) and stresses τ ′ n and θ n are updated in each iteration based on the finalised quantities in the previous iteration, which is explained in Appendix A.

Global system of equations
Given a fixed damage field (s) resulted from Equation ( 48), the hydro-mechanical model is solved for (u n+1 , θ n+1 ) using Equations ( 70) and (71) in a monolithic manner.This system is formed as The global tangent stiffness matrix is reached using Equations ( 70) and (71), and the above system is reduced to balancing internal and external forces and fluxes in each step by going through an iteration loop to solve for incremental displacements and Kirchhoff pore pressures.
The right-hand side of the above equation forms the residuals of the forces and fluxes.The finalised answer of each step is obtained by summing up all the incremental values from iterations as below.
The Finite Element matrix forms to solve the hydro-mechanical system in Equation ( 74) and the fracture evolution model in Equation (48) are mentioned in Appendix B for the three fields boundary value problem by approximating u h , θ h , and s h .

Automatic mesh refinement
Since the diffusive fracture is the result of an energy minimisation process between strain effective bulk and surface energies in the finiteelement-based phase-field method, the computational mesh needs to be fine enough to be able to capture the correct amount of energy to be released, forming the fracture over a finite band. 51Therefore, a local mesh refinement algorithm based on h-refinement approach 52 is implemented in the vicinity of the crack tip to eliminate the mesh dependency of the phase-field modelling approach.According to recommendations provided by, 53 the maximum acceptable size of the elements, the surface energy density functional γ(s, ∇s) of which is growing, is set to l 0 /2.
where, H e i and s i are size and damage parameter of the element, and s lim is the critical limit for damage parameter which is chosen as s lim = 0.7 in this study for a reliable calculation of surface energy and stresses around the crack-tip to ensure that mesh refinement has been performed before the rate of surface energy density γ(s, ∇s) growth in the element exceeds unity.To ensure the correctness of mesh refinement, a simple tensile test on a square plate with an initial notch, as presented by Miehe et al. 36 (see Fig. 3a), is simulated to model the fracture propagation from the tip of the initial notch.The change of damage and surface energy density at the crack-tip element for the case of having mesh refinement is compared to a reference fine-mesh structure in Fig. 3b to show the effectivity of the automatic refinement method in capturing the crack-tip.The reference mesh is a fine-structured mesh where the element size in the crack-path is not larger than l 0 /2, and no refinement N. Sarmadi and M. Mousavi Nezhad is applied.A pre-refinement stage has been done to set the maximum element size equal to l 0 when the damage parameter falls under s = 0.9 to increase the quality of the mesh structure.The mesh refinement is performed by splitting the critical elements in the vicinity of the crack-tip by adding an extra node on the longest boundary of the element.A schematic of the refined mesh structure around the tip of the phase-field fracture is depicted in Fig. 3c.
The second refinement criterion is set to deal with stress singularities caused by fracture propagation which results in the numerical errors for displacements and pore fluid pressures.In an updated Lagrangian formulation, the computational mesh evolves with the deformations, and the Jacobian (J) in some elements may become negative, causing the instability of analysis. 54The developed mesh refinement technique must be capable of detecting the critical elements, in which the quality of mesh is poor, and to perform appropriate refinement to guarantee the stability of analysis.The local error estimator is utilised in the iterative solving algorithm to identify the critical elements in the mesh structure.The error estimator rests on an energy norm error criterion for element m in the iteration i is introduced here depending on second Piola-Kirchhoff stress and Green-Lagrange strain as When the displacement field tends to become divergent ( , the total overall error estimator for the total number of elements (Ne) is calculated as The required mesh size for the element m, where singularity happens, can be then estimated based on a simple principle which introduced in Ref. 52 as where, 0 < ϑ < 1 is a priori convergence rate that is recommended as ϑ = 0.5 for the best result, 52 H (i) m is the current mesh size of the element m, and ϱ depends on the local error and is given as This strategy can predict the crack-path and perform the refinement when needed to increase the efficiency of the framework.

Numerical algorithm
Having defined the matrix forms of the global system of the equations, the numerical algorithm is designed and shown in Fig. 4. The hydro-mechanical solver works in an iterative manner for each step (t n+1 = t + Δt) as detailed below.e. Update the solutions: f. Update the element effective stresses, pore pressure, and fluid velocities.g.Check the maximum acceptable errors on the current mesh and perform mesh refinement (if needed) to decrease errors associated to poor mesh quality.h.Check for convergence in residuals: The presented numerical model in this study has been implemented in a MATLAB-based code to model hydraulic fractures in twodimensional plane strain configurations.The details of matrix forms are presented in Appendix B for further information.The developed finite element weak forms and computational framework can be extended to the 3D case easily by developing the relevant matrix forms and developing a 3D mesh structure coupled with a similar mesh refinement technique for tracking the hydraulic fracture over the Employing a staggered approach for coupling the hydro-mechanical model with the fracture model ensures the stability of calculations, 18 noticing that all the simulations in this study have been done using the implicit Newton-Raphson by setting the time-integration parameter Υ = 1.The linear system of equations, Equation (74), is solved using the built-in function "mldivide" in MATLAB.Solving the hydro-mechanical model, Equation ( 74) is mandatory after each cycle of mesh refinement.

Numerical examples
We first assess the ability of the developed computational code in finite strain analysis of the hydro-mechanical behaviour of porous media and the workability of the Darcy-Forchheimer equation introduced to the model.Then, we focus on simulating well-established conceptual scenarios which have been previously simulated by 15,18,19 and on comparing the numerical results to experimental data to examine the validity of our computational framework for modelling of the fluid-driven fracture propagation.

Finite strain hydro-mechanical model
Mandel-Cryer effect, 55 which has been studied experimentally by 56 and numerically by, 57 is modelled to verify the correctness of the hydro-mechanical solver.Geometry and boundary conditions of the model are depicted schematically in Fig. 5a.The model presents a saturated porous medium compressed by two impermeable jacks in vertical direction in a way that the pressurised fluid can be drained from the side boundaries.According to the symmetric plane strain condition, a domain of dimensions a = b = 1 m, made of hyperelastic material with a bulk modulus of K s = 23 MPa, poisson's ratio of 0, and a permeability of 5.08 × 10 − 7 m/s (same as the material properties in Ref. 57) is chosen.Applying symmetry boundary conditions, the uniformly distributed load p 0 is applied on the top boundary and kept constant for the whole duration of the analysis, while the bottom and left boundaries are normally fixed.
The temporal evolution of the horizontal displacements at the right boundary (x = a) and the isochrones of the excess pore fluid pressure along x-axis in various time snapshots are shown in Fig. 5c.The numerical results are compared to the analytical solutions for the case of p 0 = 1 MPa.The analytical solution is based on the infinitesimal strain theory, where the Jacobian of deformation gradient is assumed to be unity, while our analysis is built upon finite strains, and the effect of volumetric changes has been considered.In the nonlinear hydromechanical model developed in this study, the deformations affect the element's stiffness, see Appendix A on constructing the fourth-order elasticity tensor, and the elements that have undergone larger deformations show a softer response.Therefore, the right boundary does not uniformly deform in horizontal direction, and the horizontal displacement of the right boundary increases from bottom to top, since the external loading is applied on the top boundary, in the early stages of consolidation.However, the instant deformation of the right boundary after applying the load (p 0 ) is uniform and exactly matches the analytical solution because the stiffness is not affected by deformations at the beginning of the consolidation process, see Fig. 5b.The overall values of the horizontal displacements of the right boundary reasonably match the elastic analytical solution, ensuring the validity of implementations.Consideration of the finite strains causes capturing the gradual elemental volumetric changes, which delays the rate of the dissipation of the excess pore fluid pressure caused by the compaction of the solid skeleton.This effect becomes more significant when the external loads increase, as can be seen in the plotted results of the time history of pore pressure dissipation in Fig. 6a.By increasing the ratio of the applied load to the bulk modulus of the solid (i.e., p 0 /K s ), the dissipation of the pore fluid pressure is delayed because of considerable amounts of the volumetric strains, as shown in Fig. 6b.This effect has been addressed in Ref. 48, where Mandel-Cryer problem is analysed in 3D.For the smallest value of p 0 /K s , the effect of volumetric strains is negligible, so the results are in a great compatibility with the analytical solution 55 and the linear elastic analyses by 57 based on small strain assumption.

Nonlinear flow in the cracks
In this section, we analyse the fluid flow in a simplified 2D model of a pre-cracked cylindrical sandstone specimen of dimensions 100 mm by 50 mm that was experimentally investigated by Zhang and Nemcik in Ref. 58.A hydraulic head gradient is imposed on the bottom and the top boundaries of the specimen that is laterally confined by σ 3 in a triaxial testing equipment.The mechanical properties (i.e., E and ν) stated in Table 1 are taken from 58 and used for the mechanical properties of the sandstone in this study.For the case of linear flow, that is matched with the experimental results of mated samples in Ref. 58, the hydraulic aperture is considered as equal to 30 μm.The linear Darcy model is calibrated with the experimental data in Ref. 58 to obtain the permeability of the bulk sandstone.The variation of flow rates associated with different values of the hydraulic head gradient presented in Fig. 7a indicates that the developed model captures the reduction of the flow rate due to the increase in the confining stresses σ 3 by correcting the permeability of the cracked region based on the evolving hydraulic Fig. 6.The effect of considering finite strains instead of small strains 55 in the hydro-mechanical response of Mandel's problem; (a) shows the excess pore pressure dissipation in time for an element at location (0,0) for different cases of external loading (p 0 /K s ), and (b) depicts the volumetric strain of the same element in time.

Table 1
Material properties of the granite samples subjected to hydraulic fracturing, borrowed from . 61

Verification of the model with analytical KGD model
The analytical KGD model, introduced by Khristianovic and developed by Geertsma and De Klerk, 4,5 provides closed-form solutions for the length, width, and the internal fluid pressure of a propagating hydraulic fracture in 2D plane strain configuration.In this model, the fracture geometry is considered as an ellipse, and its width is maximum at the centre.For a constant injection rate Q into the hydraulic fracture, the closed-form solutions are formulated as (83) where σ min is the in-situ stress perpendicular to the direction of crack propagation, and other influential parameters are the material elastic properties (E, ν), and the viscosity μ f of a Newtonian fluid.To validate the results of our numerical model, the boundary value problem with an initial notch in the middle of the domain, shown in Fig. 8a, is considered similar to the model used by, 59 and the appropriate boundary conditions due to symmetry are applied to model half of the domain as elaborated in Fig. 8b.The material properties used in this example are chosen similar to the values in Ref. 59, where E = 20 GPa, ν = 0.25, K IC = 1 IC /E.The in-situ stresses (σ max , σ min ) are set to are set to 2 MPa and 1 MPa, respectively, and the injection rate is Q = 0.1 L/s/m.Because of the inherent sensitivity of the phase-field fracture model to the length-scale parameter l 0 24 , the numerical simulations are performed for different values of l 0 , and the results are plotted in Fig. 9. Numerical results show a better compatibility with the analytical solutions as the length-scale parameter (l 0 ) decreases.Setting l 0 as equal to 0.02 m, which is the closest representative to the sharp fracture, results in a sharp decrease of the fluid pressure inside the hydraulic fracture after the fracture initiation from the tip of the initial notch.The captured crack-width (w) for the smallest value chosen for l 0 shows a better agreement with the analytical solution.This conclusion is in line with the study conducted by, 60 who concluded that l 0 →0 revives the exact sharp fracture solution in elastic media, noticing that the KGD model is built upon the assumption of having sharp discontinuities.

The effect of injection rate
We simulated the hydraulic fracture propagation through a twodimensional plane strain granite sample which is subjected to the insitu stresses (σ min , σ max ) and a fluid injection into a circular area in the middle of the sample, so-called the wellbore.Material properties used in the simulations are presented in Table 1, which are chosen from the experimental data in Ref. 61.According to the fundamentals of the linear elastic fracture mechanics, direction of the fracture propagation is the same as the direction of the maximum principal stresses. 62Our results, shown in Fig. 10a, confirm this phenomenon via predicting the hydraulic fracture pattern developed in the direction of σ max .
The experimental study 61 has been conducted over several granite samples by imposing various values of far-field stresses and fluid injection rates into the wellbore.According to Fig. 10b, our numerical results are reasonably in agreement with the experimental observations for the case of far-field stresses σ min = 8 MPa and σ max = 12 MPa and two different fluid injection rates (Q inj ).It is seen that increasing the fluid injection rate causes an increase in the peak value of the fracture fluid pressure (so-called breakdown pressure) in both experimental and numerical results.The temporal changes of the crack-width under different fluid injection rates are depicted in Fig. 10c, which illustrates that increasing the fluid injection rate also results in an increase in the maximum fracture opening and accelerates the speed of propagation, which is seen in the reduction of the required time t cr for fracture to split the sample.As discussed previously, drained bulk modulus (K) can affect the poroelastic behaviour of rock mass subjected to hydraulic fracturing, and its range may vary between 0.2 K s and 0.5 K s , depending on the type of rock. 28To illustrate this effect, we simulated the same BVP for the N. Sarmadi and M. Mousavi Nezhad fluid injection rate of Q inj = 10 mL/min for various values of K. Fig. 10d depicts the effect of drained bulk modulus of the rock mass for granite samples having a constant K s = 30 GPa.A higher value of K increases the required breakdown pressure significantly, due to a higher capability of the rock matrix to undertake volumetric expansion implied by fluid pressure and accelerates the development of hydraulic fracture throughout the medium.

Hydraulic fracture propagation and the nonlinear effect of flow in the fracture
The ability of the developed model to simulate the fluid driven fractures in porous media is examined by simulating a two-dimensional plane-strain model of the hydraulic fracture propagation caused by injection of the pressurised water into a wellbore embedded in a cubic sandstone sample which has been originally discussed in Ref. 63 and later simulated using phase-field theory in Ref. 15. Geometry and boundary conditions of the model are shown in Fig. 11.Values of the numerical and material parameters used for simulations are taken from 15 and stated in Table 2.In the numerical analysis presented in Ref. 15, a Dirichlet-type boundary condition was applied on the left edge to model the linear increase in the fluid pore pressure, causing the fracture to propagate.In contrast to the Dirichlet-type increasing boundary condition, a Neumann type boundary condition describes the sudden fluid pressure drop at the fracture tip due to the volumetric expansion during the fracture evolution. 19Therefore, a Neumann-type boundary condition with the fluid injection rate of 6.25 mL/s is imposed on the middle of the left edge of the simulation domain as shown in Fig. 10, while an initial pore fluid pressure of and confining pressure of 8.5 MPa are assigned to the domain.Due to symmetry, the drainage boundary is only applied on the boundaries (1), (2), and (4), and impermeable condition is applied to the boundary (3).Extra compressive external stresses σ y = − 4 MPa and σ z = − 1 MPa are applied on the boundaries (1) and (2), while the boundaries (3) and ( 4) are fixed in y and z directions, respectively.
Fig. 12a and c depict the variation of the fluid pressure and the damage parameter with time at four check points along the fracture path which are located at the distances (in meters) of d = 0.02, 0.05, 0.08, and 0.17 from the left boundary.Injection of the fluid into the system increases the fluid pressure to the point of breakdown when the fracture propagation starts.It is followed by a sudden drop in the fluid pressure at the fracture tip due to volumetric expansion, the evolution of pore spaces, and the increase of permeability in the damaged solid skeleton.The values of predicted fluid pressure are compared to those presented by 15 in Fig. 12d for three check points, where a good agreement is seen between the results.The discrepancy between the results is due to the inverse pressure effect at the fracture tip, 14,63,64 which has not been addressed in Ref. 15.As illustrated in Fig. 12b, the fluid pressure in the porous zone on the right-hand-side of the fracture tip undergoes a decrease just before fracturing, which is implied by the high fluid pressure gradient around the fracture tip.This effect is more notable when the fluid pressure inside the fracture is significantly higher than the fluid pressure in the pore zone adjacent to the crack-tip.By propagation of the hydraulic fracture, the fracture fluid pressure decreases, while the surrounding regions remain pressurised.
The value of the non-Darcy coefficient β, see Equation ( 67), which theoretically depends on the width and roughness of the fractures, 65 is set to a non-zero value in the simulations to consider the nonlinearity in the fluid flow inside the fracture.As there has been no measurned value for this parameter for the senario we investigated in this section, a sensitivity analysis is performed using a range of values for β, and the effects of which on the fluid pressure inside the fracture are shown in Fig. 13a.The activation of the Darcy-Forchheimer model and consideration of non-Darcy flow effect are the result of having the Reynold's number (Re p ) higher than the linear Darcy-flow limit, which has been described in section 3.3.A realisation of the Reynold's number in the hydraulic fracture is shown in Fig. 13b, and its distribution among the fully damaged elements is depicted in Fig. 13c.The hydraulic gradients generated between the pressurised regions and the fracture result in formation of "flow-back" and "leak-off" flow regimes.The flow-back regime is formed along the hydraulic fracture, where the fluid flows from the pressurised porous medium into the fracture, and a leak-off regime is formed around the fracture tip, when the fluid leaks into the low-pressure porous zone located in the vicinity of the fracture, see Fig. 14.It is known that flow-back of the fracture fluid has significant effects on fracture closure and capability to predict the flow-back rate helps gaining knowledge on the production rate of hydraulic fracturing operations.It is understood that the difference between the amount of fracture fluid pressure in Darcy and non-Darcy cases is higher when the fracture length is short and decreases as the length grows (i.e., the fracture propagates).This effect is a consequence of having a larger velocity field inside the shorter fracture compared to that of the longer fracture while the injection rate is kept constant.For comparison purposes, the fluid flow regimes around the fracture for the Darcy case (β = 0) and non-Darcy case (β = 3E+4) are depicted in Fig. 14.It is found out that non-Darcy flow in the fracture significantly influences the extent of the flow-back region and the amount of flow-back rate.Based on the results in Fig. 15, it is concluded that considering the Darcy-Forchheimer equation, capturing the effect of nonlinear inertial flow inside the fracture because of high Reynold's number of the flow, accelerates the formation of the flow-back region and increases the flow-back rate significantly.

Conclusion
We developed a computational framework to simulate hydraulic fracture propagation in a geometrically nonlinear setting.The framework has been built upon a series of mathematical formulations derived based on the energy minimisation of a saturated porous system and solved to calculate deformations, fluid pressure, and the damage state of the system in a robust staggered manner on a continuous finite element mesh.The updated Lagrangian finite element formulation has been solved using an iterative Newton-Raphson method to consider both  stiffness degradation and porosity evolution due to fracture propagation and large deformations in the medium.A comprehensive mathematical description of the energy functionals and their spectral decomposition of the solid bulk energy for compressible and nearly incompressible hyperelastic materials has been provided, which can be used for the nonlinear finite element implementations.Procedure for finite element implementation of the mathematical formulations has been discussed.To increase the efficiency of the computational framework for modelling phase-field hydraulic fracture, a h-refinement technique has been implemented to refine a coarse mesh structure as the fracture propagates.Including the mesh refinement in the updated Lagrangian framework helps to avoid any possible mesh distortion and instability of nonlinear analysis.Our numerical simulations confirm that the developed model predict the propagation path of the hydraulic fractures, the fluid flow within a fractured system, the fluid leak-off into the porous solid, and the flow-back from the stimulated regions into the hydraulic fracture.The results of this study confirm that the consideration of nonlinear flow inside the fractures contributes to developing a lower fluid pressure gradient across the fracture walls, which can significantly affect the amount of fluid to be exchanged at the fracture boundaries.It has been concluded that in some scenarios, a simplified linear fracture flow model underestimates the amount of stimulated flow-back by almost 60%.The numerical results were compared to experimental data of hydraulic fracture propagation in granite samples subjected to various fluid injection rates, and a strong similarity was found in the fracture fluid pressure variation between our simulations and the experimental observations, showing the validity and capability of the developed framework.The variation of fluid injection rate applied into the wellbore and the drained bulk modulus of the rock mass were recognised to be influential on the breakdown pressure, the cracking speed, and the fracture openness, which were highlighted in this study.The matrix form of the weak form of the fluid mass equilibrium in a time-integrated format in the following.
The matrix [B θ ] is the derivative of shape functions corresponding to Kirchhoff pore pressure with respect to the current spatial coordinates in the last iteration.71) and is written as: Finally, the finite element matrix form of the crack evolution problem in Equation ( 46) in a total Lagrangean format is written as where, [B s ] is the derivative of shape functions with respect to the undeformed material coordinates. [

Fig. 1 .
Fig. 1.(a) Is the motion of a fracturing porous system from the reference to the current configuration; (b) is the phase field damage parameter function to represent the crack boundary.

Fig. 2 .
Fig. 2. The boundary value problem for the hydro-mechanical process in the reference and current configurations.

1 . 2 .
Setting the boundary conditions and forming external forces R EXT t+Δt and fluxes H EXT t+Δt .Solve for u n+1 and θ n+1 iteratively, assuming a fixed damaged state of the body.a. Form the internal forces F (i− 1) n+1 and fluxes Q (i− 1) n+1 , right-hand side of Equation (74).b.Update the Forchheimer number F (i− 1) v based on the updated fluid velocities.c.Update the global stiffness matrix using Equations (70) and (71).d.Solve for Δu (I) n+1 and Δθ (I) n+1 using the global system mentioned in Equation (74).

Fig. 3 .
Fig. 3. (a) Is the boundary value problem to model a simple fracture propagation from the initial notch, as presented in Ref. 36; (b) shows the evolution of damage and surface energy at the crack-tip element in the mesh structure with and without automatic mesh refinement; (c) depicts a schematic of the arrangement of nodes and element sizes in the vicinity of the phase-field crack-tip.

Fig. 4 .
Fig. 4. Numerical algorithm flowchart, solving for three-fields problem in a nonlinear finite element framework, including an iterative mesh refinement scheme.

Fig. 5 .
Fig. 5. (a) Shows a schematic of the symmetric Mandel's problem; (b) and (c) compare the numerical results to the analytical solutions for the horizontal displacement of the right boundary and the excess pore pressure along the x-axis with respect to time.

− 1 N
. Sarmadi and M. Mousavi Nezhad aperture of the cracked zone.The mismatch between numerical results and experimental data for low hydraulic head gradient (1 MPa) in Darcy-flow regime can be related to the fluid loss in the tests, as the flow rates are generally low.The non-Darcy behaviour of the flow within cracks is then investigated by increasing the hydraulic aperture to 0.1 mm for non-mated samples in Ref. 58, causing the generation of a larger amount of flow rates within the fracture.The coefficient β in the Forchheimer equation is adjusted to the model based on the roughness of the fracture walls in non-mated samples, reported in Ref. 58.Fig. 7b shows that the proposed numerical model, benefitting from the mixed Darcy-Forchheimer equation, can reproduce the experimentally measured flow rate-hydraulic gradient curves for the flow within the non-Darcy flow regime.Our results confirm that the flow rate and the nonlinearity of the flow behaviour decrease with the increase of the confining stresses σ 3 , which is in line with the experimental observations of Ref. 58.

Fig. 7 .
Fig. 7. Modelling linear and nonlinear flow regimes within a pre-defined fracture, shown in the schematics.(a) Shows the variation of the flow rate with respect to hydraulic head gradient in Darcian flow regime; while (b) shows the nonlinear behaviour of flow inside the fracture with larger hydraulic aperture equal to 0.1 mm.Numerical results are depicted in solid lines and are compared to experimental data (marker points) presented by.58

Fig. 8 .Fig. 9 .
Fig. 8. (a) Is the general configuration of the fluid injection into an initial notch, and (b) shows the boundary value problem for the model verification with regards to the KGD model; solid line represents the sealed boundary and dashed lines are drainage boundaries.

Fig. 10 .
Fig. 10.(a) Shows the geometry of the modelled granite sample (dimensions in mm) and the direction of the hydraulic fracture propagation.(b) Compares the numerical results of the fracture fluid pressure time histories to the experimental data provided in Ref. 61 for the case of σ min = 8 MPa and σ max = 12 MPa.(c) Shows the effect of increasing the injection rate on the fracture opening and the length of cracking time (t cr ).(d) Elaborates the effect of drained bulk modulus of the rock (K) on the breakdown pressure, on t cr , and on the variation of fracture fluid pressure.

Fig. 11 .
Fig. 11.Schematic of the plane-strain hydraulic fracture simulation subjected to the flow injection in Example 5.5, according to Refs.15,63.

Fig. 12 .
Fig. 12.(a) and (c) show the evolution of fluid pressure and the damage along the crack path from the injection point to the opposite outer boundary.(b) Is the magnified part of (a) since the moment crack starts propagation.(d) Compares the numerical results of the fluid pressure to the data provided by 15 with respect to the total fracture length formed.

Fig. 13 .
Fig. 13.(a) Shows the change of the fracture fluid pressure with respect to the length of the hydraulic fracture formed for different values of β; (b) depicts a realisation of Reynold's number (Re) inside the hydraulic fracture in a time snapshot, and (c) shows the distribution of the captured Re in the crack region (fully damaged elements).

Fig. 14 .
Fig. 14.Comparison between the flow fields predicted by the Darcy and non-Darcy laws: (a) and (c) show the flow-back and leak-off around the fracture channel for the Darcy flow case (β = 0) when the total fracture length is 0.30 m and 0.45 m respectively; (b) and (d) are related to the nonlinear flow case, when the fracture length is 0.30 m and 0.45 m respectively.

Fig. 15 .
Fig. 15.Consideration of the nonlinear flow inside the fracture (by choosing a nonzero β) intensifies the flow-back into the fracture.(a) and (b) depict the development of the length of the flow-back region and the rate of flow-back, respectively, with respect to length of the fracture formed.

Table 2
15terial and numerical parameters in the hydraulic fracture simulation, borrowed from .15 B s ] = N. Sarmadi and M. Mousavi Nezhad