Numerical analysis of geomechanical behavior of fractures and faults in a deformable porous medium

In this study, geomechanical behavior of fractures and faults in rocks as a saturated deformable discrete fracture-porous medium has been evaluated using coupled fluid flow-geomechanics numerical modeling. The purpose of this paper is to observe and evaluate the effects of fractures and faults on the pore pressures during fluid flow through reservoirs. This issue involves solving the equations that have been derived from the Biot consolidation theory such as fluid mass balance equation, Darcy law and momentum balance equations. Govern coupled equations were solved using the standard Galerkin finite element method for continuous porous medium. Elements called "zero-thickness elements” were also used to discretize the fault as a discontinuous part of the porous medium. In compared with the previous and similar methods, the method introduced in this paper, made modifications in either the choice of the element and the method of solving the governing equations. The main advantage of this paper is providing clear precise formulations of the double node zero-thickness element in hydro-mechanical modeling of fractures and faults. Verification of the proposed process and models presented in this paper were done by providing three index problems which their analytical and numerical solutions are available. The results of our model provide a good agreement to these reference solutions which indicates the accuracy of the method presented in this paper.


Introduction
Based on the reports, fossil fuels are one of the dominant energy sources until 2035 (Azin et al. 2022a). One of the important challenges regarding to fluid flow simulation in the porous medium of oil reservoirs is the existence of fractures and faults that can exhibit different flow behavior in the lifetime of a reservoir. These behaviors are caused by two domains of mechanics and fluid flow, they are generally known as hydro-mechanical behaviors (Azin et al. 2022b). Therefore, in order to numerical modeling of such behavior of faults, it is necessary to use coupled methods that connect fluid flow to geomechanics. In this regard, the existence of a coupled model will be necessary to study the hydromechanical behavior of the reservoir rock and the evaluation of their faults during the production or injection of fluid into the reservoir. It should be noted that the governing equations of this problem have been developed based on the theory of poroelasticity of the Biot (1941). In previous studies, researchers have used several methods to model faults using finite element-based methods such as extended finite element method (Prevost and Sukumar 2016), control volume-finite element method (Zhang et al. 2016), finite volume-finite difference (Rutqvist et al. 2013), etc. For more information on the application of different methods in this area, the reader is referred to the study of Jing (2003). Deb (2010) also introduced the application of the finite element method in geomechanics of rock fractures and joints in his book.
One of the most effective methods to model the geomechanical behavior of fractures and faults is the use of interface elements which could be zero-thickness elements, solid elements or combination of solid elements and ubiquitousjoints oriented as weak planes. Cappa and Rutqvist (2011) examined and compared these types of elements in their coupled modeling using software. In this study, the zerothickness element is used to model the fractures and faults as a discontinuous part of a continuous medium which meshed with the finite element method. This element was introduced by Goodman et al. (1968). This zero-thickness element was simple and considered only the absolute displacement. Ghaboussi et al. (1973) developed this type of element to include relative displacement for slipping or debonding surfaces. After him, many researchers used and developed this element (Cerfontaine et al. 2015;Garolera Vinent et al. 2015;Carol 2004, 2010). This element has the ability to describe the slip surfaces as well as the opening of the discontinuities. In this paper, Ghaboussi's element type formulation with some modifications in number of element nodes and application of different elements in the calculations of displacement and pressure is used to model the hydromechanical behavior of fractures and faults inside deformable porous medium. Among the articles related to this article is the study of researchers such as Carol 2010, 2008a, b), Watanabe et al. (2012), Alonso et al. (2013) and Guiducci et al. (2002Guiducci et al. ( , 2003. Guiducci et al. (2002) proposed a method for the coupled hydro-mechanical cracks behavior. Their model is validated by comparison between experimental data and numerical predictions. They stated that the coupled model can predict the fracture behavior with higher accuracy (Guiducci et al. 2002). Segura and Carol (2008) used the finite element method to develop a coupled hydro-mechanical formulation. This model can be used for fractured geo-materials. Their formulation is applicable for both developing and preexisting discontinuities. They compared their results with other numerical solution in the literature. Moreover, they conducted a sensitivity analysis on the effects of mechanical parameters of the discontinuity (Segura and Carol 2008a, b). Watanabe et al. (2012) developed interface element for representing pre-existing cracks. They used coupled hydromechanical problems with the use of finite element method in discrete porous-fractures. They verified their model by examining various numerical examples (Watanabe et al. 2012). Alonso et al. (2013) also presented a formulation of thermo-hydro-mechanical problems in different scales. Energy and mass balance equations were derived in equilibrium state for air and water phases (Alonso et al. 2013). Benedetto et al. (2018) proposed a new technique for analyzing the fracture processes. They combined interface elements and virtual element method. They also used the classical zero-thickness interface element for simulating the stress-fracture processes. Various numerical results were employed by them to show the capability of their proposed model (Benedetto et al. 2018). Urpi et al. (2020) proposed a model for presenting a shear failure in mechanical discontinuities such as fractures, faults and bedding planes. Their model was validated using an analytical solution and field data (Urpi et al. 2020). Wang (2020) reported detailed studies about the analysis of damage and fractures in his book. Six different methods with the finite element basis are introduced in this book. More information are available in Wang (2020).
In this study, hydro-mechanical behavior of rock discontinuities such as fractures and faults is investigated using coupled fluid flow-geomechanical model. Moreover, the equations are solved for a continuous part of a faulted porous medium using the standard finite element Galerkin method. Since the conventional method of finite element is inefficient in discontinuous issues, it is necessary to use other methods to model the discontinuities in mediums such as fractures and faults. In the proposed model, conventional isoparametric finite elements and zero-thickness finite elements are used to continuous and discontinuous parts of model environment, respectively. In order to achieve a more accurate solutions, each of these elements is selected according to the displacement and pore pressure, which are the main unknown parameters in governing equations. The main priority and advantage of this paper can be considered as the presentation of clear and precise hydro-mechanical formulation of zero-thickness interface elements in fractures and faults modeling. In this regard, three problems were selected for verification. The numerical examples presented in this paper are such that the capacity and ability of the stated method are well represented in fault modeling.

Governing equations
The governing equations of coupled fluid flow and geomechanical modeling in isothermal medium are mass equilibrium equation and momentum equilibrium, which are written for liquid and solid phases. It should be noted that these equations have been developed based on Biot theory of poroelasticity (Biot 1941). In order to solve these fundamental equations, construction equations and assumptions are also used. In the following, the equations are presented.

Momentum balance equation
Considering small displacements, the linear equilibrium equation for a porous medium, as well as kinematic equations and construction relations for a solid phase, is written as follows (Watanabe et al. 2012): where ′ is an effective stress tensor, and represents Biot coefficient of the porous medium. m is a vector that causes fluid pressures to affect only on the principal stresses. u is solid displacement vector, and D e is elastic matrix. e denotes strain, and g is the gravitational acceleration of the earth. is a combination of fluid density and solid density calculated as follows: where , f and s , respectively, represent solid phase(rock) porosity, density of the fluid (water) and density of the solid, respectively. With the assumption that the strain tensor is symmetric, kinematics is defined by Eq. (1). In nonlinear cases, Eq. (3) is expressed as incremental form, while in linear elastic conditions, the stress-strain relation is defined by elastic matrix ( D e ) as follows (Ugural 2003):

Fluid mass balance equation
The mass balance equation for the water phase is written as follows (Lewis et al. 1998): where K s and K l are the bulk modulus for solid and fluid phases, respectively. The expression ( − K s + K l ) in above equation is equal to the inverse of Biot modulus (M). v s denotes the velocity of the solid phase, and q f is fluid flux. Here, it is assumed that the flow of fluid follows the Darcy's law, so: where k , p and g represent the solid phase permeability, fluid pressure in the porous medium and gravity acceleration of the earth, respectively. f and ρ indicate fluid viscosity and density, respectively. These two parameters are assumed to be constant due to the conditions of the model.

Discretization of governing equations
Since the purpose of this paper is to investigate the geomechanical behavior of faults (as discontinues medium) inside the reservoir (as continuous medium), Therefore, the equations for each continuous and discontinues mediums should be separately developed and discrete, and finally, using the fully coupled approach, the equations are solved through a computer code. The main unknown of these equation is the displacement and pressure. In this study, the standard Galerkin finite element method is used to solve the governing equations for both continuous and discontinuous domains. In this regard, 4-nodes elements are used for pore pressure and 9-nodes elements for displacement parameters. Also, the interface element (element with zero thickness) is used to solve and discretize the fault in modeling. In this regard, 4-nodes interface elements are used for stresses and displacement parameters and 6-nodes interface elements for pore pressure. Figure 1 shows the elements used in this study. In this figure, ξ and η represent the axes of local coordinates in each element.
In the process of finite element method estimation, the displacements and pore pressure are expressed as vectors of nodal values of the displacement and the pore pressure in the limited points of space. In order to ensure the continuity of variables, shape functions (N) are used to express variations of variables within a single element. In this study, the shape functions used for the pressure and displacement variables are different. These two types of functions are depicted with N u and N p for displacement and pressure shape functions, respectively. We have: The stress and strain in the Cartesian coordinate system and for the continuous medium are expressed in the following form: The above relations for the discontinuous fault medium are changed as follows:

Discretization of fluid mass balance equation
By applying the standard Galerkin method and using the principle of virtual work, as well as applying Eqs. (5), (1) will discretized as follows: Fig. 1 The square elements used to interpolate the variables. a 4-nodes element for pressure b 9-nodes element for displacement and their correspond zero-thickness elements for fault where K , Q and F u represent the elastic stiffness matrix, the coupling matrix and the vector of the nodal forces, respectively, which are defined as follows: In the above relations, B is strain-displacement matrix which calculated as follows (for two-dimensional continuous medium) (Lewis et al. 1998): L is the derivative operator. In the discontinuous fault medium, the matrix B is determined as follows: And the fluid pressure shape function ( N p ) matrix is also determined for 4-nodes element and 4-node interface element as follows: In continuous medium and for square elements (matrix) In discontinuous medium and for interface elements (fault) Also, the elastic matrix D m e is calculated with the assumption of plane strain conditions and in two-dimensional mode as follows: where E and represent the Young's modulus and Poisson's ratio, respectively. Regarding the discontinuous medium, this matrix is defined for a zero-thickness interface element in plane strain conditions as follows: where k s and k n represent normal and shear stiffness of fault, respectively. The characteristics of the interface element are determined by the conditions governing the corresponding physical surface. Of course, it should be noted that during numerical modeling to ensure continuity of displacement,k s (23) D f e = k s 0 0 k n and k n values are considered to be relatively large for these surfaces, which can be described mathematically as parameters of the penalty method.
It should be noted that due to the possibility of angularity of the faults and therefore the interface elements in compared with main network of model, the rotational matrix of R is defined as follows: where In this case, the first and second expressions of Eq. (14) are rewritten as follows: Of course, due to the fact that the half of the nodes of the interface elements are placed on each other, and all of the components are positive in the N m p matrix, the values of coupling matrix (Q) are calculated twice as much as the real value, Therefore, the value of 0.5 should multiplied to the total Q-matrix. Also, the right expressions of Eq. (14), which represent external forces, are removed after the final assembly of the matrices in common nodes of the square elements of the continuous medium and the fault interface elements. So, their calculations in the interface elements are not necessary and can be considered as zero value.
Instead of the tensor I in Eq. (1), the matrix m is used in discrete equations. This matrix is determined in the twodimensional mode in the calculations of continuous reservoir domain in the form of m m = [1, 1, 0] T and in the discontinuous fault domain in the form of m f = [0, 1] T Finally, the discrete equations of fluid mass equilibrium will be as follows: In continuous medium In discontinuous medium

Discretization of momentum balance equation
Here, the standard Galerkin method has applied to Eq. (6), and by use of Eq. (7) and other relations expressed in Eqs. (8) and (9), following series of linear equations are derived: where H, S and F p represent the flow stiffness matrix, flow capacity matrix and flow sources vector, respectively, which defined as below: The above equations are valid in continuous part of reservoir, while in the same way as in the section of discretization of fluid mass balance equation, here also the discrete equations for the discontinuous part are slightly different from the continuous part of the reservoir. In the following, these differences are described in detail.
For the fluid flow along the fault, longitudinal transmissivity term (κ) and also cubic law for fault permeability should be considered in Eq. (31) (Ng and Small 1997). According to this law, the permeability of a fault with aperture of b is: Also, due to the multiplication of the two matrixes of the shape functions in this equation, the value of 0.25 is multiplied in this matrix. It is due to the reason mentioned in the previous section. Therefore, the corrected matrix for the discontinuous fault medium will be follows: And also, if the fluid flow perpendicular to the fault is applied into equations, the components of the matrix H should be calculated by following equation: where ∼ κ t is transversal conductivity coefficient of fault, and N f P * is defined as: It should be noted that the value of the Biot modulus (M) in Eq. (32) for the faulted medium is determined by the ratio of the fault aperture to the fluid compressibility (Segura and Carol 2008b). Finally, calculating the matrices of Eqs. (14) and (30) in two continuous and discontinuous mediums will form a series of equations that will solve the coupled fluid flow-geomechanics problems in a saturated compressible medium with discontinuities: This series of equations can be solved using the initial boundary conditions defined in each problem as well as the time discretization. In this study, a trapezoidal method with a (34) k frac = b 2 12 value of θ = 0.6 has been used for time discretization of above equations.
It should be noted that since the value of the term of the transmissibility of the fault in Eq. (35) depends on the value of its aperture openings in each stage, the system of Eqs. (38) is completely nonlinear. There are at least three methods to deal with this issue. The first is using a simple iterative method. This means that the values of the solution of the last step are used to evaluate and compute the coefficients in a new time step. The second method is the direct iterative method, which is also referred as fixed point method. This method, at each time step, uses the unknown values in each last repeat to evaluate and compute the coefficients in the next iteration, and the convergence is obtained when the error between the two successive repetitions is less than the threshold (Dlala et al. 2007). The third method is the Newton-Raphson method, which its approach is very similar to the second method, except that the convergence rate is higher (Saitz 1999). In this study, the Newton-Raphson method has been used to linearize the equations.

Verification and model application
In order to validate the method presented in this study, three index problems from various papers have been solved with this method. The first problem is from Watanabe et al. (Watanabe et al. 2012), in which the effect of different orientation of a single fracture on displacement jump has been evaluated. The second problem is a numerical example from Ng and Small (Ng and Small 1997), in which the hydromechanical behavior of a vertical fracture in the porous medium is investigated. Many researchers, such as Segura and Carol (Segura and Carol 2008a, b), have used this problem to validate their model and have solved this problem. The third problem is to investigate the effect of a single horizontal fault on the discharge of a porous reservoir. This is introduced in the study of Guiducci et al. (Guiducci, et al. 2003). It should be noted that the validity of all the results has been evaluated by comparing the analytical responses and the responses reported by other researchers. Noted that the model code is written using MATLAB software.

First problem: fractured rock sample under uniaxial load
In this case, the rock sample is in small dimensions, with discontinuities and under uniaxial pressure. The purpose of this study is to investigate the effect of the discontinuity of the rock sample on the displacement functions as well as relative displacement. In this paper, this problem is used to validate the proposed method in an uncoupled state and also without consideration of fluid flow, using existing analytical responses. The analytical solution to this problem is presented in Watanabe et al. (2012). The data used in this paper is similar to the data used in Watanabe et al. (2012). These data are reported in Table 1. Figure 2 also shows the geometry of the problem and its boundary conditions. Since fluid flow has not been considered in this problem, only 9-node square elements and 6-node interface elements (Fig. 1) are used to spatial discretization of the problem environment. As shown in Fig. 2, the fracture is assumed to be in different orientations within the rock sample. In this study, the angles 0, 15, 30, 45 and 60 degrees are examined. Figure 3 shows the numerical results obtained from the model presented in this paper as well as the analytic results of the problem for the vertical displacement profile along the line AA 'for all angles mentioned above. As shown in this figure, numerical results show very good agreement with analytical results. Also the results show that the highest displacement jump occurs when α = 60°.

Second problem: consolidation of saturated porous medium with vertical discontinuity
In this case, a square block is considered as a porous medium with vertical discontinuity. The purpose of this issue is to study the consolidation of such medium under plane strain condition and also considering the linear elastic behavior of the entire problem environment. This problem has also been solved in Ranjbar et al. (2020).The purpose of expressing this problem is to express the ability of the proposed method in modeling fluid flow in a compressible porous medium with vertical discontinuity. Figure 4 shows the boundary conditions of displacement and pressure as well as the related model mesh with the help of the elements presented in this paper.
As shown in Fig. 4, 9-nodes and 4-nodes elements for the continuous part, as well as 6-nodes and 4-nodes interface elements with zero thickness for the discontinuous part, have been used so that the common points of these elements are adapted to each other. Table 2 shows the values of the data used for this problem.
Conditions of this problem are such that the increase in reservoir pore pressure caused due to upper boundary applied load. The boundary of the discharge of the fluid is top boundary, so the pore pressure value for this boundary is considered zero. In this problem, the fracture opening is variable, and the permeability values are calculated according to the cubic law.
With the start of the process of applying the upper load and also the discharge of the fluid from the upper boundary, a decrease in the excess pore pressure from the upper boundary to the lower boundary could be observed, which results in effective stress reduction. Due to such behavior of pore pressure, the values of discontinuity opening will be vary in the upper boundary and the lower bound, so that, as time elapses, the fracture closing begins with the upper boundary and progresses gradually toward the lower boundary.
The various excess fluid pressures obtained using the model presented in this paper are shown in Fig. 5 for different time coefficients. The relationship of time coefficient (T) is introduced in the study of Ng and Small (Eq. 39) (Ng and Small 1997). The results are very consistent with the analytical results. Figure 6 schematically shows the closure of the discontinuity of the problem with a magnification of 100 at T = 0.9. Figure 7 also shows the where T, c v , t and d are time coefficient, coefficient of consolidation, time and drainage path, respectively.

Third problem: the effect of fluid flow in a horizontal fault
The issue of simulating geomechanics and fluid flow in a medium with a horizontal fault has been evaluated and studied. Many researchers have used this problem to validate their models (Segura and Carol 2008a;Guiducci, et al. 2003). Used this problem to study the effect of a single horizontal discontinuity on the discharge of a reservoir. Also, Noorishad et al. (1982), Bower and Zyvoloski (1997) and Watanabe et al. (2012) have used this problem to study the fluid injection inside a discontinuity that is surrounded by a low permeability rock matrix. In this paper, this problem is also used to identify the geomechanical behavior of a single fault inside the reservoir due to fluid injection as well as validating the model presented throughout the paper. It should be noted that the data used in this article is taken from the study of Watanabe et al. (2012) and Bower and Zyvoloski (1997). The semi-analytical solution of this problem is introduced in Wijesinghe's study (Wijesinghe 1986). The domain of this problem consists of a block of rock with a low permeability and a length of 25 m, divided by a horizontal fault from the middle to the two sub-black. Because of the weight of the overlaying layers uniform stress equal to 50 MPa applied from top boundary of model. It is assumed that an injection well and a production well are located on the left and right side of the reservoir, respectively. The pore pressure of the fluid in the initial state is 11 MPa. In the early stages of the simulation, the fluid is injected at a pressure of 11.9 MPa through a well that is located on the left side of the reservoir into a fault which it's aperture opening equal to 0.01 mm. In this case, we will see an increase in the pressure inside the reservoir. Figure 8 shows the geometry of this issue, along with the elements used. Table 3 lists the data used. Here, as with the two previous issues, 9-node and 4-node elements for the continuous part of reservoir and 6-node and 4-node interface elements with zero thickness for the discontinuous part are used.
In Figs. 9 and 10, respectively, simulation results regarding the fluid pressure along the fault and the fault aperture opening rate are compared with the results presented in Watanabe et al. (Watanabe et al. 2012) and the semi-analytical results. As can be seen in these figures, the results show very good agreements with semi-analytical results. Also, by comparing the results of prepared model and the model presented in Watanabe et al. (Watanabe et al. 2012), it can be seen that the simulation results are approximately equal to each other.

Conclusions
The focus of this paper is to study the geomechanical behavior of faults in the reservoir due to fluid flow. In this regard, the process of coupled geomechanics and fluid flow simulation in deformable discontinues porous medium has been investigated. The presented numerical modeling is based on the finite element method. In this study, square elements (4 nodes and 9 nodes) and zero-dimensional elements (4 nodes and 6 nodes) were used for continuous and discontinuities parts of problems, respectively. The method and model presented in this article have been evaluated using analytical and numerical examples from other researchers in this field. Three index problems were selected as the basis issues for verification. The results of the model show a very good agreement with analytical and other reported results, which indicates the accuracy and capability of the described method in coupled fluid flow-geomechanics simulations.

Funding
The authors received no finantial support for the research, authorship, and/or publication of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.