An Integrally Embedded Discrete Fracture Model for Flow Simulation in Anisotropic Formations

: The embedded discrete fracture model (EDFM), among different flow simulation models, achieves a good balance between efficiency and accuracy. In the EDFM, micro-scale fractures that cannot be characterized individually need to be homogenized into the matrix, which may bring anisotropy into the matrix. However, the simplified matrix–fracture fluid exchange assumption makes it difficult for EDFM to address the anisotropic flow. In this paper, an integrally embedded discrete fracture model (iEDFM) suitable for anisotropic formations is proposed. Structured mesh is employed for the anisotropic matrix, and the fracture element, which consists of a group of connected fractures, is integrally embedded in the matrix grid. An analytic pressure distribution is derived for the point source in anisotropic formation expressed by permeability tensor, and applied to the matrix–fracture transmissibility calculation. Two case studies were conducted and compared with the analytic solution or fine grid result to demonstrate the advantage and applicability of iEDFM to address anisotropic formation. In addition, a two-phase flow example with a reported dataset was studied to analyze the effect of the matrix anisotropy on the simulation result, which also showed the feasibility of iEDFM to address anisotropic formation with complex fracture networks.


Introduction
With the increasing demand for energy, the development of reservoirs bearing essential resources, such as oil and gas, is of great importance. In many exploitation practices, the target reservoirs are naturally fractured [1,2]. During the development of the reservoirs, numerical simulation is a reasonable and affordable tool to manage production and evaluate economic feasibility. A number of numerical models for fractured reservoir have been established in the past several decades. The dual porosity model (DPM), with the assumption of fracture and matrix uniformity, was proposed as one of the earliest methods [3][4][5]. In addition to the subsequent multi-continuum methods, such as dual-porosity dual-permeability (DPDK) model [6,7], multiple interacting continua (MINC) model [8,9], subdomain model [10], triple porosity dual-permeability (TPDK) model [11,12] and multi-porosity model [13], the DPM provides an efficient approach to simulate micro-scale fractures. However, it is also limited by the loss of detailed information relating to discovered macro-scale factures, such as geometry and location [14]. The discrete fracture model (DFM) uses unstructured grids to discretize the discovered macro-scale fractures explicitly, which is more physics-based and accurate [15][16][17][18][19][20], but is limited because of its complexity in gridding and computational cost [19]. Therefore, improvements of the DFM have mainly been focused on calculation efficiency. For example, by applying the fracture crossflow equilibrium ∂ ∂t ϕS β ρ β = −∇· ρ β v β + q β (1) where β represents different components: β = g stands for gas, β = w stands for water, and β = o stands for oil; ϕ is the effective porosity; S β is the saturation of β; ρ β is the density of β; q β is the sink/source term; and v β is the velocity of β, which can be defined by Darcy's law: k xx k xy k xz k xy k yy k yz k xz k yx k zz where K is an absolute permeability tensor; k r β is the relative permeability to β; µ β is the viscosity of β; P β is the pressure of β under reservoir conditions; g is gravitational acceleration; D is the depth from the reference datum. For a certain porous medium, the permeability and the velocity under a certain pressure gradient are determined, but the component of the permeability tensor and velocity vector varies with different coordinate systems. If the coordinate axis is rotated to an appropriate direction where k xy = k yz = k xz = 0, the coordinate axis direction becomes the principal axis direction of the permeability. The permeability tensor is diagonal as follows: For example, tensor diagonalization for two-dimensional cases can be obtained as: Unless otherwise stated, the coordinate axis direction in this paper is considered rotated to the principal axis direction of the permeability, where k xy = k yz = k xz = 0. Therefore, we can use k x , k y , k z to represent the absolute permeability in three directions. In addition, the embedding algorithm when the two sets of axes do not coincide is also introduced, which can be applied in simulators that support permeability tensor.
As implemented numerically, Equation (1) is written into discretized form with an integral finite-difference or control-volume scheme in the space domain, and a finite-difference, backward, first-order scheme in the time domain [11]: where n is the previous time level; n + 1 is the current time level; ∆t is the time step size; V i is the volume of element i; η i contains the set of neighboring elements j to which element i is directly connected; Q β,i is the sink/source term at element i; F β,ij is the flow term between element i and j, which can be described by a discrete version of Darcy's law, given by: F β,ij = λ β,ij+1/2 T i j P β,j − ρ β,ij+1/2 gD j − P β,i − ρ β,ij+1/2 gD i (8) Energies 2020, 13, 3070 4 of 21 where D i is the depth from the reference datum to the center of element i; i j + 1/2 is a proper averaging or weighting of properties; λ β,ij+1/2 is the mobility term, defined as: λ β,ij+1/2 = ρ β k rβ µ β ij+1/2 (9) where T ij is transmissibility, which is the harmonic average of T i (element i to common interface) and T j (element j to common interface): where A ij is the common interface area between elements i and j; d i and d j are the feature distances; k i and k j are the absolute permeability along each element center to the interface. The specific physical meaning of the parameters varies for different connections, which will be discussed in the next section.

Integrally Embedded Discrete Fracture Model
The integrally embedded discrete fracture model (iEDFM) is a pre-processor with the inputs of reservoir features and fracture geometries, and outputs of element and connection information. It is implemented on EDFM with a different gridding method. In this work, a new transmissibility calculation method for matrix-fracture connection in iEDFM suitable for anisotropic formation is presented, based on a more realistic pressure distribution assumption in the vicinity of fractures inside the anisotropic matrix element domain. The transmissibility can be determined before simulation, so it can be easily implemented in various simulators.
In EDFM/iEDFM, Cartesian grids are used for the matrix. Every fracture element consists of one or several connected fractures. Each fracture element is connected with its corresponding matrix element. Thus, there mainly exist 3 kinds of connections: matrix-matrix (M-M), fracture-fracture (F-F), and matrix-fracture (M-F). Due to the existence of wells, well-matrix/fracture is also an important connection. The focus of this work is on the flow between fractures and the matrix, so the vertical wells are simplified as sink points.
In Figure 1, we illustrate the procedure to discretize fracture elements with a 2D case with 4 matrix blocks and 2 intersecting fractures. Six fracture elements will be generated in the EDFM, as shown in Figure 1a. Both the boundary of the matrix blocks and intersection of the two fractures discretize the fracture elements. Meanwhile, in iEDFM, the fractures can be embedded integrally-either discretizing by only matrix block boundaries (Figure 1b) or taking the whole intersecting fracture group as one element (Figure 1c) is permitted. As a result, the number of fracture elements, as well as the number of F-F and F-M connections, can be reduced. In iEDFM(II), the entire set of the intersecting fracture group is regarded as a single element, which is also an equipotential body, thus, the fluid velocity through this fracture group is infinite. It is obviously an ideal approximation, only suitable for a situation where the fracture permeability is far greater than the matrix. Therefore, iEDFM(I) is more recommended for more realistic situations.
Energies 2020, 13, 3070 5 of 21 fracture group as one element (Figure 1c) is permitted. As a result, the number of fracture elements, as well as the number of F-F and F-M connections, can be reduced. In iEDFM(II), the entire set of the intersecting fracture group is regarded as a single element, which is also an equipotential body, thus, the fluid velocity through this fracture group is infinite. It is obviously an ideal approximation, only suitable for a situation where the fracture permeability is far greater than the matrix. Therefore, iEDFM(I) is more recommended for more realistic situations. The pressure difference between adjacent fracture elements is relatively small due to the high conductivity in the fracture. Therefore, the coarser discretization method of iEDFM can reduce calculation costs and improve the convergence [30].
The transmissibility of the three kinds of connections can be calculated referring to Equations (10) and (11) but with different physical meaning of the parameters: For M-M connections, the permeability difference between the two connected elements is relatively small, and the pressure inside each element domain nearly has a linear distribution (shown by the green line in Figure 2a), which approximately meets the assumption of one-dimensional Darcy flow. Thus, the parameters give clear physical meanings: the feature distance d i and d j are the distances between the center of each matrix element to the common interface; and the absolute permeability k i and k j are equal to the matrix permeability k x , k y or k z because the Cartesian matrix block is usually built in the direction of the coordinate system. The pressure difference between adjacent fracture elements is relatively small due to the high conductivity in the fracture. Therefore, the coarser discretization method of iEDFM can reduce calculation costs and improve the convergence [30].
The transmissibility of the three kinds of connections can be calculated referring to Equations (10) and (11) but with different physical meaning of the parameters: For M-M connections, the permeability difference between the two connected elements is relatively small, and the pressure inside each element domain nearly has a linear distribution (shown by the green line in Figure 2a), which approximately meets the assumption of one-dimensional Darcy flow. Thus, the parameters give clear physical meanings: the feature distance i d and j d are the distances between the center of each matrix element to the common interface; and the absolute permeability i k and j k are equal to the matrix permeability x k , y k or z k because the Cartesian matrix block is usually built in the direction of the coordinate system.   Figure 2b). Thus, if we still assume the pressure has a linear distribution (as shown by the green line in Figure 2b) to use the one-dimensional Darcy flow equation, non-negligible errors may be introduced when calculating M F T − . To avoid this, a common idea is to further subdivide the matrix meshes in the vicinity of the fracture to approximate the cone pressure distribution, which is why the DFM needs denser matrix meshes around the fracture.
EDFM chooses not to add extra matrix elements to maintain efficiency, and still follows the assumption of a linear distribution of pressure, which may introduce errors, especially when the For F-F connections, the physical meanings of parameters are also clear for the same reason: d i and d j are the shortest distances from the center of each fracture element along the fracture direction to the common interface; k i and k j are fracture permeability k F .
For M-F connections, T F is much larger than T M . Thus, according to Equation (10), T M−F is approximately equal to T M . Because there is a large difference between the permeability of the matrix and the fracture, the pressure inside the matrix element domain in the vicinity of the fracture nearly has a cone distribution (as shown by red line in Figure 2b). Thus, if we still assume the pressure has a linear distribution (as shown by the green line in Figure 2b) to use the one-dimensional Darcy flow equation, non-negligible errors may be introduced when calculating T M−F . To avoid this, a common Energies 2020, 13, 3070 6 of 21 idea is to further subdivide the matrix meshes in the vicinity of the fracture to approximate the cone pressure distribution, which is why the DFM needs denser matrix meshes around the fracture.
EDFM chooses not to add extra matrix elements to maintain efficiency, and still follows the assumption of a linear distribution of pressure, which may introduce errors, especially when the matrix mesh is coarsely divided or the fracture needs to be subdivided due to its complex geometry [27,30].
In iEDFM, neither the matrix nor the fracture is required to be subdivided. Furthermore, T M−F is solved based on a more realistic pressure distribution derived by a semi-analytic algorithm for anisotropic formation, which will be discussed in the next section.

Algorithm for the Embedding of Fractures into Anisotropic Formation
The complexity of fracture geometry and anisotropy of the matrix can strongly influence or even dominate the matrix-fracture fluid exchange. Therefore, calculation methods for T M−F are the key point to embed fractures into the anisotropic formation.
To calculate T M−F , an analytic solution of pressure distribution around a point source in an anisotropic formation is introduced. Then, we assume that there are a series of point sources distributed along the fractures. With the flow strength of each point source calculated, the semi-analytic pressure distribution inside the anisotropic matrix domain can be obtained by integrating the analytic solution of all the point sources. Finally, the pressure distribution is used to compute the value of T M−F .
In this section, the algorithm for 2D cases (gravity ignored) in which the coordinate axis direction coincides with the principal axis direction of the permeability are both considered.

Analytic Point-Source Solution in Anisotropic Formation
As shown in Figure 3a, when the coordinate axis direction coincides with the principal axis direction of the permeability, k xy = k yx = 0. The flow strength of the point source can be obtained as (Appendix A presents the detailed derivation): where h is the thickness of the reservoir.   The steady-state equation of a point source is: Energies 2020, 13, 3070 7 of 21 By introducing an anisotropy coefficient β = k x /k y , we can transform the problem of anisotropic formation into an equivalent form for isotropic formation as: where: The inner boundary condition is: Thus, the solution of the steady-state equation can be obtained as: As shown in Figure 3b, when an angle θ k 0 exists between the coordinate axis direction (xoy) and the principal axis direction of the permeability (ζoη), k xy = k yx 0. By the coordinate transformation procedures (Appendix B presents the detailed derivation), the point-source solution can be obtained as:

Calculation Methods for Matrix-Fracture Transmissibility
We assume there are a series of point sources S i x S i , y S i (i = 1 to N, normally N > 50) distributed along the fractures to calculate the pressure distribution inside the anisotropic matrix domain.
Based on the superposition principle of potential, the pressure of any point X(x, y) inside the matrix domain can be obtained from the integral of Equations (17) or (18): where q S i is the flow strength of point-source S i ; k e is the equivalent permeability; and r e,iX is the equivalent distance between this point X to point source S i .
If the coordinate axis direction coincides with the principal axis direction of the permeability, we have k xy = k yx = 0 and the Equations (20) and (21) can be written as: Energies 2020, 13, 3070 As shown in Figure 4b, a series of reference points F j (j = 1 to N, the same number as for point sources) on the surface of the fracture can be used to determine the unknown parameters q S i and C in Equation (19).
where , e ij r is the equivalent distance between the point source Si and the reference point j F .
(a) (b) If we define: where i ξ can be solved out from the linear system Equation (24), then, the pressure of any point X ( ) , x y inside the matrix domain can be written as: In addition, we have: 1 where M V is the volume of this matrix element; x dV is an element volume at point X.
Combining Equations (26) and (27) and the definition equation of transmissibility (Equation (8)), the M-F transmissibility and the pressure anywhere inside the matrix domain can be calculated by: where M ∈ means that the point is inside this matrix domain and: 1 , From Equations (28) (22) and (23) or Equations (20) and (21). If these reference points are on the surface of a same fracture element, they could be assumed to have the same pressure P F . Then, Equation (19) can be established for every reference point, forming an N-dimensional linear equation system: where r e,ij is the equivalent distance between the point source S i and the reference point F j . If we define: where ξ i can be solved out from the linear system Equation (24), then, the pressure of any point X (x, y) inside the matrix domain can be written as: In addition, we have: where V M is the volume of this matrix element; dV x is an element volume at point X.
Combining Equations (26) and (27) and the definition equation of transmissibility (Equation (8)), the M-F transmissibility and the pressure anywhere inside the matrix domain can be calculated by: Energies 2020, 13, 3070 where ∈ M means that the point is inside this matrix domain and: From Equations (28)-(30), we note that T M−F is only related to the properties of the reservoir and can be determined at the step of pre-processing before the numerical simulation starts.
Regardless of whether the coordinate axis direction coincides with the principal axis direction of the permeability, T M−F can be calculated with k e and r e,iX in Equations (22) and (23) or Equations (20) and (21).

Case 1: Single-Phase Flow in Anisotropic Formation with Vertical Fractured Well
In a rectangular anisotropic formation of single-phase fluid, the analytic solution of a vertical fractured well at a fixed production rate under a quasi-steady state can be derived.
In this paper, numerical results of a vertical fractured well case are presented and compared with the analytic solution. Various fracture azimuths, anisotropy coefficients, and fracture permeabilities are considered to illustrate the accuracy of iEDFM.

Analytic Solution of Quasi-Steady-State Flow
As shown in Figure 5, the dimensions of the quasi-three-dimensional rectangular reservoir are x e × y e × h, the half-length of the vertical fracture is x f , and the aperture is w f . A vertical well with constant production q is located at the midpoint of the fracture (called a vertical fracture well). The coordinate axis direction coincides with the principal axis direction of the permeability, and the fracture azimuth is θ.

Case 1: Single-Phase Flow in Anisotropic Formation with Vertical Fractured Well
In a rectangular anisotropic formation of single-phase fluid, the analytic solution of a vertical fractured well at a fixed production rate under a quasi-steady state can be derived.
In this paper, numerical results of a vertical fractured well case are presented and compared with the analytic solution. Various fracture azimuths, anisotropy coefficients, and fracture permeabilities are considered to illustrate the accuracy of iEDFM.

Analytic Solution of Quasi-Steady-State Flow
As shown in Figure 5  If the formation is isotropic, the analytic solution of the entire reservoir pressure field can be expressed as [35,39,40]:  If the formation is isotropic, the analytic solution of the entire reservoir pressure field can be expressed as [35,39,40]: where p avg (t) represents the average pressure of the entire rectangular reservoir at time t; p w f represents the pressure of the vertical well.
In Equation (31), G L (g wr ; g w ) is a two-point function to characterize the effect of the geometric features, such as the size of the calculation area and the distribution of fractures; the first variable g wr represents the coordinates of the pressure measurement point; the second variable g w represents the coordinates of the vertical fracture midpoint, that is also the vertical well position. Gringarten [41] proved that the pressure drop at 0.732 half-length of a fracture is equal to the pressure drop of the entire fracture, so the point x w + 0.732x f cos θ, y w + 0.732x f sin θ is taken as the pressure measurement point. f c f D is a function to characterize the effect of the fracture permeability. The modified coefficient c f D reflects the influence of relative fracture permeability. The greater the fracture permeability relative to the matrix, the larger the value of c f D .
The detailed calculations for the parameters in Equation (31) are given in Appendix C.
To compare the analytic solution with the numerical results of iEDFM, a dimensionless production index J D is defined as [35]: Spivey and Lee [37] equate anisotropic formation with isotropic formation by coordinate transformation. Parameter changes caused by this coordinate transformation are shown in Table 1. Table 1. Parameter changes from anisotropic formation to equivalent isotropic formation.

Anisotropic Formation Equivalent Isotropic Formation
After coordinate transformation, the dimensionless production index J D is written as [35]:

Comparison between Numerical and Analytic Solutions
A numerical case of iEDFM is established referring to Figure 5 with water production at a fixed rate in a micro-compressible rectangular reservoir. We calculate the production index at the pseudo-steady state from the numerical result, and compare it with that calculated from Equation (33).
The basic parameters of the numerical example are shown in Table 2. The remaining parameters can be calculated from different anisotropy coefficients, fracture azimuths, and modified coefficients (fracture permeability) used in this case.

Reservoir
Size   (a) (b)      Figure 8 shows the pseudo-steady-state pressure distribution with fracture azimuth θ = π/3 and modified coefficient c f D = 500. Cases under different anisotropy coefficients are presented. In Figure 8, the injection times are 5 × 10 5 s (Figure 8a) and 1 × 10 6 s ( Figure 8b); both times represent a pseudo-steady state.
As shown in Figures 6 and 7, the production index J D increases with c f D , namely, the fracture permeability. When the fracture permeability increases to a certain extent, J D remains flat.  When the fracture azimuth θ = π/4, the numerical solutions of β = 1/4 and β = 4 coincide, and the numerical solutions of β = 1/2 and β = 2 also coincide. It can be seen that the entire reservoir has symmetry at this time, and the numerical solution simulates this symmetry.
When fracture azimuth θ = π/3, the numerical solution of the production index J D with β = 4 is greater than that with β = 1/4, and the numerical solution with β = 2 is greater than that with β = 1/2. Because the fracture is biased towards the dominant principal direction of matrix permeability (that is, the y-direction when β > 1), the flow in that direction would be further strengthened and the production would decrease.
In this numerical case, the iEDFM simulation results agree well with the analytic solution in the anisotropic rectangular formation for a single-phase quasi-steady-state problem with a fixed-rate production vertical fracture well, which demonstrates the reliability of iEDFM.

Case 2: Two-Phase Flow in Anisotropic Formation with Two Crossed Fractures
For the iEDFM of anisotropic formations, matrix-fracture transmissibility is the key issue, as illustrated in Section 2.3.2. In order to investigate the effects of anisotropy on matrix-fracture transmissibility, a water flooding case with a set of intersecting fractures was simulated. In this case, matrix permeability in the y-direction was set to be three times the permeability in the x-direction (β = √ 3/3). A fixed rate water injection well (as the blue arrow shows) and a fixed pressure producing well (as the red arrow shows) are implied in this case. The sketch and meshes for iEDFM are shown in Figure 9. The parameters for case 2 are shown in Table 3. Table 3. Parameters for case 2.

Reservoir
Size transmissibility, a water flooding case with a set of intersecting fractures was simulated. In this case, matrix permeability in the y-direction was set to be three times the permeability in the x-direction . A fixed rate water injection well (as the blue arrow shows) and a fixed pressure producing well (as the red arrow shows) are implied in this case. The sketch and meshes for iEDFM are shown in Figure 9. The parameters for case 2 are shown in Table 3.   In order to further demonstrate the reliability of the iEDFM for anisotropic problems, a comparison with the fine grid result was carried out. In fine grid simulation, the matrix mesh was set as 540 540 × and fractures were considered as discretized elements with higher permeability. The same parameters and well settings were used in both iEDFM and fine grid simulations. As shown in Figure 10, the oil production result of iEDFM agrees well with the fine grid result, which shows the reliability and accuracy of the iEDFM for anisotropic formations. In addition, the calculation time of iEDFM was about 12 s, which demonstrates its advantage of calculation efficiency compared with the 133 min required for the fine grid simulation. In order to further demonstrate the reliability of the iEDFM for anisotropic problems, a comparison with the fine grid result was carried out. In fine grid simulation, the matrix mesh was set as 540 × 540 and fractures were considered as discretized elements with higher permeability. The same parameters and well settings were used in both iEDFM and fine grid simulations. As shown in Figure 10, the oil production result of iEDFM agrees well with the fine grid result, which shows the reliability and accuracy of the iEDFM for anisotropic formations. In addition, the calculation time of iEDFM was about 12 s, which demonstrates its advantage of calculation efficiency compared with the 133 min required for the fine grid simulation.   Figure 11a,b. At this moment, since the displacement front has just reached the two fractures, the oil saturation in the matrix after the displacement front is similar in both cases, but the flow in the fractures appears to be more prominent when anisotropy is not considered in the M-F transmissibility. In order to investigate the effect of anisotropy on M-F transmissibility, case 2 was simulated under two different circumstances. The anisotropy was considered when calculating the flow between M-M connections under both circumstances, while cases with and without consideration of anisotropy when calculating M-F transmissibility were compared.
The oil saturation distributions at 1 × 10 7 s are compared in Figure 11a,b. At this moment, since the displacement front has just reached the two fractures, the oil saturation in the matrix after the displacement front is similar in both cases, but the flow in the fractures appears to be more prominent when anisotropy is not considered in the M-F transmissibility. The oil saturation distributions at 7 1 10 × s are compared in Figure 11a,b. At this moment, since the displacement front has just reached the two fractures, the oil saturation in the matrix after the displacement front is similar in both cases, but the flow in the fractures appears to be more prominent when anisotropy is not considered in the M-F transmissibility.  The oil saturation distributions at 5 × 10 7 s are shown in Figure 11c,d. At this moment, the profiles of oil saturation are rather different between both cases. Fluid flows faster in fractures, which causes the oil saturation around the ends of the fractures to be smaller when M-F transmissibility is calculated without the consideration of anisotropy. The displacement front also progresses further when anisotropy is not considered. These results indicate that if anisotropy is not considered in the calculation of M-F transmissibility, fluid exchange between the matrix and the fracture will be exaggerated, which leads to erroneous simulation results. Figure 12 shows the comparison of oil rate and oil production between cases with and without consideration of anisotropy when calculating M-F transmissibility. Both oil rate and oil production appear to be smaller with anisotropy is considered, also demonstrating that ignoring the anisotropy when calculating M-F transmissibility will exaggerate fluid exchange between the matrix and the fracture.
Compared with the fine grid result, case 2 further demonstrates the reliability of iEDFM in two-phase flow in anisotropic formations, and also presents the obvious effect of anisotropy on matrix-fracture transmissibility. exaggerated, which leads to erroneous simulation results. Figure 12 shows the comparison of oil rate and oil production between cases with and without consideration of anisotropy when calculating M-F transmissibility. Both oil rate and oil production appear to be smaller with anisotropy is considered, also demonstrating that ignoring the anisotropy when calculating M-F transmissibility will exaggerate fluid exchange between the matrix and the fracture.
(a) (b) Figure 12. Comparison of (a) oil rate and (b) oil production between cases with and without consideration of anisotropy when matrix-fracture transmissibility is calculated.
Compared with the fine grid result, case 2 further demonstrates the reliability of iEDFM in twophase flow in anisotropic formations, and also presents the obvious effect of anisotropy on matrixfracture transmissibility.

Model Application
As shown in Figure 13, a quasi-three-dimensional reservoir is considered. The fracture distribution and fracture permeability are adapted from geological survey data of a real reservoir [41,42]. The reservoir thickness is 100 m; the reservoir boundary is enclosed; the equivalent permeability of the matrix x y k k k = is 10 −14 m 2 ; the fracture permeability is 10 −7 m 2 ; the initial oil saturation is 0.8; the injection well is located at the center of the reservoir and injects water at a constant rate of 4 × 10 4 m 3 /day; and four production wells are located at the corners of the reservoir producing at initial pressure. Figure 12. Comparison of (a) oil rate and (b) oil production between cases with and without consideration of anisotropy when matrix-fracture transmissibility is calculated.

Model Application
As shown in Figure 13, a quasi-three-dimensional reservoir is considered. The fracture distribution and fracture permeability are adapted from geological survey data of a real reservoir [41,42]. The reservoir thickness is 100 m; the reservoir boundary is enclosed; the equivalent permeability of the matrix k = k x k y is 10 −14 m 2 ; the fracture permeability is 10 −7 m 2 ; the initial oil saturation is 0.8; the injection well is located at the center of the reservoir and injects water at a constant rate of 4 × 10 4 m 3 /day; and four production wells are located at the corners of the reservoir producing at initial pressure. The matrix of the reservoir is gridded as 25 × 25 × 1 blocks. The fracture network is discretized to 107 fracture elements, of which nine are intersected fracture elements and 22 are non-planer fracture elements.
After 7000 days of injection, the oil rate curve of production well-1, well-2, well-3, well-4 are shown in Figure 14a  The matrix of the reservoir is gridded as 25 × 25 × 1 blocks. The fracture network is discretized to 107 fracture elements, of which nine are intersected fracture elements and 22 are non-planer fracture elements.
After 7000 days of injection, the oil rate curve of production well-1, well-2, well-3, well-4 are shown in Figure 14a to 107 fracture elements, of which nine are intersected fracture elements and 22 are non-planer fracture elements.
After 7000 days of injection, the oil rate curve of production well-1, well-2, well-3, well-4 are shown in Figure 14a (c) (d) Figure 14. Oil production rate of each production well under different anisotropy coefficients. Figure 14. Oil production rate of each production well under different anisotropy coefficients.
It can be seen from Figure 14 that under different anisotropy coefficients, the breakthrough times of the production wells differ. Although the equivalent permeability remains the same, the oil production efficiency and final oil distribution are different, indicating that the anisotropy of the reservoir greatly affects the fluid flow.
The matrix is homogeneous, the injection well is located at the midpoint, and the production wells are symmetrically distributed. Therefore, if we ignore the fracture network, the oil distribution under β = 1/4 and β = 4 should have a symmetrical relationship. However, it can be seen from Figure 11 that when the fracture network is considered, the oil distribution under β = 1/4 and β = 4 is clearly different. In addition, as the average remaining oil saturation shows, the recovery ratio under β = 4 is much lower than that of the isotropic case and the recovery ratio under β = 1/4 is slightly higher than that of the isotropic case. is clearly different. In addition, as the average remaining oil saturation shows, the recovery ratio under 4 β = is much lower than that of the isotropic case and the recovery ratio under 1 4 β = is slightly higher than that of the isotropic case. These differences demonstrate that the existence of the fracture network can increase the influence of matrix anisotropy. Specifically, the injection water easily breaks into the fracture network, then tends to flow in the direction of higher permeability. If there are production wells in this direction, it is easy to form a breakthrough channel (such as the water channel along the top of Figure 15c). This improves the initial oil recovery, but makes it more difficult to recover the remaining oil. In contrast, if the direction from the fracture end to the nearest production well is not in the direction of lower permeability (such as the situation in Figure 15b), the initial oil recovery is decreased but the total recovery will increase. This characteristic could be significant in the wellsite selection procedure during development of an anisotropic reservoir.
This example shows the applicability of the iEDFM in a more complex situation, and also demonstrates its potential applicability in real field studies. These differences demonstrate that the existence of the fracture network can increase the influence of matrix anisotropy. Specifically, the injection water easily breaks into the fracture network, then tends to flow in the direction of higher permeability. If there are production wells in this direction, it is easy to form a breakthrough channel (such as the water channel along the top of Figure 15c). This improves the initial oil recovery, but makes it more difficult to recover the remaining oil. In contrast, if the direction from the fracture end to the nearest production well is not in the direction of lower permeability (such as the situation in Figure 15b), the initial oil recovery is decreased but the total recovery will increase. This characteristic could be significant in the wellsite selection procedure during development of an anisotropic reservoir.
This example shows the applicability of the iEDFM in a more complex situation, and also demonstrates its potential applicability in real field studies.

Conclusions
In this paper, an integrally embedded discrete fracture model (iEDFM) suitable for anisotropic formations is proposed. In this model, a structured mesh is employed for the anisotropic matrix system, and a fracture element consisting of several fracture segments or a group of connected fractures is integrally embedded into the matrix grid. A method to calculate the transmissibility between the matrix and the fracture for anisotropic formations is presented.
An anisotropic rectangular reservoir with a fixed-rate production vertical fracture well, a single-phase quasi-steady-state problem, is used to demonstrate the applicability and accuracy of the iEDFM. Different circumstances with various fracture azimuths, anisotropy coefficients, and fracture permeabilities are considered. The simulation results agree well with the analytic solutions. The effect of anisotropy on M-F transmissibility is also investigated. The advantage of the iEDFM is demonstrated by comparison with a fine grid model. Furthermore, a water flooding example with a modified reported dataset is conducted to analyze the effect of the matrix anisotropy, which also shows the feasibility of the iEDFM to address anisotropic formation with complex fracture networks, and demonstrates potential application of the iEDFM in real field studies.
Author Contributions: R.S. and Y.D. conceived of the presented idea. R.S. developed the theory and performed the verification case 1 and the application case. D.W. performed the verification case 2. Y.D. and Y.-S.W. guided the data analysis and supervised the work. This paper was written by R.S. and D.W., and modified by Y.D. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
As shown in Figure 3a, the velocity in the r direction can be obtained as: where k r is the directional permeability. We also have:

Appendix B
As shown in Figure 3b, the point-source solution in the ζoη coordinate system can be written as: From the geometric relationship, we have: ζ = x cos θ k + y sin θ k , η = −x sin θ k + y cos θ k (A8) Thus, the velocity component in the x-direction can be expressed as: In the xoy coordinate system, the velocity component in the x-direction can also be expressed as: By comparing Equations (A9) and (A10), we have: After a similar derivation, we have: k yy = k ζ sin 2 θ k + k η cos 2 θ k (A13) Bringing Equations (A11)-(A13) into Equation (A7), we can obtain: cosh u n y 1 +cosh u n y 2 u 2 nn sinh(u n y e /x e ) sin u n x 1 cos u n x x e − cosh u n y 3 +cosh u n y 4 u 2 nn sinh(u n y e /x e ) sin u n x 2 cos u n x