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

Wu, Dawei; Di, Yuan; Wu, Yu-Shu

doi:https://doi.org/10.1155/2022/7531948

Geofluids

On this page

Abstract Introduction Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Special Issue

Multiscale and Multiphysics Approaches to Fluids Flow in Unconventional Reservoirs 2021

View this Special Issue

Research Article | Open Access

Volume 2022 | Article ID 7531948 | https://doi.org/10.1155/2022/7531948

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

Dawei Wu,^1,2Yuan Di ,^1,2and Yu-Shu Wu³

Academic Editor: Stefan Iglauer

Received08 Jul 2021

Revised16 Nov 2021

Accepted18 Jan 2022

Published23 Feb 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.

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. Figure 1 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 (F₁) and across (F₂) 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 (F_m-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].

(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 Figure 2. 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].

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:where represents fluid phases, is the porosity, is the saturation of phase , is the density of phase , is the velocity of phase , is the flux term of phase from wells, and is the flux term of phase between matrix and fracture elements.

The velocity of phase is defined by Darcy’s law:where is the absolute permeability, is the relative permeability of phase , is the viscosity of phase , is the pressure of phase , is the gravitational acceleration, 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:where subscript denotes the values of element , superscript represents the current time level, superscript represents the previous time level, is the flow term for phase between element and element where and are the same type of element (matrix or fracture), is the flow term for phase through the matrix-fracture connection of element and element , and is the flux term of phase from wells.

The flow terms and are expressed as follows:where subscripts and denote proper averages of properties at the interface, and and are the transmissibilities of the connections.

The transmissibility of a connection is the harmonic average of two half-transmissibilities:where and are the half-transmissibilities of element 1 and element 2, and are the absolute permeabilities, and and 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:where is the fracture surface area, is the matrix permeability, is the fracture permeability, is the fracture center distance from the interface, which equals half of the fracture aperture, and 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 is the distance from fracture and 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 Figure 3, the fracture element has two projections on the boundary of matrix element , which have the area of and , respectively. The fracture element is connected to its neighboring matrix element and two nonneighboring matrix elements and . 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.

For matrix-matrix connections, the projected areas of fractures are eliminated from the interface area. The modified interface area between matrix and matrix is given as follows:where is the original interface area between matrix and matrix , and is the fracture projection area along the x-direction. 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. 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.

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.

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 can be derived as the superposition of all the pressure drops caused by point sinks:where is the height of the strata, is the number of point sinks, is the flow rate of sink , is the distance between point and sink , and is a constant related to the fracture pressure.

When fractures are assumed to be equipotential, selecting several reference points on the fracture surface forms a linear equation system:where is the distance between reference point and sink .

The linear equation system can be rewritten as follows:where is defined as follows:

Solving the linear equation system gives the exact expression of pressure at point :

The transmissibility between the fracture element and the neighboring matrix element can be derived as follows:where is the volume of matrix .

Similarly, the transmissibility between the fracture element and the nonneighboring matrix element can be given as follows:where is the volume of matrix , is the area of the interface, and 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:where is the initial fracture normal stiffness, is the zero-stress fracture aperture, is the effective normal stress, 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:

For a fracture under shear stress, the shear modulus is as follows:where is the stiffness number, is the stiffness exponent, is the shear stress, is the peak shear stress, is the failure ratio, and 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:where is the characteristic length of the fracture.

The vertical shear displacement (dilation) uses an empirical model given as follows [33]:where is the peak horizontal shear displacement, is the peak vertical shear displacement, and 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:

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]:where 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:where is the density of the fluid-solid mixture, is the density of rock skeleton, is the displacement vector of rock skeleton, is the total stress tensor, and 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, is the average pressure of matrix blocks, is the average pressure of fractures, is the Biot coefficient of matrix, is the Biot coefficient of fracture, and is the identity matrix.

Equation (25) is discretized using finite element method (FEM), which gives the following:where is the nodal stiffness matrix, is the strain-nodal displacement matrix, is the elastic stiffness matrix, is the nodal displacement increment vector, is the volumetric stiffness vector, and are the average pressure increment in matrix and fractures, and is the external loading increment vector.

The establishment of the elastic stiffness matrix 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 is the compliance matrix of the equivalent jointed rock mass, and are the compliance matrixes for the rock matrix and fracture, is a coefficient matrix, and is the volume fraction of the rock matrix.

The compliance matrix in the global coordinate system can be obtained from coordinate transformation:where is the coordinate transformation matrix related to fracture orientation.

The elastic stiffness matrix is the inversion of the compliance matrix :

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]:where is the residual porosity, is the zero-stress porosity, is the mean effective stress, is the zero-stress permeability, and and 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 is the zero-stress fracture porosity and 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.

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.

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. The water flow across the fracture is underestimated in the EDFM, while the pIEDFM precisely recreates the water saturation profile with minor differences.

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.

The reservoir is initially saturated with oil with the density of 800.0 kg/m³ and the viscosity of 0.45 cP. Water with the density of 1000.0 kg/m³ and the viscosity of 1.0 cP is injected from the bottom left corner. The injection rate is 0.1 m³/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.

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.

(a)

(b)

(c)

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

(a)

(b)

(c)

(a)

(b)

(a)

(b)

(c)

A mesh sensitivity analysis is performed to examine the performance of the pIEDFM at different grid densities. The norm is introduced to represent the overall error of oil saturation:where is the number of sample points and is the upscaled fine grid oil saturation. Figure 16 compares the 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.

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.

(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 Figure 20. 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 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.

The reservoir is initially saturated with oil with the density of 800.0 kg/m³ 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/m³, 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 23 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.

(a)

(b)

(c)

(d)

(e)

(f)

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.

(a)

(b)

(a)

(b)

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

The mean effective stress profiles of the pIEDFM and the EDFM simulations after 20000 days of production are compared in Figure 27. 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.

(a)

(b)

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.

(a)

(b)

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.

Data Availability

All data used can be found in our manuscript.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this study.

Acknowledgments

This project was funded by the Joint Fund for Enterprise Innovation and Development of NSFC (Grant no. U19B6003-02-06) and the National Natural Science Foundation of China (Grant no. 51674010).

References

G. I. Barenblatt, I. P. Zheltov, and I. N. Kochina, “Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata],” Journal of Applied Mathematics and Mechanics, vol. 24, no. 5, pp. 1286–1303, 1960.
View at: Publisher Site | Google Scholar
J. E. Warren and P. J. Root, “The behavior of naturally fractured reservoirs,” Society of Petroleum Engineers Journal, vol. 3, no. 3, pp. 245–255, 1963.
View at: Publisher Site | Google Scholar
H. Kazemi, “Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution,” Society of Petroleum Engineers Journal, vol. 9, no. 4, pp. 451–462, 1969.
View at: Publisher Site | Google Scholar
H. Kazemi H, L. S. Merril, and K. L. Porterfield, “Numerical simulation of water-oil flow in naturally fractured,” Society of Petroleum Engineers Journal, vol. 16, no. 6, pp. 317–326, 1976.
View at: Publisher Site | Google Scholar
J. R. Gilman, “An efficient finite-difference method for simulating phase segregation in the matrix blocks in double-porosity reservoirs,” SPE Reservoir Engineering, vol. 1, no. 4, pp. 403–413, 1986.
View at: Publisher Site | Google Scholar
J. R. Gilman and H. Kazemi, “Improved calculations for viscous and gravity displacement in matrix blocks in dual-porosity simulators,” Journal of Petroleum Technology, vol. 40, no. 1, pp. 60–70, 1988.
View at: Publisher Site | Google Scholar
Y.-S. Wu, Y. Di, Z. Kang, and P. Fakcharoenphol, “A multiple-continuum model for simulating single-phase and multiphase flow in naturally fractured vuggy reservoirs,” Journal of Petroleum Science and Engineering, vol. 78, no. 1, pp. 13–22, 2011.
View at: Publisher Site | Google Scholar
Y. S. Wu and K. Pruess, “A multiple-porosity method for simulation of naturally fractured petroleum reservoirs,” SPE Reservoir Engineering, vol. 3, no. 1, pp. 327–336, 1988.
View at: Publisher Site | Google Scholar
S. Geiger-Boschung, S. Matthai, and J. Niessner, “Black-oil simulations for three-component, three-phase flow in fractured porous media,” SPE Journal, vol. 14, no. 2, pp. 338–354, 2007.
View at: Publisher Site | Google Scholar
P. Panfili, A. Cominelli, and R. Colin, “Efficient and effective field scale simulation of hydraulic fractured wells: methodology and application,” in SPE Reservoir Characterisation and Simulation Conference and Exhibition, UAE, Abu Dhabi, 2015.
View at: Publisher Site | Google Scholar
A. Moinfar, W. Narr, and M. H. Hui, “Comparison of discrete-fracture and dual-permeability models for multiphase flow in naturally fractured reservoirs,” in SPE Reservoir Simulation Symposium, Woodlands, Texas, USA, 2011.
View at: Publisher Site | Google Scholar
J. Noorishad and M. Mehran, “An upstream finite element method for solution of transient transport equation in fractured porous media,” Water Resources Research, vol. 18, no. 3, pp. 588–596, 1982.
View at: Publisher Site | Google Scholar
M. Karimi-Fard and A. Firoozabadi, “Numerical simulation of water injection in fractured media using the discrete-fracture model and the Galerkin method,” SPE Reservoir Evaluation and Engineering, vol. 6, no. 2, pp. 117–126, 2003.
View at: Publisher Site | Google Scholar
T. H. Sandve, I. Berre, and J. M. Nordbotten, “An efficient multi-point flux approximation method for discrete fracture-matrix simulations,” Journal of Computational Physics, vol. 231, no. 9, pp. 3784–3800, 2012.
View at: Publisher Site | Google Scholar
S. K. Matthai, A. Mezentsev, and M. Belayneh, “Finite element-node-centered finite-volume two-phase-flow experiments with fractured rock represented by unstructured hybrid-element meshes,” SPE Reservoir Evaluation and Engineering, vol. 10, no. 4, pp. 740–756, 2007.
View at: Publisher Site | Google Scholar
O. M. Olorode, C. M. Freeman, and G. J. Moridis, “High-resolution numerical modeling of complex and irregular fracture patterns in shale gas and tight gas reservoirs,” SPE Reservoir Evaluation and Engineering, vol. 16, no. 4, pp. 443–455, 2013.
View at: Publisher Site | Google Scholar
O. Al-Hinai, G. Singh, and G. Pencheva, “Modeling multiphase flow with nonplanar fractures,” in SPE Reservoir Simulation Symposium, Woodlands, Texas, USA, 2013.
View at: Publisher Site | Google Scholar
A. Moinfar, A. Varavei, and K. Sepehrnoori, “Development of an efficient embedded discrete fracture model for 3D compositional reservoir simulation in fractured reservoirs,” SPE Journal, vol. 19, no. 2, pp. 289–303, 2014.
View at: Publisher Site | Google Scholar
L. Li and S. H. Lee, “Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media,” SPE Reservoir Evaluation and Engineering, vol. 11, no. 4, pp. 750–758, 2008.
View at: Publisher Site | Google Scholar
Y. Xu, J. S. Cavalcante Filho, and W. Yu, “Discrete-fracture modeling of complex hydraulic-fracture geometries in reservoir simulators,” SPE Reservoir Evaluation and Engineering, vol. 20, no. 2, pp. 403–422, 2017.
View at: Publisher Site | Google Scholar
W. Yu, Y. Xu, M. Liu, K. Wu, and K. Sepehrnoori, “Simulation of shale gas transport and production with complex fractures using embedded discrete fracture model,” AIChE Journal, vol. 64, no. 6, pp. 2251–2264, 2018.
View at: Publisher Site | Google Scholar
A. Moinfar, A. Varavei, and K. Sepehrnoori, “Development of a coupled dual continuum and discrete fracture model for the simulation of unconventional reservoirs,” in SPE Reservoir Simulation Symposium, Woodlands, Texas, USA, 2013.
View at: Publisher Site | Google Scholar
W. Li, Z. Dong, and G. Lei, “Integrating embedded discrete fracture and dual-porosity, dual-permeability methods to simulate fluid flow in shale oil reservoirs,” Energies, vol. 10, no. 10, p. 1471, 2017.
View at: Publisher Site | Google Scholar
A. Moinfar, K. Sepehrnoori, and R. T. Johns, “Coupled geomechanics and flow simulation for an embedded discrete fracture model,” in SPE Reservoir Simulation Symposium, Woodlands, Texas, USA, 2013.
View at: Publisher Site | Google Scholar
J. Jiang and R. M. Younis, “An improved projection-based embedded discrete fracture model (pEDFM) for multiphase flow in fractured reservoirs,” Advances in Water Resources, vol. 109, pp. 267–289, 2017.
View at: Publisher Site | Google Scholar
M. Ţene, S. B. M. Bosma, and M. S. Al Kobaisi, “Projection-based embedded discrete fracture model (pEDFM),” Advances in Water Resources, vol. 105, pp. 205–216, 2017.
View at: Google Scholar
O. Olorode, B. Wang, and H. U. Rashid, “Three-dimensional projection-based embedded discrete-fracture model for compositional simulation of fractured reservoirs,” SPE Journal, vol. 25, no. 4, pp. 2143–2161, 2020.
View at: Publisher Site | Google Scholar
X. Rao, L. Cheng, R. Cao, P. Jia, H. Liu, and X. Du, “A modified projection-based embedded discrete fracture model (pEDFM) for practical and accurate numerical simulation of fractured reservoir,” Journal of Petroleum Science and Engineering, vol. 187, Article ID 106852, 2020.
View at: Publisher Site | Google Scholar
B. Yan, L. Mi, Z. Chai, Y. Wang, and J. E. Killough, “An enhanced discrete fracture network model for multiphase flow in fractured reservoirs,” Journal of Petroleum Science and Engineering, vol. 161, pp. 667–682, 2018.
View at: Publisher Site | Google Scholar
R. Shao and Y. Di, “An integrally embedded discrete fracture model with a semi-analytic transmissibility calculation method,” Energies, vol. 11, no. 12, p. 3491, 2018.
View at: Publisher Site | Google Scholar
R. Shao, Y. Di, D. Wu, and Y.-S. Wu, “An integrally embedded discrete fracture model for flow simulation in anisotropic formations,” Energies, vol. 13, no. 12, p. 3070, 2020.
View at: Publisher Site | Google Scholar
S. C. Bandis, A. C. Lumsden, and N. R. Barton, “Fundamentals of rock joint deformation,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 20, no. 6, pp. 249–268, 1983.
View at: Publisher Site | Google Scholar
P. Asadollahi and F. Tonon, “Constitutive model for rock fractures: revisiting Barton's empirical model,” Engineering Geology, vol. 113, no. 1-4, pp. 11–32, 2010.
View at: Publisher Site | Google Scholar
R. Olsson and N. Barton, “An improved model for hydromechanical coupling during shearing of rock joints,” International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 3, pp. 317–329, 2001.
View at: Publisher Site | Google Scholar
M. Bai, “On equivalence of dual-porosity poroelastic parameters,” Journal of Geophysical Research: Solid Earth, vol. 104, no. B5, pp. 10461–10466, 1999.
View at: Publisher Site | Google Scholar
J. Rutqvist, Y.-S. Wu, C.-F. Tsang, and G. Bodvarsson, “A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 39, no. 4, pp. 429–442, 2002.
View at: Publisher Site | Google Scholar
M. Karimi-Fard, L. J. Durlofsky, and K. Aziz, “An efficient discrete-fracture model applicable for general-purpose reservoir simulators,” SPE Journal, vol. 9, no. 2, pp. 227–236, 2004.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2022 Dawei Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

369

Downloads

679

Citations