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.

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 [35]. 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 -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 [810]. 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 [1215]. 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 [1921]. 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 [3032]. 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 with spatial dimension , composed by a porous medium domain and an internal pressurized fracture domain (see Figure 1). Dirichlet and Neumann boundary conditions are imposed on and , respectively.

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 is called the fracture phase field variable, and it can be used to describe the regularized fracture topology.

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 .

The phase field governing equations are derived using the variational principle [35].

Next, we can introduce to convert integrals over the fracture domain to the whole domain . where is an arbitrary function.

Similarly, is used to convert integrals over the porous medium domain to the whole domain .

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 Figure 3(a)).

Figure 3(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

2.2. Geomechanical Governing Equations

We assume that the strains are small, giving the infinitesimal strain tensor as

The constitutive equation, based on standard linear elasticity for isotropic solids, is where the Lamé coefficients and are given by and they are expressed in terms of Young’s modulus 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 where is Biot’s coefficient. Note that when the phase field variable , which represents the state for the medium, this constitutive equation Equation (9) converges to the conventional effective stress equation .

Finally, the linear momentum balance is expressed as where the additional term 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 pressure and the pore pressure are related to , the fluid pressure acting in the whole domain, as follows

Since the fracture width is much less than its length, i.e., , incompressible single phase fluid flow within a fracture is governed by Reynolds’ lubrication equation

Equation (12) is formulated in terms of the fracture surface tangential gradient operator which eliminates the contribution of the normal component that is perpendicular to the fracture surfaces. is the injection flow rate, is the leak-off rate, and is the fluid’s dynamic viscosity.

The fracture opening width is estimated by where is the local mesh size and 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 where the additional term improves the numerical stability. is the reservoir source/sink term. , , and are Biot’s modulus, the drained bulk modulus, and the permeability tensor, respectively. The relationships for these parameters are given by where is the matrix porosity. and 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. Find and for all time with fully coupled governing equations as follows:

This covers the stress equilibrium, the fluid mass balance, and the fracture regularization. The unknown leak-off term 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:

Multiplying Equation (16) by the weighting functions, integrating over , and performing integration by parts. Those equations then become the following:

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 is immediately eliminated by adding Equations (20) and (21) together. This yields

Equations (18), (19), and (22) constitute a coupled bilinear system. In our implementation, , , and 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 displacements , nodal phase field , and nodal pressure , respectively where denotes the matrix of shape functions and and are the row vector of shape functions. For 2D elements with nodes, they are given by (assume the same shape functions are adopted for , , and without loss of generality)

In a 2D settings, the strain-displacement matrix and the differential matrix , have the following components:

Substituting Equation (23) to the weak form Equation (18) and approximating each cells contribution by Gauss quadrature rule, the linear system is obtained:

The tangent stiffness matrix , the coupling matrix , and the external force vector are given by where is the elasticity tensor, and for 2D.

Similarly, the matrix-vector form of the phase field Equation (19) is obtained as follows:

Finally, the matrix-vector form of the fluid flow Equation (22) is given by where the

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:

The numerical implementation is presented in Algorithm 1. represents the iteration levels. , represent unknown solutions in current iteration, while , are known solutions from previous iteration.

1 Set phase field Dirichlet boundary conditions for pre-existing fracture networks;
2 Set initial conditions;
3 Solve the phase field formulation Equation (33) for ;
4 At the time ;
5 repeat
6  Update the fracture outward normal using Equation (5);
7  Update the fracture width using Equation (13);
8  Solve the pressure formulation Equation (32) for using , , , , and ;
9  Solve the displacement formulation Equation (31) for using and ;
10 until;
11 increment the time .

According to our numerical tests, the convergence of pressure is harder to reach than displacement due to strong nonlinear of fluid flow. Therefore, the global error of pressure 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 at its top surface. Material parameters were chosen as follows: Young’s modulus , Poisson’s ratio , matrix bulk modulus , fluid bulk modulus , Biot’s coefficient , porosity , permeability , and fluid viscosity .

The pore pressure distributions at time 20s, 200 s, and 300 s are shown in Figure 5. 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.

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.

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 as and Poisson’s ratio was specified as . Biot’s coefficient was set to , so that the poroelastic effect was neglected. The preexisting fracture was represented by the phase field with a fixed length scale , as illustrated in Figure 7. The constant fluid pressure was chosen as .

For comparison, two orthogonal horizontal stress components and are plotted for both the analytical solution and our model at the points along the 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 , the total stress . 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.

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 is 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 . 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 , Poisson’s ratio , matrix bulk modulus , fluid bulk modulus , Biot’s coefficient , porosity , permeability , the initial pore pressure , and fluid viscosity . For simplicity’s sake, the homogeneous Dirichlet/Neumann boundary conditions (i.e. ) 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.

The distributions of the pressure field at different times are shown in Figure 10. 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 [4547].

Figures 11 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.

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 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 fractures

Although 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.


Stress Equilibrium Governing Equation

The governing equation is derived from the total energy functional. To account for the fracturing fluid pressure in the energy functional, there is an additional term which can be interpreted as the external force work in increasing fracture volume. where is the jump displacement between fracture surfaces. Substituting the fracture normal Equation (5) into the energy functional Equation (A.1) and regularizing the fracture domain using the phase field gradient :

Based on Biot’s theory, the pore pressure energy functional considering the phase field is given by

The total potential energy functional is obtained by combining the pressure energy functional with the elastic energy functional and external energy functional:

The functional derivative of the total potential energy (A.4) with respect to results

With the Neumann boundary conditions on ,

Since is an arbitrary function, applying the fundamental lemma of calculus of variations, we obtain the stress equilibrium equations Equation (10).

Data Availability

No data were used in producing this manuscript.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Jinzhou Zhao contributed to the conceptualization, supervision, and review. Qing Yin contributed to the methodology, software acquisition, writing, and validation. John McLennan contributed to the conceptualization and writing. Yongming Li contributed to the methodology and acquistion of resources. Yu Peng and Xiyu Chen contributed to the acquisition of software and visualization. Cheng Chang, Weiyang Xie, and Zhongyi Zhu contributed to the review and resources and funding acquisition.


This study was supported by the National Natural Science Foundation of China (U19A2043), China Postdoctoral Science Foundation (2018M640935, 2020M683360, and 2020T130548), Science and Technology Cooperation Project of the CNPC-SWPU Innovation Alliance, and Scientific Research Starting Project of Southwest Petroleum University (SWPU) (No. 2018QHZ004).