Study on a Dual Embedded Discrete Fracture Model for
Fluid Flow in Fractured Porous Media

Simulation of fluid flow in the fractured porous media is very important and challenging. Researchers have developed some models for fractured porous media. With the development of related research in recent years, the prospect of embedded discrete fracture model (EDFM) is more and more bright. However, since the size of the fractures in the actual reservoir varies greatly, a very fine grid should be used which leads to a huge burden to the computing resources. To address this challenge, in the present paper, an upscaling based model is proposed. In this model, the flow in large-scale fractures is directly described by the EDFM while that in the small-scale fractures is upscaled through local simulation by EDFM. The EDFM is used to simulate the largeand small-scale fractures independently two times, so the new model is called dual embedded discrete fracture model (D-EDFM). In this paper, the detailed implementation process of D-EDFM is introduced and, through test cases, it is found the proposed model is a feasible method to simulate the flow in fractured porous media.


Introduction
Fractured porous media can be found in many fields, such as subsurface water flow [Barenblatt, Zheltov and Kochina (1960)], reservoir development process ], CO2 sequestration [Li, Li, Yuan et al. (2019)], geothermal exploitation [Kushnir, Heap and Baud (2018)] etc. The fluid flow simulation in the fractured porous media has become one of the research hotspots in recent decades. The key point is to find an accurate and efficient way to describe the fractures. Researchers  some models. The typical ones that can clearly depict the fractures' geometry structure are discrete fracture model (DFM) and embedded discrete fracture model (EDFM).
The discrete fracture model [Karimi-Fard, Durlofsky and Aziz (2004) ;Garipov, Karimi-Fard and Tchelepi (2016)] explicitly describes the structures of fracture system, and it is considered as the most accurate methods to describe fluid flow and heat transfer in fractured porous media. However, the shortcoming of DFM is very obvious. Since the grid of matrix and fractures are related/correlated with each other, unstructured grids should be used and matrix grid refinement is necessary near the fractures. The grids number increases significantly as the fracture number increases, which leads to large computational cost and low computational efficiency. As a result, for complex threedimensional fracture network in engineering, simulation through DFM is almost impossible [Chen, Clauser, Marquart et al. (2018)].
The embedded discrete fracture model was first put forward by Lee et al. [Lee, Lough and Jensen (2001)] based on DFM. In EDFM, the matrix grid and fracture grid are nearly independent. The matrix can be discretized through structured grids, and the fractures are embedded in the matrix grid, thus avoiding the large number of unstructured grids needed in DFM. Subsequently, Moinfar et al. [Moinfar, Varavei, Sepehrnoori et al. (2014)] proposed an improved EDFM. The improved model can simulate inclined fractures in any direction, which means the method is feasible to describe the complexity and heterogeneity of fractured reservoirs. Tene et al. [Tene, Bosma, Al Kobaisi et al. (2017)] found that if the fracture coincides with the grids interfaces when the EDFM is used to simulate fractures, the results obtained are not accurate. Therefore, Tene et al. put forward projection-based embedded discrete fracture model (p-EDFM) to solve this problem. Recently, EDFM is used in many research fields [Shakiba and Sepehrnoori (2015); Xu (2015); Ren, Jiang and Younis (2016);Yu, Tripoppoom, Sepehrnoori et al. (2018)] showing good performance.
In the actual oil, gas and geothermal reservoir, the fracture networks are usually very complex and the size of fractures varies greatly. The computation burden could be unacceptable to depict all the fractures even by EDFM. One of the solutions [Efendiev, Lee, Li et al. (2015); Bosma, Hajibeygi, Tene et al. (2017) Li, Gao, Han et al. (2020)] to deal with this problem is to upscale all the fractures as permeable matrix based on DFM or EDFM. The upscaling process is called "flow based simulation method" [Durlofsky (2003); Long, Remer, Wilson et al. (1982);Bogdanov, Mourzenko, Thovert et al. (2003)] which is contrary to the "analytical method" [Ahmed Elfeel and Geiger (2012)]. The "flow based simulation method" calculates the equivalent permeability based on local simulation with particular boundary conditions. Koudina et al. [Koudina, Garcia, Thovert et al. (1998)] studied the equivalent permeability for three dimensional fractures with matrix permeability being ignored. Since the flow through matrix may have influence, by DFM, Bogdanov et al. [Bogdanov, Mourzenko, Thovert et al. (2003)] calculated the effective permeability considering the flow through matrix as well as the mass exchange of matrix and fractures. As it is stated above, for complex fractures network, the DFM needs large amount of grids which can lead to low computational efficiency. To speed up the upscaling process, Fumagalli et al. [Fumagalli, Zonca and Formaggia (2017)] and Dong et al. [Dong, Li, Lei et al. (2019)] used EDFM to calculate transmissibility and equivalent permeability respectively. The model accuracy can be very low, if all large-and small-scale fractures are upscaled. Researchers have presented the idea that "upscale the small-scale fractures while explicitly depict the large-scale fractures", which is known as a kind of hierarchical model [Pluimers (2015); Hajibeygi, Karvounis and Jenny (2011); Akkutlu, Efendiev and Vasilyeva (2016)]. In the research by Lee et al. [Lee, Lough and Jensen (2001)], the small-scale fractures were upscaled by analytical method while the large ones were depicted by EDFM. Li et al. [Li and Lee (2008)] classified the fractures into small, medium and long fractures. To save the upscaling time consumption, the small ones were upscaled through analytical method. To improve the accuracy, the medium ones were homogenized through boundary element method. Same as Lee, the large ones were explicitly depicted by EDFM.
From the above reviews, we can find that the upscaling methods established from flow based simulation method still have some shortcomings, for instance, when simulating the flow in the real large-scale fractured porous media, the current computer power cannot be afforded. Therefore, motivated by the flow based simulation method, we proposed a new upscaling method, named dual embedded discrete fracture model (D-EDFM). In the proposed model, the fluid flow of large-scale fractures are depicted directly through EDFM, while that of small-scale fractures are upscaled through local simulation by EDFM. The EDFM is independently used to simulate the large-and small-scale fractures for two times.
The structure of this paper is organized as follows. The model description and implementation processes of D-EDFM are introduced to the Section 2. The accuracy and robustness of the D-EDFM is tested for Section 3. The conclusions are presented in the Section 4.

Model description
In the actual fractured porous media, the fracture network is very complex and the sizes of fractures vary greatly. To describe the fluid flow and heat transfer of large-and small-scale fractures, a dual embedded discrete fracture model is presented in this paper. Since the proposed model is based on the traditional EDFM, the EDFM is briefly introduced first and the D-EDFM is subsequently introduced in detail.

Embedded discrete fracture model (EDFM)
The concept of EDFM is shown in Fig. 1. In EDFM, the matrix and fractures are considered to be two separate media coupled with the flow between them [Li, Gao, Han et al. (2020)]. The flow through the matrix and fractures is governed by mass conservation equation and Darcy's law shown as Eqs. (1)-(3). The matrix and the fractures are coupled with the flow exchanged between matrix and fractures given in Eq. (4). The coupling between different fractures is governed by flow between them described by Eq. (5).

Mass conservation equation for the matrix:
Mass conservation equation for the fractures: Darcy's law: Exchange flow Q frÀm between the matrix and the fractures: Exchange flow between different fractures: In the above equations, the superscripts m and fr represent the matrix and the fracture respectively, and subscript f is for fluid. ρ f is the density of fluid, kg/m 3 . c t is coefficient of fluid compressibility, 1/Pa. f m is porosity of matrix. P m is pressure in the matrix, Pa. Q w is the well source term. f fr is porosity of fractures. k i is permeability of media i, m 2 . z is depth of matrix, m. CI is the connectivity index between the matrix and fractures (CI=ΘA/d). Θ is the mean mobility of the fluid, defined as the fraction of permeability and viscosity. A represents area of matrix grid occupied by the fracture segment, m 2 . d represents the average distance from all points in the matrix grid to the fracture segment, m, which is related to the geometry relationship between fracture segment and matrix grid [Durlofsky (2003)]. TI is the transmissivity between different fractures.
b fr represents fracture aperture, m. k fr represents permeability in the fractures, m 2 . μ f represents dynamic viscosity of working fluid, Pa·s. L i represents length of the fracture segment, m.
Based on the grid in EDFM shown in Fig. 1, the above equations can be numerically solved through finite volume method (FVM) and finite element method (FEM). The FVM is used in this paper due to its good performance on mass conservation and convenience in implementation, more details of which can be found in Li et al. [Li, Gao, Han et al. (2020)].

Dual embedded discrete fracture model (D-EDFM)
The concept of D-EDFM is shown in Fig. 2. The fractures are divided into small-scale fractures (blue lines) and large-scale fractures (red lines). The distinction between largeand small-scale fractures is related to the size of the matrix grid. This method to classify the fracture is similar to the method proposed by Lee et al. [Lee, Jensen and Lough (1999)]. If the fracture length is much larger than the matrix cell size (L f /L g >> 1), such a fracture is defined as a large fracture, otherwise it is a small fracture. The large-scale fractures are explicitly depicted through EDFM with coarse grids. Based on particular boundary condition, the small-scale fractures are upscaled as matrix with anisotropic permeability tensor by off-line local simulation, which is known as "flow based simulation method" [Ahmed Elfeel and Geiger (2012)].
The obtained effective permeability tensor is very crucial to D-EDFM which directly influences the accuracy of D-EDFM. Since its value is calculated based on particular boundary condition, the type of the boundary condition is very important. "Sealed side boundary condition" and "Linear boundary condition" are the most used ones in related references. In this paper, to find the proper boundary condition for upscaling procedure in D-EDFM, both of the two types of boundary conditions are adopted and compared. Based on the boundary conditions, the detail implementation process to obtain the effective permeability is given as follows: Step 1: For a cell containing small-scale fractures (shown in Fig. 3), a constant pressure gradient on the x-direction is given. The faces on the y-direction are set to no flow boundary. Length of the cell in x-direction is L x , and in y-direction is L y .
(1)-(5) by EDFM and obtain the pressure distribution in the cell.
Step 3: Calculate the flow rate Q x through a line perpendicular to the x-direction. Any line is feasible, since the total flow along x-direction is a constant.
Step 4: Calculate the equivalent permeability K eff x of the x-direction flow according to the Eq. (7): Step 5: Switch the boundary conditions in the x-direction and y-direction. Then, calculate the K eff y through the similar process given in Steps 2-4.

Figure 3: Schematic diagram of sealed side boundary condition
Step 6: Combine K eff x and K eff y , the effective permeability tensor is obtained shown as Eq. (8) Fig. 4 gives the schematic diagram of linear boundary condition. Similar as "sealed side boundary condition method", a fixed pressure gradient is given in one direction. Instead of sealed side boundary shown in Fig. 3, the linear boundary is given in the other direction (See Fig. 4).

Obtain the effective permeability based on linear boundary
The detail implementation procedures to upscale the small-scale fractures as anisotropic matrix are introduced as follows: Step 1: Apply a constant pressure gradient for the cell shown in Fig. 4 on the x-direction.
A linear boundary condition shown in Eq. (9) is given in the upper and lower boundaries.
(1)-(5) by EDFM and obtain the pressure distribution in the cell.
Step 3: Calculate the average flux and average pressure gradient in the cell shown in Fig. 4.

Figure 4: Schematic diagram of linear boundary condition
Step 4: Based on Darcy's law, where, the superscript x means the pressure gradient shown in Fig. 4 on the x-direction. rP x x and u x x are the value of pressure gradient and flux at any parallel faces on the x-direction respectively.
Step 5: Switch the boundary conditions given in Step 1, two similar equations can be obtained shown as Eqs. (12) and (13).
where, the superscript y means the pressure gradient given on the y-direction.
3 Equations and mathematical expressions In order to validate and evaluate the D-EDFM, two numerical examples are given in this section. Example 1 is given to compare the two types of boundary conditions for upscaling, which has 3 large-scale fractures and 14 small-scale fractures (See Fig. 5(a)).
To test the feasibility of D-EDFM for complex fracture network (See Fig. 5(b)), Example 2, which has 10 large-scale fractures and 500 small-scale fractures, is adopted.

Comparisons of the two types of boundary conditions for upscaling (Example 1)
The geometry of the fractured porous media of Example 1 is given in Fig. 5(a). It is assumed that fluid properties and reservoir parameters do not change with pressure. The left and right boundary conditions are Dirichlet type, which is 1×10 6 Pa and 1×10 5 Pa respectively. The type of the boundary conditions for the up and down boundaries is Neumann boundary condition (no flow boundary). Tab. 1 shows the basic parameters in the model.
To compare the two types of boundary conditions for upscaling given in Section 2.2, with coarse grids (40×40), the pressure field in steady state of Example 1 is calculated through D-EDFM based on "sealed side boundary upscaling" and "linear boundary upscaling" respectively. To get the benchmark solution, Example 1 is also simulated by EDFM with fine grids (100×100).
The results of D-EDFM with coarse grids and the results of EDFM with fine grids are shown in Figs. 6(a), 6(b) and 6(c), respectively. The white lines in Figs. 6(a) and 6(b) are those smallscale fractures which have been upscaled as matrix with anisotropic permeability tensor.
It can be shown in Fig. 6 that the results of D-EDFM based on coarse grids agree well with the result of EDFM based on fine grids. And both sealed side boundary upscaling and linear boundary upscaling are feasible for the D-EDFM. To quantitatively compare the accuracy of two upscaling methods, the pressure distribution of the line from the lower left corner to the top right corner is plotted (Y=X).
From Fig. 7, it is found that in the case of this paper, the average relative error and maximum relative error of D-EDFM based on sealed side boundary upscaling are 0.30% and 1.73% respectively. While those of D-EDFM based on linear boundary upscaling are 1.60% and  4.32% respectively. The accuracy of both upscaling methods is acceptable to engineering, but the sealed side boundary upscaling method has better accuracy and is easy in implementation. Therefore, the sealed side boundary upscaling is recommended in D-EDFM.

Test of D-EDFM for complex fractures network (Example 2)
Through the above simple example (Example 1), a preliminary conclusion can be obtained: the D-EDFM has good performance for simple fractures network based on sealed side boundary upscaling. In this section, the D-EDFM is tested on a much more complex fractures network (Example 2) which is given in Fig. 5(b). In Example 2, the physical parameters of the model are identical to the parameters in Example 1. The difference is that Example 2 is an unsteady problem which has an initial field of 0 Pa. The simulation results of different time are given in Figs. 8-11. In this example, the white lines also represent small fractures which have been upscaled as matrix with anisotropic permeability tensor.
From the pressure field distributions, it can be seen that results of D-EDFM under coarse grids agree well with those obtained by EDFM under fine grids. To further verify this, the pressure on three straight lines (Y=10 m, Y=50 m and Y=80 m) are compared in Fig. 12.   (1) Through the comparisons between results calculated by D-EDFM and traditional EDFM, the proposed D-EDFM shows good accuracy and robustness. The advantage of D-EDFM is reflected on the good accuracy of its results with much coarser grids, the number of which is, in the given examples, 16% of that by traditional EDFM. This will greatly improve the simulation efficiency.
(2) The accuracy of D-EDFM depends on the upscaling process of small-scale fractures. This paper compares the sealed side boundary upscaling and linear boundary upscaling. It is found that both the two upscaling processes can offer accurate effective permeability tensor. The sealed side boundary upscaling for D-EDFM is recommended in this paper, since it has better performance and is easier to be implemented.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.