# Iteratively Coupled Flow and Geomechanics in Fractured Poroelastic Reservoirs: A Phase Field Fracture Model

Authors: Jinzhou Zhao; Qing Yin; John McLennan; Yongming Li; Yu Peng; Xiyu Chen; Cheng Chang; Weiyang Xie; Zhongyi Zhu
Journal: Geofluids
 (2021)

---

## Abstract

Fluid-solid coupling in fractured reservoirs plays a critical role for optimizing and managing in energy and geophysical engineering. Computational difficulties associated with sharp fracture models motivate phase field fracture modeling. However, for geomechanical problems, the fully coupled hydromechanical modeling with the phase field framework is still under development. In this work, we propose a fluid-solid fully coupled model, in which discrete fractures are regularized by the phase field. Specifically, this model takes into account the complex coupled interaction of Darcy-Biot-type fluid flow in poroelastic media, Reynolds lubrication governing flow inside fractures, mass exchange between fractures and matrix, and the subsequent geomechanical response of the solid. An iterative coupling method is developed to solve this multifield problem efficiently. We present numerical studies that demonstrate the performance of our model.


---

## Body

## 1. Introduction

In subsurface energy extraction engineering disciplines, there is a growing interest in developing a clear understanding of coupled processes between rock deformation and fluid flow in a deformable porous media. These coupled processes are relevant for modeling reservoir stimulation, waste disposal, geothermal energy recovery, and CO2 sequestration [1]. This fluid-solid coupled problem is especially significant in naturally fractured or stress-sensitive reservoirs, due to the evolution of fracture aperture and permeability with time [2]. However, documentation of effective modeling and simulation of a complex discrete fracture network (DFN) in geological media is difficult due to the presence of discontinuities across discrete fractures. Furthermore, accounting for different governing flow laws in the matrix and in DFN-dominated systems has also not been well studied.The importance and complications of DFN flow performance have driven the development of increasingly more efficient and accurate models that are able to deal with coupling of flow and geomechanical response for fractured reservoirs. The dual porosity/permeability (DP) approach has been typically used for the last few decades [3–5]. This dual continuum concept was conceived on the premise that a fracture network with high conductivity acts as an interconnected system of secondary porosity within a preexisting porous matrix, accounting for the primary porosity. This approach is not applicable for large-scale fracturing since the fracture systems are dense, according to the inherent assumptions of the model. Moreover, the shape factors that need to introduced in this framework to define the flow transfer between fractures and matrix are not straightforward to determine.For fluid-solid coupled problems in a porous medium, at a minimum, it is necessary to explicitly represent large-scale fractures [6, 7]. As a result, the discrete fracture model (DFM) was developed. This approach is attractive for large-scale and discrete fracture systems. In the DFM approach, lower n−1-dimensional or equal-dimensional representation of fractures can be used and the fracture mesh should be small enough to coincide with the matrix at the interface [8–10]. Therefore, locally refined and unstructured grids are required to conform to the explicit fracture geometry [11]. This mesh refinement and conformity significantly increases the computational cost and becomes impractical for a large number of fractures.To overcome the requirement for complex meshing associated with DFM, other techniques based on extended finite element method (XFEM) for embedded fractures in a nonconforming mesh have been developed [12–15]. In this method, a set of enrichment functions is used to approximate the near-tip behavior of fractures [16]. XFEM does not require a fracture to be aligned with element boundaries and can achieve good precision without requiring a highly refined mesh. However, different enrichment functions need to be used depending on the different fracture patterns inside an element. In other words, different fracture patterns need to be recognized and treated explicitly by XFEM. This is cumbersome especially in three dimensions, and this also leads to significant computational cost [17].The phase field approach is based on the variational formulation proposed by Francfort and Marigo [18]. It has become numerically attractive in recent years due to its ability to deal with complex fractures and easy implementation. The phase field approximation regularizes a sharp fracture surface by a diffusive or smeared-out fracture zones governed by a scalar auxiliary variable in a fixed-grid topology. This straightforward modeling provides a significant advantage over the discrete fracture description, whose numerical implementation requires explicit (e.g., DP methods) or implicit (e.g., DFN and XFEM methods) handling of the discontinuities.Seminal extensions have been proposed to capture the effect of pressurized fractures in porous media [19–21]. The phase field fracture model and an existing reservoir simulator have been coupled [22]. In that case, simulations cannot match the benchmarks since an approximated permeability function is used rather than Reynold’s lubrication theory. A fully coupled model has been proposed by others to describe fluid flow in a fracture using the Reynolds lubrication equation with the surrounding poroelastic medium following the Biot poroelasticity model [23]. However, this model is limited to a single fracture since there is a requirement for two separate meshes to track both fractures and fluid flow.In practical implementations, there are two main numerical schemes. Different monolithic algorithms have been proposed to solve the nonconvex and nonlinear phase field formulation for pressurized fractures [24, 25]. That framework has been improved and extended to coupled problems by Lee et al. [26, 27]. An alternative thermodynamically consistent framework of phase field, based on operator splits [28, 29], has been widely used because its implementation is very simple and robust. In that scheme, the fracture problem has been decoupled as several Euler governing equations and it can be implemented in a straightforward manner to couple multifield problems [30–32]. Ambati et al. [33] further improved this scheme by transferring it to a linear problem; this enables a significant reduction in computational cost. Generally, recent theoretical advances and practical application of phase field methods demonstrate that it has remarkable potential for tackling complex fracture network problems with multifield coupling.In the present paper, we further develop a fluid-solid fully coupled phase field model of a poroelastic DFN-matrix system. This model accounts for Darcy-Biot-type fluid transport in porous media while Reynolds lubrication is used to govern flow in discrete fractures. To circumvent the issues of tracking fluid flow in the fractures, the lower-dimensional lubrication equation in the fracture domain has been converted (regularized) to the whole computational domain by using a phase field variable. Then, an arbitrary distribution of fracture networks can be considered. As a result, cumbersome configurations, such as having two separate mesh systems, are also avoided. Under this unified framework, the flow transfer between fractures and porous media has been intrinsically modeled and no additional tuning parameters such as the shape factors or Carter’s leak-off coefficient are required. Therefore, there is no need for additional techniques to deal with this. We then incorporate an extended fixed stress iterative coupling algorithm for fluid-filled fracture networks. Multifield problems can be readily considered using this framework due to the decoupled governing equations. To validate the phase field approach adopted here, we compared our model predictions with two classical analytical solutions. Along the way, we assess the performance of our model by simulating fully coupled flow in a complex discrete fracture network.This paper focuses on flow and geomechanics fully coupled modeling in fractured poroelastic reservoirs; therefore, the simulation of fracture propagation is out of scope. The paper is organized as follows. In Modeling, the displacement-pressure-phase field modeling equations and the extended fixed stress iterative coupling implementation are described. Model verification and numerical tests are presented in Numerical Tests. Finally, a summary and some conclusions are presented in Conclusion.


## 2. Modeling

We consider an arbitrary poroelastic bodyΩ⊂ℝd with spatial dimension d∈1,3, composed by a porous medium domain Ω\Γ and an internal pressurized fracture domain Γ (see Figure 1). Dirichlet and Neumann boundary conditions are imposed on ∂DΩ and ∂NΩ, respectively.Figure 1
Schematic representation of a poroelastic domain with an internal pressurized fracture domain.


### 2.1. Phase Field Fracture Model

The physical processes in the fractured domainΓ and the porous medium domain Ω\Γ are quite different. However, it is cumbersome to compute these two types of domains separately, especially for complex fracture networks. To circumvent the difficulty of tracking fracture surfaces, the phase field model is introduced to redefine all the subdomains over the whole domain.The scalarϕx∈0,1 is called the fracture phase field variable, and it can be used to describe the regularized fracture topology.
(1)ϕ≔0onΩ\Γ,1onΓ.It is worth noting that the length scale parameterℓ∈ℝ+ determines the width of the smooth approximation of the fracture (see Figure 2). Therefore, the sharp fracture is represented by a regularized or diffusive fracture. The length scale controlled area is called the “diffusive fracture zone,” which can be interpreted as a microcrack zone. Based on the theory of Γ-convergence [34], it can be proved that Γℓ converges to Γ as ℓ⟶0.Figure 2
Regularization of the discontinuity fracture domain by the phase field variableϕ.The phase field governing equations are derived using the variational principle [35].
(2)ϕ−ℓ2Δϕ=0,∇ϕ·n=0on∂NΩϕ.Next, we can introduce∇ϕ to convert integrals over the fracture domain Γ to the whole domain Ω.
(3)∫ΓxdΓ=∫Ω∇ϕfxdV,where fx is an arbitrary function.Similarly,1−ϕ2 is used to convert integrals over the porous medium domain Ω\Γ to the whole domain Ω.
(4)∫Ω∖ΓfxdV=∫Ω1−ϕ2fxdV.It is better to err on the side of caution when obtaining the fracture outward normal since it is not always the same as the mesh’s normal, especially for a curved fracture (see Figure3(a)).Figure 3
Fracture outward normal. (a) Mesh and fracture outward normals on the fracture surface. (b) Negative phase field gradient vector.
(a)(b)Figure3(b) shows that the direction of the fracture outward normal is consistent with the gradient of the phase field variable. Consequently, the fracture outward normal can be defined by
(5)nΓ≈−∇ϕ∇ϕ.


### 2.2. Geomechanical Governing Equations

We assume that the strains are small, giving the infinitesimal strain tensor as(6)εu=∇su≔12∇u+∇Tu.The constitutive equation, based on standard linear elasticity for isotropic solids, is(7)σeε=2με+λtrεI,where the Lamé coefficients μ and λ are given by
(8)μ=E21+ν,λ=Eν1+ν1−2ν,and they are expressed in terms of Young’s modulus E and Poisson’s ratio ν.Following the rock mechanics conventions, both the tensile stress and fluid pressure are considered positive in compression. For the whole domain with the regularized fractures, the constitutive equation for the total poroelastic stress tensor is formed as(9)σ=1−ϕ2σe+αpI,where α is Biot’s coefficient. Note that when the phase field variable ϕ⟶0, which represents the state for the medium, this constitutive equation Equation (9) converges to the conventional effective stress equation σ=σe+αpI.Finally, the linear momentum balance is expressed as(10)∇·σ+f−p∇ϕ=0,σ·n=t¯on∂NΩu,where the additional term p∇ϕ can be interpreted as fluid pressure forces that act on the fracture surfaces to increase/decrease the fracture volume. See the detailed derivation of the governing equation in appendix.


### 2.3. Fluid Flow Governing Equations

Assuming that fluid pressure is continuous from fracture to reservoir in the poroelastic medium, the fracture pressurepf and the pore pressure pr are related to p, the fluid pressure acting in the whole domain, as follows
(11)p=prinΩ\Γ,p=pfonΓ.Since the fracture width is much less than its length, i.e.,w<<L, incompressible single phase fluid flow within a fracture is governed by Reynolds’ lubrication equation
(12)ẇ−∇Γ·w312μf∇Γp+qL=qinjinΓ.Equation (12) is formulated in terms of the fracture surface tangential gradient operator ∇Γ=∇−nΓnΓ·∇ which eliminates the contribution of the normal component that is perpendicular to the fracture surfaces. qinj is the injection flow rate, qL is the leak-off rate, and μf is the fluid’s dynamic viscosity.The fracture opening width is estimated by(13)w=nΓ·εnΓhe,where he is the local mesh size and nΓ can be computed from Equation (5).Inspired by Biot’s theory and the fixed stress iterative approach [36, 37], the mass balance equation for incompressible single phase flow in the poroelastic medium is expressed as
(14)1M+α2Kdrṗ−∇·kμf∇p=qrs−α∇·u̇+α2KdrṗinΩ\Γ,where the additional term α2/Kdrṗ improves the numerical stability. qrs is the reservoir source/sink term. M, Kdr, and k are Biot’s modulus, the drained bulk modulus, and the permeability tensor, respectively. The relationships for these parameters are given by
(15)Kdr=E1−ν1+ν1−2ν,1M=α−ϕporKm+ϕporKf,where ϕpor is the matrix porosity. Km and Kf are the matrix and fluid bulk modulus, respectively.


### 2.4. Summary of Governing Equations

Finally, this highly nonlinear problem has been split into a set of linear governing equations and can be solved by an iterative coupling method. Findu and p for all time t with fully coupled governing equations as follows:
(16)∇·σ+f−p∇ϕ=0inΩ,ϕ−ℓ2Δϕ=0inΩ,ẇ−∇Γ·w312μf∇Γp+qL=qinjinΓ,1M+α2Kdrṗ−∇·kμf∇p=qrs−α∇·u̇+α2KdrṗinΩ∖Γ,σ·n=t¯on∂NΩu,∇ϕ·n=0on∂NΩϕ,kμf∇pk+1·nΓ=0on∂NΩp,kμf∇pk+1·nΓ=qLonΓ.This covers the stress equilibrium, the fluid mass balance, and the fracture regularization. The unknown leak-off termqL has been implicitly modeled as the fluid flux from the fracture into the reservoir. Therefore, there is no need for legacy fluid leak-off models, such as Carter’s formulation [38].Another implementation challenge in Equation (16) is that the lubrication equation is defined over the fracture domain Γ, and the mass balance equation in the poroelastic medium is defined over the porous medium domain Ω\Γ. Two separate mesh systems, for the fracture domain and the porous domain, respectively, can be used to solve these two domains separately. However, a cumbersome fracture tracking technique is required for this case. To circumvent this issue, the phase field calculus can be used (Equations (3) and (4)). This will be explained in the following section.


### 2.5. Weak Formulation

Assign three weighting functions:(17)φu∈φu∈H1Ω∣φu=0on∂DΩu,φϕ∈φϕ∈H1Ω∣φϕ=0on∂DΩϕ,φp∈φp∈H1Ω∣φp=0on∂DΩp.Multiplying Equation (16) by the weighting functions, integrating over Ω, and performing integration by parts. Those equations then become the following:
(18)∫Ω∇φu:σdV=∫Ωφu·fdV−∫Ωφu·p∇ϕdV+∫∂NΩuφu·t¯dS,(19)∫ΩφϕϕdV+ℓ2∫Ω∇φϕ·∇ϕdV=0,(20)∫Γ∇Γφp·w312μf∇ΓpdΓ=∫Γφpqinj−qL−ẇdΓ,(21)∫Ω∖Γφp1M+α2KdrṗdV+∫Ω∖Γ∇φpkμf∇pdV=∫Ω∖ΓφpqrsdV−α∫Ω∖Γφp∇·u̇dV+α2Kdr∫Ω∖ΓφpṗdV+∫ΓφpqLdΓ.The integration of the fracture flow is perfomed over the fractured domainΓ while the integration of the reservoir flow is performed over the unfractured domain Ω/Γ. The consequence is that the aforementioned domain treatments—Equations (3) and (4)—can be used for fluid flow integrals. In addition, the unknown leak-off term qL is immediately eliminated by adding Equations (20) and (21) together. This yields
(22)∫Ωφp1M+1−ϕ2α2KdrṗdV+∫Ω∇φpkμf∇pdV+∫Ω∇ϕ∇Γφp·w312μf∇ΓpdV=∫Ω1−ϕ2φpqrsdV−α∫Ω1−ϕ2φp∇·u̇dV+1M∫Ω1−1−ϕ2φpṗdV+α2Kdr∫Ω1−ϕ2φpṗdV−∫Ω∇ϕφpẇdV+∫Ω∇ϕφpqinjdV.Equations (18), (19), and (22) constitute a coupled bilinear system. In our implementation, u, ϕ, and p are discretized with quadrilateral finite elements in space and a backward Euler scheme in time.


### 2.6. Numerical Implementation

#### 2.6.1. Finite Element Discretization

The finite element approximations of the solid displacement field the fracture phase field and the fluid pressure field are interpolated in terms of the nodal displacementsu¯, nodal phase field ϕ¯, and nodal pressure p¯, respectively
(23)u≈uhx=Nuu¯,ε≈Buu¯,ϕ≈ϕhx=Nϕϕ¯,∇ϕ≈Bϕϕ¯,p≈phx=Npp¯,∇p≈Bpp¯.where Nu denotes the matrix of shape functions and Nϕ and Np are the row vector of shape functions. For 2D elements with n nodes, they are given by (assume the same shape functions Nn are adopted for u, ϕ, and p without loss of generality)
(24)Nu=N10N20⋯Nn00N10N2⋯0NnNϕ=Np=N1N2⋯Nn.,In a 2D settings, the strain-displacement matrixBu and the differential matrix Bϕ, Bp have the following components:
(25)Bu=N1,x0N2,x0⋯Nn,x00N1,y0N2,y⋯0Nn,yN1,yN1,xN2,yN2,x⋯Nn,yNn,xBϕ=Bp=N1,xN1,x⋯N1,xN1,yN1,y⋯N1,y.,Substituting Equation (23) to the weak form Equation (18) and approximating each cells contribution by Gauss quadrature rule, the linear system is obtained:
(26)Kuuu¯=Fu−Kupp¯.The tangent stiffness matrixKuu, the coupling matrix Kup, and the external force vector Fu are given by
(27)Kuu=∫Ω1−Nϕϕ¯2BuTDBudV,Kup=∫Ω1−Nϕϕ¯2αBuTmNpdV,Fu=∫∂ΩuNuTt¯dS+∫ΩρNuTfdV,where D is the elasticity tensor, and m=110T for 2D.Similarly, the matrix-vector form of the phase field Equation (19) is obtained as follows:
(28)Mϕ+Kϕϕ¯=0Mϕ=∫ΩNϕTNϕdV,Kϕ=∫Ωℓ2BϕTBϕdV.Finally, the matrix-vector form of the fluid flow Equation (22) is given by
(29)Mprp¯̇+Kpp+Kpfp¯=Q−Kpuu¯+Hprp,¯̇where the
(30)Mpr=∫Ω1M+1−Nϕϕ¯2α2KdrNpTNpdV,Kpp=∫ΩkμfBpTBpdV,Kpf=∫ΩBϕϕ¯w312μfBpTBpdV,Q=∫ΩNpT1−Nϕϕ¯2qrsdV−∫ΩNpTBϕϕ¯ẇdV+∫ΩNpTBϕϕ¯qinjdV,Kpu=∫Ω1−Nϕϕ¯2αNpTmTBudV,Hpr=∫Ω1M1−1−Nϕϕ¯2NpTNpdV+∫Ωα2Kdr1−Nϕϕ¯2NpTNpdV.


#### 2.6.2. Iterative Coupling Algorithm

Finally, an extended fixed stress iterative coupling method for fracture networks has been introduced to solve three matrix-vector systems Equations (26), (28), and (29). The central idea of this algorithm is to hold the other two variables constant to yield three extremely simple linear subsystems:
(31)Kuuu¯k+1=Fu−Kupp¯k+1,(32)Mprp¯k+1+Kpp+Kpfp¯k+1=Q−Kpuu¯k+Hprp¯k,(33)Mϕ+Kϕϕ¯=0.The numerical implementation is presented in Algorithm1. k represents the iteration levels. uk+1, pk+1 represent unknown solutions in current iteration, while uk, pk are known solutions from previous iteration.Algorithm 1: Extended fixed stress iterative coupling method.
1 Set phase field Dirichlet boundary conditionsϕ=1 for pre-existing fracture networks;2 Set initial conditions;3 Solve the phase field formulation Equation (33) for ϕ;4 At the timetn;5repeat6  Update the fracture outward normalnΓk using Equation (5);7  Update the fracture widthwk using Equation (13);8  Solve the pressure formulation Equation (32) for pk+1 using ϕ, pk, uk, nΓk, and wk;9  Solve the displacement formulation Equation (31) for uk+1 using ϕ and pk+1;10untilRp<ϵtol;11 increment the timetn⟶tn+1.According to our numerical tests, the convergence of pressurep is harder to reach than displacement u due to strong nonlinear of fluid flow. Therefore, the global L2 error of pressure Rp=pk+1−pk is sufficient to examine both the convergence of the total coupling system. In general, the fixed stress iterative coupling algorithm is very robust (unconditional stable) and efficient (the number of iterations is less than 10).


## 2.1. Phase Field Fracture Model

The physical processes in the fractured domainΓ and the porous medium domain Ω\Γ are quite different. However, it is cumbersome to compute these two types of domains separately, especially for complex fracture networks. To circumvent the difficulty of tracking fracture surfaces, the phase field model is introduced to redefine all the subdomains over the whole domain.The scalarϕx∈0,1 is called the fracture phase field variable, and it can be used to describe the regularized fracture topology.
(1)ϕ≔0onΩ\Γ,1onΓ.It is worth noting that the length scale parameterℓ∈ℝ+ determines the width of the smooth approximation of the fracture (see Figure 2). Therefore, the sharp fracture is represented by a regularized or diffusive fracture. The length scale controlled area is called the “diffusive fracture zone,” which can be interpreted as a microcrack zone. Based on the theory of Γ-convergence [34], it can be proved that Γℓ converges to Γ as ℓ⟶0.Figure 2
Regularization of the discontinuity fracture domain by the phase field variableϕ.The phase field governing equations are derived using the variational principle [35].
(2)ϕ−ℓ2Δϕ=0,∇ϕ·n=0on∂NΩϕ.Next, we can introduce∇ϕ to convert integrals over the fracture domain Γ to the whole domain Ω.
(3)∫ΓxdΓ=∫Ω∇ϕfxdV,where fx is an arbitrary function.Similarly,1−ϕ2 is used to convert integrals over the porous medium domain Ω\Γ to the whole domain Ω.
(4)∫Ω∖ΓfxdV=∫Ω1−ϕ2fxdV.It is better to err on the side of caution when obtaining the fracture outward normal since it is not always the same as the mesh’s normal, especially for a curved fracture (see Figure3(a)).Figure 3
Fracture outward normal. (a) Mesh and fracture outward normals on the fracture surface. (b) Negative phase field gradient vector.
(a)(b)Figure3(b) shows that the direction of the fracture outward normal is consistent with the gradient of the phase field variable. Consequently, the fracture outward normal can be defined by
(5)nΓ≈−∇ϕ∇ϕ.


## 2.2. Geomechanical Governing Equations

We assume that the strains are small, giving the infinitesimal strain tensor as(6)εu=∇su≔12∇u+∇Tu.The constitutive equation, based on standard linear elasticity for isotropic solids, is(7)σeε=2με+λtrεI,where the Lamé coefficients μ and λ are given by
(8)μ=E21+ν,λ=Eν1+ν1−2ν,and they are expressed in terms of Young’s modulus E and Poisson’s ratio ν.Following the rock mechanics conventions, both the tensile stress and fluid pressure are considered positive in compression. For the whole domain with the regularized fractures, the constitutive equation for the total poroelastic stress tensor is formed as(9)σ=1−ϕ2σe+αpI,where α is Biot’s coefficient. Note that when the phase field variable ϕ⟶0, which represents the state for the medium, this constitutive equation Equation (9) converges to the conventional effective stress equation σ=σe+αpI.Finally, the linear momentum balance is expressed as(10)∇·σ+f−p∇ϕ=0,σ·n=t¯on∂NΩu,where the additional term p∇ϕ can be interpreted as fluid pressure forces that act on the fracture surfaces to increase/decrease the fracture volume. See the detailed derivation of the governing equation in appendix.


## 2.3. Fluid Flow Governing Equations

Assuming that fluid pressure is continuous from fracture to reservoir in the poroelastic medium, the fracture pressurepf and the pore pressure pr are related to p, the fluid pressure acting in the whole domain, as follows
(11)p=prinΩ\Γ,p=pfonΓ.Since the fracture width is much less than its length, i.e.,w<<L, incompressible single phase fluid flow within a fracture is governed by Reynolds’ lubrication equation
(12)ẇ−∇Γ·w312μf∇Γp+qL=qinjinΓ.Equation (12) is formulated in terms of the fracture surface tangential gradient operator ∇Γ=∇−nΓnΓ·∇ which eliminates the contribution of the normal component that is perpendicular to the fracture surfaces. qinj is the injection flow rate, qL is the leak-off rate, and μf is the fluid’s dynamic viscosity.The fracture opening width is estimated by(13)w=nΓ·εnΓhe,where he is the local mesh size and nΓ can be computed from Equation (5).Inspired by Biot’s theory and the fixed stress iterative approach [36, 37], the mass balance equation for incompressible single phase flow in the poroelastic medium is expressed as
(14)1M+α2Kdrṗ−∇·kμf∇p=qrs−α∇·u̇+α2KdrṗinΩ\Γ,where the additional term α2/Kdrṗ improves the numerical stability. qrs is the reservoir source/sink term. M, Kdr, and k are Biot’s modulus, the drained bulk modulus, and the permeability tensor, respectively. The relationships for these parameters are given by
(15)Kdr=E1−ν1+ν1−2ν,1M=α−ϕporKm+ϕporKf,where ϕpor is the matrix porosity. Km and Kf are the matrix and fluid bulk modulus, respectively.


## 2.4. Summary of Governing Equations

Finally, this highly nonlinear problem has been split into a set of linear governing equations and can be solved by an iterative coupling method. Findu and p for all time t with fully coupled governing equations as follows:
(16)∇·σ+f−p∇ϕ=0inΩ,ϕ−ℓ2Δϕ=0inΩ,ẇ−∇Γ·w312μf∇Γp+qL=qinjinΓ,1M+α2Kdrṗ−∇·kμf∇p=qrs−α∇·u̇+α2KdrṗinΩ∖Γ,σ·n=t¯on∂NΩu,∇ϕ·n=0on∂NΩϕ,kμf∇pk+1·nΓ=0on∂NΩp,kμf∇pk+1·nΓ=qLonΓ.This covers the stress equilibrium, the fluid mass balance, and the fracture regularization. The unknown leak-off termqL has been implicitly modeled as the fluid flux from the fracture into the reservoir. Therefore, there is no need for legacy fluid leak-off models, such as Carter’s formulation [38].Another implementation challenge in Equation (16) is that the lubrication equation is defined over the fracture domain Γ, and the mass balance equation in the poroelastic medium is defined over the porous medium domain Ω\Γ. Two separate mesh systems, for the fracture domain and the porous domain, respectively, can be used to solve these two domains separately. However, a cumbersome fracture tracking technique is required for this case. To circumvent this issue, the phase field calculus can be used (Equations (3) and (4)). This will be explained in the following section.


## 2.5. Weak Formulation

Assign three weighting functions:(17)φu∈φu∈H1Ω∣φu=0on∂DΩu,φϕ∈φϕ∈H1Ω∣φϕ=0on∂DΩϕ,φp∈φp∈H1Ω∣φp=0on∂DΩp.Multiplying Equation (16) by the weighting functions, integrating over Ω, and performing integration by parts. Those equations then become the following:
(18)∫Ω∇φu:σdV=∫Ωφu·fdV−∫Ωφu·p∇ϕdV+∫∂NΩuφu·t¯dS,(19)∫ΩφϕϕdV+ℓ2∫Ω∇φϕ·∇ϕdV=0,(20)∫Γ∇Γφp·w312μf∇ΓpdΓ=∫Γφpqinj−qL−ẇdΓ,(21)∫Ω∖Γφp1M+α2KdrṗdV+∫Ω∖Γ∇φpkμf∇pdV=∫Ω∖ΓφpqrsdV−α∫Ω∖Γφp∇·u̇dV+α2Kdr∫Ω∖ΓφpṗdV+∫ΓφpqLdΓ.The integration of the fracture flow is perfomed over the fractured domainΓ while the integration of the reservoir flow is performed over the unfractured domain Ω/Γ. The consequence is that the aforementioned domain treatments—Equations (3) and (4)—can be used for fluid flow integrals. In addition, the unknown leak-off term qL is immediately eliminated by adding Equations (20) and (21) together. This yields
(22)∫Ωφp1M+1−ϕ2α2KdrṗdV+∫Ω∇φpkμf∇pdV+∫Ω∇ϕ∇Γφp·w312μf∇ΓpdV=∫Ω1−ϕ2φpqrsdV−α∫Ω1−ϕ2φp∇·u̇dV+1M∫Ω1−1−ϕ2φpṗdV+α2Kdr∫Ω1−ϕ2φpṗdV−∫Ω∇ϕφpẇdV+∫Ω∇ϕφpqinjdV.Equations (18), (19), and (22) constitute a coupled bilinear system. In our implementation, u, ϕ, and p are discretized with quadrilateral finite elements in space and a backward Euler scheme in time.


## 2.6. Numerical Implementation

### 2.6.1. Finite Element Discretization

The finite element approximations of the solid displacement field the fracture phase field and the fluid pressure field are interpolated in terms of the nodal displacementsu¯, nodal phase field ϕ¯, and nodal pressure p¯, respectively
(23)u≈uhx=Nuu¯,ε≈Buu¯,ϕ≈ϕhx=Nϕϕ¯,∇ϕ≈Bϕϕ¯,p≈phx=Npp¯,∇p≈Bpp¯.where Nu denotes the matrix of shape functions and Nϕ and Np are the row vector of shape functions. For 2D elements with n nodes, they are given by (assume the same shape functions Nn are adopted for u, ϕ, and p without loss of generality)
(24)Nu=N10N20⋯Nn00N10N2⋯0NnNϕ=Np=N1N2⋯Nn.,In a 2D settings, the strain-displacement matrixBu and the differential matrix Bϕ, Bp have the following components:
(25)Bu=N1,x0N2,x0⋯Nn,x00N1,y0N2,y⋯0Nn,yN1,yN1,xN2,yN2,x⋯Nn,yNn,xBϕ=Bp=N1,xN1,x⋯N1,xN1,yN1,y⋯N1,y.,Substituting Equation (23) to the weak form Equation (18) and approximating each cells contribution by Gauss quadrature rule, the linear system is obtained:
(26)Kuuu¯=Fu−Kupp¯.The tangent stiffness matrixKuu, the coupling matrix Kup, and the external force vector Fu are given by
(27)Kuu=∫Ω1−Nϕϕ¯2BuTDBudV,Kup=∫Ω1−Nϕϕ¯2αBuTmNpdV,Fu=∫∂ΩuNuTt¯dS+∫ΩρNuTfdV,where D is the elasticity tensor, and m=110T for 2D.Similarly, the matrix-vector form of the phase field Equation (19) is obtained as follows:
(28)Mϕ+Kϕϕ¯=0Mϕ=∫ΩNϕTNϕdV,Kϕ=∫Ωℓ2BϕTBϕdV.Finally, the matrix-vector form of the fluid flow Equation (22) is given by
(29)Mprp¯̇+Kpp+Kpfp¯=Q−Kpuu¯+Hprp,¯̇where the
(30)Mpr=∫Ω1M+1−Nϕϕ¯2α2KdrNpTNpdV,Kpp=∫ΩkμfBpTBpdV,Kpf=∫ΩBϕϕ¯w312μfBpTBpdV,Q=∫ΩNpT1−Nϕϕ¯2qrsdV−∫ΩNpTBϕϕ¯ẇdV+∫ΩNpTBϕϕ¯qinjdV,Kpu=∫Ω1−Nϕϕ¯2αNpTmTBudV,Hpr=∫Ω1M1−1−Nϕϕ¯2NpTNpdV+∫Ωα2Kdr1−Nϕϕ¯2NpTNpdV.


### 2.6.2. Iterative Coupling Algorithm

Finally, an extended fixed stress iterative coupling method for fracture networks has been introduced to solve three matrix-vector systems Equations (26), (28), and (29). The central idea of this algorithm is to hold the other two variables constant to yield three extremely simple linear subsystems:
(31)Kuuu¯k+1=Fu−Kupp¯k+1,(32)Mprp¯k+1+Kpp+Kpfp¯k+1=Q−Kpuu¯k+Hprp¯k,(33)Mϕ+Kϕϕ¯=0.The numerical implementation is presented in Algorithm1. k represents the iteration levels. uk+1, pk+1 represent unknown solutions in current iteration, while uk, pk are known solutions from previous iteration.Algorithm 1: Extended fixed stress iterative coupling method.
1 Set phase field Dirichlet boundary conditionsϕ=1 for pre-existing fracture networks;2 Set initial conditions;3 Solve the phase field formulation Equation (33) for ϕ;4 At the timetn;5repeat6  Update the fracture outward normalnΓk using Equation (5);7  Update the fracture widthwk using Equation (13);8  Solve the pressure formulation Equation (32) for pk+1 using ϕ, pk, uk, nΓk, and wk;9  Solve the displacement formulation Equation (31) for uk+1 using ϕ and pk+1;10untilRp<ϵtol;11 increment the timetn⟶tn+1.According to our numerical tests, the convergence of pressurep is harder to reach than displacement u due to strong nonlinear of fluid flow. Therefore, the global L2 error of pressure Rp=pk+1−pk is sufficient to examine both the convergence of the total coupling system. In general, the fixed stress iterative coupling algorithm is very robust (unconditional stable) and efficient (the number of iterations is less than 10).


## 2.6.1. Finite Element Discretization

The finite element approximations of the solid displacement field the fracture phase field and the fluid pressure field are interpolated in terms of the nodal displacementsu¯, nodal phase field ϕ¯, and nodal pressure p¯, respectively
(23)u≈uhx=Nuu¯,ε≈Buu¯,ϕ≈ϕhx=Nϕϕ¯,∇ϕ≈Bϕϕ¯,p≈phx=Npp¯,∇p≈Bpp¯.where Nu denotes the matrix of shape functions and Nϕ and Np are the row vector of shape functions. For 2D elements with n nodes, they are given by (assume the same shape functions Nn are adopted for u, ϕ, and p without loss of generality)
(24)Nu=N10N20⋯Nn00N10N2⋯0NnNϕ=Np=N1N2⋯Nn.,In a 2D settings, the strain-displacement matrixBu and the differential matrix Bϕ, Bp have the following components:
(25)Bu=N1,x0N2,x0⋯Nn,x00N1,y0N2,y⋯0Nn,yN1,yN1,xN2,yN2,x⋯Nn,yNn,xBϕ=Bp=N1,xN1,x⋯N1,xN1,yN1,y⋯N1,y.,Substituting Equation (23) to the weak form Equation (18) and approximating each cells contribution by Gauss quadrature rule, the linear system is obtained:
(26)Kuuu¯=Fu−Kupp¯.The tangent stiffness matrixKuu, the coupling matrix Kup, and the external force vector Fu are given by
(27)Kuu=∫Ω1−Nϕϕ¯2BuTDBudV,Kup=∫Ω1−Nϕϕ¯2αBuTmNpdV,Fu=∫∂ΩuNuTt¯dS+∫ΩρNuTfdV,where D is the elasticity tensor, and m=110T for 2D.Similarly, the matrix-vector form of the phase field Equation (19) is obtained as follows:
(28)Mϕ+Kϕϕ¯=0Mϕ=∫ΩNϕTNϕdV,Kϕ=∫Ωℓ2BϕTBϕdV.Finally, the matrix-vector form of the fluid flow Equation (22) is given by
(29)Mprp¯̇+Kpp+Kpfp¯=Q−Kpuu¯+Hprp,¯̇where the
(30)Mpr=∫Ω1M+1−Nϕϕ¯2α2KdrNpTNpdV,Kpp=∫ΩkμfBpTBpdV,Kpf=∫ΩBϕϕ¯w312μfBpTBpdV,Q=∫ΩNpT1−Nϕϕ¯2qrsdV−∫ΩNpTBϕϕ¯ẇdV+∫ΩNpTBϕϕ¯qinjdV,Kpu=∫Ω1−Nϕϕ¯2αNpTmTBudV,Hpr=∫Ω1M1−1−Nϕϕ¯2NpTNpdV+∫Ωα2Kdr1−Nϕϕ¯2NpTNpdV.


## 2.6.2. Iterative Coupling Algorithm

Finally, an extended fixed stress iterative coupling method for fracture networks has been introduced to solve three matrix-vector systems Equations (26), (28), and (29). The central idea of this algorithm is to hold the other two variables constant to yield three extremely simple linear subsystems:
(31)Kuuu¯k+1=Fu−Kupp¯k+1,(32)Mprp¯k+1+Kpp+Kpfp¯k+1=Q−Kpuu¯k+Hprp¯k,(33)Mϕ+Kϕϕ¯=0.The numerical implementation is presented in Algorithm1. k represents the iteration levels. uk+1, pk+1 represent unknown solutions in current iteration, while uk, pk are known solutions from previous iteration.Algorithm 1: Extended fixed stress iterative coupling method.
1 Set phase field Dirichlet boundary conditionsϕ=1 for pre-existing fracture networks;2 Set initial conditions;3 Solve the phase field formulation Equation (33) for ϕ;4 At the timetn;5repeat6  Update the fracture outward normalnΓk using Equation (5);7  Update the fracture widthwk using Equation (13);8  Solve the pressure formulation Equation (32) for pk+1 using ϕ, pk, uk, nΓk, and wk;9  Solve the displacement formulation Equation (31) for uk+1 using ϕ and pk+1;10untilRp<ϵtol;11 increment the timetn⟶tn+1.According to our numerical tests, the convergence of pressurep is harder to reach than displacement u due to strong nonlinear of fluid flow. Therefore, the global L2 error of pressure Rp=pk+1−pk is sufficient to examine both the convergence of the total coupling system. In general, the fixed stress iterative coupling algorithm is very robust (unconditional stable) and efficient (the number of iterations is less than 10).


## 3. Numerical Tests

In this section, two verification benchmarks (i.e., Terzaghi’s and Sneddon’s problems) are considered to ensure that the model developments are constrained by the appropriate underlying physics and mathematics. Then, in order to assess the performance of our model, a complex discrete fracture network flow case is presented. The numerical formulation was implemented using C++.


### 3.1. Terzaghi’s Consolidation Problem

The simplest validation benchmark for poroelasticity is Terzaghi’s one-dimensional consolidation problem [39], as shown in Figure 4. All of the boundaries are impermeable except the top one. The sample is loaded by a constant vertical compressive force t¯ at its top surface. Material parameters were chosen as follows: Young’s modulus E=20GPa, Poisson’s ratio ν=0.2, matrix bulk modulus Km=5GPa, fluid bulk modulus Kf=1GPa, Biot’s coefficient α=0.6, porosity ϕpor=0.3, permeability k=2.8e−12m2, and fluid viscosity μf=100mPa·s.Figure 4
Schematic of Terzaghi’s consolidation problem domain and boundary conditions.The pore pressure distributions at time 20s, 200 s, and 300 s are shown in Figure5. The decrease in pore pressure is accompanied by a gradual increase of displacement. In other words, during the consolidation process, the sample developed from an undrained to a fully drained state.Figure 5
Pore pressure 2D surfaces at timet=20s, 200 s, and 300 s for Terzaghi’s problem.Our model was compared with the analytical solutions which were summarized by [40]. The results are plotted in Figure 6, where analytical solutions are lines while the numerical solutions are shown as circles. The pore pressure at different times, as determined by the numerical solutions, are in excellent agreement with the analytical solutions.Figure 6
Dimensionless pore pressure-height curve for for Terzaghi’s problem.


### 3.2. Sneddon’s Pressurized Fracture Problem

Sneddon’s pressurized fracture problem is a classical benchmark for a single pressurized, plane strain and static fracture in an infinite elastic medium. The analytical solution [41] of this problem is based on the assumption of a uniform pressure, acting on constant length fracture surfaces. Our model was validated by comparing the numerical results with the analytical solution.Since the analytical solution is for an infinite domain, we considered a large domain size of 40 m relative to the preexisting fracture length of 4 m, to eliminate boundary effects. Young’s modulus was chosen asE=20GPa and Poisson’s ratio was specified as v=0.2. Biot’s coefficient was set to α=0, so that the poroelastic effect was neglected. The preexisting fracture was represented by the phase field with a fixed length scale ℓ=0.2m, as illustrated in Figure 7. The constant fluid pressure was chosen as p=2MPa.Figure 7
Sneddon’s pressurized fracture. (a) Internally pressurized plane-strain fracture is centered in the 2D domain. (b) Phase field contour surface in a 2D mesh (40m×40m). ϕ=1 represents the preexisting fracture and ϕ=0 represents the elastic medium. Zero displacement Dirichlet boundary conditions are applied for all the domain boundaries. However, the numerical simulation is insensitive to other types of displacement boundary conditions since the computational domain is large relative to the fracture length.
(a)(b)For comparison, two orthogonal horizontal stress componentsσx and σy are plotted for both the analytical solution and our model at the points along the y axis (see Figure 8(a)). Our results matched the analytical solution well except in the diffusive fracture zone (i.e., near the fracture surface). The stresses of the phase field model degraded to zero near the fracture surface because both stress and pressure are smeared. This phenomenon is often called strain softening in elastic plastic fracture mechanics (which is ignored in linear elastic fracture mechanics). From a mathematical point of view, when the phase field variable ϕ⟶1, the total stress σ=1−ϕ2σe+αpI⟶0. This continuum feature brings a bridge between the strong discontinuous fracture and the continuum rock matrix and avoids dealing with complicated boundary conditions along fracture surfaces.Figure 8
Stress-distance curves for Sneddon’s problem. (a) Comparison between numerical and analytical results. (b) Effect of the length scale on the stress-distance curves. The stress component along they axis σy is plotted.
(a)(b)The diffusive fracture zone is commonly interpreted as a fracture process zone in which micro-cracks evolved. We observed that the length scale controls the width of this fracture process zone. The effect of different length scalesℓ is shown in Figure 8(b). In other words, a larger length scale value leads to a wider diffusive fracture zone. Consequently, the peak stress value increased as the length scale decreased. When the length scale approaches zero, the stresses increase, i.e., the sharp fracture will be recovered. By means of this analysis, it is implied that the length scale ℓ might be considered a material parameter since it influences the stress peak value. This issue has also been discussed in [42].


### 3.3. Discrete Fracture Network Flow

A final evaluation shows a two-dimension simulation of a simplified discrete fracture network problem. The domain size is10m×10m and includes a complex discrete fracture network consisting of several connected or disconnected fractures (see Figure 9). The injection point is located at the center of the domain and the injection rate was chosen as 2.5e−4m2/s. Although depletion tests can be dealt with straightforwardly in an identical manner, we defer interested readers to the study of geomechanical response during depletion [43, 44]. Material parameters were chosen as follows: Young’s modulus E=17GPa, Poisson’s ratio ν=0.2, matrix bulk modulus Km=36GPa, fluid bulk modulus Kf=2.2GPa, Biot’s coefficient α=0.8, porosity ϕpor=0.2, permeability k=9e−15m2, the initial pore pressure p=1MPa, and fluid viscosity μf=0.89mPa·s. For simplicity’s sake, the homogeneous Dirichlet/Neumann boundary conditions (i.e. u=0,∇p=0) are applied to all domain boundaries. In order to assess the impact of solid-fluid coupling, all of the preexisting fractures are originally inactive, i.e., the initial widths were set to zero.Figure 9
A 2D domain with a discrete fracture network. The whole domain is discretized into a structured quadratic mesh. The mesh around fractures is locally refined. The orange spot at the center0,0 is the injection point. Discrete fractures are represented by line segments. The coordinates of two endpoints xa,ya,xb,yb for each fracture line segment are predefined: ((-2, 0), (2, 0)), ((-2.8, 0.625), (-1.25, 0.625)), ((-3.125, 0.3125), (-4.0625, -0.625)), ((2.8125, 0.3125), (4.375, -1.25)), ((1.5625,-0.0390625),(1.5625,-0.625)), ((-1.25, -0.625),(0, -1.25)), ((-1.71875, 0.585938),(-1.71875, 0.0390625)), ((1.60156, -0.625),(3.71094, -0.625)), and ((0, 1.25),(0.625, 0.625)).The distributions of the pressure field at different times are shown in Figure10. As can be seen, the complex discrete fracture network leads to a diffuse but heterogeneous distribution of the pressure field. Another notable result is that the contribution of the disconnected fractures was very small compared to the interconnected fractures. In other words, the activated and interconnected fractures became the main fluid flow paths. This is caused by the “stress shadow” effect; that is, the interconnected open fractures exert additional stresses on the surrounding matrix and adjacent disconnected fractures. This stress shadow significantly impacts the opening of the disconnected fractures. This interaction (shadowing) can not be observed in a pure fluid diffusion process where this kind of coupling is neglected. Moreover, some researchers have also shown that the isolated fractures can not only be constricted but also be activated when considering shear dilation effects [45–47].Figure 10
Pressure contours at different times. The time is shown at the top of each subplot. (a–d) Subfigures indicate the progress that the interconnected fractures become the diffusion-dominated area due to the coupling effect.
(a)(b)(c)(d)Figures11 and 12 demonstrate that the phase field model can easily address a discontinuous displacement field across fracture surfaces. In these figures, the displacement contour is asymmetrical due to the impact of the fracture networks. By virtue of the gradually increasing pressure, the poroelastic displacement tends to expand the fracture surfaces.Figure 11
ux displacement contour at different times. The time is presented at the top of each subplot.
(a)(b)(c)(d)Figure 12
uy displacement contour at different times. The time is presented at the top of each subplot.
(a)(b)(c)(d)


## 3.1. Terzaghi’s Consolidation Problem

The simplest validation benchmark for poroelasticity is Terzaghi’s one-dimensional consolidation problem [39], as shown in Figure 4. All of the boundaries are impermeable except the top one. The sample is loaded by a constant vertical compressive force t¯ at its top surface. Material parameters were chosen as follows: Young’s modulus E=20GPa, Poisson’s ratio ν=0.2, matrix bulk modulus Km=5GPa, fluid bulk modulus Kf=1GPa, Biot’s coefficient α=0.6, porosity ϕpor=0.3, permeability k=2.8e−12m2, and fluid viscosity μf=100mPa·s.Figure 4
Schematic of Terzaghi’s consolidation problem domain and boundary conditions.The pore pressure distributions at time 20s, 200 s, and 300 s are shown in Figure5. The decrease in pore pressure is accompanied by a gradual increase of displacement. In other words, during the consolidation process, the sample developed from an undrained to a fully drained state.Figure 5
Pore pressure 2D surfaces at timet=20s, 200 s, and 300 s for Terzaghi’s problem.Our model was compared with the analytical solutions which were summarized by [40]. The results are plotted in Figure 6, where analytical solutions are lines while the numerical solutions are shown as circles. The pore pressure at different times, as determined by the numerical solutions, are in excellent agreement with the analytical solutions.Figure 6
Dimensionless pore pressure-height curve for for Terzaghi’s problem.


## 3.2. Sneddon’s Pressurized Fracture Problem

Sneddon’s pressurized fracture problem is a classical benchmark for a single pressurized, plane strain and static fracture in an infinite elastic medium. The analytical solution [41] of this problem is based on the assumption of a uniform pressure, acting on constant length fracture surfaces. Our model was validated by comparing the numerical results with the analytical solution.Since the analytical solution is for an infinite domain, we considered a large domain size of 40 m relative to the preexisting fracture length of 4 m, to eliminate boundary effects. Young’s modulus was chosen asE=20GPa and Poisson’s ratio was specified as v=0.2. Biot’s coefficient was set to α=0, so that the poroelastic effect was neglected. The preexisting fracture was represented by the phase field with a fixed length scale ℓ=0.2m, as illustrated in Figure 7. The constant fluid pressure was chosen as p=2MPa.Figure 7
Sneddon’s pressurized fracture. (a) Internally pressurized plane-strain fracture is centered in the 2D domain. (b) Phase field contour surface in a 2D mesh (40m×40m). ϕ=1 represents the preexisting fracture and ϕ=0 represents the elastic medium. Zero displacement Dirichlet boundary conditions are applied for all the domain boundaries. However, the numerical simulation is insensitive to other types of displacement boundary conditions since the computational domain is large relative to the fracture length.
(a)(b)For comparison, two orthogonal horizontal stress componentsσx and σy are plotted for both the analytical solution and our model at the points along the y axis (see Figure 8(a)). Our results matched the analytical solution well except in the diffusive fracture zone (i.e., near the fracture surface). The stresses of the phase field model degraded to zero near the fracture surface because both stress and pressure are smeared. This phenomenon is often called strain softening in elastic plastic fracture mechanics (which is ignored in linear elastic fracture mechanics). From a mathematical point of view, when the phase field variable ϕ⟶1, the total stress σ=1−ϕ2σe+αpI⟶0. This continuum feature brings a bridge between the strong discontinuous fracture and the continuum rock matrix and avoids dealing with complicated boundary conditions along fracture surfaces.Figure 8
Stress-distance curves for Sneddon’s problem. (a) Comparison between numerical and analytical results. (b) Effect of the length scale on the stress-distance curves. The stress component along they axis σy is plotted.
(a)(b)The diffusive fracture zone is commonly interpreted as a fracture process zone in which micro-cracks evolved. We observed that the length scale controls the width of this fracture process zone. The effect of different length scalesℓ is shown in Figure 8(b). In other words, a larger length scale value leads to a wider diffusive fracture zone. Consequently, the peak stress value increased as the length scale decreased. When the length scale approaches zero, the stresses increase, i.e., the sharp fracture will be recovered. By means of this analysis, it is implied that the length scale ℓ might be considered a material parameter since it influences the stress peak value. This issue has also been discussed in [42].


## 3.3. Discrete Fracture Network Flow

A final evaluation shows a two-dimension simulation of a simplified discrete fracture network problem. The domain size is10m×10m and includes a complex discrete fracture network consisting of several connected or disconnected fractures (see Figure 9). The injection point is located at the center of the domain and the injection rate was chosen as 2.5e−4m2/s. Although depletion tests can be dealt with straightforwardly in an identical manner, we defer interested readers to the study of geomechanical response during depletion [43, 44]. Material parameters were chosen as follows: Young’s modulus E=17GPa, Poisson’s ratio ν=0.2, matrix bulk modulus Km=36GPa, fluid bulk modulus Kf=2.2GPa, Biot’s coefficient α=0.8, porosity ϕpor=0.2, permeability k=9e−15m2, the initial pore pressure p=1MPa, and fluid viscosity μf=0.89mPa·s. For simplicity’s sake, the homogeneous Dirichlet/Neumann boundary conditions (i.e. u=0,∇p=0) are applied to all domain boundaries. In order to assess the impact of solid-fluid coupling, all of the preexisting fractures are originally inactive, i.e., the initial widths were set to zero.Figure 9
A 2D domain with a discrete fracture network. The whole domain is discretized into a structured quadratic mesh. The mesh around fractures is locally refined. The orange spot at the center0,0 is the injection point. Discrete fractures are represented by line segments. The coordinates of two endpoints xa,ya,xb,yb for each fracture line segment are predefined: ((-2, 0), (2, 0)), ((-2.8, 0.625), (-1.25, 0.625)), ((-3.125, 0.3125), (-4.0625, -0.625)), ((2.8125, 0.3125), (4.375, -1.25)), ((1.5625,-0.0390625),(1.5625,-0.625)), ((-1.25, -0.625),(0, -1.25)), ((-1.71875, 0.585938),(-1.71875, 0.0390625)), ((1.60156, -0.625),(3.71094, -0.625)), and ((0, 1.25),(0.625, 0.625)).The distributions of the pressure field at different times are shown in Figure10. As can be seen, the complex discrete fracture network leads to a diffuse but heterogeneous distribution of the pressure field. Another notable result is that the contribution of the disconnected fractures was very small compared to the interconnected fractures. In other words, the activated and interconnected fractures became the main fluid flow paths. This is caused by the “stress shadow” effect; that is, the interconnected open fractures exert additional stresses on the surrounding matrix and adjacent disconnected fractures. This stress shadow significantly impacts the opening of the disconnected fractures. This interaction (shadowing) can not be observed in a pure fluid diffusion process where this kind of coupling is neglected. Moreover, some researchers have also shown that the isolated fractures can not only be constricted but also be activated when considering shear dilation effects [45–47].Figure 10
Pressure contours at different times. The time is shown at the top of each subplot. (a–d) Subfigures indicate the progress that the interconnected fractures become the diffusion-dominated area due to the coupling effect.
(a)(b)(c)(d)Figures11 and 12 demonstrate that the phase field model can easily address a discontinuous displacement field across fracture surfaces. In these figures, the displacement contour is asymmetrical due to the impact of the fracture networks. By virtue of the gradually increasing pressure, the poroelastic displacement tends to expand the fracture surfaces.Figure 11
ux displacement contour at different times. The time is presented at the top of each subplot.
(a)(b)(c)(d)Figure 12
uy displacement contour at different times. The time is presented at the top of each subplot.
(a)(b)(c)(d)


## 4. Conclusion

In this paper, we presented a fully coupled model in a poroelastic, fractured reservoir linked to phase field modeling of discrete fracture networks.(i)
This model captures the complex coupled interaction of Darcy-Biot-type fluid flow in poroelastic media, Reynolds lubrication governing flow inside fracture networks, mass exchange between fractures and matrix, and the subsequent geomechanical response of the medium(ii)
This approach is capable of modeling arbitrary number of fractures. The fractures can be added trivially through setting the initial phase field valueϕ=1 for the corresponding points(iii)
The proposed iterative coupling algorithm provides a simple and efficient implementation to subdivide the nonlinear coupled problem without any complex numerical techniques or consequences(iv)
The numerical model has been validated against classical benchmarks and the effect of the length scale on the diffusive fracture zone has been analyzed(v)
The numerical test of discrete fractures demonstrates that the contribution of disconnected fractures was dramatically smaller than that of the interconnected fractures due to the poor fracture opening. This suggests that the well stimulation treatments, e.g., hydraulic fracturing and acidizing, should endeavor to connect nature fracturesAlthough some good results were obtained, there are many issues worthy of further investigation. The computational cost of the phase field model is notoriously expensive. This highly computational cost caused by fine mesh around fractures can be alleviated using an adaptive mesh refinement technique. Furthermore, since the coupling between the fluid flow inside fractures and the rock deformation strongly dependents on the cubic law (i.e., the volumetric flow rate is proportional to the cube of the fracture width), the fracture width computation needs to be improved through the development of a more robust and accurate algorithm.



---

*Source: 6235441-2021-12-20.xml*