An Integrally Embedded Discrete Fracture Model with a Semi-Analytic Transmissibility Calculation Method

The embedded discrete fracture model (EDFM) combines the advantages of previous numerical models for fractured reservoirs, achieving a good balance between calculation cost and simulation accuracy. In this work, an integrally embedded discrete fracture model (iEDFM) is introduced to further improve the simulation accuracy and expand the application of the model. The iEDFM has a new gridding method that can arbitrarily grid the fractures according to the requirements rather than finely subdividing fracture elements. Then, with a more precise pressure distribution assumption inside the matrix blocks, we are able to obtain a semi-analytic calculation method of matrix-fracture transmissibility applied to iEDFM. Several case studies were conducted to demonstrate the advantage of iEDFM and its applicability for intersecting and nonplanar fractured reservoirs, and a 3D case with a modified dataset from a reported seismic survey could be used to demonstrate the potential application of the iEDFM in real field studies.


Introduction
Fractured reservoirs are commonly found all over the world.In many geoscience applications, such as petroleum extraction, the target formations are fractured [1,2].In these formations, matrix rock is crossed by several fractures at multiple length scales, behaving as hydraulic conductors.
In order to evaluate the economic feasibility and to manage production, numerical simulation tools have to be used.However, when modeling flow in fractured reservoirs, the high heterogeneity caused by complex fractures and complicated matrix-fracture fluid exchange will cause inefficiency and inaccuracy [3][4][5].
Research in this area has been advanced significantly over the past several decades.The dual porosity model (DPM) [6][7][8] is one of the earliest methods for modelling fractured systems, and is still widely used in the petrol industry because of its simplicity and practicability.Since then, many other methods based on a similar concept have been developed to expand the application of DPM, such as the dual-porosity dual-permeability (DPDK) model [9,10], multiple interacting continua (MINC) model [11,12], subdomain model [13], triple-porosity dual-permeability (TPDK) model [14,15], and multi-porosity model [16].These multi-continuum methods provide an efficient approach to simulate micro scale fractures.However, the assumption of fracture uniformity is limited due to losing detailed information (such as geometry and location) of the discovered macro-scale factures.
A more accurate and physics-based approach was proposed to discretize the discovered macro-scale fractures explicitly, which is called the discrete fracture model (DFM).Unstructured grids are commonly used in DFM [17][18][19][20][21][22].Based on DFM, some complex fractures, such as intersecting fractures and nonplanar fractures, can be represented through appropriate gridding [22][23][24][25].However, the use of unstructured fine grids in real field studies is still limited because of its complexity in gridding and computational cost [26,27].
The embedded discrete fracture model (EDFM) has been developed as a compromise.EDFM borrows the dual-porosity concept from DPM, using traditional Cartesian gridding for the matrix to keep the efficiency, but also incorporates the effect of each fracture explicitly as with DFM to account for the complexity and heterogeneity of reservoirs.The concept was first proposed by Lee [28].The fracture element within its parent matrix block (element) is represented by a control volume and is connected to the parent matrix block and other fracture element.This concept was implemented by Li and Lee [29] to vertical fracture cases and implemented by Moinfar [26] to non-vertical fracture cases, in which the fractures have arbitrary dip and strike angles.Xu [30] and Yu [31] implemented EDFM to some more complex fracture-networks, such as nonplanar shape and variable aperture.Li [32] combined EDFM with DPDK for reservoirs with different scale fractures.
However, all these EDFMs are based on a simplified matrix-fracture fluid exchange assumption, which leads to inaccuracy in some cases and also limits the application of the model, as it is necessary to finely subdivide fracture elements.Matei [33] proposed a projection-based EDFM (pEDFM) where the fracture and matrix grids are independently defined.The pEDFM is proposed to deal with highly conductive fractures and flow barriers, thus it makes little improvement on computational efficiency, but still provides a useful method to improve gridding and transmissibility calculation upon the classic EDFM.
In this work, an integrally embedded discrete fracture model (iEDFM) is introduced to improve simulation accuracy and expand the applications of the model.The iEDFM has a new gridding method that can arbitrarily grid the fractures according to the requirements, and then embed them integrally in matrix blocks.A semi-analytic matrix-fracture transmissibility calculation method is applied to iEDFM with a more precise pressure distribution assumption inside the matrix blocks.
This paper is organized as follows.First of all, the basic mathematical method is introduced.Second, the improvements of iEDFM upon EDFM are described, where the gridding method and the semi-analytic calculation method of transmissibility are the most important.Subsequently, we demonstrate the applicability and advantages of iEDFM through a single-phase case and two flooding tests.Finally, by using a modified 3D case with a reported dataset from the seismic survey, we are able to demonstrate the potential applications of the iEDFM in real field studies.

Basic Mathematical Method
Models in this paper (including the fine-grid model, EDFM and iEDFM) have been implemented into a multiphase multidimensional black-oil reservoir simulator, named MSFLOW [34].The basic mathematical model of this multiphase multidimensional black-oil method is introduced as follows.
In an isothermal system containing three mass components, three mass-balance equations are needed to describe flow in fracture elements and matrix blocks.For the flow of phase β (β = g for gas, β = w for water, and β = o for oil), the mass-balance equation is given by: where the velocity of phase β is defined by Darcy's law: where S β is the saturation of phase β; k is the absolute permeability of the formation; k r β is relative permeability to phase β; µ β is the viscosity of phase β; P β is the pressure of phase β; ρ β is the density of phase β under reservoir conditions; g is gravitational acceleration; φ is the effective porosity; q β is the sink/source term of phase β; and D is depth from a reference datum.As implemented numerically, Equation ( 1) is discretized in space with an integral finite-difference or control-volume scheme for the fracture elements and matrix blocks.And in time, it's discretized with a backward, first-order, finite-difference scheme.The discrete equations are as follows: where n is the previous time level; n + 1 is the current time level; t is time step size; V i is the volume of element i (matrix or fracture element); η i contains the set of neighboring elements j (matrix or fracture element) to which element i is directly connected; F β,i j is the flow term for phase β between element i and j; and Q βi is the sink/source term at element i of Phase β.
The flow term F β,i j in Equation ( 3) is described by a discrete version of Darcy's law, given by: where λ β,ij+1/2 is the mobility term to phase β, defined as: where ij + 1/2 is a proper averaging or weighting of properties at the interface between element i and j.The flow potential term is defined as: where D i is the depth from a reference datum to the center of element i, and T ij is transmissibility.
If the integral finite-difference scheme [14] is used, the transmissibility will be calculated as: where A ij is the common interface area between elements i and j; d i is the distance from the center of element i to the interface; and k i j+1/2 is an averaged (such as harmonically weighted) absolute permeability along the connection between elements i and j.

Embedded Discrete Fracture Model
EDFM creates fracture elements connected with corresponding matrix blocks (each represents a matrix element) to account for the mass transfer between matrix and fractures.Once a fracture penetrates a matrix block, an additional element is created to represent the fracture segment in the physical domain (Figure 1a).Each individual fracture is discretized into several fracture elements by the fracture intersections (Figure 1b) and the matrix block boundaries (Figure 1c).Thus, there exist three kinds of connections: matrix-matrix (M-M), fracture-fracture (F-F), and matrix-fracture (M-F).The transmissibility of each connection can be calculated referring to Equation (7).
The parameters for the M-M connection give clear physical meanings, and so the transmissibility can be easily obtained.For the F-F connection, a simplified approximation from Karimi-Fard [23] is used.The two-point flux approximation scheme is: where F k is the absolute permeability of fracture, c A is the common interface area for these two fracture elements, and 1 F d and 2 F d are the average distances from fracture elements 1 and 2 to the common interface.
Transmissibility of M-F depends on the matrix permeability and fracture geometry.When a fracture segment fully penetrates a matrix block, EDFM assumes a uniform pressure gradient in the matrix element, and that pressure gradient is normal to the fracture plane, creating a linear pressure distribution assumption.Then, the M-F transmissibility referring to the equation: where F A is the area of the fracture element on one side, M k is the absolute permeability of matrix (when using the harmonically weighted average permeability, the huge fracture permeability can be ignored), and F M d − is the average normal distance from matrix to fracture, which is calculated as: where V is the volume of the matrix element, dV is the volume element of matrix, and n x is the distance from the volume element to the fracture plane.If the fracture does not fully penetrate the matrix element, most of the EDFMs make the same assumption as Li [29] that the transmissibility is proportional to the area of the fracture element inside the matrix element, which actually further simplifies the previous linear pressure distribution assumption.Thus, there exist three kinds of connections: matrix-matrix (M-M), fracture-fracture (F-F), and matrix-fracture (M-F).The transmissibility of each connection can be calculated referring to Equation (7).
The parameters for the M-M connection give clear physical meanings, and so the transmissibility can be easily obtained.For the F-F connection, a simplified approximation from Karimi-Fard [23] is used.The two-point flux approximation scheme is: where k F is the absolute permeability of fracture, A c is the common interface area for these two fracture elements, and d F 1 and d F 2 are the average distances from fracture elements 1 and 2 to the common interface.Transmissibility of M-F depends on the matrix permeability and fracture geometry.When a fracture segment fully penetrates a matrix block, EDFM assumes a uniform pressure gradient in the matrix element, and that pressure gradient is normal to the fracture plane, creating a linear pressure distribution assumption.Then, the M-F transmissibility referring to the equation: where A F is the area of the fracture element on one side, k M is the absolute permeability of matrix (when using the harmonically weighted average permeability, the huge fracture permeability can be ignored), and d M−F is the average normal distance from matrix to fracture, which is calculated as: where V is the volume of the matrix element, dV is the volume element of matrix, and x n is the distance from the volume element to the fracture plane.If the fracture does not fully penetrate the matrix element, most of the EDFMs make the same assumption as Li [29] that the transmissibility is proportional to the area of the fracture element inside the matrix element, which actually further simplifies the previous linear pressure distribution assumption.

The Improvement of iEDFM
The integrally embedded discrete fracture model (iEDFM) is implemented on EDFM with a new gridding method, a semi-analytic matrix-fracture transmissibility equation from a more realistic pressure distribution assumption inside the matrix element, which would improve simulation accuracy and expand the scope of application.An iEDFM preprocessor was developed with the inputs of reservoir features and fracture geometries.The gridding and transmissibility calculation processes are conducted in this preprocessor.
In Figure 2, we illustrate the procedure to add fracture elements with a 2D case with 4 matrix blocks and 2 fractures.Figure 2a shows that in EDFM, 6 fracture elements have to be added with 14 F-F connections and 6 F-M connections because of the matrix element boundaries and fracture intersections.However, in iEDFM, the fractures can be embedded integrally-either discretizing by matrix element boundaries (Figure 2b) or taking the intersecting fracture group as one element (Figure 2c) is permitted.As a result, the number of fracture elements can be reduced to 3 or 1, and the number of F-F and F-M connections can be reduced to 2 or 0, 3 or 3.

The Improvement of iEDFM
The integrally embedded discrete fracture model (iEDFM) is implemented on EDFM with a new gridding method, a semi-analytic matrix-fracture transmissibility equation from a more realistic pressure distribution assumption inside the matrix element, which would improve simulation accuracy and expand the scope of application.An iEDFM preprocessor was developed with the inputs of reservoir features and fracture geometries.The gridding and transmissibility calculation processes are conducted in this preprocessor.
In Figure 2, we illustrate the procedure to add fracture elements with a 2D case with 4 matrix blocks and 2 fractures.Figure 2a shows that in EDFM, 6 fracture elements have to be added with 14 F-F connections and 6 F-M connections because of the matrix element boundaries and fracture intersections.However, in iEDFM, the fractures can be embedded integrally-either discretizing by matrix element boundaries (Figure 2b) or taking the intersecting fracture group as one element (Figure 2c) is permitted.As a result, the number of fracture elements can be reduced to 3 or 1, and the number of F-F and F-M connections can be reduced to 2 or 0, 3 or 3.With the fracture added, we determined the calculation method of transmissibility in iEDFM.For M-M and F-F connections, Equations ( 7) and ( 8) are applicable.However, other than using Equation ( 10), iEDFM has a new semi-analytic transmissibility equation for M-F connections.
When calculating M-F transmissibility in iEDFM, the analytic solution of pressure distribution around the fractures and superposition principle of potential are applied.The detailed method will be explained later.
In general, EDFM has four weaknesses which have been overcome by iEDFM: 1.The pressure difference between adjacent fracture pieces is relatively small due to the high conductivity in the fracture.Therefore, such fine gridding for fracture in EDFM may bring unnecessary calculation costs and difficulty in convergence; 2. When the matrix block is coarse or the embedded fractures are more complicated, the linear assumption in EDFM is too rough (showed in Case 1); 3.If there are more than one fracture pieces inside one matrix element or the fracture has complex geometries, the pressure distribution inside the matrix element will no longer be available, which will limit the application of EDFM (shown in Case 1(c) and Case 3); 4.Only the fracture piece and its background matrix block are used for the transmissibility calculation, while the global fracture network is not taken into consideration.For example, when calculating the transmissibility of F5-M3 in Figure 2a, according to the linear pressure distribution assumption in EDFM, the pressure drop at the red point is only influenced by F5.This is inaccurate because there is also F6 nearby, which will also cause a pressure drop at the red point (showed in Case 1(b)).With the fracture added, we determined the calculation method of transmissibility in iEDFM.For M-M and F-F connections, Equations ( 7) and ( 8) are applicable.However, other than using Equation ( 10), iEDFM has a new semi-analytic transmissibility equation for M-F connections.
When calculating M-F transmissibility in iEDFM, the analytic solution of pressure distribution around the fractures and superposition principle of potential are applied.The detailed method will be explained later.
In general, EDFM has four weaknesses which have been overcome by iEDFM: 1.
The pressure difference between adjacent fracture pieces is relatively small due to the high conductivity in the fracture.Therefore, such fine gridding for fracture in EDFM may bring unnecessary calculation costs and difficulty in convergence; 2.
When the matrix block is coarse or the embedded fractures are more complicated, the linear assumption in EDFM is too rough (showed in Case 1); 3.
If there are more than one fracture pieces inside one matrix element or the fracture has complex geometries, the pressure distribution inside the matrix element will no longer be available, which will limit the application of EDFM (shown in Case 1(c) and Case 3); 4.
Only the fracture piece and its background matrix block are used for the transmissibility calculation, while the global fracture network is not taken into consideration.For example, when calculating the transmissibility of F5-M3 in Figure 2a, according to the linear pressure distribution assumption in EDFM, the pressure drop at the red point is only influenced by F5.This is inaccurate because there is also F6 nearby, which will also cause a pressure drop at the red point (showed in Case 1(b)).

Calculation of M-F Transmissibility in iEDFM
The M-F transmissibility calculation is based on the assumption of pressure distribution near the fracture.The linear pressure distribution assumed in EDFM is rough, especially when the embedded fractures are non-penetrating or complicated.In iEDFM, fractures can be considered as a bunch of point sinks.Around each point sink, an analytic pressure distribution formula exists.Then, the superposition principle of potential allows us to obtain the semi-analytic pressure distribution near the fracture.With this pressure distribution, the M-F transmissibility equation is obtained.
Given a 2D gridded fractured reservoir as an example (Figure 3a), we could take these two intersecting fractures as one element (iEDFM gridding method II) and consider the incompressible steady state single-phase flow in this reservoir.

Calculation of M-F Transmissibility in iEDFM
The M-F transmissibility calculation is based on the assumption of pressure distribution near the fracture.The linear pressure distribution assumed in EDFM is rough, especially when the embedded fractures are non-penetrating or complicated.In iEDFM, fractures can be considered as a bunch of point sinks.Around each point sink, an analytic pressure distribution formula exists.Then, the superposition principle of potential allows us to obtain the semi-analytic pressure distribution near the fracture.With this pressure distribution, the M-F transmissibility equation is obtained.
Given a 2D gridded fractured reservoir as an example (Figure 3a), we could take these two intersecting fractures as one element (iEDFM gridding method II) and consider the incompressible steady state single-phase flow in this reservoir.The pressure of the fracture element is PF (for the iEDFM gridding method I, a fracture piece in each matrix block has a different pressure, but the formula derivation is also similar, as follows).The pressure of the fracture element is P F (for the iEDFM gridding method I, a fracture piece in each matrix block has a different pressure, but the formula derivation is also similar, as follows).
We use point sink S i (i = 1 to N) to replace whole intersecting fractures.If the point sinks are dense enough (e.g., more than 50 sinks inside one matrix block), the pressure distribution near the fracture could be imitated near these point sinks.Point F j (j = 1 to N) is on the surface of the fracture, thus having the same fracture pressure P F .The top view of the dotted frame part in Figure 3a is presented in Figure 3b.
We define a potential: when the flow reaches the steady state, we have the analytic potential distribution formula of a single point sink from the integral of the plane radial flow equation: Then, the N-dimensional linear equations are obtained from Equation ( 13) and the superposition principle of potential: where q S i is the flow rate of the point sink S i , h is the height of the reservoir, r ij is the distance between the point sink S i and the point F j , and C and c are constant numbers.We define: where ξ i can be solved out from Equation (14).

The Semi-Analytic Calculation Method
Here, we consider the transmissibility between a specific matrix element (for example, the red matrix block in Figure 3a, the top view shown in Figure 3c and the fracture element.The average pressure of the whole matrix element (block) is P M .
Similar to Equation ( 14), the potential of point X near the fracture inside this specific matrix element can be determined by: Thus, we obtain: Combining Equations ( 16) and ( 17) and the definition of transmissibility (Equation ( 4)), the M-F transmissibility and the pressure anywhere inside the matrix block can be calculated by: Energies 2018, 11, 3491 8 of 20 where: where V M is the volume of the specific matrix element, ∈ M means that the point is inside this matrix block, dV x is an element volume of this matrix element, r ix is the distance between the point sink S i and dV x .T M−F is only related to the properties of the reservoir and can be determined at the step of pre-processing before the simulation starts.
For the 3D situation, referring to Equations ( 6) and ( 12), the analytic potential distribution formula of a single point sink is: After a similar derivation, the-F transmissibility and the pressure distribution can be written as: where,

Verifications and Applications
In the following simulation studies, we present four cases to demonstrate the applicability of iEDFM.
First, a single-phase case is considered to demonstrate the improvement of iEDFM upon EDFM.Then, in the second case, we demonstrate the accuracy of our model by comparing the flow rates and saturation profiles with the fine-grid model through a flooding test.Both gridding methods I and II are used in this case.Then, a nonplanar fractures case is presented to show the applicability of iEDFM for a complex geometry situation.At last, a 3D field case demonstrates the potential application of iEDFM in real field studies.
Most of the reservoir properties and operation parameters were kept the same in all the cases, as shown in Table 1.The calculation method of M-F transmissibility is based on the assumption of pressure distribution pattern in the vicinity of fractures.EDFM assumes a uniform pressure gradient in the matrix element, while iEDFM uses a semi-analytic method to calculate the pressure distribution.
In this case, we consider an ideal reservoir with two intersecting fractures (Figure 4), using simulation results from the fine-grid model to verify the effectiveness of iEDFM over EDFM.EDFM or iEDFM simulation is not conducted in this case.Only the pressure distribution and M-F transmissibility for red block (a)/(b)/(c) are calculated through the pre-processing procedure of iEDFM and EDFM, comparing with the simulation result of the whole reservoir of the fine-grid model.Most of the reservoir properties and operation parameters were kept the same in all the cases, as shown in Table 1.The calculation method of M-F transmissibility is based on the assumption of pressure distribution pattern in the vicinity of fractures.EDFM assumes a uniform pressure gradient in the matrix element, while iEDFM uses a semi-analytic method to calculate the pressure distribution.
In this case, we consider an ideal reservoir with two intersecting fractures (Figure 4), using simulation results from the fine-grid model to verify the effectiveness of iEDFM over EDFM.EDFM or iEDFM simulation is not conducted in this case.Only the pressure distribution and M-F transmissibility for red block (a)/(b)/(c) are calculated through the pre-processing procedure of iEDFM and EDFM, comparing with the simulation result of the whole reservoir of the fine-grid model.Single-phase fluid (water) is pumping out from the constant pressure vertical fractures simulated in the fine-grid model discretized by 900 × 900 fine gridblocks horizontally.The element dimensions are non-uniform in the x and y directions to accommodate refinement around fractures.The widths of the fracture elements and their adjacent matrix elements were equal to the fracture aperture.Table 2 supplements some parameters in addition to Table 1.Single-phase fluid (water) is pumping out from the constant pressure vertical fractures simulated in the fine-grid model discretized by 900 × 900 fine gridblocks horizontally.The element dimensions are non-uniform in the x and y directions to accommodate refinement around fractures.The widths of the fracture elements and their adjacent matrix elements were equal to the fracture aperture.Table 2 supplements some parameters in addition to Table 1.This pumping pressure represents the pressure of the fracture element, and the average pressure of the matrix inside the red blocks represents the pressure of the matrix element.The fluid exchange can be obtained from simulation results of the fine-grid model.Therefore, the equivalent M-F transmissibility by the fine-grid model can be calculated through Equation ( 4).
Figure 5b shows that the transmissibility calculation methods in EDFM and iEDFM are more accurate when the fracture fully penetrates the matrix block.It can be seen from Figure 5a that when the fracture does not fully penetrate, the EDFM method will cause some error which cannot be ignored, which is in agreement with the estimation by Xu [30], while the iEDFM method can still be relatively accurate in this situation.In Figure 5c, we can see that the EDFM method produces more obvious errors in the complex situation with fracture intersection, whereas the iEDFM method is still relatively accurate.This pumping pressure represents the pressure of the fracture element, and the average pressure of the matrix inside the red blocks represents the pressure of the matrix element.The fluid exchange can be obtained from simulation results of the fine-grid model.Therefore, the equivalent M-F transmissibility by the fine-grid model can be calculated through Equation ( 4).
Figure 5b shows that the transmissibility calculation methods in EDFM and iEDFM are more accurate when the fracture fully penetrates the matrix block.It can be seen from Figure 5a that when the fracture does not fully penetrate, the EDFM method will cause some error which cannot be ignored, which is in agreement with the estimation by Xu [30], while the iEDFM method can still be relatively accurate in this situation.In Figure 5c, we can see that the EDFM method produces more obvious errors in the complex situation with fracture intersection, whereas the iEDFM method is still relatively accurate.The errors and differences are mainly due to different assumptions regarding the pressure distribution in the vicinity of fractures.In Figure 6, the pressure profiles in block (a) (b) and (c) after 0.1s of pumping are presented.Three kinds of pressure profiles are considered here: pressure distribution calculated by the fine-grid method, and the pressure distribution assumed by EDFM / iEDFM.
As shown in Figure 6, the results of the fine-grid model are quite similar to iEDFM, with negligible difference.However, the results of EDFM indicate that when the matrix block is not fully penetrated by the fracture (Figure 6a) or the embedded fractures are complicated (Figure 6c), the linear distribution of EDFM's assumption can be rough.Even if the fracture penetrates the block (b), the linear distribution still does not exactly reflect the real situation, because of the effect of the fracture segment outside this matrix block, which EDFM does not take into consideration but iEDFM does.The errors and differences are mainly due to different assumptions regarding the pressure distribution in the vicinity of fractures.In Figure 6, the pressure profiles in block (a) (b) and (c) after 0.1s of pumping are presented.Three kinds of pressure profiles are considered here: pressure distribution calculated by the fine-grid method, and the pressure distribution assumed by EDFM/iEDFM.
As shown in Figure 6, the results of the fine-grid model are quite similar to iEDFM, with negligible difference.However, the results of EDFM indicate that when the matrix block is not fully penetrated by the fracture (Figure 6a) or the embedded fractures are complicated (Figure 6c), the linear distribution of EDFM's assumption can be rough.Even if the fracture penetrates the block (b), the linear distribution still does not exactly reflect the real situation, because of the effect of the fracture segment outside this matrix block, which EDFM does not take into consideration but iEDFM does.In summary, iEDFM has a more accurate M-F transmissibility calculation method based on a more realistic pressure distribution assumption.In addition, iEDFM can assume a specific pressure field for any complex embedding situation.If pressure-related physical properties, such as adsorption analysis and diffusion effects, need to be considered, iEDFM is able to show a greater applicability.

Case 2: Intersecting-Fractured Reservoir Flooding Test
Figure 7 shows a 2D fractured reservoir containing three intersecting vertical fractures.This case is a displacement of oil by water, applied in the fine-grid model and the iEDFM model.Table 3 supplements some parameters of this case in addition to Table 1, which will also be used in the following cases.

Parameter
Value Unit Initial reservoir pressure 1 × 10 6  Pa Producer pressure 1 × 10 6 Pa Initial water saturation 20% - In summary, iEDFM has a more accurate M-F transmissibility calculation method based on a more realistic pressure distribution assumption.In addition, iEDFM can assume a specific pressure field for any complex embedding situation.If pressure-related physical properties, such as adsorption analysis and diffusion effects, need to be considered, iEDFM is able to show a greater applicability.

Case 2: Intersecting-Fractured Reservoir Flooding Test
Figure 7 shows a 2D fractured reservoir containing three intersecting vertical fractures.This case is a displacement of oil by water, applied in the fine-grid model and the iEDFM model.Table 3 supplements some parameters of this case in addition to Table 1, which will also be used in the following cases.For the fine-grid simulation, the grid is 900 × 900 elements in the x and y directions.For the iEDFM simulation, two sets of uniform matrix grids (30 × 30/10 × 10) are used.Gridding methods I and II are also compared in this case (gridding method I is applied if not mentioned, as below).The oil rate and the profiles of oil saturation after 115 days of water injection (8.64 × 10 −5 m 3 /day) calculated by iEDFM and fine-grid model are presented in Figures 8 and 9.For the fine-grid simulation, the grid is 900 × 900 elements in the x y directions.For the iEDFM simulation, two sets of uniform matrix grids (30 30/10 × 10) are used.Gridding methods I and II are also compared in this case (gridding method I is applied if not mentioned, as below).The oil rate and the profiles of oil saturation after 115 days of water injection (8.64 × 10 −5 m 3 /day) calculated by iEDFM and fine-grid model are presented in Figures 8 and 9.For the fine-grid simulation, the grid is 900 × 900 elements in the x and y directions.For the iEDFM simulation, two sets of uniform matrix grids (30 × 30/10 × 10) are used.Gridding methods I and II are also compared in this case (gridding method I is applied if not mentioned, as below).The oil rate and the profiles of oil saturation after 115 days of water injection (8.64 × 10 −5 m 3 /day) calculated by iEDFM and fine-grid model are presented in Figures 8 and 9.   Figure 8 compares the oil rates, confirming the accuracy of the iEDFM approach.A good agreement exists in both 30 × 30/10 × 10.The curves of gridding methods I and II almost coincide, indicating that the results of these two methods are not much different when the permeability of fracture is much higher than that of the matrix, respectively.A curve of the same reservoir without any fracture is also present on this figure to show the influence of the existence of fractures.The effect of phase behavior on oil saturations is also pronounced in this case, as Figure 9 shows.The computational times for iEDFM (30 × 30), iEDFM (10 × 10) and the fine-grid model are 9.1 s, 4.3 s and 6.7 h, respectively, which indicates a high efficiency of iEDFM.

Case 3: Nonplanar-Fractured Reservoir Flooding Test
Recent advances in fracture-diagnostic tools and fracture-propagation models make it necessary to model fractures with complex geometries in reservoir-simulation studies.A nonplanar shape is one of the most common complex geometries [30].
Fractures tend to grow in the direction perpendicular to the minimum horizontal stress.In some cases, the preferred direction of fracture may not keep the same, which often leads to a nonplanar shape.
The methodology introduced in this study can be directly applied in a nonplanar fractures case.As Figure 10 shows, the location of the assumed points Si and Fi will reflect the geometry of the fracture.As a result, the semi-analytic method can bring iEDFM enough flexibility in modeling nonplanar fractures.Figure 8 compares the oil production rates, confirming the accuracy of the iEDFM approach.A good agreement exists in both 30 × 30/10 × 10.The curves of gridding methods I and II almost coincide, indicating that the results of these two methods are not much different when the permeability of fracture is much higher than that of the matrix, respectively.A curve of the same reservoir without any fracture is also present on this figure to show the influence of the existence of fractures.The effect of phase behavior on oil saturations is also pronounced in this case, as Figure 9 shows.The computational times for iEDFM (30 × 30), iEDFM (10 × 10) and the fine-grid model are 9.1 s, 4.3 s and 6.7 h, respectively, which indicates a high efficiency of iEDFM.

Case 3: Nonplanar-Fractured Reservoir Flooding Test
Recent advances in fracture-diagnostic tools and fracture-propagation models make it necessary to model fractures with complex geometries in reservoir-simulation studies.A nonplanar shape is one of the most common complex geometries [30].
Fractures tend to grow in the direction perpendicular to the minimum horizontal stress.In some cases, the preferred direction of fracture may not keep the same, which often leads to a nonplanar shape.
The methodology introduced in this study can be directly applied in a nonplanar fractures case.As Figure 10 shows, the location of the assumed points S i and F i will reflect the geometry of the fracture.As a result, the semi-analytic method can bring iEDFM enough flexibility in modeling nonplanar fractures.We present an ideal case as shown in Figure 11 in which the fracture is formed as two connected arcs seen in the top view and is vertical the z direction.The same parameters as Case 2 are applied here, and a fine-grid model similar to Case 2 is built.For the iEDFM simulation, gridding method I and a uniform matrix grid (30 × 30 with one layer in z direction) are used.The oil rate curves are presented in Figure 12 and the oil saturation profiles after 115 days of injection are shown in Figure 13, where a good agreement demonstrates the accuracy of the iEDFM in modeling the nonplanar fractures reservoir.The curves of the same reservoir without any fracture and with a planar fracture with the same starting and ending location are also present on this figure to show the influence of the existence of nonplanar fractures.We present an ideal case as shown in Figure 11 in which the fracture is formed as two connected arcs seen in the top view and is vertical in the z direction.The same parameters as Case 2 are applied here, and a fine-grid model similar to Case 2 is built.For the iEDFM simulation, gridding method I and a uniform matrix grid (30 × 30 with one layer in z direction) are used.We present an ideal case as shown in Figure 11 in which the fracture is formed as two connected arcs seen in the top view and is vertical in the z direction.The same parameters as Case 2 are applied here, and a fine-grid model similar to Case 2 is built.For the iEDFM simulation, gridding method I and a uniform matrix grid (30 × 30 with one layer in z direction) are used.The oil rate curves are presented in Figure 12 and the oil saturation profiles after 115 days of injection are shown in Figure 13, where a good agreement demonstrates the accuracy of the iEDFM in modeling the nonplanar fractures reservoir.The curves of the same reservoir without any fracture and with a planar fracture with the same starting and ending location are also present on this figure to show the influence of the existence of nonplanar fractures.The oil rate curves are presented in Figure 12 and the oil saturation profiles after 115 days of injection are shown in Figure 13, where a good agreement demonstrates the accuracy of the iEDFM in modeling the nonplanar fractures reservoir.The curves of the same reservoir without any fracture and with a planar fracture with the same starting and ending location are also present on this figure to show the influence of the existence of nonplanar fractures.Actually, as Figures 3b and 10 show, different locations point Si and Fi are able to reflect any geometry of the fractures, such as fractures with variable apertures and vuggy-fractures, which means that iEDFM is naturally suitable for any other complex geometry cases beside nonplanar fractures.

Case 4: 3D Case with a Modified Dataset of a Real Field
As mentioned previously, iEDFM can be used as a general procedure in both 2D and 3D cases.In real field applications, the reservoir may have multiple layers, and the height of the fractures can be smaller than the reservoir height.Therefore, a 3D simulation example is presented to show the application of iEDFM in a typical field study.
In this case, the geological model is modified from a reported dataset which is interpreted from a 3D seismic survey [35,36].Some irregular, sparsely distributed large-scale fractures are present as main fractures (black planes in Figure 14a).Some stochastic medium-scale fractures (blue planes in Figure 14a) are added to test iEDFM's applicability in a complex fracture-network situation.The main fractures and the wells are assumed to extend throughout the entire depth of the reservoir, while the Actually, as Figures 3b and 10 show, different locations of point Si and Fi are able to reflect any geometry of the fractures, such as fractures with variable apertures and vuggy-fractures, which means that iEDFM is naturally suitable for any other complex geometry cases beside nonplanar fractures.

Case 4: 3D Case with a Modified Dataset of a Real Field
As mentioned previously, iEDFM can be used as a general procedure in both 2D and 3D cases.In real field applications, the reservoir may have multiple layers, and the height of the fractures can be smaller than the reservoir height.Therefore, a 3D simulation example is presented to show the application of iEDFM in a typical field study.
In this case, the geological model is modified from a reported dataset which is interpreted from a 3D seismic survey [35,36].Some irregular, sparsely distributed large-scale fractures are present as main fractures (black planes in Figure 14a).Some stochastic medium-scale fractures (blue planes in Figure 14a) are added to test iEDFM's applicability in a complex fracture-network situation.The main fractures and the wells are assumed to extend throughout the entire depth of the reservoir, while the Actually, as Figures 3b and 10 show, different locations of point S i and F i are able to reflect any geometry of the fractures, such as fractures with variable apertures and vuggy-fractures, which means that iEDFM is naturally suitable for any other complex geometry cases beside nonplanar fractures.

Case 4: 3D Case with a Modified Dataset of a Real Field
As mentioned previously, iEDFM can be used as a general procedure in both 2D and 3D cases.In real field applications, the reservoir may have multiple layers, and the height of the fractures can be smaller than the reservoir height.Therefore, a 3D simulation example is presented to show the application of iEDFM in a typical field study.
In this case, the geological model is modified from a reported dataset which is interpreted from a 3D seismic survey [35,36].Some irregular, sparsely distributed large-scale fractures are present as main fractures (black planes in Figure 14a).Some stochastic medium-scale fractures (blue planes in Figure 14a) are added to test iEDFM's applicability in a complex fracture-network situation.The main fractures and the wells are assumed to extend throughout the entire depth of the reservoir, while the stochastic medium-scale fractures are of different heights.Figure 14b shows the matrix grids (25 × 25 × 5) and the position of fractures and wells.An injection well (4 × 10 4 m 3 /day) is placed in the center, while producers are placed in four corners of the reservoir.Again, the reservoir without any fracture is also simulated a comparison.The oil rate curves are presented in Figure 15.The curves of wells 1-4 in the reservoir without fractures are coincident because of the symmetry.The pressure profiles of the top and bottom layers after 7300 days of injection are shown in Figure 16.As we can see, the reservoir shows strong heterogeneity due to the existence of fractures, and the water's breakthrough is advanced, which reduces the recovery efficiency.As medium-scale fractures are added, these phenomena become more pronounced.Because of gravity, the average oil saturation of the first layer is higher than that of the bottom layer (0.6952 of Figure 16c   The oil rate curves are presented in Figure 15.The curves of wells 1-4 in the reservoir without fractures are coincident because of the symmetry.The pressure profiles of the top and bottom layers after 7300 days of injection are shown in Figure 16.As we can see, the reservoir shows strong heterogeneity due to the existence of fractures, and the water's breakthrough is advanced, which reduces the recovery efficiency.As medium-scale fractures are added, these phenomena become more pronounced.Because of gravity, the average oil saturation of the first layer is higher than that of the bottom layer (0.6952 of Figure 16c    (e) (f) The case study showed the applicability of the iEDFM in reservoirs with complex fracture networks.The influence of different scales of fractures can be modeled appropriately by iEDFM.The 3D multiphase simulation example demonstrates the potential application of the iEDFM in real field studies.

Conclusions and Future Work
In this study, we developed a new approach called the integrally embedded discrete fracture model (iEDFM).This approach, for the first time, avoids the limitations of the need to subdivide fracture elements in EDFM.
In iEDFM, we can arbitrarily grid fractures according to the requirements, and then embed them integrally in matrix blocks.As a precise pressure distribution assumption inside the matrix blocks is introduced, we can obtain a semi-analytic calculation method of matrix-fracture transmissibility.As a result, the simulation accuracy is improved and the application is also expanded to fractures with complex geometries.Several cases have been presented to support these conclusions.The potential application of the iEDFM in real field studies has also been testified through a 3D case.
Applying the iEDFM to real field study and guiding production is our ultimate goal.Thus, heterogeneous and more complex actual reservoir examples with actual production data will be considered in our on-going project.The case study showed the applicability of the iEDFM in reservoirs with complex fracture networks.The influence of different scales of fractures can be modeled appropriately by iEDFM.The 3D multiphase simulation example demonstrates the potential application of the iEDFM in real field studies.

Conclusions and Future Work
In this study, we developed a new approach called the integrally embedded discrete fracture model (iEDFM).This approach, for the first time, avoids the limitations of the need to subdivide fracture elements in EDFM.
In iEDFM, we can arbitrarily grid fractures according to the requirements, and then embed them integrally in matrix blocks.As a precise pressure distribution assumption inside the matrix blocks is introduced, we can obtain a semi-analytic calculation method of matrix-fracture transmissibility.As a result, the simulation accuracy is improved and the application is also expanded to fractures with complex geometries.Several cases have been presented to support these conclusions.The potential application of the iEDFM in real field studies has also been testified through a 3D case.
Applying the iEDFM to real field study and guiding production is our ultimate goal.Thus, heterogeneous and more complex actual reservoir examples with actual production data will be considered in our on-going project.

Figure 1 .
Figure 1.Explanation of the embedded discrete fracture model (EDFM) gridding [26]: (a) a fracture segment is embedded in a matrix block; (b) two fracture planes intersect in a matrix block; and (c) two fracture segments embedded in neighboring matrix blocks.

Figure 1 .
Figure 1.Explanation of the embedded discrete fracture model (EDFM) gridding [26]: (a) a fracture segment is embedded in a matrix block; (b) two fracture planes intersect in a matrix block; and (c) two fracture segments embedded in neighboring matrix blocks.

Figure 2 .
Figure 2. Explanation of the integrally embedded discrete fracture model (iEDFM) gridding: (a) EDFM gridding method; (b) iEDFM gridding method I-discretizing by matrix element boundaries; and (c) iEDFM gridding method II-taking the intersecting fracture group as one element.

Figure 2 .
Figure 2. Explanation of the integrally embedded discrete fracture model (iEDFM) gridding: (a) EDFM gridding method; (b) iEDFM gridding method I-discretizing by matrix element boundaries; and (c) iEDFM gridding method II-taking the intersecting fracture group as one element.

Figure 3 .
Figure 3. Illustration showing the calculation of M-F transmissibility in iEDFM: (a) a 2D example of the fractured reservoir; (b) the top view of the dotted frame part; and (c) the top view of the red matrix block.

Figure 3 .
Figure 3. Illustration showing the calculation of M-F transmissibility in iEDFM: (a) a 2D example of the fractured reservoir; (b) the top view of the dotted frame part; and (c) the top view of the red matrix block.

Figure 4 .
Figure 4.The fractured reservoir in Case 1.The 9 × 9 grids showed a possible gridding choice for iEDFM or EDFM.

Figure 4 .
Figure 4.The fractured reservoir in Case 1.The 9 × 9 grids showed a possible gridding choice for iEDFM or EDFM.

Figure 5 .
Figure 5. TM-F of the fine-grid model (equivalent)/iEDFM/EDFM for block (a)/(b)/(c).The equivalent transmissibility calculated by the fine-grid model will change over time until the flow approaches a steady state, while the ones calculated through the pre-processing procedure of iEDFM/EDFM stay the same.

Figure 5 .
Figure 5. T M-F of the fine-grid model (equivalent)/iEDFM/EDFM for block (a)/(b)/(c).The equivalent transmissibility calculated by the fine-grid model will change over time until the flow approaches a steady state, while the ones calculated through the pre-processing procedure of iEDFM/EDFM stay the same.

Figure 7 .
Figure 7.The fractured reservoir in Case 2. A water injector is placed in one corner of the reservoir and a producer is located in the opposite corner.

Figure 8 .
Figure 8.Comparison of oil production rate by the iEDFM/fine-grid model in Case 2.

Figure 7 .
Figure 7.The fractured reservoir in Case 2. A water injector is placed in one corner of the reservoir and a producer is located in the opposite corner.

Energies 2018 , 20 Figure 7 .
Figure 7.The fractured reservoir in Case 2. A water injector is placed in one corner of the reservoir and a producer is located in the opposite corner.

Figure 8 .
Figure 8.Comparison of oil production rate by the iEDFM/fine-grid model in Case 2.Figure 8. Comparison of oil production rate by the iEDFM/fine-grid model in Case 2.

Figure 8 .
Figure 8.Comparison of oil production rate by the iEDFM/fine-grid model in Case 2.Figure 8. Comparison of oil production rate by the iEDFM/fine-grid model in Case 2.

Figure 10 .
Figure 10.Schematic of equivalent point sinks and other parameters used in the nonplanar fracture case corresponding to Figure 3b.

Figure 11 .
Figure 11.The nonplanar-fractured reservoir in Case 3. A water injector (8.64 × 10 −5 m 3 /day) is located in one corner of the reservoir, and a producer is placed in the opposite corner.

Figure 10 .
Figure 10.Schematic of equivalent point sinks and other parameters used in the nonplanar fracture case corresponding to Figure 3b.

Energies 2018 , 20 Figure 10 .
Figure 10.Schematic of equivalent point sinks and other parameters used in the nonplanar fracture case corresponding to Figure 3b.

Figure 11 .
Figure 11.The nonplanar-fractured reservoir in Case 3. A water injector (8.64 × 10 −5 m 3 /day) is located in one corner of the reservoir, and a producer is placed in the opposite corner.

Figure 11 .
Figure 11.The nonplanar-fractured reservoir in Case 3. A water injector (8.64 × 10 −5 m 3 /day) is located in one corner of the reservoir, and a producer is placed in the opposite corner.

Figure 12 .Figure 13 .
Figure 12.Comparison of oil production rate by the iEDFM/fine-grid model in Case 3.

Energies 2018 ,Figure 14 .
Figure 14.(a) Three-dimensional view of main fractures and medium-scale fractures; (b) The matrix grids of the reservoir and the position of each fracture and each well.

Figure 14 .
Figure 14.(a) Three-dimensional view of main fractures and medium-scale fractures; (b) The matrix grids of the reservoir and the position of each fracture and each well.

Figure 11 .
Figure 11.Comparison of oil production rate for Case 4.

Figure 16 .
Figure 16.Profiles of oil saturation after 7300 days' injection for Case 4 of: (a) the top layer of the nonfracture reservoir; (b) the bottom layer of the non-fracture reservoir; (c) the top layer of the reservoir with only main fractures; (d) the bottom layer of the reservoir with only main fractures; (e) the top layer of the reservoir with all fractures; and (f) the bottom layer of the reservoir with all fracture.

Figure 16 .
Figure 16.Profiles of oil saturation after 7300 days' injection for Case 4 of: (a) the top layer of the non-fracture reservoir; (b) the bottom layer of the non-fracture reservoir; (c) the top layer of the reservoir with only main fractures; (d) the bottom layer of the reservoir with only main fractures; (e) the top layer of the reservoir with all fractures; and (f) the bottom layer of the reservoir with all fracture.

Table 1 .
Basic reservoir and fluid parameters in simulations for case 1~4.

Table 1 .
Basic reservoir and fluid parameters in simulations for case 1~4.

Table 2 .
Some parameters in simulations for Case 1.

Table 2 .
Some parameters in simulations for Case 1.

Table 3 .
Some parameters in simulations for Cases 2-4.

Table 3 .
Some parameters in simulations for Cases 2-4.