# A New Projection-Based Integrally Embedded Discrete Fracture Model and Its Application in Coupled Flow and Geomechanical Simulation for Fractured Reservoirs

Authors: Dawei Wu; Yuan Di; Yu-Shu Wu
Journal: Geofluids
 (2022)

---

## Abstract

The embedded discrete fracture model (EDFM) has been popular for the modeling of fractured reservoirs due to its flexibility and efficiency while maintaining the complex geometry of fracture networks. Though the EDFM has been validated for single-phase flow simulations, some recent cases show that the EDFM might result in apparent errors in multiphase flow situations. The projection-based embedded discrete fracture model (pEDFM) and the integrally embedded discrete fracture model (IEDFM) are two recently developed methods, which intend to improve the accuracy of the EDFM. In this study, a projection-based integrally embedded discrete fracture model (pIEDFM) is proposed, which combines the advantages of the pEDFM and the IEDFM. Similar to the pEDFM, the pIEDFM uses a kind of additional connections between fracture and nonneighboring matrix cells to obtain more physically authentic velocity fields. As a significant improvement, a semi-analytical cone-shaped pressure distribution that follows the IEDFM is adopted in the pIEDFM to capture the sharp pressure change near the fracture surfaces. Comparisons with benchmark results and explicit-fracture fine grid simulation results show that the pIEDFM provides accurate solutions using a moderate amount of grids. The proposed pIEDFM is also applied to coupled flow and geomechanical simulation for fractured reservoirs. Comparison of our coupled simulation results with that of the EDFM shows that the pIEDFM is applicable for the coupled simulation, and the different methods for matrix-fracture transmissibility between the pIEDFM and the EDFM may lead to deviations in stress fields predicted by geomechanical modeling, which eventually affects the oil production, water cut, and oil saturation profiles.


---

## Body

## 1. Introduction

Numerical simulation approaches for fractured reservoirs have drawn great attention in past decades. Due to the significant permeability difference between the less permeable rock matrix and the highly conductive fractures, modeling multiphase flow in fractured media has become challenging. Barenblatt et al. [1] first proposed a concept of dual porosity to describe the seepage process in fractured rocks. Warren and Root [2] extended the concept of dual porosity and developed the dual-porosity model. Kazemi et al. [3, 4] introduced the dual-porosity model into petroleum engineering and applied the method in the modeling of fractured reservoirs. Later, a series of numerical approaches were developed as extensions of the dual-porosity model, including the dual-porosity dual-permeability model [5, 6], the triple-porosity dual-permeability model [7], and the multiple interacting continuum model [8]. All these models can be classified as dual-continuum models. Dual-continuum models provide an efficient way of simulating fractured reservoirs. However, the geometries of the fractures are lost in the assumption of dual-porosity models, and the real fracture network cannot be represented in dual-continuum models [9]. Panfili et al. [10] showed that the homogenization used in dual-continuum models could introduce unphysical fracture flows between disconnected areas. Moinfar et al. [11] investigated the examples of reservoirs with complex fracture networks where dual-continuum models cannot provide precise solutions.A discrete fracture model (DFM) [12–15] was developed to solve the limitations of dual-continuum models. In the discrete fracture model, every fracture is treated explicitly and individually, providing a more physical and realistic representation of fractures, especially with complex geometries [16, 17]. However, to adapt to the complex fracture networks, unstructured matrix grids are commonly used in the DFM, which causes difficulties in griding. In fracture networks composed of microscale fractures, the generated unstructured grids are usually in large numbers, which causes high computational cost, making the DFM impractical in actual field studies [18].The embedded discrete fracture model (EDFM) was proposed by Li and Lee [19] as an alternative to the DFM and dual-continuum models. In the EDFM, the Cartesian grids are used in matrix discretization to keep the efficiency of the method. The fractures are discretized explicitly as control volumes, also known as fracture elements. The fracture elements are embedded into their parent matrix grids through matrix-fracture connections. The EDFM has been successfully implemented in vertical fracture cases [19], nonvertical fracture cases [18], and nonplanar fracture cases [20, 21]. Some authors [22, 23] combined the EDFM with dual-continuum methods to model reservoirs with multiscale fractures. Moinfar et al. [24] applied the EDFM in coupled flow and geomechanical simulations of fractured reservoirs.Though the EDFM has been validated in various studies of its accuracy for single-phase flow simulations, it may result in apparent errors in some cases for multiphase simulations. Figure1 shows a matrix block M intersected by a fracture. The two matrix fractions of M separated by the fracture are marked as A and B. As illustrated in Figure 1(a), in the realistic situation, an incoming water flow that moves towards a fracture in the water flooding process first enters fraction B and then is split into two parts that move along (F1) and across (F2) the fracture, respectively. However, in the EDFM, the two fractions A and B are combined as an intact matrix block instead of being considered as two individual grids. Thus, the water flow across the fracture from fraction A to fraction B cannot be exhibited. Instead, the averaged flow between the fracture and the matrix block M is evaluated (Fm-f), as shown in Figure 1(b). The water from fractions A and B can flow towards the fracture simultaneously, which increases the water flux that moves along the fracture and results in an “unphysical flux split” in the EDFM. Therefore, the water flux along the fracture is overestimated, and the water flux across the fracture is underestimated [25].Figure 1
Flow process for water displacement across a fracture in (a) realistic scenario and (b) EDFM.
(a)(b)To solve this problem, a projection-based embedded discrete fracture model (pEDFM) is proposed by Tene and others. [26]. Jiang and Younis proved [25] that a more physical flux split could be achieved in the pEDFM, thus fixing the erroneous water displacement process predicted by the EDFM. Olorode et al. [27] extended the pEDFM to three-dimensional cases and investigated 3D compositional modeling with the pEDFM. Rao et al. [28] modified the original pEDFM and developed a micro-translation algorithm to help select projection-face combinations.Another limitation of the EDFM is the oversimplified assumption for pressure distribution in the matrix domain adjacent to fracture. In the EDFM, an approximate linear pressure distribution assumption is used around the fractures. However, a cone-shaped distribution of pressure is usually generated due to the large difference in permeability between matrix and fracture, as shown in Figure2. The oversimplified assumption in the EDFM may lead to errors in calculating the transmissibilities of matrix-fracture connections [29]. An integrally embedded discrete fracture model (IEDFM) [30] has been proposed to improve the transmissibility calculation of the EDFM. In the IEDFM, the transmissibilities of matrix-fracture connections are derived semi-analytically, which obtains the more realistic pressure distribution near the fracture surfaces and improves the accuracy of modeling flow in fractured reservoirs. The IEDFM is later extended to the modeling of anisotropic fractured reservoirs [31].Figure 2
Illustration of linear and cone-shaped pressure distributions in fractured media.A projection-based integrally embedded discrete fracture model (pIEDFM) is proposed in this study, which combines the advantages of the pEDFM and the IEDFM. Similar to the pEDFM, additional matrix-fracture connections are added in the pIEDFM between a fracture element and the nonneighboring matrix elements along the fracture projection directions. The transmissibilities of neighboring and nonneighboring matrix-fracture connections in the pIEDFM are derived semi-analytically using the methods in the IEDFM. The accuracy of the pIEDFM is validated by benchmark results and fine grid simulation results. The proposed pIEDFM is then applied in coupled flow and geomechanical simulation for fractured reservoirs. The applicability of the proposed numerical method is examined.


## 2. Governing Equations for Fractured Reservoir Simulation

### 2.1. Mass Conservation Equations

The mass conservation equations in fractured media are given as follows:(1)∂∂tϕSβρβ+∇⋅ρβvβ=qβW+qβfm,where β represents fluid phases, ϕ is the porosity, Sβ is the saturation of phase β, ρβ is the density of phase β, vβ is the velocity of phase β, qβW is the flux term of phase β from wells, and qβfm is the flux term of phase β between matrix and fracture elements.The velocity of phaseβ is defined by Darcy’s law:(2)vβ=−kkrβμβ∇Pβ−ρβg∇Z=−kkrβμβ∇Ψβ,where k is the absolute permeability, krβ is the relative permeability of phase β, μβ is the viscosity of phase β, Pβ is the pressure of phase β, g is the gravitational acceleration, Z is the depth, and Ψβ is the flow potential of phase β.Equation (1) is discretized using the control-volume finite difference scheme in space and first-order scheme in time, which gives the following:(3)ϕSβρβin+1−ϕSβρβinViΔt=∑jFβ,ijn+1+∑kFβ,ikfm,n+1+Qβ,iW,n+1,where subscript i denotes the values of element i, superscript n+1 represents the current time level, superscript n represents the previous time level, Fβ,ij is the flow term for phase β between element i and element j where i and j are the same type of element (matrix or fracture), Fβ,ikfm is the flow term for phase β through the matrix-fracture connection of element i and element k, and Qβ,iW is the flux term of phase β from wells.The flow termsFβ,ij and Fβ,ikfm are expressed as follows:(4)Fβ,ij=Tijρβkrβμβij+1/2Ψβ,j−Ψβ,i,Fβ,ikfm=Tikρβkrβμβik+1/2Ψβ,k−Ψβ,i,where subscripts ij+1/2 and ik+1/2 denote proper averages of properties at the interface, and Tij and Tik are the transmissibilities of the connections.The transmissibility of a connection is the harmonic average of two half-transmissibilities:(5)T12=T1T2T1+T2,T1=A12k1d1,T2=A12k2d2.where T1 and T2 are the half-transmissibilities of element 1 and element 2, k1 and k2 are the absolute permeabilities, and d1 and d2 are the distances from the centers of elements to the interface.


### 2.2. EDFM

The reservoir matrix is discretized using the Cartesian grids, and additional fracture elements are added to represent the fracture control volumes. Thus, there are three kinds of connections between elements in the EDFM: matrix-matrix, fracture-fracture, and matrix-fracture connections. For matrix-matrix and fracture-fracture connections, the transmissibilities can be derived geometrically from the two-point flux approximation (TPFA) using (5).For matrix-fracture connection in 2D reservoir cases, the half-transmissibilities of matrix and fracture can be derived as follows:(6)Tm=2Afkmdfm,Tf=2Afkfdf,where Af is the fracture surface area, km is the matrix permeability, kf is the fracture permeability, df is the fracture center distance from the interface, which equals half of the fracture aperture, and dfm is the equivalent distance between matrix and fracture elements.Using the approximation of linear pressure distribution around fractures, the equivalent distance can be given as follows:(7)dfm=∫rifdV​V,where rif is the distance from fracture and V is the volume of the matrix.


## 2.1. Mass Conservation Equations

The mass conservation equations in fractured media are given as follows:(1)∂∂tϕSβρβ+∇⋅ρβvβ=qβW+qβfm,where β represents fluid phases, ϕ is the porosity, Sβ is the saturation of phase β, ρβ is the density of phase β, vβ is the velocity of phase β, qβW is the flux term of phase β from wells, and qβfm is the flux term of phase β between matrix and fracture elements.The velocity of phaseβ is defined by Darcy’s law:(2)vβ=−kkrβμβ∇Pβ−ρβg∇Z=−kkrβμβ∇Ψβ,where k is the absolute permeability, krβ is the relative permeability of phase β, μβ is the viscosity of phase β, Pβ is the pressure of phase β, g is the gravitational acceleration, Z is the depth, and Ψβ is the flow potential of phase β.Equation (1) is discretized using the control-volume finite difference scheme in space and first-order scheme in time, which gives the following:(3)ϕSβρβin+1−ϕSβρβinViΔt=∑jFβ,ijn+1+∑kFβ,ikfm,n+1+Qβ,iW,n+1,where subscript i denotes the values of element i, superscript n+1 represents the current time level, superscript n represents the previous time level, Fβ,ij is the flow term for phase β between element i and element j where i and j are the same type of element (matrix or fracture), Fβ,ikfm is the flow term for phase β through the matrix-fracture connection of element i and element k, and Qβ,iW is the flux term of phase β from wells.The flow termsFβ,ij and Fβ,ikfm are expressed as follows:(4)Fβ,ij=Tijρβkrβμβij+1/2Ψβ,j−Ψβ,i,Fβ,ikfm=Tikρβkrβμβik+1/2Ψβ,k−Ψβ,i,where subscripts ij+1/2 and ik+1/2 denote proper averages of properties at the interface, and Tij and Tik are the transmissibilities of the connections.The transmissibility of a connection is the harmonic average of two half-transmissibilities:(5)T12=T1T2T1+T2,T1=A12k1d1,T2=A12k2d2.where T1 and T2 are the half-transmissibilities of element 1 and element 2, k1 and k2 are the absolute permeabilities, and d1 and d2 are the distances from the centers of elements to the interface.


## 2.2. EDFM

The reservoir matrix is discretized using the Cartesian grids, and additional fracture elements are added to represent the fracture control volumes. Thus, there are three kinds of connections between elements in the EDFM: matrix-matrix, fracture-fracture, and matrix-fracture connections. For matrix-matrix and fracture-fracture connections, the transmissibilities can be derived geometrically from the two-point flux approximation (TPFA) using (5).For matrix-fracture connection in 2D reservoir cases, the half-transmissibilities of matrix and fracture can be derived as follows:(6)Tm=2Afkmdfm,Tf=2Afkfdf,where Af is the fracture surface area, km is the matrix permeability, kf is the fracture permeability, df is the fracture center distance from the interface, which equals half of the fracture aperture, and dfm is the equivalent distance between matrix and fracture elements.Using the approximation of linear pressure distribution around fractures, the equivalent distance can be given as follows:(7)dfm=∫rifdV​V,where rif is the distance from fracture and V is the volume of the matrix.


## 3. Projection-Based Integrally Embedded Discrete Fracture Model

In the proposed pIEDFM, the connection relationship establishment follows the rules of the pEDFM, where additional connections are introduced between the fracture element and the nonneighboring matrix element along the fracture projection directions. As shown in Figure3, the fracture element f has two projections on the boundary of matrix element i, which have the area of AmXP and AmYP, respectively. The fracture element f is connected to its neighboring matrix element i and two nonneighboring matrix elements j and k. The criterion of selecting the matrix faces of fracture projections follows the work of Jiang and Younis [25]. The matrix faces that are closer to the fracture center are selected as the projected faces in each dimension.Figure 3
Illustration of connections in the pEDFM.For matrix-matrix connections, the projected areas of fractures are eliminated from the interface area. The modified interface area between matrixi and matrix j is given as follows:(8)A¯ij=Aij−AmXP,where Aij is the original interface area between matrix i and matrix j, and AmXP is the fracture projection area along the x-direction. When a fracture cuts through the matrix elements, the matrix-matrix connection will be removed.Figure4 shows an example of the connection establishment in the pIEDFM and the EDFM. There are four matrix elements marked as M1, M2, M3, and M4. The two fractures are discretized into several fracture segment elements by the matrix block boundaries. The fracture segment elements represented by blue lines with red dots are marked as F1, F2, F3, F4, F5, and F6. In the EDFM, there are 4 matrix-matrix connections, 5 fracture-fracture connections, and 6 matrix-fracture connections. However, in the pIEDFM, 3 matrix-matrix connections, 5 fracture-fracture connections, and 14 matrix-fracture connections are included. The number of matrix-fracture connections increases in the pIEDFM due to the additional connections between the fracture segments and the nonneighboring matrix elements.Figure 4
Schematic for connections in the pIEDFM and the EDFM.In the pIEDFM, the matrix-matrix and fracture-fracture connection transmissibilities can be directly derived using (5). The calculation formulations of matrix-fracture connection transmissibilities are derived on the basis of the IEDFM, where the fractures are considered as series of point sinks and the transmissibilities are solved semi-analytically. The pressure distribution inside a matrix can be derived from the point sinks that form the fracture, thus reproducing the cone-shaped distribution of pressure to improve the simulation accuracy of fluid exchange between matrix and fracture. In the pIEDFM, the transmissibility between the fracture and the neighboring matrixes and the transmissibility between the fracture and the nonneighboring matrixes are calculated separately.Figure5 presents a schematic for the calculation of matrix-fracture transmissibilities in the pIEDFM. In the vicinity of a fracture, the pressure of point X can be derived as the superposition of all the pressure drops caused by point sinks:(9)PXx,y=∑i=1NSqSiμ2πhkmlnriX+CP,where h is the height of the strata, NS is the number of point sinks, qSi is the flow rate of sink Si, riX is the distance between point X and sink Si, and CP is a constant related to the fracture pressure.Figure 5
Schematic for the calculation of matrix-fracture transmissibility.When fractures are assumed to be equipotential, selecting several reference pointsFj on the fracture surface forms a linear equation system:(10)PF=∑i=1NSqSiμ2πhkmlnrij+CPj=1,…,NS,where rij is the distance between reference point Fj and sink Si.The linear equation system can be rewritten as follows:(11)∑i=1NSξi2πhlnrij=1j=1,…,NS,where ξi is defined as follows:(12)ξi=μqSikmPF−CP.Solving the linear equation system gives the exact expression of pressure at pointX:(13)PXx,y=∑i=1NSξilnriX2πhPF−CP+CP.The transmissibility between the fracture elementf and the neighboring matrix element m can be derived as follows:(14)Tfm=kmλ−1∑ξiλ=1Vm∫∫∑i=1NSξilnriX2πhdVX,where Vm is the volume of matrix m.Similarly, the transmissibility between the fracture elementf and the nonneighboring matrix element mP can be given as follows:(15)TfmP=kmAmPλP−1A∑ξiλP=1VmP∫∫∑i=1NSξilnriXP2πhdVXP,where VmP is the volume of matrix mP, A is the area of the interface, and AmP is the area of the fracture projection on the interface.


## 4. Modeling of Geomechanics

### 4.1. Deformation of Fractures

The deformation behavior of fracture is strongly stress-dependent and nonlinear. An empirical model based on experimental laboratory data is used to calculate the fracture moduli [32]. For a fracture under normal stress, Young’s modulus is as follows:(16)EF=kniw01−σn/Δumaxkni+σn2,(17)kni=−7.15+1.75JRC+0.02JCSw0,Δumax=−0.1032−0.0074JRC+0.02JCSw0−0.251,where kni is the initial fracture normal stiffness, w0 is the zero-stress fracture aperture, σn is the effective normal stress, Δumax is the maximum normal closure of fracture, JRC is the joint roughness coefficient, and JCS is the joint compressive strength.The fracture aperture change under normal stress is given as follows:(18)Δw=σnkni+σn/Δumax.For a fracture under shear stress, the shear modulus is as follows:(19)GF=Kjσnnj1−τRfτpeak2,(20)Kj=−17.19+3.86JRCτpeak=σntanJRClog10JCSσn+φr,where Kj is the stiffness number, nj is the stiffness exponent, τ is the shear stress, τpeak is the peak shear stress, Rf is the failure ratio, and φr is the residue friction angle.Deformation of a fracture under shear stress is the combination of horizontal shear displacement and vertical shear displacement, also known as fracture dilation. The horizontal shear displacement of a fracture is given as follows:(21)δh=τLGF,where L is the characteristic length of the fracture.The vertical shear displacement (dilation) uses an empirical model given as follows [33]:(22)δv=13tanJRClogJCSσnδh2δhδh,peak−1δh<δh,peak,∫δh,peakδhtanJRClogJCSσnδh,peakδh0.381dδh+δv,peakδh≥δh,peak,δh,peak=0.0077L0.45σnJCS0.34cosJRClogJCSσn,δv,peak=δh,peaktan13MDJRClogJCSσn,where δh,peak is the peak horizontal shear displacement, δv,peak is the peak vertical shear displacement, and MD is a damage coefficient.For a fracture under normal and shear stress, both normal aperture change and fracture dilation contribute to the overall fracture aperture change. Thus, the fracture aperture under stress is given as follows:(23)wm=w0−Δw+δv.The fracture aperture in (23) is the average point-to-point distance between two fracture surfaces, which is defined as the “mechanical” aperture. However, actual fractures have rough walls and variable aperture, and the mechanical aperture is not appropriate in calculating the hydraulic conductivity of the fracture. The “hydraulic” aperture is determined by flow analysis and is better for describing the fracture conducting capacity. An empirical relationship between hydraulic aperture and mechanical aperture can be given as follows [34]:(24)wh=wm2JRC2.5,δh≤0.75δh,peak,wm0.5JRCmob,δh≥δh,peak,JRCmob=arctanτ/σn−φrlogJCS/σn,where JRCmob is the mobilized joint roughness coefficient.Generally, the hydraulic aperture of a fracture is smaller than the mechanical aperture due to the roughness and tortuosity of fracture walls.


### 4.2. Mechanical Equilibrium Equation

The governing equation of geomechanics, also known as the mechanical equilibrium equation, is obtained from the momentum conservation law:(25)ρd2usdt2=∂σ∂x+ρf,where ρ=1−ϕρs+ϕ∑ρβSβ is the density of the fluid-solid mixture, ρs is the density of rock skeleton, us is the displacement vector of rock skeleton, σ is the total stress tensor, and f is the body force, which is gravity in this work. In static analysis, the dynamic term in the left-hand side of (25) can be omitted.In a fractured porous media, matrix and fracture are considered as two separate porous spaces that contain fluid, and the dual-porosity effective stress law is given as follows [35]:(26)σ′=σ+αMPMI+αFPFI,where σ' is the effective stress tensor, PM is the average pressure of matrix blocks, PF is the average pressure of fractures, αM is the Biot coefficient of matrix, αF is the Biot coefficient of fracture, and I is the identity matrix.Equation (25) is discretized using finite element method (FEM), which gives the following:(27)KΔuN−αMKvΔPM−αFKvΔPF=ΔF,where K=∫VBTDBdV is the nodal stiffness matrix, B is the strain-nodal displacement matrix, D is the elastic stiffness matrix, ΔuN is the nodal displacement increment vector, Kv is the volumetric stiffness vector, ΔPM and ΔPF are the average pressure increment in matrix and fractures, and ΔF is the external loading increment vector.The establishment of the elastic stiffness matrixD of a fractured porous media uses the equivalent continuum approach. In the local coordinate system of fracture, the compliance matrix can be written as follows:(28)C=γPCMP+1−γCF,where C is the compliance matrix of the equivalent jointed rock mass, CM and CF are the compliance matrixes for the rock matrix and fracture, P is a coefficient matrix, and γ is the volume fraction of the rock matrix.The compliance matrix in the global coordinate systemCglobal can be obtained from coordinate transformation:(29)Cglobal=TTCT,where T is the coordinate transformation matrix related to fracture orientation.The elastic stiffness matrixD is the inversion of the compliance matrix Cglobal:(30)D=Cglobal−1.


### 4.3. Coupling Strategy

An iterative coupling strategy is employed. The flow and geomechanical modules are invoked sequentially in every iterative loop. In each iteration, the flow simulation is first performed. The pressure and saturation information from flow simulation results are then be passed to the geomechanical module to calculate the stress, strain, and displacement results. The changes in hydraulic properties, such as porosity and permeability, are updated before the flow simulation of the next iterative step. The update for hydraulic properties of the matrix and fractures is performed separately. For matrix, empirical relationships are used [36]:(31)ϕ=ϕr+ϕ0−ϕre−aσm′k=k0ecϕ/ϕ0−1,where ϕr is the residual porosity, ϕ0 is the zero-stress porosity, σm′ is the mean effective stress, k0 is the zero-stress permeability, and a and c are the update parameters.In addition to the porosity and permeability updating functions, the capillary pressure in rock matrix is updated by the Leverett function:(32)Pc=Pc0k0/ϕ0k/ϕ.For fracture, the porosity and permeability are updated from the aperture change:(33)ϕF=ww0ϕF0.(34)kF=kF0ww02,where ϕF0 is the zero-stress fracture porosity and kF0 is the zero-stress fracture permeability.The fracture stiffness is also updated using (16) and (19) before the next iteration loop. The iteration process stops when the convergence criteria of both flow and geomechanical modules are reached. Then, the simulation process of the next time step starts.


## 4.1. Deformation of Fractures

The deformation behavior of fracture is strongly stress-dependent and nonlinear. An empirical model based on experimental laboratory data is used to calculate the fracture moduli [32]. For a fracture under normal stress, Young’s modulus is as follows:(16)EF=kniw01−σn/Δumaxkni+σn2,(17)kni=−7.15+1.75JRC+0.02JCSw0,Δumax=−0.1032−0.0074JRC+0.02JCSw0−0.251,where kni is the initial fracture normal stiffness, w0 is the zero-stress fracture aperture, σn is the effective normal stress, Δumax is the maximum normal closure of fracture, JRC is the joint roughness coefficient, and JCS is the joint compressive strength.The fracture aperture change under normal stress is given as follows:(18)Δw=σnkni+σn/Δumax.For a fracture under shear stress, the shear modulus is as follows:(19)GF=Kjσnnj1−τRfτpeak2,(20)Kj=−17.19+3.86JRCτpeak=σntanJRClog10JCSσn+φr,where Kj is the stiffness number, nj is the stiffness exponent, τ is the shear stress, τpeak is the peak shear stress, Rf is the failure ratio, and φr is the residue friction angle.Deformation of a fracture under shear stress is the combination of horizontal shear displacement and vertical shear displacement, also known as fracture dilation. The horizontal shear displacement of a fracture is given as follows:(21)δh=τLGF,where L is the characteristic length of the fracture.The vertical shear displacement (dilation) uses an empirical model given as follows [33]:(22)δv=13tanJRClogJCSσnδh2δhδh,peak−1δh<δh,peak,∫δh,peakδhtanJRClogJCSσnδh,peakδh0.381dδh+δv,peakδh≥δh,peak,δh,peak=0.0077L0.45σnJCS0.34cosJRClogJCSσn,δv,peak=δh,peaktan13MDJRClogJCSσn,where δh,peak is the peak horizontal shear displacement, δv,peak is the peak vertical shear displacement, and MD is a damage coefficient.For a fracture under normal and shear stress, both normal aperture change and fracture dilation contribute to the overall fracture aperture change. Thus, the fracture aperture under stress is given as follows:(23)wm=w0−Δw+δv.The fracture aperture in (23) is the average point-to-point distance between two fracture surfaces, which is defined as the “mechanical” aperture. However, actual fractures have rough walls and variable aperture, and the mechanical aperture is not appropriate in calculating the hydraulic conductivity of the fracture. The “hydraulic” aperture is determined by flow analysis and is better for describing the fracture conducting capacity. An empirical relationship between hydraulic aperture and mechanical aperture can be given as follows [34]:(24)wh=wm2JRC2.5,δh≤0.75δh,peak,wm0.5JRCmob,δh≥δh,peak,JRCmob=arctanτ/σn−φrlogJCS/σn,where JRCmob is the mobilized joint roughness coefficient.Generally, the hydraulic aperture of a fracture is smaller than the mechanical aperture due to the roughness and tortuosity of fracture walls.


## 4.2. Mechanical Equilibrium Equation

The governing equation of geomechanics, also known as the mechanical equilibrium equation, is obtained from the momentum conservation law:(25)ρd2usdt2=∂σ∂x+ρf,where ρ=1−ϕρs+ϕ∑ρβSβ is the density of the fluid-solid mixture, ρs is the density of rock skeleton, us is the displacement vector of rock skeleton, σ is the total stress tensor, and f is the body force, which is gravity in this work. In static analysis, the dynamic term in the left-hand side of (25) can be omitted.In a fractured porous media, matrix and fracture are considered as two separate porous spaces that contain fluid, and the dual-porosity effective stress law is given as follows [35]:(26)σ′=σ+αMPMI+αFPFI,where σ' is the effective stress tensor, PM is the average pressure of matrix blocks, PF is the average pressure of fractures, αM is the Biot coefficient of matrix, αF is the Biot coefficient of fracture, and I is the identity matrix.Equation (25) is discretized using finite element method (FEM), which gives the following:(27)KΔuN−αMKvΔPM−αFKvΔPF=ΔF,where K=∫VBTDBdV is the nodal stiffness matrix, B is the strain-nodal displacement matrix, D is the elastic stiffness matrix, ΔuN is the nodal displacement increment vector, Kv is the volumetric stiffness vector, ΔPM and ΔPF are the average pressure increment in matrix and fractures, and ΔF is the external loading increment vector.The establishment of the elastic stiffness matrixD of a fractured porous media uses the equivalent continuum approach. In the local coordinate system of fracture, the compliance matrix can be written as follows:(28)C=γPCMP+1−γCF,where C is the compliance matrix of the equivalent jointed rock mass, CM and CF are the compliance matrixes for the rock matrix and fracture, P is a coefficient matrix, and γ is the volume fraction of the rock matrix.The compliance matrix in the global coordinate systemCglobal can be obtained from coordinate transformation:(29)Cglobal=TTCT,where T is the coordinate transformation matrix related to fracture orientation.The elastic stiffness matrixD is the inversion of the compliance matrix Cglobal:(30)D=Cglobal−1.


## 4.3. Coupling Strategy

An iterative coupling strategy is employed. The flow and geomechanical modules are invoked sequentially in every iterative loop. In each iteration, the flow simulation is first performed. The pressure and saturation information from flow simulation results are then be passed to the geomechanical module to calculate the stress, strain, and displacement results. The changes in hydraulic properties, such as porosity and permeability, are updated before the flow simulation of the next iterative step. The update for hydraulic properties of the matrix and fractures is performed separately. For matrix, empirical relationships are used [36]:(31)ϕ=ϕr+ϕ0−ϕre−aσm′k=k0ecϕ/ϕ0−1,where ϕr is the residual porosity, ϕ0 is the zero-stress porosity, σm′ is the mean effective stress, k0 is the zero-stress permeability, and a and c are the update parameters.In addition to the porosity and permeability updating functions, the capillary pressure in rock matrix is updated by the Leverett function:(32)Pc=Pc0k0/ϕ0k/ϕ.For fracture, the porosity and permeability are updated from the aperture change:(33)ϕF=ww0ϕF0.(34)kF=kF0ww02,where ϕF0 is the zero-stress fracture porosity and kF0 is the zero-stress fracture permeability.The fracture stiffness is also updated using (16) and (19) before the next iteration loop. The iteration process stops when the convergence criteria of both flow and geomechanical modules are reached. Then, the simulation process of the next time step starts.


## 5. Validation of the pIEDFM

### 5.1. Validation with Benchmark Results

The simulation results of the pIEDFM are compared with the benchmark results presented in the work of Karimi-Fard et al. [37]. As shown in Figure 6, a simple fractured reservoir with horizontal and vertical fractures is considered.Figure 6
Configuration of the reservoir.The porosity of the matrix is 0.20, and the permeability is 1 mD. The fracture aperture is 0.1 mm. The block is initially saturated with oil. Water is injected from the bottom left corner. The mixture of oil and water is produced from the top right corner. The viscosities of oil and water are 0.45 cP and 1.0 cP. The relative permeability curves in both rock matrix and fractures are linear. Capillary pressure is neglected in both matrix and fractures. The simulation is performed until 2 PV of water is injected.Figure7 shows a good agreement of the cumulative oil produced between the pIEDFM results and the fine grid results from the work of Karimi-Fard et al. Figure 8 shows the water saturation profiles across the reservoir after 0.1 PV, 0.3 PV, and 0.5 PV of water injection. For comparison, the results from the EDFM simulations are presented as well. It is indicated that a better agreement is reached between the pIEDFM results and the benchmark results. The water flow across the fracture is underestimated in the EDFM, while the pIEDFM precisely recreates the water saturation profile with minor differences.Figure 7
Comparison of cumulative oil produced.Figure 8
Water saturation profiles after 0.1 PV, 0.3 PV, and 0.5 PV of water injection.


### 5.2. Validation with Fine Grid Simulation

The proposed pIEDFM is further validated by a reservoir case with one fracture, and the performances of the pIEDFM and the traditional EDFM are compared. The reservoir is 0.15 m× 0.15 m × 1.0 m in size and discretized by 30 × 30 Cartesian grids in the pIEDFM and EDFM, as shown in Figure 9. The porosity of the matrix is 0.20. The permeability of the matrix is 1 mD. The fracture aperture is 0.5 mm. The initial stress of the reservoir is 10.0 MPa in the x-direction and y-direction and 25.0 MPa in the z-direction.Figure 9
Mesh of the reservoir with one fracture.The reservoir is initially saturated with oil with the density of 800.0 kg/m3 and the viscosity of 0.45 cP. Water with the density of 1000.0 kg/m3 and the viscosity of 1.0 cP is injected from the bottom left corner. The injection rate is 0.1 m3/day. The mixture of oil and water is produced from the top right corner with fixed bottom hole pressure of 10.0 MPa. Linear relative permeability curves are used in both matrix and fractures, and capillary pressure is neglected. The initial pressure in the reservoir is 10.0 MPa. The simulation is run for 1800 seconds.The pIEDFM and EDFM simulation results are compared to the fine grid explicit-fracture results, which are assumed to be exact. In the fine grid explicit-fracture simulation, the reservoir is discretized by 300× 300 Cartesian grids. The grid size equals the aperture of the fracture, and the fracture is treated explicitly as a series of highly permeable grids. The oil saturation profile after 1800 seconds of injection is shown in Figure 10 and used as the reference results in the following discussions.Figure 10
Oil saturation profile for the fine grid simulation after 1800 seconds of injection.In Figure11, the oil saturation profiles of the pIEDFM and EDFM after 1800 seconds of injection are compared with the fine grid oil saturation profile. To support the discussion, an upscaling is performed on the fine grid oil saturation profile. It can be found that the oil saturation profile of the pEDFM shows better agreement with the fine grid results than the EDFM. Figure 12 presents the profiles of absolute oil saturation error of the pIEDFM and EDFM simulation results against the fine grid result, which also indicates that the pIEDFM outperforms the EDFM in recreating the realistic oil saturation distribution. In the EDFM results, the oil saturation around the center of the fracture is higher than the fine grid results, which shows that the water displacement across the fracture is underestimated. Besides, the oil saturation around the upper edge of the fracture is lower than the fine grid results, which shows that the water displacement along the fracture is overestimated.Figure 11
Oil saturation profiles for (a) upscaled fine grid simulation, (b) pIEDFM, and (c) EDFM after 1800 seconds of injection.
(a)(b)(c)Figure 12
Oil saturation error profiles for (a) pIEDFM and (b) EDFM after 1800 seconds of injection.
(a)(b)The simulation scenarios are repeated using finer mesh grids (60× 60). Figure 13 and Figure 14 show the oil saturation profiles and oil saturation error profiles of the pIEDFM and EDFM simulations after 1800 seconds of injection. After mesh refinement, the pIEDFM saturation profile shows a much better match with the fine grid results. However, the improvement of the EDFM results is limited, which indicates that mesh refinements cannot effectively reduce the errors caused by the unphysical flux split in the EDFM. By adding additional connections between the fractures and the nonneighboring matrix cells, the pIEDFM successfully captures the realistic flux split and reduces the errors in predicting oil saturation distributions. Figure 15 shows the flow fields of the oil phase in the fine grid, pIEDFM, and EDFM results, which demonstrates that the pIEDFM can obtain more realistic velocity fields than the EDFM around the fracture surfaces.Figure 13
Oil saturation profiles after mesh refinement for (a) upscaled fine grid simulation, (b) pIEDFM, and (c) EDFM after 1800 seconds of injection.
(a)(b)(c)Figure 14
Oil saturation error profiles after mesh refinement for (a) pIEDFM and (b) EDFM after 1800 seconds of injection.
(a)(b)Figure 15
Oil phase flow fields for (a) upscaled fine grid simulation, (b) pIEDFM, and (c) EDFM after 1800 seconds of injection.
(a)(b)(c)A mesh sensitivity analysis is performed to examine the performance of the pIEDFM at different grid densities. TheL2 norm is introduced to represent the overall error of oil saturation:(35)L2=∑i=1NSoi−Soi,fine−grid2N,where N is the number of sample points and Soi,fine−grid is the upscaled fine grid oil saturation. Figure 16 compares the L2 norm of the pIEDFM and the EDFM at different grid sizes. The results show that the error of pIEDFM decreases with the refinement of matrix grids, which demonstrates the convergence of the method. The pIEDFM has better agreement with the fine grid results than the EDFM at all grid sizes, showing the superiority of the proposed method.Figure 16
L2 norm against grid size.To demonstrate the improvement of the pIEDFM in capturing the realistic pressure fields, pressure distributions in three slices of the reservoir that pass through the lower edge of the fracture, the center of the fracture, and the upper edge of the fracture, respectively, are investigated, as shown in Figure17. Figure 18 shows the pressure distributions of the fine grid results in the three slices, which are obviously nonlinear.Figure 17
Oil pressure distribution of the fine grid simulation after 1800 seconds of injection.Figure 18
Oil pressure distributions in the slices of the fine grid pressure profile.Figure19 shows the distributions of errors in the oil pressure of the pIEDFM and EDFM results compared with the fine grid results. It can be found that the pIEDFM has significant advantages in predicting the pressure distribution than the EDFM. In the EDFM, the nonlinear pressure changes near fractures cannot be considered, which results in apparent errors near the fracture surface, while the pIEDFM better captures the sharp pressure change in the vicinity of fractures.Figure 19
Oil pressure errors in the (a) lower edge, (b) center, and (c) upper edge slices.
(a)(b)(c)


## 5.1. Validation with Benchmark Results

The simulation results of the pIEDFM are compared with the benchmark results presented in the work of Karimi-Fard et al. [37]. As shown in Figure 6, a simple fractured reservoir with horizontal and vertical fractures is considered.Figure 6
Configuration of the reservoir.The porosity of the matrix is 0.20, and the permeability is 1 mD. The fracture aperture is 0.1 mm. The block is initially saturated with oil. Water is injected from the bottom left corner. The mixture of oil and water is produced from the top right corner. The viscosities of oil and water are 0.45 cP and 1.0 cP. The relative permeability curves in both rock matrix and fractures are linear. Capillary pressure is neglected in both matrix and fractures. The simulation is performed until 2 PV of water is injected.Figure7 shows a good agreement of the cumulative oil produced between the pIEDFM results and the fine grid results from the work of Karimi-Fard et al. Figure 8 shows the water saturation profiles across the reservoir after 0.1 PV, 0.3 PV, and 0.5 PV of water injection. For comparison, the results from the EDFM simulations are presented as well. It is indicated that a better agreement is reached between the pIEDFM results and the benchmark results. The water flow across the fracture is underestimated in the EDFM, while the pIEDFM precisely recreates the water saturation profile with minor differences.Figure 7
Comparison of cumulative oil produced.Figure 8
Water saturation profiles after 0.1 PV, 0.3 PV, and 0.5 PV of water injection.


## 5.2. Validation with Fine Grid Simulation

The proposed pIEDFM is further validated by a reservoir case with one fracture, and the performances of the pIEDFM and the traditional EDFM are compared. The reservoir is 0.15 m× 0.15 m × 1.0 m in size and discretized by 30 × 30 Cartesian grids in the pIEDFM and EDFM, as shown in Figure 9. The porosity of the matrix is 0.20. The permeability of the matrix is 1 mD. The fracture aperture is 0.5 mm. The initial stress of the reservoir is 10.0 MPa in the x-direction and y-direction and 25.0 MPa in the z-direction.Figure 9
Mesh of the reservoir with one fracture.The reservoir is initially saturated with oil with the density of 800.0 kg/m3 and the viscosity of 0.45 cP. Water with the density of 1000.0 kg/m3 and the viscosity of 1.0 cP is injected from the bottom left corner. The injection rate is 0.1 m3/day. The mixture of oil and water is produced from the top right corner with fixed bottom hole pressure of 10.0 MPa. Linear relative permeability curves are used in both matrix and fractures, and capillary pressure is neglected. The initial pressure in the reservoir is 10.0 MPa. The simulation is run for 1800 seconds.The pIEDFM and EDFM simulation results are compared to the fine grid explicit-fracture results, which are assumed to be exact. In the fine grid explicit-fracture simulation, the reservoir is discretized by 300× 300 Cartesian grids. The grid size equals the aperture of the fracture, and the fracture is treated explicitly as a series of highly permeable grids. The oil saturation profile after 1800 seconds of injection is shown in Figure 10 and used as the reference results in the following discussions.Figure 10
Oil saturation profile for the fine grid simulation after 1800 seconds of injection.In Figure11, the oil saturation profiles of the pIEDFM and EDFM after 1800 seconds of injection are compared with the fine grid oil saturation profile. To support the discussion, an upscaling is performed on the fine grid oil saturation profile. It can be found that the oil saturation profile of the pEDFM shows better agreement with the fine grid results than the EDFM. Figure 12 presents the profiles of absolute oil saturation error of the pIEDFM and EDFM simulation results against the fine grid result, which also indicates that the pIEDFM outperforms the EDFM in recreating the realistic oil saturation distribution. In the EDFM results, the oil saturation around the center of the fracture is higher than the fine grid results, which shows that the water displacement across the fracture is underestimated. Besides, the oil saturation around the upper edge of the fracture is lower than the fine grid results, which shows that the water displacement along the fracture is overestimated.Figure 11
Oil saturation profiles for (a) upscaled fine grid simulation, (b) pIEDFM, and (c) EDFM after 1800 seconds of injection.
(a)(b)(c)Figure 12
Oil saturation error profiles for (a) pIEDFM and (b) EDFM after 1800 seconds of injection.
(a)(b)The simulation scenarios are repeated using finer mesh grids (60× 60). Figure 13 and Figure 14 show the oil saturation profiles and oil saturation error profiles of the pIEDFM and EDFM simulations after 1800 seconds of injection. After mesh refinement, the pIEDFM saturation profile shows a much better match with the fine grid results. However, the improvement of the EDFM results is limited, which indicates that mesh refinements cannot effectively reduce the errors caused by the unphysical flux split in the EDFM. By adding additional connections between the fractures and the nonneighboring matrix cells, the pIEDFM successfully captures the realistic flux split and reduces the errors in predicting oil saturation distributions. Figure 15 shows the flow fields of the oil phase in the fine grid, pIEDFM, and EDFM results, which demonstrates that the pIEDFM can obtain more realistic velocity fields than the EDFM around the fracture surfaces.Figure 13
Oil saturation profiles after mesh refinement for (a) upscaled fine grid simulation, (b) pIEDFM, and (c) EDFM after 1800 seconds of injection.
(a)(b)(c)Figure 14
Oil saturation error profiles after mesh refinement for (a) pIEDFM and (b) EDFM after 1800 seconds of injection.
(a)(b)Figure 15
Oil phase flow fields for (a) upscaled fine grid simulation, (b) pIEDFM, and (c) EDFM after 1800 seconds of injection.
(a)(b)(c)A mesh sensitivity analysis is performed to examine the performance of the pIEDFM at different grid densities. TheL2 norm is introduced to represent the overall error of oil saturation:(35)L2=∑i=1NSoi−Soi,fine−grid2N,where N is the number of sample points and Soi,fine−grid is the upscaled fine grid oil saturation. Figure 16 compares the L2 norm of the pIEDFM and the EDFM at different grid sizes. The results show that the error of pIEDFM decreases with the refinement of matrix grids, which demonstrates the convergence of the method. The pIEDFM has better agreement with the fine grid results than the EDFM at all grid sizes, showing the superiority of the proposed method.Figure 16
L2 norm against grid size.To demonstrate the improvement of the pIEDFM in capturing the realistic pressure fields, pressure distributions in three slices of the reservoir that pass through the lower edge of the fracture, the center of the fracture, and the upper edge of the fracture, respectively, are investigated, as shown in Figure17. Figure 18 shows the pressure distributions of the fine grid results in the three slices, which are obviously nonlinear.Figure 17
Oil pressure distribution of the fine grid simulation after 1800 seconds of injection.Figure 18
Oil pressure distributions in the slices of the fine grid pressure profile.Figure19 shows the distributions of errors in the oil pressure of the pIEDFM and EDFM results compared with the fine grid results. It can be found that the pIEDFM has significant advantages in predicting the pressure distribution than the EDFM. In the EDFM, the nonlinear pressure changes near fractures cannot be considered, which results in apparent errors near the fracture surface, while the pIEDFM better captures the sharp pressure change in the vicinity of fractures.Figure 19
Oil pressure errors in the (a) lower edge, (b) center, and (c) upper edge slices.
(a)(b)(c)


## 6. Applications of pIEDFM in Coupled Flow and Geomechanical Simulations

The proposed pIEDFM is applied in coupled flow and geomechanical reservoir simulations. The simulation scenarios are repeated using the original EDFM for comparisons of the reservoir geomechanical behaviors predicted by the pIEDFM and the EDFM.A 2D reservoir with natural fractures that penetrate the entire depth of the reservoir is investigated, and the configuration of which is shown in Figure20. The principle stress directions of the reservoir are assumed to align with the axis directions, and zero internal friction angle is considered. Thus, the average strike angle of the fractures is ±45° from the axis directions. The average length of the fractures is 10m. The reservoir is discretized by 30 × 30 Cartesian grids. The porosity of the matrix is 0.15. The zero-stress permeability of the matrix is 5 mD. The zero-stress fracture aperture is 0.3 mm. The initial stress is 6.0 MPa in the x-direction and y-direction and 18.0 MPa in the z-direction. Under the initial stress condition, the initial matrix permeability is 2.023 mD, and the initial fracture aperture is 0.222 mm.Figure 20
Structure of the fractured reservoir.The reservoir is initially saturated with oil with the density of 800.0 kg/m3 and the viscosity of 5.0 cP. Water is injected from an injection well located at the bottom left corner at the constant bottom hole pressure of 20.0 MPa. The water density is 1000.0 kg/m3, and the water viscosity is 0.9 cP. The mixture of oil and water is produced from the top right corner with fixed bottom hole pressure of 5.0 MPa. The relative permeability for fracture is assumed to be linear, and the Brooks–Corey model is used for the relative permeability of the matrix, as shown in Figure 21. The capillary pressure curve of the rock matrix is presented in Figure 22, while capillary pressure in the fractures is neglected. The initial pressure in the reservoir is 20.0 MPa. The mechanical parameters are listed in Table 1. The simulations are run for 20000 days of production.Figure 21
Relative permeability curves of the rock matrix and fracture.Figure 22
Capillary pressure of the rock matrix.Table 1
Mechanical parameters.
ParameterValueMatrix Young’s modulus,E15.0 GPaMatrix Poisson ratio,v0.25Joint roughness coefficient, JRC12.0Joint compressive strength, JCS65.0 MPaFailure ratio,Rf1.0Stiffness exponent,nj0.8Fracture residue angle,φr30.0Matrix residue porosity,ϕr0.18Matrix porosity update parameter,a0.05 MPa−1Matrix permeability update parameter,c23.0Biot coefficient of matrix,αM0.875Biot coefficient of fracture,αF0.025Figure23 compares the oil saturation profiles for the pIEDFM and EDFM simulations after 5000, 10000, and 20000 days of production. The oil saturation profiles demonstrate that the water flow across the fractures is more prominent in the pIEDFM results. Due to the limitation in capturing the proper flux split through a fracture, the EDFM is incapable of representing the accurate multiphase displacements across a fracture.Figure 23
Oil saturation profiles for the pIEDFM simulation after (a) 5000 days, (c) 10000 days, and (e) 20000 days of production and the EDFM simulation after (b) 5000 days, (d) 10000 days, and (f) 20000 days of production.
(a)(b)(c)(d)(e)(f)Figure24 shows the producing rate histories of oil and water from the production well. In both coupled and uncoupled simulations, the water producing rate predicted by the pIEDFM is larger than the EDFM, which can be attributed to the underestimation of water displacement process across fractures by the EDFM predictions. Due to the more prominent water flow across the fractures in the pIEDFM simulations, more oil is displaced from the reservoir, resulting in the higher steady-state oil rate. Figure 25 presents the cumulative oil and water production histories, showing that the difference in the water displacement process between the pIEDFM and the EDFM affects not only the results of flow simulation but also the geomechanical performance of the reservoir. At 15000 days of production, the ratio of cumulative oil production in the coupled case to the uncoupled case is 92.7% for the pIEDFM simulation and 93.9% for the EDFM simulation, and the ratio of cumulative water production is 20.8% for the pIEDFM simulation and 19.0% for the EDFM simulation.Figure 24
(a) Oil and (b) water producing rate histories.
(a)(b)Figure 25
Cumulative (a) oil and (b) water production histories.
(a)(b)Figure26 compares the water cut histories of pIEDFM and EDFM simulations. In uncoupled simulations, the water breakthrough time is day 6112 for the pIEDFM and day 5904 for the EDFM. The primary reason for the earlier water breakthrough in the EDFM simulations can be attributed to the overestimation of flow along the fractures since more than half of the fractures orient towards the production well. When geomechanics is considered, the water breakthrough time will be delayed due to the permeability reduction in the rock matrix and the closure of fractures. In coupled simulations, the water breakthrough time is day 10926 for the pIEDFM and day 11157 for the EDFM. In contrast, the predicted water breakthrough time is earlier in the pIEDFM simulation, which indicates that the EDFM simulation is more affected by the coupling effect of geomechanics.Figure 26
Water cut histories.The mean effective stress profiles of the pIEDFM and the EDFM simulations after 20000 days of production are compared in Figure27. The mean effective stress profiles show that the reservoir is more consolidated during the depletion in the EDFM simulations, which is also supported by the contour map presented in Figure 28.Figure 27
Mean effective stress profiles for (a) pIEDFM and (b) EDFM after 20000 days of production.
(a)(b)Figure 28
Contour map of mean effective stress field after 20000 days of production.In coupled simulations, the error in the effective stress field leads to differences in permeability distribution between the pIEDFM and EDFM results. Figure29 presents the permeability profiles for the pIEDFM and EDFM simulations at 20000 days of production, which indicates that the reservoir is more prone to the permeability reduction in the EDFM simulations. Figure 30 shows the distribution of differences in the permeability field between the pIEDFM and EDFM results, and the contour maps are compared in Figure 31. In most areas across the reservoir, the permeability of the pIEDFM results is higher than the EDFM results, and the main differences in permeability appear near the production well.Figure 29
Permeability profiles for (a) pIEDFM and (b) EDFM after 20000 days of production.
(a)(b)Figure 30
Permeability difference between the pIEDFM and EDFM results.Figure 31
Contour map of permeability distribution after 20000 days of production.The apparent delay in water breakthrough time, the differences in effective stress, and permeability distributions demonstrate that the EDFM might lead to deviations in the results of coupled flow and geomechanical simulations, and the pIEDFM could effectively reduce these deviations.


## 7. Conclusions

A projection-based integrally embedded discrete fracture model is proposed. In the pIEDFM, additional connections are added between fracture elements and the nonneighboring matrix elements to provide a more realistic representation of the flux split process of water inflow across a fracture. The transmissibility of connections between the fracture element and the neighboring matrix element and the transmissibility of connections between the fracture element and the additional nonneighboring matrix element are calculated separately using a semi-analytical cone-shaped pressure distribution around the fracture surfaces. The accuracy of the pIEDFM is validated by the benchmark results and the explicit-fracture fine grid results. The performances of the pIEDFM and the EDFM are compared and discussed through a single fracture water flooding case. The pIEDFM is applied in coupled flow and geomechanical simulations, and the results are compared with that of the EDFM. The conclusions are as follows:(1)
Good agreements in the saturation profiles are reached between the pIEDFM simulations and the benchmark results, which addresses the accuracy of the pIEDFM in modeling the multiphase flow process in fractured media.(2)
Comparison of the EDFM and the pIEDFM results with the explicit-fracture fine grid simulation results shows that the pIEDFM can obtain a more physically authentic velocity field and better predict the multiphase flow process in fractured reservoirs. The oil saturation of the pIEDFM shows good agreement with the fine grid results using a moderate amount of meshes.(3)
The pIEDFM has significant advantages in predicting the pressure distribution than the EDFM. Compared with the EDFM results, the pressure errors around the fracture surfaces are obviously reduced in the pIEDFM, showing that the nonlinear pressure change in the vicinity of the fracture can be captured.(4)
Application of the pIEDFM in coupled flow and geomechanical simulations shows that the differences in predicting the water displacement process between the pIEDFM and the EDFM affect not only the results of production and saturation profiles but also the geomechanical performance of the reservoir. The biased water displacement process in the EDFM leads to deviations in the predictions of the water breakthrough time, the effective stress field, and the stress-dependent permeability distribution. Thus, it is promising to incorporate the pIEDFM in coupled flow and geomechanical simulations for fractured reservoirs.



---

*Source: 7531948-2022-02-23.xml*