Embedded Fracture Model for Coupled Flow and Geomechanics
EEmbedded Fracture Model for Coupled Flow and Geomechanics
I. Shovkun and T. Garipov and H. A. Tchelepi
Energy Resources Engineering, Stanford University, USA
Abstract:
Fluid injection and production cause changes in reservoir pressure, which result in deformationsin the subsurface. This phenomenon is particularly important in reservoirs with abundant fractures and faultsbecause the induced slip and opening of the fractures may significantly alter their hydraulic properties. Mod-eling strongly coupled poro-mechanical processes in naturally fractured reservoirs is a challenging problem.The Discrete Fracture Model (DFM) is a state-of-art method for modeling coupled flow and mechanics infractured reservoirs. This method requires constructing computational grids that comform to fractures, whichis very challenging in complex 3D settings. The objective of this study is to develop a numerical methodthat does not require gridding near fractures and can e ffi ciently model hydromechanical interactions in frac-tured reservoirs. We utilize formulations based on the Strong Discontinuity Approach (SDA) for mechanicsand Embedded Discrete Fracture Model (EDFM) for flow. We first present a mathematical formulation andemphasize the kinematic aspects of fracture slip and opening. We then introduce a series of mechanical teststhat investigate the spatial convergence of the model and compare its accuracy with the Discrete FractureModel (DFM). We finally consider a synthetic coupled case of a reservoir with several fractures and com-pare the performance of the SDA and DFM methods. Our results indicate super-linear spatial convergenceof the proposed SDA algorithm. Numerical simulations confirm the applicability of the proposed method tomodeling the coupling e ff ects in subsurface applications. Keywords : Embedded Discrete Fracture Model, Strong Discontinuity Approach, Geomechanics, ReservoirSimulation. 1 a r X i v : . [ c s . C E ] A ug Introduction
Field operations in petroleum, geothermal, and waste disposal applications frequently involve injecting andwithdrawing fluids from the subsurface. Changes in reservoir pressure in naturally fractured reservoirs canlead to reactivation of natural fractures and faults [Grasso, 1992, Zoback and Zinke, 2002, Rutledge et al.,2004, Zoback and Gorelick, 2012, Hwang et al., 2015]. Fault reactivation often leads to shear of the wellborecasing or damage of the reservoir integrity [Elf-Aquitaine, 1992, Wiprut and Zoback, 2000]. In tight rocks,fracture reactivation may cause significant changes in reservoir permeability [Min et al., 2004].Modeling fracture reactivation is possible with a variety of numerical methods. It has been shown that the Fi-nite Element Method (FEM) with discrete fracture representation (DFM) is an accurate and e ffi cient methodfor modeling discontinuities in porous media [Cappa and Rutqvist, 2011, Fu et al., 2016, Salimzadeh et al.,2018]. This method, however, requires the usage of sophisticated gridding techniques for preparing high-quality computational grids. The Extended / Generalized Finite Element (xFEM) and Phase-Field methodsallow to handle arbitrary fracture geometries and, therefore, suit for modeling fracture propagation [Duarteet al., 2000, Heister et al., 2015]. These methods introduce additional global variables to improve the repre-sentation of the deformation field [Haddad and Sepehrnoori, 2015, Shovkun and Espinoza, 2018]. This leadsto a significant computational overhead compared to conventional FE-based methods. Several approachesto modeling fractures have been developed based on non-local flomulations of continuum mechanics. Peri-dynamics and Discrete Element methods are among the most capable techniques for modeling complexgeometrical fracture configurations [Zhao et al., 2009, Ouchi et al., 2015, Sun et al., 2016]. Similarly toxFEM and Phase-Field, however, these non-local methods have very high computational cost. The Bound-ary Element and Displacement Discontinuity methods are very e ffi cient techniques; they are widely appliedto mechanical problems with fractures [Olson, 2008, Sesetty et al., 2012]. These methods allow to modelreservoirs with high numbers of fractures; however, they cannot capture heterogeneous and anisotropic rockproperties [McClure, 2012]. It has been shown that the Boundary Element methods can be combined withthe Finite Element methods to solve coupled hydro-mechanical problems [Norbeck et al., 2016].Another type of Finite Element mechanical models is based on the Strong Discontinuity Approach (SDA)[Oliver, 1996]. This method uses the “embedded” fracture representation and can be applied to modelfracture activation and propagation [Oliver et al., 1999]. The method shares similarities with xFEM butcaptures the deformation of discontinuity locally [Linder and Armero, 2007]. The nonlinear return mappinglgorithm allows SDA to reduce the number of degrees of freedom to that of a linear elastic system withoutfractures [Regueiro and Borja, 2001]. This localization allows seamless incorporation of the method intothe existing Finite Element software and renders the method very attractive from the numerical simulationstandpoint [Mosler, 2005].There are several numerical techniques to simulate fluid flow in fractured formations. Dual permeability anddual porosity models [Warren et al., 1963], the Discrete Fracture Network approach (DFN) [Cacas et al.,1990], Discrete Fracture Models [Karimi-Fard et al., 2004], and the Multiple Interacting Continua (MINC)[Narasimhan and Pruess, 1988, Pruess, 1992] methods are widely used in the industry. Several approachesbased on local permeability modification have been proposed to model flow in fracture media in an upscaleddual-porosity fashion [Gong et al., 2008, Li et al., 2015]. The application of the Finite Volume method withthe embedded fracture representation (EDFM) has been investigated for the past decade [Li et al., 2008, ? ].While substantial research exists on using SDA in solid mechanics, very few studies investigate the appli-cation of SDA to coupled flow-mechanical problems. The approach of enriching a global variable to ap-proximate a discontinuity used in SDA was also applied to model discontinuities in fluid pressure [Alfaiateet al., 2010]. The authors, however, did not consider hydraulically active fractures with flow within them.A cobbination of SDA enrichments for both displacement and flow discontinuities can be used to model thecoupled behavior of fractured media [Callari and Armero, 2002]. This approach and the aforementionedSDA model for flow [Alfaiate et al., 2010] do not allow to model highly conductive fractures.A recent advancement combined a mechanical model with embedded fracture representation and a conven-tional EDFM for flow [Deb and Jenny, 2017]. In their work, Deb and Jenny utilized the Extended FiniteVolume method with a stabilization term for fracture activation. Their model, however, is limited to a 2Dsetting and slipping fractures, and it treats the fracture jump as an independent variable. An approach thatconsists of a Finite Element mechanical model with SDA and an embedded fracture model for flow wouldallow to model reservoirs with conductive fractures e ffi ciently.In this manuscript, we present a new hybrid numerical method for modeling coupled fluid flow and geome-chanics with embedded fracture representation. We combine an FEM-based mechanical SDA model andan FVM-based flow EDFM model, which, to the best of our knowledge, has not been done before. Themechanical formulation utilizes SDA with the return-mapping algorithm [Borja, 2013]. This paper consistsof comparison cases between the proposed EDFM model and the DFM model described in [Garipov et al.,016]. We first investigate the spatial convergence of the mechanical SDA model in plane strain scenarioswith single fractures. Second, we investigate the spatial convergence of the EDFM flow model on a single-phase stationary problem. We finally compare the results produced by the EDFM and DFM models in acoupled 3D case with several fractures. In this section, we introduce a system of equations that describes coupled fluid flow and rock deformationsin a saturated porous medium with fractures. We employ linear poroelasticity equations and Darcy’s lawto describe the single-phase flow and mechanics within a compressible porous medium [Coussy, 2004].We then outline the governing equations for single-phase flow within fractures and corresponding fracturesdeformations.
The mass conservation of a single-phase fluid in a porous medium follows the continuity equation ∂ (cid:16) ρ f φ (cid:17) ∂ t = ∇ · (cid:16) ρ f v (cid:17) + Q FM , (1)where ρ f is the fluid density, φ is the rock porosity, v is the fluid superficial velocity, Q FM is the massflux from the fracture into the matrix, and t is the time. Assuming that in an isothermal system a sligtlycompressible fluid d ρ f d p M = c f ρ f , (2)complies with the Darcy’s law v = − k M µ f (cid:16) ∇ p M − ρ f g (cid:17) (3)and that the poroelastic solid is linear and isotropic d φ = b d ( ∇ · u ) + K (1 − b ) ( b − φ ) d p M , (4)q. (1) becomes [Coussy, 2010]:1 M ∂ p M ∂ t + b ∂ ( ∇ · u ) ∂ t − ∇ · k M µ f (cid:16) ∇ p M − ρ f g (cid:17) = q MF , (5)where u is the solid matrix displacement, b is the rock Biot coe ffi cient, k M is the matrix permeability, µ f is the fluid viscosity, ρ f is the fluid mass density, M = φ c f + / K (1 − b ) ( b − φ ) is the Biot modulus, K is the drained bulk modulus of the rock, g is the gravity acceleration constant, and q MF = Q MF /ρ F is thevolumetric matrix-fracture flux. The mechanical equilibrium of the system is described by the quasi-static Cauchy equation: ∇ · S = − ρ b g (6)where S is the total stress tensor and ρ b is the bulk mass density of the saturated rock. We use the followingisotropic relation between total and e ff ective stresses [Zoback, 2013]: S = σ + bp M I (7)where σ is the e ff ective stress tensor, p is pore pressure, b is the scalar Biot coe ffi cient, and I is the second-order unit tensor. We further assume a linear relationship between the stress and strain tensors: σ = (cid:131) : ε , (8)where (cid:131) is the fourth-order linear elastic sti ff ness tensor. The governing equation of the fluid flow in fractures is very similar to that in the matrix. We neglect, however,the coupling terms that describe the changes in matrix aperture permeability due to pressure perturbations: c f ∂ p F ∂ t − ∇ · k F µ f (cid:16) ∇ p F − ρ f g (cid:17) = q FF − q MF (9)here c f is the fluid compressibility, p F is the fluid pressure in the fracture, F is the fracture permeability, q FF is the flux between intersecting fractures. The computation of the terms q FF and q MF is discussed inSection 3.1. Ω Γ n ux [u] ux[u]u nc u c (a) (b) (c)u ux Figure 1: (a) A body Ω that contains a discontinuity Γ with a normal n . (b) Displacement field exhibitsa jump across the fracture. (c) Decomposition of the displacement field into a conforming (regular) andnon-conforming (singular) parts. Kinematics:
We further consider a domain Ω separated by a two-dimensional surface Γ that represents afracture (Fig. 1a). Since fractures are discontinuities, the displacement field u is discontinuous (Fig. 1b) andcan be decomposed as [Simo and Oliver, 1994] u = u + [ u ] H Γ , (10)where [ u ] is the jump vector and H Γ is the Heaviside step function that is equal to zero on one side of Γ andequal to one on the other side. SDA modifies Eq. 10 by adding and subtracting a continuous scalar function f u = u + [ u ] f (cid:124) (cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32) (cid:125) u c + [ u ] ( H Γ − f ) (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) u nc (11)so that u c is the conformal part of displacement, u nc is the non-conformal part of displacement [Simo andOliver, 1994] as shown in Fig. 1c. The form of the level-set (or ramp) function f depends on the discretizationand is discussed elsewhere [Oliver et al., 2003, Foster et al., 2007]. An expression for the small strain ε canbe obtained by taking a symmetric gradient of Eq 11: ε = ∇ s u = ∇ s u c − ( [ u ] ⊗ ∇ f ) s + ( H Γ − f ) ∇ s [ u ] + ( [ u ] ⊗ n ) s δ Γ (12)here δ Γ is the Dirac delta-function and n is the unit normal vector to the fracture surface. We neglectthe second term in the right-hand side of Eq. 12 due to the jump-discontinuity assumption, so that Eq. 12becomes [Mosler, 2005] ε = ∇ s u c − ( [ u ] ⊗ ∇ f ) s + ( [ u ] ⊗ n ) s δ Γ (13) Plasticity formulation for SDA:
We incorporate SDA into the non-associated plasticity framework de-scribed by the system [Oliver, 1996, Mosler, 2005, Foster et al., 2007]:˙ σ = (cid:131) : (cid:0) ˙ ε − ˙ ε p (cid:1) (14a)˙ ε p = λ ∂ G ∂ σ (14b)˙ q = λ H ∂ G ∂ q (14c) F ( σ , q ) = ε p is the plastic strain, λ is the plastic multiplier, q is the stress-like internal variable, G and F arethe plastic potential and yield function, respectively, and H is the softerning / hardening parameter. In system(14), Eq. 14a is the stress-strain relation. Eq. 14b described the plastic strain evolution. Eq. 14c is thesoftening / hardening law, and Eq. 14d is the yield surface. We emphasize that the system (14) contains thee ff ective stress σ as opposed to the total stress S [Callari et al., 2010].The governing system of equation for SDA can be obtained by substituting Eq. 13 into system (14), assumingthat the stress σ is bounded and that the distributions of the plastic multiplier λ = λ δ δ Γ and hardeningmodulus H = H Γ δ Γ are singular [Simo et al., 1993, Mosler and Meschke, 2003]:˙ σ = (cid:131) : (cid:2) ∇ s ˙ u c − ([ ˙ u ] ⊗ ∇ f ) s (cid:3) (15a)[ ˙ u ] = λ δ ∂ G ∂ t (15b)˙ q = − λ δ H δ ∂ G ∂ q (15c) F ( t , q ) = In this section, we specify the flow rule and plastic potential for three possible states of a fracture: slip,opening, and stick. lip:
During slip, the jump vector in the fracture is nearly-parallel to its surface. We use the following flowrule and plastic potential to impose friction in the fracture: F ( t , q ) = t τ − µ t n − q (16a) G ( t , q ) = t τ − θ t n − q (16b)where t τ = n · σ · τ is the tangent traction, t n = n · σ · n is the normal traction, µ is the tangent of the frictionangle, and θ is the tangent of the dilation angle. Eq. 16b transforms Eq. 15b into[ ˙ u ] n = θλ [ ˙ u ] τ = λ (17)where [ ˙ u ] n and [ ˙ u ] τ are the normal and tangential components of [ ˙ u ]. Opening:
For an open fracture, we propose to use three components of the jump function [ u ] as inde-pendent local variables. We select the flow rule and plastic potential to enforce the absence of all tractioncomponents independently: F ( t , q ) = G ( t , q ) = t n − qF ( t , q ) = G ( t , q ) = t τ − qF ( t , q ) = G ( t , q ) = t τ − q (18)Note that in the case of an opening fracture the plastic potential G and flow rule F are equal (associativeplasticity) and have three components. This choice of the plastic potential transforms Eq. 15b to the followingform: [ ˙ u ] τ = λ [ ˙ u ] τ = λ [ ˙ u ] n = λ (19) Stick:
The case of a closed non-slipping fracture (stick case) also requires special treatment. Althoughduring loading, a rock with a fracture behaves elastically until the yield conditions are met, it is important toenforce the zero-normal jump and no-slip conditions during the unloading of a fracture:[ u ] n =
0[ ˙ u ] τ = Discretization
Our numerical model for embedded fractures utilizes a hybrid formulation with first-order Galerkin dis-placement approximation, Petrov-Galerkin for the enhanced strain approximations, and piecewise-constantfinite-volume representation of pressure. The distribution of variables in a cell is shown in Fig. 3. Due tothe choice of shape functions, the conforming displacement is defined at grid nodes. As discussed in Sec-tion 2.2.2, the displacement jump [ u ] on the fracture is assumed to be discontinuous and is approximated asa piecewise-constant function. We define the displacement jump at the cell centers, although, as discussed inSection 3.2, its location is irrelevant. Fluid pressure in the rock matrix and fracture is defined in the matrixand fracture control volume centers, respectively. [u],p M u c p F Figure 2: Location of variables in the coupled formulation. The conforming displacement u c is defined incell vertices. The jump [ u ] and the fluid pressure in the rock matrix p M are defined at the cell centers. Thefluid pressure p M in the fracture is defined at the center of the fracture control volume.In the following sections, we provide the relevant details regarding the numerical treatment of fractures. Wefirst discuss fluid flow and then provide the details of discretizing the mechanics equilibrium equations. The EDFM approach for approximating fluid flow flow consists of approximating the matrix-fracture flow q MF and fracture-fracture flow q FF terms in Eq. 5 and 9. We opt to select a simplified approach first intro-duced in [Li et al., 2008]. More accurate schemes are available as discussed in Section 5.3.We compute the total mass flux between a reservoir control volume and a fracture control volume as: Q MF = A n · k R · n d ρ f µ f ( p M − p F ) , (21)where the coe ffi cient d is the average distance between the points in the reservoir cell and the fracture cells.When a fracture element with its center in M F bisects a cell with the volume V and the center in M , thenhe initial volume splits into shape I with the center in M and volume V and shape II with the center in M and volume V (see Fig. 3). Obviously, the sum of the volumes of the shapes I and II is equal to the originalcell volume V + V = V . Then the coe ffi cient d can be evaluated as d = V V || M − M F || + V V || M − M F || (22) (a) (b) DFM reservoir cellM M M M F DFM fracture CVEDFM fracture CVF F F F F F EDFM reservoir cellFracture CV
Figure 3: Fluid domain discretization. (a) Control volumes (CV) for EDFM and DFM models. (b) Schemat-ics for fracture-fracture transmissibility calculation.We next discuss the flux term that arises at the intersection of two fractures. We assume that the flux betweentwo intersecting fracture elements F and F is proportional to the pressure di ff erence (see Fig.3b): Q = T ρ f µ ( p − p ) (23)Hereafter, we denote the transmissibility between two segments i and j by T ij . To compute the transmissi-bility T between two intersecting fracture elements F and F , we split these elements by the intersectionline into the segments F , F and F , F , respectively Fig. 3b). Then the transmissibility T is obtainedwith sum of the corresponding transmissibilities: T = T + T + T + T , (24)where the transmissibilities between the fracture segments (in the right-hand side of Eq. 24) are computedwith the star-delta transformation [Karimi-Fard et al., 2003]. .2 Mechanics We now discuss the discretization of the mechanical system that consists of Eq. 6 and the system (15). Usingthe enhanced assumed strain approach, we apply finite-element approximations to Eq. 6 [Simo and Rifai,1990]: (cid:90) Ω e ∇ s η : S d Ω = (cid:90) Ω e ρ g d Ω , (25a) (cid:90) Ω e γ : S d Ω = , (25b)where η denotes a continuous test function, γ is the variation of the enhanced strain, and Ω e is an elementthat contains a fracture. Eq. 6 transforms into a system of two equations due to singular strain across thefracture. To solve the coupled system (6) and (5), we use the contact integral in Eq. 25b as an additionalconstraint. This expression must be modified for the case of piecewise constant jump values [Borja, 2000]and can be replaced by the local stress-continuity condition [Regueiro and Borja, 1999, Callari et al., 2010]:1 V e (cid:90) Ω e S · n d Ω = t + p f n (26)Eq. 26 relates the average total stress S in an element to the fracture traction t and fracture pressure p F .Hereafter, we use the following notation ( • ) : = V e (cid:90) Ω e ( • ) d Ω . Therefore, in order to complete the system(15), we use the volume average of Eq. 15a˙ σ = (cid:131) : (cid:104) ∇ s ˙ u c − (cid:16) [ ˙ u ] ⊗ ∇ f (cid:17) s (cid:105) , (27)and the following relationship between the average e ff ective stress and fracture traction: t + p F n = σ · n + bp M n . (28)We emphasize that the average e ff ective stress σ and jump [ u ] are functions of the matrix and fracturepressures and have the corresponding terms in linearization as discussed in Section 3.2.2. Note the absenceof the overline above the jump [ u ] as [ u ] = (cid:2) u (cid:3) due to the piecewise-constant jump assumption.In our formulation, we solve the local system of equations (15), (27), and (28) for all the elements containingfractures using the Newton-Raphson method. We use slip rate [ ˙ u ] , average stress tensor σ , and the Lagrangeultiplier λ as primary unknowns. Nodal displacement u c , matrix and fracture pressure p M and p F are fixedduring the nonlinear iterations. In order to obtain quadratic convergence for the system (25), a consistentstress linearization must be performed [Borja, 2013]. Eq. 27 and Eq. 15b-d constitute a complete system, which is used to find the jump [ u ] and traction vector t from the strain ε . After computing these quantities, the e ff ective stress σ in Gauss quadrature points can beobtained with Eq. 15a. The e ff ective-stress derivatives can be computed from the following linearization:d σ = (cid:131) : (cid:34) I + (cid:32) ∂ [ u ] ∂ ε ⊗ ∇ f (cid:33) s (cid:35) d ε + (cid:131) : (cid:32) ∂ [ u ] ∂ p f ⊗ ∇ f (cid:33) s d p f + (cid:131) : (cid:32) ∂ [ u ] ∂ p m ⊗ ∇ f (cid:33) s d p m . (29)The derivatives ∂ [ u ] ∂ p f and ∂ [ u ] ∂ p m can be obtained from the system (15). To evaluate these derivatives and obtainquadratic convergence of the Newton scheme, we use an approach based on the inverse theorem and auto-matic di ff erentiation as discussed elsewhere [Garipov et al., 2018]. Using Eq. 7, we write the linearizationfor the total stress tensor:d S = (cid:131) : (cid:34) I + (cid:32) ∂ [ u ] ∂ ε ⊗ ∇ f (cid:33) s (cid:35) d ε + (cid:131) : (cid:32) ∂ [ u ] ∂ p f ⊗ ∇ f (cid:33) s d p f + (cid:131) : (cid:34)(cid:32) ∂ [ u ] ∂ p m ⊗ ∇ f (cid:33) s + b I (cid:35) d p m (30)We emphasize that the linearization of the total stress with respect to the matrix pressure p m is not isotropicdue to the term (cid:16) ∂ [ u ] ∂ p m ⊗ ∇ f (cid:17) s . In this section we validate presented mechanical model with plane-strain analytical solutions and compareits spatial convergence with that of the DFM model. The simulation cases presented in this section arecomputed in a domain with a single cell in the third (z) direction since the available analytical solutionsare derived based on the plain-strain assumption. The domain size was assumed ten times larger than thefracture length to mimic an infinite plane. Further we use the term EDFM (embedded discrete fracturemodel) to refer to both the mechanical and flow models. In the context of mechanics, the terms EDFM andSDA are interchangeable. .1 Mechanical tests
First, we consider an inclined fracture slipping under compressive loading conditions. Here, a single fracturewith length l =
10 m and striking at α = ◦ is slipping due to the applied compressive load σ x =
10 MPa( σ y = u τ ] = (cid:16) − ν (cid:17) E σ x sin α (cos α − µ sin α ) (cid:113) l − ( x − l ) (31)where E is the rock Young’s modulus and ν is the rock Poisson’s ratio.Second, we consider an open fracture test. In this case, an x-aligned fracture ( α =
0) opens due to the appliedtensile load σ y =
10 MPa ( σ x = u n ] = − νµ σ y (cid:113) l − ( x − l ) (32)where µ is the rock shear modulus. All the geometrical and mechanical parameters used in the simulationsare listed in Table 1. α
10 m σ x200 m m σ y (a) (b) (c) αα Figure 4: (a) Simulation domain. (b) A fracture conforming to the grid. (c) A fracture not conforming tothe grid.roperty ValueDomain size [m ] 200 × × E [MPa] 10 Matrix Poisson’s ratio ν [-] 0.25Friction coe ffi cient µ [-] 0.6Dilation coe ffi cient θ [-] 0.0Initial cohesive strength q [MPa] 0.0Hardening parameter H [MPa] 0.0Table 1: Geometrical and material parameters used in the mechanical simulations. We first present modeling results obtained on a domain with the grid conforming to the fracture (Fig. 4b).We performed several simulations on hexahedral grids with various numbers of cell elements containing thefracture n F =
4, 8, 16, and 32. Identical cases with the same grid element sizes were simulated with theDFM model. The tangential jump on the fracture surfaces given by EDFM and DFM models at various gridrefinement levels, as well as the analytical solution, are shown in Fig. 5a-b. All slip values are normalizedby the maximum of Eq. 31. The coordinates are normalized by the fracture length l . The solution obtainedfrom the EDFM model overestimates the jump value and converges monotonically to the analytical solution(Fig. 5a). In contrast, the solution obtained from the DFM model underestimates the jump and also convergesmonotonically to the analytical solution (Fig. 5b).The L -error in the solutions produced by the EDFM and DFM models as functions of the grid elementsize is shown in Fig. 5b. Both errors given by DFM and EDFM models reduce while refining the gridat approximately the same rate. In the same figure, we provided linear and quadratic trends to illustrate theorder of spatial convergence of the models. The errors given by DFM and EDFM models exhibit super-linearconvergence. .0 0.2 0.4 0.6 0.8 1.0Normalized coordinate0.00.20.40.60.81.01.21.41.61.8 N o r m a li z e d t a n g e n t j u m p r e f i n e m e n t r e f i n e m e n t (a) DFM analyticalh=2.500 mh=1.250 mh=0.625 m 0.0 0.2 0.4 0.6 0.8 1.0Normalized coordinate0.00.20.40.60.81.01.21.41.61.8 N o r m a li z e d t a n g e n t j u m p (b) EDFM analyticalh=2.500 mh=1.250 mh=0.625 mh=0.312 m 0.3 0.6 1.0 3.0 6.0Element Size [m]10 L E rr o r E D F M DFM L i n e a r Q u a d r a t i c (c) Convergence Figure 5: Spatial convergence exhibited by the EDFM and DFM models for the shear problem. (a) Tangen-tial jump as a function of coordinate obtained at various refinement levels. Note that DFM model underes-timates the jump, whereas EDFM overestimates it. (b) L -error in the tangential jump value as a functionof the element size h. Both models manifest super-linear convergence. The error in the tangential jumpobtained with the DFM model is lower than that obtained with the EDFM model.For the second test problem, we also performed a convergence analysis of the DFM and EDFM models onCartesian grids with the numbers of grid comprising the fracture n F =
2, 4, 8, 16, and 32. The fractureaperture given by the EDFM and DFM models at various grid refinement levels, as well as the analyticalsolution, are shown in Fig. 6a-b. Both models underestimate the aperture value and converge monotonicallyto the analytical solution.The L -errors in fracture aperture as functions of the grid element size are shown in Fig. 6c. Same as in theprevious case, the results obtained from the DFM model manifest super-linear convergence. In contrast, theEDFM convergence in this case is di ff erent: at coarse grids the L error e as a function the of element size h has the slope less steep than the linear trend; at fine grids the slope of the convergence curve is super-linear( e ∝ h . ). Another notable di ff erence with the case of a sliding fracture is that the numerical error given bythe EDFM model at coarse grids (less than eight elements per fracture) is less than that of the DFM model.At fined grids (more than eight elements per fracture) the DFM model predicts the fracture aperture moreaccurately than the EDFM model. .0 0.2 0.4 0.6 0.8 1.0Normalized coordinate0.00.20.40.60.81.0 N o r m a li z e d a p e r t u r e (a) DFM analyticalh=2.500 mh=1.250 mh=0.312 m 0.0 0.2 0.4 0.6 0.8 1.0Normalized coordinate0.00.20.40.60.81.0 N o r m a li z e d a p e r t u r e (b) EDFM analyticalh=2.500 mh=1.250 mh=0.312 m 10 Element Size [m]10 L e rr o r L i n e a r Q u a d r a t i c EDFMDFM (c) Convergence
Figure 6: Spatial convergence manifested by the EDFM and DFM models for the opening problem. (a-b)Fracture aperture as a function of the coordinate obtained with the DFM (a) and EDFM (b) models. (c)L -error in the fracture aperture as a function of the grid element size. DFM model manifests super-linearconvergence. EDFM model exhibits convergence that is worse than linear trend on coarse grids and linearconvergence on fine grids. In this section, we present EDFM simulation results for the cases of slipping and opening fractures with theproblem setup identical to that in Section 4.1.1. The only di ff erence with the previously-presented simula-tions is the grid, which does not conform to the fracture (see Fig. 4c).For the case of a slipping fracture we performed numerical simulations on Cartesian grids with the elementsizes h =
2. m, 1.25 m, and 0.625 m. In these simulations, we also varied the fracture strike angle α = ◦ , 10 ◦ , 20 ◦ , 30 ◦ , and 40 ◦ in order to eliminate the e ff ect of fracture orientation. The spatial convergenceof the EDFM model for the case of a slipping fracture on non-conforming mesh is summarized in Fig. 7.Fig. 7a-e show the slip along the fracture normalized by the maximum slip given by Eq. 31 as a functionof a normalized coordinate along the fracture for particular fracture orientations. Fig. 7a-e indicate that thefracture slip given by EDFM on fine grids ( h = h = α . Fig. 7f shows the L -error in the numerical solutionas a function of the grid element size at various fracture strike angles. The error in the numerical solutiondecreases monotonically while refining the mesh for each fracture orientation. The average polynomial fit ofthe data in Fig. 7f yields e ∝ h . , which indicates super-linear convergence. .0 0.5 1.0Normalized Coordinate0.00.51.01.52.02.5 N o r m a li z e d J u m p (a) = 5 °h=2.500 mh=1.250 mh=0.625 manalytical 0.0 0.5 1.0Normalized Coordinate0.00.51.01.52.02.5 N o r m a li z e d J u m p (b) = 10 °h=2.500 mh=1.250 mh=0.625 manalytical 0.0 0.5 1.0Normalized Coordinate0.00.51.01.52.02.5 N o r m a li z e d J u m p (c) = 20 °h=2.500 mh=1.250 mh=0.625 manalytical0.0 0.5 1.0Normalized Coordinate0.00.51.01.52.02.5 N o r m a li z e d J u m p (d) = 30 °h=2.500 mh=1.250 mh=0.625 manalytical 0.0 0.5 1.0Normalized Coordinate0.00.51.01.52.02.5 N o r m a li z e d J u m p (e) = 40 °h=2.500 mh=1.250 mh=0.625 manalytical 0.6 1 1.5 2 2.5Element size h [m]10 L - E rr o r L i n e a r t r e n d Q u a d r a t i c t r e n d avg. fit1.491 (f) Spatial Convergence Figure 7: Spatial convergence of the EDFM mode on non-conforming mesh for the fracture slip problem.(a-e) Normalized fracture slip as a function of normalized distance along fracture at various refinement levelfor various fracture strike angles α . (f) L -error in tangential jump as a function of mesh element size forvarious fracture strike angles α .We now present the results of a convergence study for the case of an opening fracture on non-conforminggrids. Fig. 8a-e show the normalized fracture aperture as a function of the normalized coordinate for variousgrid element sizes and orientations. Similarly to the case of a slipping fracture, the EDFM model capturespoorly the shape of the aperture profile on coarse grids (approximately 3-5 elements per fracture). Upongrid refinement, however, the numerical results converge to the analytical solution for each case of gridorientation.he spatial convergence of the EDFM model on non-conforming grids for the Sneddon’s problem is shownin Fig. 8f. The spatial convergence is monotonic, nonlinear, and di ff ers for various orientation cases. Onaverage, however, the error decreases super-linearly with a decrease in the mesh element size as estimatedfrom the polynomial fit of the data in Fig. 8f e ∝ h . .Overall, the proposed EDFM mechanical model is accurate and converges to the analytical solution. Weobserve super-linear convergence behavior in all simulation cases on both conforming and non-conforminggrids. The DFM model also converges super-linearly and is more accurate than the EDFM model on con-forming grids. N o r m a li z e d A p e r t u r e (a) = 5 °h=2.500 mh=1.250 mh=0.625 manalytical 0.00 0.25 0.50 0.75 1.00Normalized Coordinate0.00.20.40.60.81.01.21.41.6 N o r m a li z e d A p e r t u r e (b) = 10 °h=2.500 mh=1.250 mh=0.625 manalytical 0.00 0.25 0.50 0.75 1.00Normalized Coordinate0.00.20.40.60.81.01.21.41.6 N o r m a li z e d A p e r t u r e (c) = 20 °h=2.500 mh=1.250 mh=0.625 manalytical0.00 0.25 0.50 0.75 1.00Normalized Coordinate0.00.20.40.60.81.01.21.41.6 N o r m a li z e d A p e r t u r e (d) = 30 °h=2.500 mh=1.250 mh=0.625 manalytical 0.00 0.25 0.50 0.75 1.00Normalized Coordinate0.00.20.40.60.81.01.21.41.6 N o r m a li z e d A p e r t u r e (e) = 40 °h=2.500 mh=1.250 mh=0.625 manalytical 0.6 1 2 3 4 5Element size h [m]10 L - E rr o r Q u a d r a t i c t r e n d L i n e a r t r e n d avg. fit1.511 (f) Spatial Convergence ° ° ° ° ° Figure 8: Spatial convergence of the EDFM mode on non-conforming mesh for the fracture opening prob-lem. (a) Normalized fracture aperture as a function of normalized distance along fracture. (b) L -error infracture aperture as a function of mesh element size. .2 Fluid flow test In this section, we validate the proposed approach for the flow problem. The 2D-domain geometry is thesame as in the previous tests (200 × ×
10 m ). A single fracture striking at 140 ◦ to the x-axis is locatedin the domain center. The flow parameters of the problem are listed in Table 2. Four wellbores are placedin the corners of the domain as shown in Fig. 9. The wellbore in the left upper corner of the domain isinjecting fluid at constant bottom-hole pressure 15 MPa. The rest three wells are production wells with fixedbottom-hole pressure 5 MPa. The initial reservoir pressure is 10 MPa. We apply no-flow condition on allboundaries. o
200 m m
160 m
Figure 9: Domain schematics for the flow test problem. The red circle indicates an injection well. The threeblue circles indicates production wells. The red straight line shows the fracture.To validate the EDFM flow model and investigate its spatial convergence, we compare the pressure in thefracture with that obtained with the DFM model on a fine grid ( h = h = ] 200 × × ◦ ] 140Fracture conductivity [mD · m] 20Matrix permeability [mD] 10Rock porosity [-] 0.2Initial reservoir pressure [MPa] 10Table 2: Geometrical and material parameters used in the flow test problem.The simulation results for the flow test are shown in Fig. 10. Fig. 10a shows the fracture pressure as afunction of the coordinate within the fracture at 10, 20, 40 days, and the steady-state solution. The fracturepressure is normalized by the initial reservoir pressure, and the coordinate is normalized by the fracturelength l . The pressure in the fracture decreases with time due to the higher number of producing wells thaninjection wells. The results in Fig. 10a were obtained from the DFM and EDFM models on the grids withthe characteristic element size h = .0 0.2 0.4 0.6 0.8 1.0Normalized Coordinate0.70.80.91.01.1 N o r m a li z e d p r e ss u r e
10 days20 days40 dayssteady state (a) Profiles at various time steps
DFMEDFM 0.0 0.2 0.4 0.6 0.8 1.0Normalized Coordinate0.700.750.800.850.90 N o r m a li z e d p r e ss u r e (b) Steady-state EDFM profiles referenceh = 4.0h = 2.0h = 1.0h = 0.50.0 0.2 0.4 0.6 0.8 1.0Normalized Coordinate0.700.750.800.850.90 N o r m a li z e d p r e ss u r e (c) Steady-state DFM profiles referenceh = 4.0h = 2.0h = 1.0h = 0.5 0.6 1 2 3 4Element size [m]12345 L - E rr o r [ ] (d) Spatial convergence EDFMDFMLinear trend
Figure 10: Simulation results for the flow problem. (a) Normalized fracture pressure as a function ofcoordinate along the fracture at various time instants given by DFM and EDFM models. (b)-(c) Steady-state fracture pressure given by the EDFM (b) and DFM (c) models at various grid resolutions. (d) Spatialconvergence of the steady-state solutions: both DFM and EDFM models exhibit linear convergence.
In this section we present a comparison of the DFM and EDFM models in a coupled 3D case. Fig. 11 showsa reservoir with nine variously-oriented fractures and five wellbores. The domain dimensions are the sameas in the previous examples. The initial reservoir pressure is 10 MPa. Four wells located in the four cornersof the domain produce fluid with the bottom-hole pressure 10 MPa. Another well at the center of the domaininjects fluid with constant rate 50 m / day and is connected to the longest fracture. The fracture geometricalproperties are listed in Table 3. /dayPres=10 MPa Figure 11: Domain geometry for a coupled problem. The domain comprises nine fractures and five wells.The producing wells are shown in blue. The injection well is connected to the long fracture in the center ofthe domain and is shown in red.The reservoir is subjected to anisotropic stresses as follows: the maximum total horizontal stress S Hmax = S x =
30 MPa, the minimum total horizontal stress S hmin = S y =
24 MPa, and the total vertical stress S z = =
0) results in various tangent traction t in di ff erent fractures. Fractures S Hmax . Therefore,they have the lowest initial tangent traction t τ ≈
10 MPa. Fractures ◦ to the direction of the maximum horizontal stress S Hmax , which causes a medium initial tangenttraction t τ ≈
26 MPa. Fractures ◦ , which results in the initial tangent traction to t τ ≈ t τ ≈
80 MPa because of the non-vertical orientation (dip = ◦ ) and the strike aligned closely to themaximum horizontal stress S Hmax . All the fractures have constant conductivity 20 mD · m.We applied the presented DFM and EDFM models to simulate 60 days of injection. The computationaldomain for the DFM simulation was discretized with an unstructured grid with 59032 tetrahedrons. The gridis fine at the fracture faces (2.2 m) and coarse at the outer domain boundaries (10 m). The mesh has fourells in the vertical direction near the fractures. The EDFM simulations utilized a Cartesian grid with 40804hexahedrons. The element size in the horizontal plane is 2 m. Same as the DFM grid, the EDFM grid hasfour cells in the vertical (z) direction.Fracture ◦ ] Dip [ ◦ ]1 120 30 902 and 3 80 30 904 and 5 60 80 906 and 7 60 80 808 and 9 60 30 80Table 3: Fracture geometrical parameters used in the coupled problem. S z = 70 MPaS y = 24 MPa v e r t i c a l d i p o v e r t i c a l v e r t i c a l o T a n g e n t T r a c t i o n [ M P a ] o S x = 30 MPa v e r t i c a l xyz v e r t i c a l d i p o d i p o Figure 12: Initial stresses and tangential tractions on fracture surfaces. Vertical fractures t t ≈
10 MPa. Vertical fractures t t ≈
25 MPa. Fractures ◦ and have a high initial tangent traction t t ≈
70 MPa.Critically-oriented fractures ◦ and have the highest tangent traction t t ≈
80 MPa.The process of induced fracture reactivation is illustrated in Fig. 13. At approximately 15 days of injectionhe longest Fracture ≈
30% area) of Fracture iscrete Fractures Embedded Fractures(a)(b)(c) Stick 27 days31 days33 days1 24 58 6 79 xyz
Figure 13: Development of the fracture states predicted by the DFM (left column) and EDFM (right column)methods. The figure shows which fracture elements are inactive (stick) or are active and exhibit slip oropening at various time steps. (a) At 27 days of injection the longest fracture a) 34 days OpeningSlipStick1 2 1 24 58 6 934 x yz Figure 14: Development of the fracture states predicted by the DFM (left column) and EDFM (right column)methods. (a) Fracture ≈ ff erence in the tangential jump given by the DFM and EDFMmodels are due to slight di ff erences in fluid pressure and because the EDFM model overestimates fractureslip, whereas the DFM model underestimates it as exemplified in previous sections.
40 50 60 70 80Coordinate [m]012345678 T a n g e n t j u m p [ mm ]
27 days34 days40 days (a) Fracture
DFMEDFM 0 20 40 60Coordinate [m]012345678 T a n g e n t j u m p [ mm ]
34 days40 days (b) Fracture T a n g e n t j u m p [ mm ]
34 days40 days (c) Fracture
Figure 15: Fracture slip values as a functions of the coordinate at various time steps in (a) Fracture
Discussion
This section evaluates the performance of the proposed EDFM model and discusses its limitations. We alsoprovide an insight into the modeling cases where using EDFM or DFM is more preferable.
All the examples with slipping fractures shown in Section 4 indicate that EDFM overestimates the slip.Moreover, the results evidence that the tangential jumps given by EDFM are always higher than those com-puted with the DFM model. This is presumably due to the locking e ff ect in the proximity of the fracturetip caused by the continuous slip interpolation in the DFM model [Borja, 2013]. The near-tip solution canpotentially be improved by assuming non-constant jump within an element [Linder and Armero, 2007]. Ad-ditionally, using EDFM on non-conforming grids results in non-smooth displacement profiles, which, in turn,may cause non-physical stress and traction values in active fractures. This trait of the numerical model maybe undesirable in stress-sensitive reservoirs, especially if fracture permeability is computed as a function ofstress [Rutqvist, 2015, Shovkun and Espinoza, 2017]. In addition, this also may cause convergence issuesin fully-coupled models. Therefore, we suggest to use the DFM approach in cases when a smooth jump isrequired.The issue with non-realistic stresses in the SDA models can be circumvented by enforcing the continu-ity of the jump [Deb and Jenny, 2017]. Requiring it, however, results in a non-locality in the numericaltreatment of the jump. An alternative way to circumvent the deficiencies associated with computing themechanics-dependent fracture permeability with the presented SDA model is to assign permeability as afunction of fracture aperture [Zimmerman et al., 1991, Shovkun and Espinoza, 2019]. Since fracture aper-ture is widely used to calculate the permeability of open-mode fractures with the lubrication theory [Leeet al., 2017, Shovkun and Espinoza, 2019] and given the capability to handle arbitrary fracture geometries,SDA is potentially suitable for modeling open-mode fracture propagation. A hybrid formulation that makesuse of both SDA and DFM methods for various types of fractures can potentially be an optimal approach forcoupled hydro-mechanical modelling of fractured reservoirs. .2 SDA formulation for intersecting fractures In this section we outline the basic idea behind the mechanical treatment of fracture intersection in our SDAmodel. When more than one embedded fracture are located in a grid element, the contact mechanics isdescribed with the following system: ˙ σ = (cid:131) : ∇ s ˙ u c − N F (cid:88) i = ([ ˙ u ] i ⊗ ∇ f i ) s (33a)[ ˙ u ] i = λ δ i ∂ G i ∂ t i , i = ... N F (33b)˙ q i = − λ δ i H δ i ∂ G i ∂ q i , i = ... N F (33c) λ δ i F i ( t i , q i ) = , i = ... N F (33d)where N F is the number of embedded fractures in the cell, [ u ] i are the fracture jump vectors, λ i are fractureplastic multipliers, q i are the fracture internal state variables, H δ i are the fracture hardening moduli, t i is thefracture traction, and F i and G i are the flow rules plastic potentials, respectively. The system (33) is thedirect analogue to (15). All other equations and the solution procedure remain the same as in the case of asingle fracture. Therefore, adding a fracture to the system consists of (1) imposing independent constraints(33)b-c on variables related to a particular fracture and (2) using the superposition of jumps to determine thestress (Eq. 33a). Thus, coupling several EDFM fractures within a cell occurs via the e ff ective stress σ .This strategy has two drawbacks. First, formulation (33) permits the maximum of two embedded frac-tures per grid element. Applying this approach to more than two fractures per element results in an over-determined system.Second, each combination of fracture states (e.g. the first fracture is slipping and the second fracture isopening) requires a special treatment. This treatment consists in relaxing the constraints (16)- (20) to avoidthe over-determination of the system (33). Future work is needed to develop a robust procedure for treatingfracture intersection with various statuses. The approach to modeling embedded fractures discussed in Section 3.1, is a simplified way of describingfluid flow. While providing a relatively accurate approximation for single-phase fluids for highly-permeableractures, the traditional EDFM approach has two important limitations. First, standard EDFM cannot sim-ulate fractures with permeability lower than that of the rock matrix [Flemisch et al., 2018]. Second, thestandard EDFM cannot accurately capture cross-fracture flow, sometimes referred to as fracture sweeping[Karvounis, 2013]. This e ff ect becomes particularly significant when solving transport equations. Two so-lutions have been suggested to remove these limitations. Karvounis in proposed to modify the computationof the matrix-fracture flow q MF terms based on the direction of the cell fluxes [Karvounis, 2013]. ProjectionEDFM (pEDFM) is another approach that modifies the number of connections in an approximation pattern[T¸ ene et al., 2017]. The latter approach can be easily implemented in the connection-based flow simulatoras a part of the preprocessing stage [Karimi-Fard et al., 2004], In this paper we presented a numerical model that uses the Strong Discontinuity Approach (SDA) for ge-omechanics and the Embedded Fracture (EDFM) approach for fluid flow. This study validates the use of acoupled EDFM model in relatively complex geological settings.A series of mechanical tests was performed to compare the performance of the Discrete Fracture Model(DFM) and the proposed EDFM model. Both DFM and EDFM manifest asymptotic super-linear conver-gence on conforming grids when modeling slipping and opening fractures. On non-conforming grids, theEDFM model manifests a range of convergence behaviors depending on the fracture orientation with respectto the computational grids. On average, however, the EDFM model retains the super-linear convergencetrend for both opening and slipping fractures.Numerical simulation evidences discontinuous slip and aperture within EDFM fractures. On very coarsegrids, the absence of jump continuity may result in non-physical jump profiles. Based on the simulations weconclude that EDFM overestimates fracture slip, whereas DFM underestimates it.We presented a simplified EDFM approach to simulate single-phase fluid flow in fractured porous media.The presented model provides accurate results that are very close to those given by the DFM model. Wedemonstrated that both methods have similar linear convergence behavior.Finally, we applied the DFM and EDFM methods to simulate a complex 3D coupled hydro-mechanicalproblem. We considered a fluid injection scenario that resulted in natural fractures activation induced byxcessive pressure and stress alterations.Our results indicate that EDFM and DFM methods provide qualitatively similar results. Based on multiplenumerical tests, we conclude that EDFM can be considered an alternative approach to model highly fracturedreservoirs.
Acknowledgments
We thank the Reservoir Simulation Industrial A ffi liates Consortium at Stanford University (SUPRI-B) forfunding this research. References [Alfaiate et al., 2010] Alfaiate, J., Moonen, P., Sluys, L., and Carmeliet, J. (2010). On the use of strongdiscontinuity formulations for the modeling of preferential moisture uptake in fractured porous media.Computer methods in applied mechanics and engineering, 199(45-48):2828–2839.[Borja, 2000] Borja, R. I. (2000). A finite element model for strain localization analysis of strongly discon-tinuous fields based on standard galerkin approximation. Computer Methods in Applied Mechanics andEngineering, 190(11-12):1529–1549.[Borja, 2013] Borja, R. I. (2013). Plasticity, volume 2. Springer.[Cacas et al., 1990] Cacas, M.-C., Ledoux, E., Marsily, G. d., Tillie, B., Barbreau, A., Durand, E., Feuga, B.,and Peaudecerf, P. (1990). Modeling fracture flow with a stochastic discrete fracture network: calibrationand validation: 1. the flow model. Water Resources Research, 26(3):479–489.[Callari and Armero, 2002] Callari, C. and Armero, F. (2002). Finite element methods for the analysisof strong discontinuities in coupled poro-plastic media. Computer Methods in Applied Mechanics andEngineering, 191(39-40):4371–4400.[Callari et al., 2010] Callari, C., Armero, F., and Abati, A. (2010). Strong discontinuities in partially satu-rated poroplastic solids. Computer Methods in Applied Mechanics and Engineering, 199(23-24):1513–1535.Cappa and Rutqvist, 2011] Cappa, F. and Rutqvist, J. (2011). Modeling of coupled deformation and per-meability evolution during fault reactivation induced by deep underground injection of co2. InternationalJournal of Greenhouse Gas Control, 5(2):336–346.[Coussy, 2004] Coussy, O. (2004). Poromechanics. Wiley.[Coussy, 2010] Coussy, O. (2010). Mechanics and Physics of Porous Solids. Wiley.[Deb and Jenny, 2017] Deb, R. and Jenny, P. (2017). Modeling of shear failure in fractured reservoirs witha porous matrix. Computational Geosciences, 21(5-6):1119–1134.[Duarte et al., 2000] Duarte, C. A., Babuˇska, I., and Oden, J. T. (2000). Generalized finite element methodsfor three-dimensional structural mechanics problems. Computers & Structures, 77(2):215–232.[Elf-Aquitaine, 1992] Elf-Aquitaine, F. (1992). Monitoring of subsidence and induced seismicity in thel-ael gas field (france): the consequences on gas production and field operation. Engineering Geology,32:123–135.[Flemisch et al., 2018] Flemisch, B., Berre, I., Boon, W., Fumagalli, A., Schwenck, N., Scotti, A., Stefans-son, I., and Tatomir, A. (2018). Benchmarks for single-phase flow in fractured porous media. Advancesin Water Resources, 111:239–258.[Foster et al., 2007] Foster, C., Borja, R., and Regueiro, R. (2007). Embedded strong discontinuity finiteelements for fractured geomaterials with variable friction. International Journal for Numerical Methodsin Engineering, 72(5):549–581.[Fu et al., 2016] Fu, P., Hao, Y., Walsh, S. D. C., and Carrigan, C. R. (2016). Thermal drawdown-inducedflow channeling in fractured geothermal reservoirs. Rock Mechanics and Rock Engineering, 49(3):1001–1024.[Garipov et al., 2016] Garipov, T., Karimi-Fard, M., and Tchelepi, H. (2016). Discrete fracture model forcoupled flow and geomechanics. Computational Geosciences, 20(1):149–160.[Garipov et al., 2018] Garipov, T., Tomin, P., Rin, R., Voskov, D., and Tchelepi, H. (2018). Unified thermo-compositional-mechanical framework for reservoir simulation. Computational Geosciences, 22(4):1039–1057.Gong et al., 2008] Gong, B., Karimi-Fard, M., Durlofsky, L. J., et al. (2008). Upscaling discrete fracturecharacterizations to dual-porosity, dual-permeability models for e ffi cient simulation of flow with stronggravitational e ff ects. SPE Journal, 13(01):58–67.[Grasso, 1992] Grasso, J.-R. (1992). Mechanics of seismic instabilities induced by the recovery of hydro-carbons. Pure and Applied Geophysics, 139(3-4):507–534.[Haddad and Sepehrnoori, 2015] Haddad, M. and Sepehrnoori, K. (2015). Integration of xfem and czmto model 3d multiple-stage hydraulic fracturing in quasi-brittle shale formations: Solution-dependentpropagation direction. In Proceedings of the AADE National Technical Conference and Exhibition,AADE2015, San Antonio, Texas, 8-9 April 2015.[Hansbo and Hansbo, 2002] Hansbo, A. and Hansbo, P. (2002). An unfitted finite element method, basedon nitsches method, for elliptic interface problems. Computer methods in applied mechanics andengineering, 191(47-48):5537–5552.[Heister et al., 2015] Heister, T., Wheeler, M. F., and Wick, T. (2015). A primal-dual active set methodand predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach.Computer Methods in Applied Mechanics and Engineering, 290:466–495.[Hwang et al., 2015] Hwang, J., Bryant, E. C., Sharma, M. M., et al. (2015). Stress reorientation in water-flooded reservoirs. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.[Karimi-Fard et al., 2004] Karimi-Fard, M., Durlofsky, L., and Aziz, K. (2004). An e ffi cient discrete-fracture model applicable for general-purpose reservoir simulators. SPE Journal, 9(2):227–236.[Karimi-Fard et al., 2003] Karimi-Fard, M., Durlofsky, L. J., Aziz, K., et al. (2003). An e ffi cient dis-crete fracture model applicable for general purpose reservoir simulators. In SPE Reservoir SimulationSymposium. Society of Petroleum Engineers.[Karvounis, 2013] Karvounis, D. C. (2013). Simulations of enhanced geothermal systems with an adaptivehierarchical fracture representation. PhD thesis, ETH Zurich.[Lee et al., 2017] Lee, S., Wheeler, M. F., and Wick, T. (2017). Iterative coupling of flow, geomechanicsand adaptive phase-field fracture including level-set crack width approaches. Journal of Computationaland Applied Mathematics, 314:40–60.Li et al., 2015] Li, J., Lei, Z., Qin, G., and Gong, B. (2015). E ff ective local-global upscaling of fracturedreservoirs under discrete fractured discretization. Energies, 8(9):10178–10197.[Li et al., 2008] Li, L., Lee, S. H., et al. (2008). E ffi cient field-scale simulation of black oil in a naturallyfractured reservoir through discrete fracture networks and homogenized media. SPE Reservoir Evaluation& Engineering, 11(04):750–758.[Linder and Armero, 2007] Linder, C. and Armero, F. (2007). Finite elements with embedded strong discon-tinuities for the modeling of failure in solids. International Journal for Numerical Methods in Engineering,72(12):1391–1433.[Liu and Borja, 2010] Liu, F. and Borja, R. I. (2010). Stabilized low-order finite elements for frictional con-tact with the extended finite element method. Computer Methods in Applied Mechanics and Engineering,199(37-40):2456–2471.[McClure, 2012] McClure, M. W. (2012). Modeling and characterization of hydraulic stimulation andinduced seismicity in geothermal and shale gas reservoirs. PhD thesis, Stanford University Stanford,California.[Min et al., 2004] Min, K.-B., Rutqvist, J., Tsang, C.-F., and Jing, L. (2004). Stress-dependent permeabil-ity of fractured rock masses: a numerical study. International Journal of Rock Mechanics and MiningSciences, 41(7):1191–1210.[Mosler, 2005] Mosler, J. (2005). A novel algorithmic framework for the numerical implementation oflocally embedded strong discontinuities. Computer Methods in Applied Mechanics and Engineering,194(45-47):4731–4757.[Mosler and Meschke, 2003] Mosler, J. and Meschke, G. (2003). 3d modelling of strong discontinuitiesin elastoplastic solids: fixed and rotating localization formulations. International Journal for NumericalMethods in Engineering, 57(11):1553–1576.[Narasimhan and Pruess, 1988] Narasimhan, T. and Pruess, K. (1988). Minc: An approach for analyzingtransport in strongly heterogeneous systems. In Groundwater Flow and Quality Modelling, pages 375–391. Springer.Norbeck et al., 2016] Norbeck, J. H., McClure, M. W., Lo, J. W., and Horne, R. N. (2016). An embeddedfracture modeling framework for simulation of hydraulic fracturing and shear stimulation. ComputationalGeosciences, 20(1):1–18.[Oliver, 1996] Oliver, J. (1996). Modelling strong discontinuities in solid mechanics via strain softeningconstitutive equations. part 1: Fundamentals. International journal for numerical methods in engineering,39(21):3575–3600.[Oliver et al., 1999] Oliver, J., Cervera, M., and Manzoli, O. (1999). Strong discontinuities and continuumplasticity models: the strong discontinuity approach. International journal of plasticity, 15(3):319–351.[Oliver et al., 2003] Oliver, J., Huespe, A., and Samaniego, E. (2003). A study on finite elements for captur-ing strong discontinuities. International journal for numerical methods in engineering, 56(14):2135–2161.[Olson, 2008] Olson, J. E. (2008). Multi-fracture propagation modeling: Applications to hydraulic fractur-ing in shales and tight gas sands. In 42nd US Rock Mechanics Symposium and 2nd US-Canada RockMechanics Symposium, ARMA 08-237.[Ouchi et al., 2015] Ouchi, H., Katiyar, A., York, J., Foster, J. T., and Sharma, M. M. (2015). A fully coupledporous flow and geomechanics model for fluid driven cracks: a peridynamics approach. ComputationalMechanics, 55(3):561–576.[Phan et al., 2003] Phan, A.-V., Napier, J., Gray, L., and Kaplan, T. (2003). Symmetric-galerkin bem sim-ulation of fracture with frictional contact. International journal for numerical methods in engineering,57(6):835–851.[Pruess, 1992] Pruess, K. (1992). Brief guide to the minc-method for modeling flow and transport in frac-tured media. Technical report, Lawrence Berkeley Lab., CA (United States).[Regueiro and Borja, 1999] Regueiro, R. A. and Borja, R. I. (1999). A finite element model of localizeddeformation in frictional materials taking a strong discontinuity approach. Finite Elements in Analysisand Design, 33(4):283–315.[Regueiro and Borja, 2001] Regueiro, R. A. and Borja, R. I. (2001). Plane strain finite element analysisof pressure sensitive plasticity with strong discontinuity. International Journal of Solids and Structures,38(21):3647–3672.Rutledge et al., 2004] Rutledge, J., Phillips, W., and Mayerhofer, M. (2004). Faulting induced by forcedfluid injection and fluid flow forced by faulting: An interpretation of hydraulic-fracture microseismicity,carthage cotton valley gas field, texas. Bulletin of the Seismological Society of America, 94(5):1817–1830.[Rutqvist, 2015] Rutqvist, J. (2015). Fractured rock stress-permeability relationships from in situ data ande ff ects of temperature and chemical-mechanical couplings. Geofluids, 15(1-2):48–66.[Salimzadeh et al., 2018] Salimzadeh, S., Paluszny, A., Nick, H. M., and Zimmerman, R. W. (2018). Athree-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal sys-tems. Geothermics, 71:212 – 224.[Sesetty et al., 2012] Sesetty, V., Ghassemi, A., et al. (2012). Simulation of hydraulic fractures and theirinteractions with natural fractures. In 46th US Rock Mechanics / Geomechanics Symposium. AmericanRock Mechanics Association.[Shovkun and Espinoza, 2017] Shovkun, I. and Espinoza, D. N. (2017). Coupled fluid flow-geomechanicssimulation in stress-sensitive coal and shale reservoirs: Impact of desorption-induced stresses, shear fail-ure, and fines migration. Fuel, 195:260–272.[Shovkun and Espinoza, 2018] Shovkun, I. and Espinoza, D. N. (2018). Geomechanical implications ofdissolution of mineralized natural fractures in shale formations. Journal of Petroleum Science andEngineering, 160:555–564.[Shovkun and Espinoza, 2019] Shovkun, I. and Espinoza, D. N. (2019). Propagation of toughness-dominated fluid-driven fractures in reactive porous media. International Journal of Rock Mechanics andMining Sciences.[Simo and Oliver, 1994] Simo, J. and Oliver, J. (1994). A new approach to the analysis and simulation ofstrain softening in solids. Fracture and damage in quasibrittle structures, pages 25–39.[Simo et al., 1993] Simo, J. C., Oliver, J., and Armero, F. (1993). An analysis of strong discontinuitiesinduced by strain-softening in rate-independent inelastic solids. Computational mechanics, 12(5):277–296.Simo and Rifai, 1990] Simo, J. C. and Rifai, M. (1990). A class of mixed assumed strain methods and themethod of incompatible modes. International journal for numerical methods in engineering, 29(8):1595–1638.[Sneddon and Elliot, 1946] Sneddon, I. and Elliot, H. (1946). The opening of a gri ffiffi