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

The embedded discrete fracture model (EDFM) has been popular for the modeling of fractured reservoirs due to its ﬂexibility and eﬃciency while maintaining the complex geometry of fracture networks. Though the EDFM has been validated for single-phase ﬂow simulations, some recent cases show that the EDFM might result in apparent errors in multiphase ﬂow 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 ﬁelds. As a signiﬁcant 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 ﬁne grid simulation results show that the pIEDFM provides accurate solutions using a moderate amount of grids. The proposed pIEDFM is also applied to coupled ﬂow 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 diﬀerent methods for matrix-fracture transmissibility between the pIEDFM and the EDFM may lead to deviations in stress ﬁelds predicted by geomechanical modeling, which eventually aﬀects the oil production, water cut, and oil saturation proﬁles.


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 dualcontinuum 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][13][14][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].
e 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. e fractures are discretized explicitly as control volumes, also known as fracture elements. e fracture elements are embedded into their parent matrix grids through matrix-fracture connections. e 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 dualcontinuum methods to model reservoirs with multiscale fractures. Moinfar et al. [24] applied the EDFM in coupled flow and geomechanical simulations of fractured reservoirs.
ough 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. Figure 1 shows a matrix block M intersected by a fracture. e 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 (F 1 ) and across (F 2 ) 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. us, 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 (F m-f ), as shown in Figure 1(b). e 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. erefore, the water flux along the fracture is overestimated, and the water flux across the fracture is underestimated [25].
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 Figure 2. e 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. e IEDFM is later extended to the modeling of anisotropic fractured reservoirs [31].
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. e transmissibilities of neighboring and nonneighboring matrix-fracture connections in the pIEDFM are derived semi-analytically using the methods in the IEDFM. e accuracy of the pIEDFM is validated by benchmark results and fine grid simulation results. e proposed pIEDFM is then applied in coupled flow and geomechanical simulation for fractured reservoirs. e applicability of the proposed numerical method is examined.

Governing Equations for Fractured
Reservoir Simulation

Mass Conservation Equations.
e mass conservation equations in fractured media are given as follows: 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. e velocity of phase β is defined by Darcy's law: where k is the absolute permeability, k rβ 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 β.

Geofluids
Equation (1) is discretized using the control-volume finite difference scheme in space and first-order scheme in time, which gives the following: 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 fm β, ik is the flow term for phase β through the matrix-fracture connection of element i and element k, and Q W β,i is the flux term of phase β from wells. e flow terms F β, ij and F fm β, ik are expressed as follows: where subscripts ij + 1/2 and ik + 1/2 denote proper averages of properties at the interface, and T ij and T ik are the transmissibilities of the connections. e transmissibility of a connection is the harmonic average of two half-transmissibilities:   Geofluids where T 1 and T 2 are the half-transmissibilities of element 1 and element 2, k 1 and k 2 are the absolute permeabilities, and d 1 and d 2 are the distances from the centers of elements to the interface.

EDFM.
e reservoir matrix is discretized using the Cartesian grids, and additional fracture elements are added to represent the fracture control volumes. us, 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: where A f is the fracture surface area, k m is the matrix permeability, k f is the fracture permeability, d f is the fracture center distance from the interface, which equals half of the fracture aperture, and d fm 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: where r if is the distance from fracture and V is the volume of the matrix.

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 Figure 3, the fracture element f has two projections on the boundary of matrix element i, which have the area of A P mX and A P mY , respectively.
e fracture element f is connected to its neighboring matrix element i and two nonneighboring matrix elements j and k. e criterion of selecting the matrix faces of fracture projections follows the work of Jiang and Younis [25]. e matrix faces that are closer to the fracture center are selected as the projected faces in each dimension.
For matrix-matrix connections, the projected areas of fractures are eliminated from the interface area. e modified interface area between matrix i and matrix j is given as follows: where A ij is the original interface area between matrix i and matrix j, and A P mX is the fracture projection area along the xdirection. When a fracture cuts through the matrix elements, the matrix-matrix connection will be removed. Figure 4 shows an example of the connection establishment in the pIEDFM and the EDFM. ere are four matrix elements marked as M1, M2, M3, and M4. e two fractures are discretized into several fracture segment elements by the matrix block boundaries. e 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 matrixfracture connections. However, in the pIEDFM, 3 matrixmatrix connections, 5 fracture-fracture connections, and 14 matrix-fracture connections are included. e number of matrix-fracture connections increases in the pIEDFM due to the additional connections between the fracture segments and the nonneighboring matrix elements.
In the pIEDFM, the matrix-matrix and fracture-fracture connection transmissibilities can be directly derived using (5). e 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. e 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. Figure 5 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: where h is the height of the strata, N S is the number of point sinks, q S i is the flow rate of sink S i , r iX is the distance between point X and sink S i , and C P is a constant related to the fracture pressure.

Geofluids
When fractures are assumed to be equipotential, selecting several reference points F j on the fracture surface forms a linear equation system: where r ij is the distance between reference point F j and sink S i . e linear equation system can be rewritten as follows: where ξ i is defined as follows: Solving the linear equation system gives the exact expression of pressure at point X: e transmissibility between the fracture element f and the neighboring matrix element m can be derived as follows: where V m is the volume of matrix m. Similarly, the transmissibility between the fracture element f and the nonneighboring matrix element m P can be given as follows: where V m P is the volume of matrix m P , A is the area of the interface, and A P m is the area of the fracture projection on the interface.

Deformation of Fractures.
e 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: where k ni is the initial fracture normal stiffness, w 0 is the zero-stress fracture aperture, σ n is the effective normal stress, Δu max is the maximum normal closure of fracture, JRC is the joint roughness coefficient, and JCS is the joint compressive strength. e fracture aperture change under normal stress is given as follows: For a fracture under shear stress, the shear modulus is as follows: where K j is the stiffness number, n j is the stiffness exponent, τ is the shear stress, τ peak is the peak shear stress, R f 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. e horizontal shear displacement of a fracture is given as follows: Geofluids 5 where L is the characteristic length of the fracture. e vertical shear displacement (dilation) uses an empirical model given as follows [33]: where δ h,peak is the peak horizontal shear displacement, δ v,peak is the peak vertical shear displacement, and M D 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.
us, the fracture aperture under stress is given as follows: e 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.
e "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]: where JRC mob 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.

Mechanical Equilibrium Equation.
e governing equation of geomechanics, also known as the mechanical equilibrium equation, is obtained from the momentum conservation law: where ρ � (1 − ϕ)ρ s + ϕ ρ β S β is the density of the fluidsolid mixture, ρ s is the density of rock skeleton, u s 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]: where σ ' is the effective stress tensor, P M is the average pressure of matrix blocks, P F 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: where K � V B T DBdV is the nodal stiffness matrix, B is the strain-nodal displacement matrix, D is the elastic stiffness matrix, Δu N is the nodal displacement increment vector, K v is the volumetric stiffness vector, ΔP M and ΔP F are the average pressure increment in matrix and fractures, and ΔF is the external loading increment vector. e establishment of the elastic stiffness matrix D 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: where C is the compliance matrix of the equivalent jointed rock mass, C M and C F are the compliance matrixes for the rock matrix and fracture, P is a coefficient matrix, and c is the volume fraction of the rock matrix. e compliance matrix in the global coordinate system C global can be obtained from coordinate transformation: where T is the coordinate transformation matrix related to fracture orientation. e elastic stiffness matrix D is the inversion of the compliance matrix C global :

Coupling
Strategy. An iterative coupling strategy is employed. e flow and geomechanical modules are invoked sequentially in every iterative loop. In each iteration, the flow simulation is first performed. e pressure and saturation 6 Geofluids  information from flow simulation results are then be passed to the geomechanical module to calculate the stress, strain, and displacement results. e changes in hydraulic properties, such as porosity and permeability, are updated before the flow simulation of the next iterative step. e update for hydraulic properties of the matrix and fractures is performed separately. For matrix, empirical relationships are used [36]: where ϕ r is the residual porosity, ϕ 0 is the zero-stress porosity, σ m ′ is the mean effective stress, k 0 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: For fracture, the porosity and permeability are updated from the aperture change: where ϕ F0 is the zero-stress fracture porosity and k F0 is the zero-stress fracture permeability. e fracture stiffness is also updated using (16) and (19) before the next iteration loop. e iteration process stops when the convergence criteria of both flow and geomechanical modules are reached. en, the simulation process of the next time step starts.

Validation with Benchmark Results.
e 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. e porosity of the matrix is 0.20, and the permeability is 1 mD. e fracture aperture is 0.1 mm. e block is initially saturated with oil. Water is injected from the bottom left corner. e mixture of oil and water is produced from the top right corner. e viscosities of oil and water are 0.45 cP and 1.0 cP. e relative permeability curves in both rock matrix and fractures are linear. Capillary pressure is neglected in both matrix and fractures. e simulation is performed until 2 PV of water is injected. Figure 7 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. e water flow across the fracture is underestimated in the EDFM, while the pIEDFM precisely recreates the water saturation profile with minor differences.

Validation with Fine Grid Simulation.
e proposed pIEDFM is further validated by a reservoir case with one fracture, and the performances of the pIEDFM and the traditional EDFM are compared. e reservoir is 0. e simulation is run for 1800 seconds. e 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. e grid size equals the aperture of the fracture, and the fracture is treated explicitly as a series of highly permeable grids. e oil saturation profile after 1800 seconds of injection is shown in Figure 10 and used as the reference results in the following discussions.
In Figure 11, 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. e 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.
A mesh sensitivity analysis is performed to examine the performance of the pIEDFM at different grid densities. e L 2 norm is introduced to represent the overall error of oil saturation: where N is the number of sample points and S i, fine−grid o is the upscaled fine grid oil saturation. Figure 16 compares the L 2 norm of the pIEDFM and the EDFM at different grid sizes.
e results show that the error of pIEDFM decreases with the refinement of matrix grids, which demonstrates the convergence of the method. e pIEDFM has better 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 Figure 17. Figure 18 shows the pressure distributions of the fine grid results in the three slices, which are obviously nonlinear. Figure 19 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.

Applications of pIEDFM in Coupled Flow and
Geomechanical Simulations e proposed pIEDFM is applied in coupled flow and geomechanical reservoir simulations. e 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 Figure 20. e principle stress directions of the reservoir are assumed to align with the axis directions, and zero internal friction angle is considered. us, the average strike angle of the fractures is ±45°from the axis directions. e average length of the fractures is 10m. e reservoir is discretized by 30 × 30 Cartesian grids. e porosity of the matrix is 0.15. e zero-stress permeability of the matrix is 5 mD. e zero-stress fracture aperture is 0.3 mm. e 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. cP. e mixture of oil and water is produced from the top right corner with fixed bottom hole pressure of 5.0 MPa. e 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. e capillary pressure curve of the rock matrix is presented in Figure 22, while capillary pressure in the fractures is neglected. e initial pressure in the reservoir is 20.0 MPa. e mechanical parameters are listed in Table 1. e simulations are run for 20000 days of production. Figure 23 compares the oil saturation profiles for the pIEDFM and EDFM simulations after 5000, 10000, and 20000 days of production. e 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 24 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 26 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. e 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. e mean effective stress profiles of the pIEDFM and the EDFM simulations after 20000 days of production are compared in Figure 27. e 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.
In coupled simulations, the error in the effective stress field leads to differences in permeability distribution between the pIEDFM and EDFM results. Figure 29 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. e 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.

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. e 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. e accuracy of the pIEDFM is validated by the benchmark results and the explicit-fracture fine grid results. e performances of the pIEDFM and the EDFM are compared and discussed through a single fracture water flooding case. e pIEDFM is applied in coupled flow and geomechanical simulations, and the results are compared with that of the EDFM. e 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. e oil saturation of the pIEDFM shows good agreement with the fine grid results using a moderate amount of meshes. (3) e 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. e 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. us, it is promising to incorporate the pIEDFM in coupled flow and geomechanical simulations for fractured reservoirs.

Data Availability
All data used can be found in our manuscript.