An Improved Embedded Discrete Fracture Model with Fracture Growth for Water-Induced Fracture Simulation in Low Permeability Reservoirs

Water injection is often used to supply formation energy in low permeability reservoirs. A lot of ﬁ eld data shows that fractures nearby the wellbores may extend hundreds of meters because of the poor injectability which is called water-induced fractures. The presence of these fractures signi ﬁ cantly changes the ﬂ ow of injected water thus a ﬀ ecting the development e ﬀ ects. However, few studies of reservoir numerical simulation methods have been presented to simulate the dynamic changes of fractures in the long-term water ﬂ ooding process. This study is based on the embedded discrete fracture model (EDFM) in which the fracture system has less dependence on the grid system. By improving the preprocessing algorithm, fracture growth can be considered in EDFM which is called the dynamic embedded discrete model (dEDFM). This new method provides an innovative idea for fracture propagation simulation by attaching new fracture elements to the original fractures once the failure criterion is satis ﬁ ed. Meanwhile, there is no need to calculate the stress ﬁ eld once the in situ stress is provided, so the computational e ﬃ ciency is greatly improved compared with the ﬂ uid-solid coupling method. Besides, dEDFM has a ﬂ exible way of handling fracture elements, so it is suitable for water-induced fracture simulation. Results show that fracture propagation speed is signi ﬁ cantly in ﬂ uenced by matrix permeability, in situ stress, and injection intensity, and if the water-induced fractures propagate faster than the original waterfront movement velocity, changing law of the water content rate rising of production wells in di ﬀ erent directions will be di ﬀ erent. The contribution of this work lies in that it provides a suitable and e ﬃ cient way for reservoir engineers to deal with water-induced fractures since the propagation of fractures can be considered when it comes to reservoir numerical simulation, and it is helpful for production data prediction and development adjustment in low permeability reservoirs.


Introduction
Low and ultralow permeability reservoirs play an important role in Ordos Basin where water injection development is one of the most commonly used methods to enhance oil recovery. Meanwhile, a lot of oilfield data such as production curves, water injection profiles, and field test data show a characteristic of directional water penetration after a period of exploitation, and water-induced fractures are proved to be widespread in these reservoirs [1,2]. In unconventional reservoirs, fractures were studied as the main channel of fluid seepage in many aspects [3][4][5][6]. Thus, traditional reservoir numerical simulation methods cannot satisfy the dynamic changes of fractures in low permeability reservoirs when it comes to water injection development. Hagoort et al. [7] first used a single-phase numerical simulator coupled with a fracture propagation pressure model to simulate fracture propagation in a symmetrical waterflood pattern; however, only the change in fracture length was studied in the mathematical model, and the pressure drop in the fracture was ignored. Gadde and Sharma [8] coupled a single-well model with the reservoir simulator to study the rate of fracture growth in which the effect of particle plugging and thermal stress were considered. Abou-Sayed and Zaki [9] focused on the influence of filter cake on the fracture propagation law. de Souza et al. [10] combined injection with fracture propagation pressure (IFPP) with a commercial reservoir simulator to evaluate the displacement efficiency of different vertical well patterns and horizontal well patterns; The Net Present Value (NPV) was also calculated considering different fracture behaviors to predict oilfield benefits. Besides, some scholars dealt with fracture propagation problems during the waterflooding process in a similar way, respectively [11][12][13]. In recent years, the discrete fracture model (DFM) was used to realize the dynamic changes of fractures [14,15] but was limited to the complexity of grid division. A dynamic fracture model framework was established in which the opening and closure of fractures were reflected by permeability change of grids [16,17]; however, the model oversimplifies the characterization of fractures. So far, analytical fracture models coupling with traditional reservoir numerical simulator were used to simulate water-induced fractures by most scholars; some innovative solutions such as dynamic DFM and dynamic fracture model were proposed in recent years, but shortcomings exist when dealing with fracture propagation.
EDFM was first proposed by Lee et al. [18][19][20], which can be seen as an improvement of the dual medium; EDFM can characterize fractures in a flexibly way because the fracture system and matrix system are independent, so it is more convenient for characterizing fractures compared with DFM in which fractures are characterized by unstructured grids. Later, dual continuum and EDFM were coupled to simulate fractures of different scales and EDFM was firstly extended to three-dimensional [21,22]. Tene et al. [23] proposed a projection-based embedded discrete fracture model (pEDFM) to deal with the condition when the fracture conductivity is less than the matrix conductivity. Jiang and Younis [24] improved pEDFM in terms of the limitation in a "naive implementation" of the original pEDFM method. Rao et al. [25] studied two-phase heat and mass transfer based on pEDFM, and Liu et al. [26] used pEDFM to appraise CO 2 sequestration capacity with complex boundary shapes in shale gas reservoirs. In traditional EDFM, fractures are considered to be static, but in recent years, EDFM was combined with the extended finite element method (XFEM) in order to simulate the propagation of fractures [27][28][29]. However, because of the introduction of the stress field, it is ineffective to simulate the long-term waterflooding process.
Our work is aimed at efficiently simulating waterinduced fractures during the waterflooding development process in ultralow permeability reservoirs; therefore, the dynamic embedded discrete fracture model (dEDFM) is established by improving the preprocessing algorithm of EDFM. In this paper, basic governing equations are introduced firstly; then, a detailed description of the process of dEDFM is presented. In the next section, the model is verified with Abaqus and an inverted-nine-point well pattern case is carried out by using dEDFM; detailed production data are analyzed simultaneously. Finally, we will give an evaluation of the dEDFM and a summary of the main results.

Methodology
2.1. Governing Equations. The model we proposed is based on the framework of EDFM in which each fracture is explicitly characterized. Similar to the dual medium model, we will give the matrix equation, the fracture equation, and the fracture-matrix mass transfer equation in the next section. The model considers the oil-water two-phase seepage flow with the assumption that the flow of each phase conforms to Darcy's law and the gravity term is ignored. Fluid and rock are slightly compressible, and capillary force can be considered in the model.
Similarly, the two-phase fluid flow in the fracture can be expressed as The additional equations are where K m represents matrix permeability; K f represents fracture permeability; k ro represents relative permeability of oil phase; k rw represents relative permeability of water phase; μ o represents oil viscosity; μ w represents water viscosity; B o represents oil formation volume factor; B w represents water formation volume factor; p o represents the pressure of oil phase; p w represents the pressure of water phase; p cow represents the oil-water capillary force; ϕ m represents matrix porosity; ϕ f represents fracture porosity; s o represents saturation of oil; s w represents saturation of water; t represents time; _ q omw represents flow strength (flow per unit time) of oil phase to matrix from the wellbore; _ q wmw represents flow strength of water phase to matrix from the wellbore; _ q of w represents flow strength of oil phase to fracture from the wellbore; _ q wf w represents flow strength of water phase to fracture from the wellbore; _ q omf represents flow strength of oil phase from fracture to matrix; and _ q wmf represents flow strength of water phase from fracture to matrix.  Lithosphere coefficient and its geometric factor. Take the oil phase as an example, the flow coefficient and geometric factor for all conditions (matrix grid and matrix grid, matrix grid and fracture element, and fracture element and fracture element) can be expressed as where λ ij represents the flow coefficient between element i and j; k roðupstreamÞ represents oil relative permeability of the element with higher pressure; μ o,i and μ o,j represent the oil viscosity of i and j element; B o,i and B o,j represent the oil formation volume factor of i and j element; G ij represents the geometric factor between element i and j; and G i and G j represent the geometric factor of i and j element. G i can be expressed in a general term as where K i represents permeability of i element; A represents the contact area; and hd i i represents the distance between two elements. For two adjacent rectangular grids or fracture elements, hd i i is the distance between the center of element i with the contact surface; when a fracture element is connected to a matrix grid, Equation (9) is used to calculate the average normal distance by [20] where S represents the area of matrix grid i. However, because of the relatively low accuracy in traditional EDFM, a novel Boundary Element Method (BEM) model with higher early accuracy is used in our model to calculate the geometric factor instead of the average normal distance [30]. Equations (10) and (11) give the final expression.
where q omf represents the flow rate of oil phase from fracture to matrix; q wmf represents the flow rate of oil phase from fracture to matrix; nf represents the number of fracture elements; and h represents matrix grid thickness. The forms of ½G ff , ½Gf, ½G, ½G fm , and ∂G/∂n are given in the appendix.

Well
Equations. The Peaceman model is used to calculate the flow rate between wellbore and matrix or between wellbore and fracture [31]. Take the oil phase as an example, when the well is connected with matrix grid, r e = 0:28 When the well is connected with fracture, where r e represents the equivalent supply radius; r w represents wellbore radius; and w f represents fracture width.

Discretization and Solution of
Equations. In our model, the reservoir is generated by the rectangular grid system and the finite volume method is used for numerical discretization. Take Equation (1) as an example, we obtain the discrete form as where TI ij can be calculated by Equations (5)-(7).

dEDFM
Method. dEDFM is the dynamic improvement of EDFM for the purpose to simulate water-induced fracture in a long-term water injection process. The main idea of dEDFM is to add dynamic fracture elements once the fracture criterion is met. Here, we will give a detailed process description of dEDFM, a single fracture driven by fluid pressure is studied in the following part: (1) Firstly, assume that the fracture propagates in the direction of the maximum principal stress. The preprocessing program includes grid division and fracture division, Figure 1 shows the schematic of a grid-fracture system in which fractures are cut by grid boundary for BEM discussed in Equations (10) and (11). To simulate fracture propagation, a slight fracture is connected to the injection well to initialize the water-induced fracture as is shown in Figure 2(a). The initial fracture should be short enough to affect the seepage field at the lowest degree (2) Next, the model is calculated according to the input parameters, if the convergence condition is satisfied after a time step, we will make a judgment whether the fracture will propagate, considering a hydraulic fracture under the force of in situ stress and internal pressure as is shown in Figure 3. According to rock mechanics, when a hydraulic tension fracture is at the point of propagating, the fluid pressure at the fracture tip should overcome the normal stress vertical to the fracture surface and the rock cohesion force; the fracture extension pressure is calculated by Equations (15) and (16) [32]. In some cases, the simplified form is used as is shown in Equation (17). If the propagate condition is satisfied, a new fracture element is connected to the initial fracture to characterize the propagation of the fracture which is called the dynamic fracture element. Figure 2(b) shows the situation when dynamic fracture elements are introduced. If not, the model will continue with the operation until Equation (15) is met or time up (3) For the propagating fracture, the red fracture section in Figure 2(b) has internal pressure higher than the  Lithosphere fracture pressure (after a certain injection time) while the green fracture sections are the dynamic fracture elements with the pressure of initial formation pressure which is much lower than the fracture pressure. During the next certain operation time, the injected water will continuously flow from the red part of fracture to the green part of fracture; meanwhile, pressure will rise until the tip pressure exceeds the fracture pressure in the whole fracture as is shown in Figure 2(c) (4) Steps (1)-(3) are the explanation of one cycle of fracture propagation; for the next propagation, steps (2)-(3) will be repeated. Figure 4 gives the framework of the whole dEDFM process. It is worthwhile to note that, when a dynamic fracture element is connected to the original fracture using EDFM, two types of fracture have the same coordinate at one point which results in the generation of "virtual units." These "virtual units" have no fracture length so the flow between fracture units will be blocked. Only by eliminating these "virtual units" can the fracture continue to propagate where p f represents fluid pressure along the fracture; l f represents fracture half-length; σ y represents normal stress vertical to the fracture surface; U represents rock surface energy; E represents elastic modulus; ν represents Poisson's ratio; s l represents rock tensile strength; σ 1 represents the maximum principal stress; σ 3 represents the minimum principal stress; and θ represents the angle between the direction of fracture and the maximum in situ stress.

Model Validation
In this section, dEDFM is compared with a single fracture propagation model simulated by Abaqus using a cohesive module. The grid number of dEDFM is 50 × 50 × 1ðx × y × zÞ, and different grid size of dEDFM is considered.
The main parameters are listed in Table 1. Assume that the fracture propagates in a fixed grid length (each time fracture will grow a length of dy). Because the change of fracture length is the main concern in water injection development, the fracture width is considered to be constant in dEDFM.  The fracture length-time curve for the fixed propagation length is verified as Figure 5 shows. Results show that fracture propagation law has a consistent trend in different conditions and the fitting effect is better when the grid size is smaller.

Numerical Examples
In this section, three kinds of main factors that affect fracture propagation including matrix permeability, in situ stress, and injection intensity will be discussed. Simultaneously, the influence of the dynamic change of fracture on production indexes is studied in an inverted nine-spot pattern. Figure 6 shows the specific well location diagram. The red points represent side wells in the horizontal direction; the green points represent side wells in the vertical direction; the orange points represent corner wells. We assume the maximum principal stress direction is horizontal and the minimum principal stress direction is vertical so that the water-induced fracture propagates horizontally in the following numerical cases.

Matrix Permeability.
In this case, we will discuss how matrix permeability will affect fracture propagate speed thus affecting the water content ratio of different production wells in the inverted nine-spot pattern as is shown in Figure 6. Three levels of matrix permeability (5 mD, 1 mD, and 0.5 mD) are chosen for numerical simulation. The model scale is 100 × 100 × 1 in a rectangular grid system, the grid size is 5 m × 5 m × 5 m. In this condition, we set an injection-production spacing of 200 m, the remaining main parameters are listed in Table 2.
The pressure field and water saturation field of dEDFM and EDFM are compared in the case of different matrix permeability, which is shown in Figures 7-12. Meanwhile, Figure 13 shows the fracture length and water content rate of the horizontal well change with time in different conditions. Results show that fracture propagates faster when the matrix permeability is lower but propagates slower or even no propagation when the matrix permeability is higher as is shown in Figure 13. There are two main reasons: Firstly, the bottom hole pressure (BHP) is high in a low permeability reservoir because of the poor injectability for the case of constant liquid injection, which is shown in Figures 7, 9, and 11; the main seepage channel is the fracture which is connected to the injection well. Secondly, mass transfer between fracture and matrix is limited when the matrix permeability is low, so injected water is more inclined to flow along with the fractures instead of vertical to the fracture. Both of the above two reasons will lead to a higher pressure at fracture tips which will finally cause a faster propagation of the fracture. Besides, the fracture propagates faster in the early stage but slower when the fracture propagates to a certain length; this is because the pressure at fracture tips is easy to rise when the fracture is short; however, as the fracture propagates, the accumulative total mass transfer is larger both because the larger fracture surface and the lower matrix pressure far from the injection well. From Figures 8, 10,   10 Lithosphere and 12, we can see that waterfront advanced velocity is anisotropic because of the propagation of fracture compared with EDFM; the water cut of the side well in the horizontal direction increases more rapidly while a little slower when it comes to the well in the vertical direction and the corner well as is shown in Figure 14. There is one more thing to emphasize that though fracture propagation speed affects water cut law to a great extent, matrix seepage is also an important factor that we can not ignore. For example, if the matrix permeability exceeds a certain value, the original waterfront will propagate faster than the fracture, and matrix seepage becomes the main factor instead of fracture seepage at this moment.

In Situ Stress.
In this case, we will discuss the influence of in situ stress on fracture propagation. Three levels of minimum principal stress in situ stress are chosen for simulation under the matrix permeability of 1 mD, and the other numerical model paraments are the same as case 1. The main parameters are listed in Table 3. Figure 15 shows the water saturation field in case of different minimum principal stress at the time point of 1 year, and Figure 16 shows the fracture length and water content rate of the side wells in the horizontal direction change with time in different conditions. Results show that fracture propagates faster than the waterfront when the minimum principal stress is small as is shown in Figure 15(a) but at a lesser degree when the minimum principal stress is larger as is shown in Figures 15(b) and 15(c); this is because that the fracture needs higher inner pressure to resist in situ stress as is discussed in Equations (15) and (16) when the minimum principal stress is large. Similarly, the fast propagation of the fracture leads to the early breakthrough of the side wells in the horizontal direction.

Injection Intensity.
In this case, we will discuss the influence of injection intensity on fracture propagation. Three levels of injection intensity are chosen for simulation under the matrix permeability of 1 mD and minimum principal stress of 26 MPa. The other numerical model parameters are the same as case 1 and case 2. The main parameters are listed in Table 4. Figure 17 shows the water saturation field in case of different injection intensity stress at the time   Figure 18 shows the fracture length and water content rate of the side wells in the horizontal direction change with time in different conditions. Results show that when the injection intensity is high, the original waterfront advance velocity and fracture propagation speed both increase which is different from case 2. When the fracture propagation speed is lower than original waterfront advance velocity as is shown in Figure 17(a), little influence is caused to the production index. However, when the fracture propagation speed exceeds original waterfront advance velocity as is shown in Figures 17(b) and 17(c), the side wells in horizontal direction breakthrough at an earlier time is shown in Figure 18.

Conclusion
In this paper, we propose a novel efficient numerical simulation method based on EDFM to simulate water-induced fracture in a long-term water injection development process, and then, the main factors that affect fracture propagation are discussed. Some main conclusions are summarized as the following: (1) dEDFM is the dynamic improvement of EDFM; fracture propagation is realized by introducing dynamic fracture elements and controlled by the discriminant module. dEDFM provides a good idea for the long-term water injection simulation when water-induced fractures should be considered and deals with tensile fracture propagation problems without calculating the stress field which greatly improves the calculation efficiency (2) Because of the simplified treatment of fractures, EDFM has less accuracy compared with DFM. In this condition, a novel BEM model is introduced to improve the early accuracy of the fracture-matrix mass transfer process  12 Lithosphere (3) Fracture propagation usually happens when the matrix permeability is small and propagates faster when the matrix permeability is smaller because of poor injectability and less fracture-matrix mass transfer. The propagation of fracture may lead to the directional water breakthrough which depends on the fracture propagate speed and waterfront seepage speed. The minimum principal stress and injection intensity are two other main factors that affect fracture propagation. The smaller of minimum principal stress and the larger of injection intensity, the faster that fracture propagates