Embedded Discrete Fracture Modeling as a Method to Upscale Permeability for Fractured Reservoirs

: Fractured reservoirs are distributed widely over the world, and describing ﬂuid ﬂow in fractures is an important and challenging topic in research. Discrete fracture modeling (DFM) and equivalent continuum modeling are two principal methods used to model ﬂuid ﬂow through fractured rocks. In this paper, a novel method, embedded discrete fracture modeling (EDFM), is developed to compute equivalent permeability in fractured reservoirs. This paper begins with an introduction on EDFM. Then, the paper describes an upscaling procedure to calculate equivalent permeability. Following this, the paper carries out a series of simulations to compare the computation cost between DFM and EDFM. In addition, the method is veriﬁed by embedded discrete fracture modeling and ﬁne grid methods, and grid-block and multiphase ﬂow are studied to prove the feasibility of the method. Finally, the upscaling procedure is applied to a three-dimensional case in order to study performance for a gas injection problem. This study is the ﬁrst to use embedded discrete fracture modeling to compute equivalent permeability for fractured reservoirs. This paper also provides a detailed comparison and discussion on embedded discrete fracture modeling and discrete fracture modeling in the context of equivalent permeability computation with a single-phase model. Most importantly, this study addresses whether this novel method can be used in multiphase ﬂow in a reservoir with fractures.


Introduction
Fractures, whether naturally occurring or hydraulically created, occur at different length scales with different densities in geological porous media.Fractures have a great impact on the quality of a porous medium.Understanding fractures is important for describing the challenges inherent to hydrocarbon flow in porous media.Two main methods are used to model fluid flow through fractured rocks.One is discrete fracture modeling (DFM), and the other is equivalent continuum modeling.
DFM explicitly describes matrix and fracture systems [1], can be used for modeling fractured reservoirs when the properties of the porous media are heterogeneous [2], and is considered one of the most accurate methods to describe fluid flow through a fracture and matrix system.However, a main drawback of this method is that meshing a complex three-dimensional (3D) fracture network is very difficult.Another drawback is the high computational cost due to a large number of grids and convergence problems.Given the improvement of computational efficiency, though, DFM is becoming more and more attractive; however, its application at the field scale is not realistic at present [3,4].
Energies 2019, 12, 812 3 of 15 single-phase and multiphase flow and study the performance for different grid-block numbers in the upscaling process.We address comparisons between the proposed upscaling approach and the EDFM method used in the oil industry for a realistic 3D case.Section 5 gathers conclusions and presents an outlook on flow-based upscaling for fractured reservoirs.

Embedded Discrete Fracture Modeling
EDFM was introduced by Li and Lee [5] and was expanded on by Moinfar et al. [25] and Tene et al. [26].Unlike DFM, the EDFM method does not require generating complex meshing and can use conventional corner-point grids, which alleviate the computational cost for complex fracture modeling.
In this method, we need to define the corner-point cell for the matrix and then embed fractures into the matrix cells (Figure 1); after that, we define intersection faces between different media, for example, matrix and fracture.Once we have defined all cells and have found the intersections between them, we can proceed to compute transmissibility, which is a key aspect of EDFM. Figure 2 illustrates the connection list of continua in the computational domain for a simple scenario.As mentioned by Li et al. [9], there are four types of connection: Type one: Conductivity between matrix grids Type two: Conductivity between fracture segments inside an individual fracture Type three: Conductivity between intersecting fracture segments Type four: Conductivity between fracture and matrix Figure 2a shows four matrix cells and two fractures, with fracture one divided into segments F1 and F2 by matrix cells and fracture two divided into segments F3 and F4 by matrix cells.Figure 2b shows the connections among different media.Black lines indicate the connection between matrix cells; green lines indicate the connection between matrix and fracture; red lines show the connection between fracture segments in an individual fracture; and the blue line is the connection of intersecting fractures.
Energies 2019, 12, x FOR PEER REVIEW 3 of 15 for single-phase and multiphase flow and study the performance for different grid-block numbers in the upscaling process.We address comparisons between the proposed upscaling approach and the EDFM method used in the oil industry for a realistic 3D case.Section 5 gathers conclusions and presents an outlook on flow-based upscaling for fractured reservoirs.

Embedded Discrete Fracture Modeling
EDFM was introduced by Li and Lee [5] and was expanded on by Moinfar et al. [25] and Tene et al. [26].Unlike DFM, the EDFM method does not require generating complex meshing and can use conventional corner-point grids, which alleviate the computational cost for complex fracture modeling.
In this method, we need to define the corner-point cell for the matrix and then embed fractures into the matrix cells (Figure 1); after that, we define intersection faces between different media, for example, matrix and fracture.Once we have defined all cells and have found the intersections between them, we can proceed to compute transmissibility, which is a key aspect of EDFM. Figure 2 illustrates the connection list of continua in the computational domain for a simple scenario.As mentioned by Li et al. [9], there are four types of connection:   for single-phase and multiphase flow and study the performance for different grid-block numbers in the upscaling process.We address comparisons between the proposed upscaling approach and the EDFM method used in the oil industry for a realistic 3D case.Section 5 gathers conclusions and presents an outlook on flow-based upscaling for fractured reservoirs.

Embedded Discrete Fracture Modeling
EDFM was introduced by Li and Lee [5] and was expanded on by Moinfar et al. [25] and Tene et al. [26].Unlike DFM, the EDFM method does not require generating complex meshing and can use conventional corner-point grids, which alleviate the computational cost for complex fracture modeling.
In this method, we need to define the corner-point cell for the matrix and then embed fractures into the matrix cells (Figure 1); after that, we define intersection faces between different media, for example, matrix and fracture.Once we have defined all cells and have found the intersections between them, we can proceed to compute transmissibility, which is a key aspect of EDFM. Figure 2 illustrates the connection list of continua in the computational domain for a simple scenario.As mentioned by Li et al. [9], there are four types of connection:

Type One: Conductivity between Matrix Grids
The transmissibility factor between matrix cells (the black lines indicated in Figure 2b) in the model depends on the matrix permeability and matrix geometry, which was given by Li et al. [9].Equation (1) refers to the transmissibility factor between matrix cells: where A mij is the contact area between matrix cells i and j. k mij is the average permeability tensor for matrix cells i and j, which can be defined as follows: where k mi , k mj are the permeability tensors for matrix cells i and j, respectively.d mij is the average normal distance between matrix cells i and j.

Type Two: Conductivity between Fracture Segments Inside an Individual Fracture
Karimi-Fard et al. [3] introduced a two-point flux approximation scheme to calculate the transmissibility factor between fracture segments in an individual fracture (F1 and F2, F3 and F4 in Figure 2b).The transmissibility factor was defined as follows: where k f i , k f j are the permeability tensors for fracture segments i and j, respectively.A c is the contact area for two segments.d segi , d segj are the distances from the centroids of fracture segments i and j to the contact face, respectively.

Type Three: Conductivity between Intersecting Fracture Segments
The calculation method for transmissibility between intersecting fracture segments (the blue lines indicated in Figure 2b) was given by Moinfar et al. [25] and is defined as follows: where L int is the length of the intersection line.k f i , k f j are the permeabilities of fractures i and j, respectively.w f i , w f j are the apertures of fractures i and j, respectively.d f i , d f j are the distances from the centroids of the segments i and j to the intersection line, respectively.

Type Four: Conductivity between Fracture and Matrix
The transmissibility factor between the matrix and fracture segment (the green lines indicated in Figure 2b) depends on the matrix permeability and fracture geometry, which was given by Li et al. [9].When a fracture segment fully penetrates a matrix cell, the matrix-fracture transmissibility factor is as follows: where A m f is the area of the fracture segment on one side.k m f is the permeability tensor, which is calculated as follows: where k f is the permeability tensor for the fracture.k m is the permeability tensor for the matrix.d f −m is the average normal distance from the matrix to the fracture.

Upscaling Procedure
As pointed out by Fumagalli [4,24], along with the increase in the number of fractures, the degree of difficulty in generating fracture grids using EDFM/DFM is rising, and the number of grids is increasing, both leading to an increase in computation cost, which limits the application of EDFM/DFM for complex fracture systems.Therefore, different strategies are required to solve this problem, such as equivalent continuum modeling.In equivalent continuum modeling, an upscaling process is needed.Durlofsky [16] gave a comprehensive review of upscaling techniques in porous media without fractures.In this study, our main goal is to upscale grid properties (focusing on permeability) for fractured reservoirs.This section shows the procedure for upscaling based on EDFM.The entire procedure presented in this section is summarized as Figure 3.
During the upscaling process, we first receive the matrix grid-block and its properties from commercial software, such as Petrel [27] or CMG [28].The information usually includes eight-point positions of each grid-block in a corner-point grid system, as well as its porosity, permeability, and saturation.We also obtain fracture information, such as fracture distribution, permeability, and porosity, from commercial software, such as Fracman [29] or the fracture module in Petrel.Using the information gained from the above two steps, we can generate a matrix grid and fracture distribution for a two-dimensional (2D) case (Figure 4).
of difficulty in generating fracture grids using EDFM/DFM is rising, and the number of grids is increasing, both leading to an increase in computation cost, which limits the application of EDFM/DFM for complex fracture systems.Therefore, different strategies are required to solve this problem, such as equivalent continuum modeling.In equivalent continuum modeling, an upscaling process is needed.Durlofsky [16] gave a comprehensive review of upscaling techniques in porous media without fractures.In this study, our main goal is to upscale grid properties (focusing on permeability) for fractured reservoirs.This section shows the procedure for upscaling based on EDFM.The entire procedure presented in this section is summarized as Figure 3.During the upscaling process, we first receive the matrix grid-block and its properties from commercial software, such as Petrel [27] or CMG [28].The information usually includes eight-point positions of each grid-block in a corner-point grid system, as well as its porosity, permeability, and Next, for all cells, we apply the following process: Choose a cell, keep the matrix and fracture geometry and property information, and divide the cell into NX × NY × NZ, which may change for different cases and should be optimized sometimes (as shown in Figure 4a,b).
Define the flow boundary for the new grid system.In this study, we use a constant pressure boundary in the flow direction and a closed boundary for the other two sides, as shown in Figure 5.
Step 1: Initialize the new grid system.Assign properties values for the matrix and fractures, such as porosity, permeability, and saturation, which are inherited from the original big model.We also need to define the connection and calculate the transmissibility among different media using Equations ( 1)- (10).Figure 6 shows an example of connection.We also define the well and production in this step.
Step 2: We calculate the flow rate and apply the Darcy flow equation to calculate the whole grid system's permeability, which combines the matrix and fracture flow capacity.
saturation.We also obtain fracture information, such as fracture distribution, permeability, and porosity, from commercial software, such as Fracman [29] or the fracture module in Petrel.Using the information gained from the above two steps, we can generate a matrix grid and fracture distribution for a two-dimensional (2D) case (Figure 4).
Next, for all cells, we apply the following process: Choose a cell, keep the matrix and fracture geometry and property information, and divide the cell into NX × NY × NZ, which may change for different cases and should be optimized sometimes (as shown in Figure 4a and 4b).
Define the flow boundary for the new grid system.In this study, we use a constant pressure boundary in the flow direction and a closed boundary for the other two sides, as shown in Figure 5.
Step 1: Initialize the new grid system.Assign properties values for the matrix and fractures, such as porosity, permeability, and saturation, which are inherited from the original big model.We also need to define the connection and calculate the transmissibility among different media using Equations ( 1)- (10).Figure 6 shows an example of connection.We also define the well and production in this step.
Step 2: We calculate the flow rate and apply the Darcy flow equation to calculate the whole grid system's permeability, which combines the matrix and fracture flow capacity.

Numerical Examples
To evaluate and validate the performance of the new method introduced in Section 3, a series of simulation cases are presented in this section.In Section 4.1, we validate the EDFM by comparison with the DFM model.In Section 4.2, we compare the computation time and performance among the DFN method, the EDFM method, and the proposed upscaling approach for a sample of fractures.Then, in Section 4.3, we validate the new upscaling approach.A series of simulation cases is studied for the grid block and multiphase in Section 4.4.In Section 4.5, we compare the proposed approach and the EDFM method for a 3D case with complex fractures.

Validating the EDFM Method
As mentioned in the introduction section, DFM is the classic method used for upscaling fracture permeability.In this section, computation time and gas rate for EDFM and DFM are compared by testing a model with 251 fractures (Figure 7a).The matrix grid number is 5000 (50 × 50 × 2) in the EDFM model (Figure 7b), while the DFM model includes 238,762 cells (Figure 7c).In the model, matrix permeability is set to 1 mD, and fracture permeability is 5000 mD.A gas well in the center of the model with a constant rate is assumed in the study.Figure 8 presents that the two methods obtain the same gas rate.The computation time for the DFM is 287 seconds and for the EDFM is 73 seconds based on the same computation environment for single-phase flow.These results indicate that EDFM performs better than the DFM method, as one of EDFM's main advantages is that it avoids the need to generate complex grids for fractures and the matrix.

Numerical Examples
To evaluate and validate the performance of the new method introduced in Section 3, a series of simulation cases are presented in this section.In Section 4.1, we validate the EDFM by comparison with the DFM model.In Section 4.2, we compare the computation time and performance among the DFN method, the EDFM method, and the proposed upscaling approach for a sample of fractures.Then, in Section 4.3, we validate the new upscaling approach.A series of simulation cases is studied for the grid block and multiphase in Section 4.4.In Section 4.5, we compare the proposed approach and the EDFM method for a 3D case with complex fractures.

Validating the EDFM Method
As mentioned in the introduction section, DFM is the classic method used for upscaling fracture permeability.In this section, computation time and gas rate for EDFM and DFM are compared by testing a model with 251 fractures (Figure 7a).The matrix grid number is 5000 (50 × 50 × 2) in the EDFM model (Figure 7b), while the DFM model includes 238,762 cells (Figure 7c).In the model, matrix permeability is set to 1 mD, and fracture permeability is 5000 mD.A gas well in the center of the model with a constant rate is assumed in the study.Figure 8 presents that the two methods obtain the same gas rate.The computation time for the DFM is 287 seconds and for the EDFM is 73 seconds based on the same computation environment for single-phase flow.These results indicate that EDFM performs better than the DFM method, as one of EDFM's main advantages is that it avoids the need to generate complex grids for fractures and the matrix.
As mentioned in the introduction section, DFM is the classic method used for upscaling fracture permeability.In this section, computation time and gas rate for EDFM and DFM are compared by testing a model with 251 fractures (Figure 7a).The matrix grid number is 5000 (50 × 50 × 2) in the EDFM model (Figure 7b), while the DFM model includes 238,762 cells (Figure 7c).In the model, matrix permeability is set to 1 mD, and fracture permeability is 5000 mD.A gas well in the center of the model with a constant rate is assumed in the study.Figure 8 presents that the two methods obtain the same gas rate.The computation time for the DFM is 287 seconds and for the EDFM is 73 seconds based on the same computation environment for single-phase flow.These results indicate that EDFM performs better than the DFM method, as one of EDFM's main advantages is that it avoids the need to generate complex grids for fractures and the matrix.

Validating the New Upscaling Approach by Comparison with the DFM Method
As shown in Figure 9, a cell model with two fractures is used to compare the EDFM and DFN model for upscaling.The EDFM includes 400 cells for the matrix and six fractures, while the fine grid model includes 862 cells.Fracture permeability is set to 5000 mD with an aperture of 0.001 ft, and matrix permeability is 0.1 mD for both models.A constant pressure boundary is used in the flow direction, and closed boundaries are set for the other two sides.The upscaling permeabilities calculated from the EDFM and the fine grid model in the X direction are 0.1235 mD and 0.1227 mD, respectively, and in the Y direction are 0.1236 mD and 0.1228 mD, respectively.The gap between the two methods is less than 1%.The computation time for the EDFM and the fine grid model is 4.2 s mD and 9.8 mD, respectively.The results indicate that using the EDFM method for upscaling can speed up the upscaling process.

Validating the New Upscaling Approach by Comparison with the DFM Method
As shown in Figure 9, a cell model with two fractures is used to compare the EDFM and DFN model for upscaling.The EDFM includes 400 cells for the matrix and six fractures, while the fine grid model includes 862 cells.Fracture permeability is set to 5000 mD with an aperture of 0.001 ft, and matrix permeability is 0.1 mD for both models.A constant pressure boundary is used in the flow direction, and closed boundaries are set for the other two sides.The upscaling permeabilities calculated from the EDFM and the fine grid model in the X direction are 0.1235 mD and 0.1227 mD, respectively, and in the Y direction are 0.1236 mD and 0.1228 mD, respectively.The gap between the two methods is less than 1%.The computation time for the EDFM and the fine grid model is 4.2 s mD and 9.8 mD, respectively.The results indicate that using the EDFM method for upscaling can speed up the upscaling process.
direction, and closed boundaries are set for the other two sides.The upscaling permeabilities calculated from the EDFM and the fine grid model in the X direction are 0.1235 mD and 0.1227 mD, respectively, and in the Y direction are 0.1236 mD and 0.1228 mD, respectively.The gap between the two methods is less than 1%.The computation time for the EDFM and the fine grid model is 4.2 s mD and 9.8 mD, respectively.The results indicate that using the EDFM method for upscaling can speed up the upscaling process.

Comparison with Other Models for Flow Results
In this section, the simulation result of the new method is compared with results from other methods, such as fine grid and EDFM.In the first case, only two cross-fractures are considered.The fracture distribution, fine model, and EDFM model are shown in Figure 10.The permeability of the matrix is 1 mD, while the fracture permeability is 5000 mD with an aperture of 0.001 ft.In the upscaling, we divide each cell into 10 × 10 cells and use a constant pressure boundary.The permeability distribution after upscaling is shown in Figure 11.

Comparison with Other Models for Flow Results
In this section, the simulation result of the new method is compared with results from other methods, such as fine grid and EDFM.In the first case, only two cross-fractures are considered.The fracture distribution, fine model, and EDFM model are shown in Figure 10.The permeability of the matrix is 1 mD, while the fracture permeability is 5000 mD with an aperture of 0.001 ft.In the upscaling, we divide each cell into 10 × 10 cells and use a constant pressure boundary.The permeability distribution after upscaling is shown in Figure 11.To validate the method, we compare the production performance among different approaches.We define a well in the center of the model, and the well intersects with fractures.As shown in Figure 12, although the middle part of the cumulative production curves has a slight difference, the three methods obtain similar performance for rate and cumulative production.Figure 13 shows the pressure profile of the different methods.The results of EDFM and the new upscaling method are highly consistent in these plots, and the fine grid model is slightly different among these methods.To validate the method, we compare the production performance among different approaches.We define a well in the center of the model, and the well intersects with fractures.As shown in Figure 12, although the middle part of the cumulative production curves has a slight difference, the three methods obtain similar performance for rate and cumulative production.Figure 13 shows the pressure profile of the different methods.The results of EDFM and the new upscaling method are highly consistent in these plots, and the fine grid model is slightly different among these methods.
120 1600 To validate the method, we compare the production performance among different approaches.We define a well in the center of the model, and the well intersects with fractures.As shown in Figure 12, although the middle part of the cumulative production curves has a slight difference, the three methods obtain similar performance for rate and cumulative production.Figure 13 shows the pressure profile of the different methods.The results of EDFM and the new upscaling method are highly consistent in these plots, and the fine grid model is slightly different among these methods.
permeability distribution in the X direction, Kx, and (b) shows the permeability distribution in the Y direction, Ky.
To validate the method, we compare the production performance among different approaches.We define a well in the center of the model, and the well intersects with fractures.As shown in Figure 12, although the middle part of the cumulative production curves has a slight difference, the three methods obtain similar performance for rate and cumulative production.Figure 13 shows the pressure profile of the different methods.The results of EDFM and the new upscaling method are highly consistent in these plots, and the fine grid model is slightly different among these methods.For the second case, the fracture system in Figure 14a is used for investigation.In this model, fracture permeability is 5000 mD, matrix permeability is 1 mD, and a well that is used to produce liquid at a constant liquid rate is located in the center.Figure 14b,c show the upscaled permeability for the X and Y directions.We find that the distribution of permeability is consistent with that of the fracture.Figure 15 presents the simulation results of oil rate and cumulative oil at the field scale.We observe that the results obtained by the proposed upscaling method and EDFM are close to each other in this study.Figure 16 also indicates that the pressure distribution obtained by the two methods is in good agreement.For the second case, the fracture system in Figure 14a is used for investigation.In this model, fracture permeability is 5000 mD, matrix permeability is 1 mD, and a well that is used to produce liquid at a constant liquid rate is located in the center.Figure 14b,c show the upscaled permeability for the X and Y directions.We find that the distribution of permeability is consistent with that of the fracture.Figure 15 presents the simulation results of oil rate and cumulative oil at the field scale.We observe that the results obtained by the proposed upscaling method and EDFM are close to each other in this study.Figure 16 also indicates that the pressure distribution obtained by the two methods is in good agreement.For the second case, the fracture system in Figure 14a is used for investigation.In this model, fracture permeability is 5000 mD, matrix permeability is 1 mD, and a well that is used to produce liquid at a constant liquid rate is located in the center.Figure 14b,c show the upscaled permeability for the X and Y directions.We find that the distribution of permeability is consistent with that of the fracture.Figure 15 presents the simulation results of oil rate and cumulative oil at the field scale.We observe that the results obtained by the proposed upscaling method and EDFM are close to each other in this study.Figure 16 also indicates that the pressure distribution obtained by the two methods is in good agreement.

Sensitivity Study for the New Approach
Grid-block: As mentioned in Section 3, in the upscaling procedure, we need to divide each objective cell into small cells and then use EDFM to model flow rate and calculate permeability.The key point here is the number of cells required to obtain an accurate result.The test domain is represented in Figure 17a, where the fractures and the 11 × 11 matrix grid-blocks are represented.The permeability of the matrix is 1 mD, while the fracture permeability is 5000 mD with an aperture of 0.001 ft.Each cell (for example, Figure 17b) would be divided into four scenarios (5 × 5, 10 × 10, 15 × 15, 20 × 20), which are shown in Figure 17c-f.Figure 18a shows the upscaling permeability in cell (5,3) obtained by different scenarios.The result indicates that the computed permeability converges toward a single value.Figure 18b shows upscaling permeability for each cell obtained by different scenarios.The result represents that the slopes of curves (20 × 20 versus n × n) tend to be one as the grid number (n) increases, which means that the calculated permeability tends to be uniform for each cell with an increase in meshes in the upscaling process.The corresponding computation costs for the four scenarios

Sensitivity Study for the New Approach
Grid-block: As mentioned in Section 3, in the upscaling procedure, we need to divide each objective cell into small cells and then use EDFM to model flow rate and calculate permeability.
The key point here is the number of cells required to obtain an accurate result.The test domain is represented in Figure 17a, where the fractures and the 11 × 11 matrix grid-blocks are represented.The permeability of the matrix is 1 mD, while the fracture permeability is 5000 mD with an aperture of 0.001 ft.Each cell (for example, Figure 17b) would be divided into four scenarios (5 × 5, 10 × 10, 15 × 15, 20 × 20), which are shown in Figure 17c-f.Figure 18a shows the upscaling permeability in cell (5,3) obtained by different scenarios.The result indicates that the computed permeability converges toward a single value.Figure 18b shows upscaling permeability for each cell obtained by different scenarios.The result represents that the slopes of curves (20 × 20 versus n × n) tend to be one as the grid number (n) increases, which means that the calculated permeability tends to be uniform for each cell with an increase in meshes in the upscaling process.The corresponding computation costs for the four scenarios (5 × 5, 10 × 10, 15 × 15, 20 × 20) are 311 s, 320 s, 340 s, 360 s, respectively.
(5,3) obtained by different scenarios.The result indicates that the computed permeability converges toward a single value.Figure 18b shows upscaling permeability for each cell obtained by different scenarios.The result represents that the slopes of curves (20 × 20 versus n × n) tend to be one as the grid number (n) increases, which means that the calculated permeability tends to be uniform for each cell with an increase in meshes in the upscaling process.The corresponding computation costs for the four scenarios (5 × 5, 10 × 10, 15 × 15, 20 × 20) are 311 s, 320 s, 340 s, 360 s, respectively.Multiphase flow: We consider multiphase flow in the model based on the model shown in Figure 13.Water is pumped at a constant rate into the injection well, whereas the production well is operating at a constant pressure.After upscaling, fracture permeability is average to the matrix, and equivalent permeability is less than 1/10 or even one percent of fracture permeability.For example, fracture permeability is 5000 mD, matrix permeability is 0.01 mD, and equivalent permeability is 0.023 mD after upscaling.Injection fluid (water) flow through fractures is faster that through the matrix with equivalent permeability.As shown in Figure 19a, the upscaled solution shows a delay of water breakthrough (water cut) compared to the EDFM method.The water and recovery obtained from the two methods are close to each other for multiphase flow over the long term in this study (Figure 19b).The global difference for water cut and recovery between two methods is 1.87% and 0.86%, respectively.Multiphase flow: We consider multiphase flow in the model based on the model shown in Figure 13.Water is pumped at a constant rate into the injection well, whereas the production well is operating at a constant pressure.After upscaling, fracture permeability is average to the matrix, and equivalent permeability is less than 1/10 or even one percent of fracture permeability.For example, fracture permeability is 5000 mD, matrix permeability is 0.01 mD, and equivalent permeability is 0.023 mD after upscaling.Injection fluid (water) flow through fractures is faster that through the matrix with equivalent permeability.As shown in Figure 19a, the upscaled solution shows a delay of water breakthrough (water cut) compared to the EDFM method.The water and recovery obtained from the two methods are close to each other for multiphase flow over the long term in this study (Figure 19b).The global difference for water cut and recovery between two methods is 1.87% and 0.86%, respectively.after upscaling.Injection fluid (water) flow through fractures is faster that through the matrix with equivalent permeability.As shown in Figure 19a, the upscaled solution shows a delay of water breakthrough (water cut) compared to the EDFM method.The water and recovery obtained from the two methods are close to each other for multiphase flow over the long term in this study (Figure 19b).The global difference for water cut and recovery between two methods is 1.87% and 0.86%, respectively.

Comparison Using a Complex Fracture System
In this section, the objective is to test the robustness and the efficiency of the new approach in a complex fracture model.A gas injection simulation is performed on a complex fracture reservoir, and simulation results obtained from the new upscaling procedure and EDFM are compared.
The fracture system in Figure 20 is employed in the study.The model includes 41 × 41 × 7 matrix cells, along with 309 fractures.The matrix porosity is 0.10, and the matrix and fracture permeabilities are 0.1 mD and 5000 mD, respectively.Gas is pumped at a constant rate into the injection well, whereas the nine production wells are operating at a constant pressure.Figure 21 presents the simulation results of the gas-oil ratio and oil recovery for each producer.The gas-oil ratio curve shows a slight difference at the end of injection, and the recovery curve shows good consistency between the two methods.A slight difference gas-oil ratio is caused by the different relative

Comparison Using a Complex Fracture System
In this section, the objective is to test the robustness and the efficiency of the new approach in a complex fracture model.A gas injection simulation is performed on a complex fracture reservoir, and simulation results obtained from the new upscaling procedure and EDFM are compared.
The fracture system in Figure 20 is employed in the study.The model includes 41 × 41 × 7 matrix cells, along with 309 fractures.The matrix porosity is 0.10, and the matrix and fracture permeabilities are 0.1 mD and 5000 mD, respectively.Gas is pumped at a constant rate into the injection well, whereas the nine production wells are operating at a constant pressure.Figure 21 presents the simulation results of the gas-oil ratio and oil recovery for each producer.The gas-oil ratio curve shows a slight difference at the end of injection, and the recovery curve shows good consistency between the two methods.A slight difference gas-oil ratio is caused by the different relative permeability between fracture and matrix.However, the computation time for the new approach and for EDFM are 175 s, and 2560 s, respectively.permeability between fracture and matrix.However, the computation time for the new approach and for EDFM are 175 s, and 2560 s, respectively.

Conclusions
In this work, we implement a novel upscaling simulation method that uses EDFM to calculate equivalent permeability for a complex fracture network.The main objective of this work was to introduce the procedure and study the feasibility of this novel method.The new upscaling approach permeability between fracture and matrix.However, the computation time for the new approach and for EDFM are 175 s, and 2560 s, respectively.

Conclusions
In this work, we implement a novel upscaling simulation method that uses EDFM to calculate equivalent permeability for a complex fracture network.The main objective of this work was to introduce the procedure and study the feasibility of this novel method.The new upscaling approach

Figure. 1 AFigure 2 .
Figure2ashows four matrix cells and two fractures, with fracture one divided into segments F1 and F2 by matrix cells and fracture two divided into segments F3 and F4 by matrix cells.Figure2bshows the connections among different media.Black lines indicate the connection between matrix cells; green lines indicate the connection between matrix and fracture; red lines show the connection between fracture segments in an individual fracture; and the blue line is the connection of intersecting fractures.
Figure2ashows four matrix cells and two fractures, with fracture one divided into segments F1 and F2 by matrix cells and fracture two divided into segments F3 and F4 by matrix cells.Figure2bshows the connections among different media.Black lines indicate the connection between matrix cells; green lines indicate the connection between matrix and fracture; red lines show the connection between fracture segments in an individual fracture; and the blue line is the connection of intersecting fractures.

Figure 2 .
Figure 2.An example of connection lists for EDFM.(a) shows the matrix grid and fracture distribution.(b) is the connection list.Reproduced from [Li et al.], publisher: 2017.

Figure 3 .
Figure 3.The fracture permeability upscaling procedure for fractured reservoirs.

Figure 3 .
Figure 3.The fracture permeability upscaling procedure for fractured reservoirs.

Figure 4 .
Figure 4.The matrix grid and fracture distribution for a two-dimensional (2D) case.(a) is a matrix system of the corner-point grid (5 × 8), along with nine fractures.(b) shows an arbitrary grid from the left figure.The bottom-right figure shows the grid system for both the matrix and fractures after meshing the grid to 5 × 7.

Figure 4 .
Figure 4.The matrix grid and fracture distribution for a two-dimensional (2D) case.(a) is a matrix system of the corner-point grid (5 × 8), along with nine fractures.(b) shows an arbitrary grid from the left figure.The bottom-right figure shows the grid system for both the matrix and fractures after meshing the grid to 5 × 7.

Figure 4 .
Figure 4.The matrix grid and fracture distribution for a two-dimensional (2D) case.(a) is a matrix system of the corner-point grid (5 × 8), along with nine fractures.(b) shows an arbitrary grid from the left figure.The bottom-right figure shows the grid system for both the matrix and fractures after meshing the grid to 5 × 7.

Figure 5 .
Figure 5.The boundary condition for a new grid system shown in Figure 4b.Blue lines are fractures, green lines represent the constant pressure boundary, and thick black lines indicate closed (no flow) boundaries.In this case, the flow direction is from right to left or from left to right.

Figure 5 .Figure 6 .
Figure 5.The boundary condition for a new grid system shown in Figure 4b.Blue lines are fractures, green lines represent the constant pressure boundary, and thick black lines indicate closed (no flow) boundaries.In this case, the flow direction is from right to left or from left to right.Energies 2019, 12, x FOR PEER REVIEW 7 of 15

Figure 6 .
Figure 6.Connection information for the EDFM method for the cells shown in the red block of Figure 5. (a) shows four matrix cells with two fractures, with fracture one divided into segments F1 and F4 by matrix cells and fracture two divided into segments F2, F3, and F5 by matrix cells.(b) shows the connection among different media.Line colors represent the same elements as in Figure 2.

Figure 7 .
Figure 7. Fracture distribution for the EDFM and DFM model used in the study.(a) shows the fracture distribution, including 250 fractures with 5000-mD permeability and 0.01-ft aperture.(b) is the EDFM model, which presents fracture explicitly.(c) is the DFM model, which includes 238,762 cells.

Figure 7 .Figure 8 .
Figure 7. Fracture distribution for the EDFM and DFM model used in the study.(a) shows the fracture distribution, including 250 fractures with 5000-mD permeability and 0.01-ft aperture.(b) is the EDFM model, which presents fracture explicitly.(c) is the DFM model, which includes 238,762 cells.Energies 2019, 12, x FOR PEER REVIEW 8 of 15

Figure 9 .
Figure 9.A conceptual model of fracture distribution and the EDFM and DFN model for

Figure 8 .
Figure 8.A comparison of gas rate (a) and cumulative gas (b) between EDFM and DFM.

Figure 9 .
Figure 9.A conceptual model of fracture distribution and the EDFM and DFN model for one cell.(a) is the fracture distribution, including six fractures with 1000-mD permeability and 0.01-ft aperture.(b) is the EDFM model, and (c) is the DFN model.

Figure 9 .
Figure 9.A conceptual model of fracture distribution and the EDFM and DFN model for one cell.(a) is the fracture distribution, including six fractures with 1000-mD permeability and 0.01-ft aperture.(b) is the EDFM model, and (c) is the DFN model.

Figure 10 .Figure 11 .
Figure 10.A simple model with two cross-fractures used to validate the new approach.(a) shows the fractures, (b) is the fine grid model, which is used to represent fractures, and (c) is the EDFM.

Figure 10 .Figure 10 .Figure 11 .
Figure 10.A simple model with two cross-fractures used to validate the new approach.(a) shows the fractures, (b) is the fine grid model, which is used to represent fractures, and (c) is the EDFM.

Figure 11 .
Figure 11.The permeability upscaling result based on the new method.(a) shows the permeability distribution in the X direction, Kx, and (b) shows the permeability distribution in the Y direction, Ky.

Figure 12 .
Figure 12.A comparison of oil rate (a) and cumulative oil (b) for the new upscaling method, fine grid model, and EDFM model in a single-phase flow case.

Figure 12 .Figure 13 .
Figure 12.A comparison of oil rate (a) and cumulative oil (b) for the new upscaling method, fine grid model, and EDFM model in a single-phase flow case.Energies 2019, 12, x FOR PEER REVIEW 10 of 15

Figure 14 .
Figure 14.The fracture distribution and oil saturation profiles.(a) shows the fracture distribution for the model used in the study.(b) and (c) present upscaled permeability for the X and Y directions, respectively.The yellow lines in the figures represent fractures.

Figure 13 .
Figure 13.A comparison of the pressure profile at the same time level for the new upscaling method, EDFM, and fine grid model in one-phase flow.The pressure figures for (a), (b), and (c) are the new upscaling method, EDFM, and fine grid model, respectively.

Figure 13 .
Figure 13.A comparison of the pressure profile at the same time level for the new upscaling method, EDFM, and fine grid model in one-phase flow.The pressure figures for (a), (b) , and (c) are the new upscaling method, EDFM, and fine grid model, respectively.

Figure 14 .
Figure 14.The fracture distribution and oil saturation profiles.(a) shows the fracture distribution for the model used in the study.(b) and (c) present upscaled permeability for the X and Y directions, respectively.The yellow lines in the figures represent fractures.

Figure 14 .
Figure 14.The fracture distribution and oil saturation profiles.(a) shows the fracture distribution for the model used in the study.(b) and (c) present upscaled permeability for the X and Y directions, respectively.The yellow lines in the figures represent fractures.

Figure 14 .Figure 15 .
Figure 14.The fracture distribution and oil saturation profiles.(a) shows the fracture distribution for the model used in the study.(b) and (c) present upscaled permeability for the X and Y directions, respectively.The yellow lines in the figures represent fractures.

Figure 16 .
Figure 16.Pressure profiles for different methods.(a) shows the pressure distribution for the new upscaling method model, and (b) and (c) show the matrix and fracture pressure distribution for EDFM.The pressure profile for the two methods is at the same time level.The blue line in the center of the model represents the producer.

Figure 17 .
Figure 17.The test domain with fractures and grid-blocks.(a) shows the domain with 14 fractures and an 11 × 11 grid-block.(b) represents a cell (5,3) with four fractures from the test domain, which would be subdivided into smaller cells.The bottom figures show four scenarios :(c) uses 5 × 5 grid-block to model one cell shown in figure (b), (d) 10 × 10 gridblock to model one cell shown in figure (b), (e) 15 × 15 grid-block to model one cell shown in figure (b), (f) 20 × 20 grid-block to model one cell shown in figure (b).Note that serial small cells are added to define the boundary condition in the study for each scenario.

Figure 17 .
Figure 17.The test domain with fractures and grid-blocks.(a) shows the domain with 14 fractures and an 11 × 11 grid-block.(b) represents a cell (5,3) with four fractures from the test domain, which would be subdivided into smaller cells.The bottom figures show four scenarios: (c) uses 5 × 5 grid-block to model one cell shown in figure (b), (d) 10 × 10 grid-block to model one cell shown in figure (b), (e) 15 × 15 grid-block to model one cell shown in figure (b), (f) 20 × 20 grid-block to model one cell shown in figure (b).Note that serial small cells are added to define the boundary condition in the study for each scenario.

Figure 18 .
Figure 18.A comparison of upscaling permeability with different grid-block numbers.(a) shows the upscaling permeability with different grid-block numbers for cell (5,3).(b) represents the upscaling permeability with different grid-block numbers for all cells.

Figure 18 .
Figure 18.A comparison of upscaling permeability with different grid-block numbers.(a) shows the upscaling permeability with different grid-block numbers for cell (5,3).(b) represents the upscaling permeability with different grid-block numbers for all cells.

Figure 19 .
Figure 19.A comparison of water cut (a) and recovery (b) for different methods in multiphase flow.

Figure 19 .
Figure 19.A comparison of water cut (a) and recovery (b) for different methods in multiphase flow.

Figure 20 .Figure 21 .
Figure 20.The fracture distribution and EDFM model used in the study.(a) is the fracture distribution, including 309 fractures with 5000-mD permeability and 0.01 ft aperture.(b) is the EDFM model, which presents fracture explicitly.

Figure 20 .
Figure 20.The fracture distribution and EDFM model used in the study.(a) is the fracture distribution, including 309 fractures with 5000-mD permeability and 0.01 ft aperture.(b) is the EDFM model, which presents fracture explicitly.

Figure 20 .Figure 21 .
Figure 20.The fracture distribution and EDFM model used in the study.(a) is the fracture distribution, including 309 fractures with 5000-mD permeability and 0.01 ft aperture.(b) is the EDFM model, which presents fracture explicitly.

Figure 21 .
Figure 21.The simulation results of the gas-oil ratio (a) and oil recovery (b) for each producer.