Streamline simulation of water-oil displacement in a heterogeneous fractured reservoir using different transfer functions

. Main parts of oil and gas reserves are stored in fractured reservoirs. Simulation of multiphase ﬂ ow in fracturedreservoirsrequiresalargeamountofcalculationsduetothecomplexity,reservoirscaleandheterogeneityoftherockproperties.Theaccuracyandspeedofthestreamlinemethodforsimulatinghydrocarbonreservoirsat ﬁ eld scale make it more applicable than conventional Eulerian simulators using ﬁ nite difference and ﬁ nite element techniques. Conventional simulators for fractured reservoirs consume a great deal of time and expense and require powerfulCPUslikesupercomputers.Thismakesthedevelopmentofafast,powerfulandprecisesimulationmethodofgreatimportance.Thepresentstudywasundertakentodevelopacomputationalcodeasastreamlinesimulatortostudywater ﬂ ooding in a two-dimensional fractured reservoir with heterogeneous permeability using the Dual Porosity-Single Permeability (DPSP) model. In this simulator, the pressure equation is solved implicitly over an EuleriangridandthenthestreamlinesaregeneratedusingPollock ’ ssemi-analyticalmethodandaretraced.Atthis point,theTime-Of-Flight(TOF) isdevelopedandthesaturationequationsaremappedandsolvedexplicitlyalong thestreamlines.Next,theresultsaretransferredbackontotheEuleriangridandthecalculationsarerepeateduntilthesimulationendtime.Infracturedreservoirs,theinteractionbetweenthematrixandfractureisincludedinthetransferfunctions.Transferfunctionsmodel uid ow and production mechanisms between the matrix and fracture. They introduce source/sink equations between the matrix and fracture and they are distributed throughoutthemedia.Inthecurrentstudy,a problemissimulatedusingstreamlinemethodandseveralimportant transfer functions. A new linear transfer function with a constant coef ﬁ cient is introduced that is based on differences in water saturation between the matrix and fracture. The simulation results were then compared and a commercial software is applied to solve the same problem. The results of the streamline simulator were compared withthoseofthecommercialsoftwareandshowedappropriateaccuracyforthenewlyintroducedtransferfunction.Theaccuracyandef ﬁ ciencyofthestreamlinesimulatorforsimulationoftwo-phase ﬂ owinfracturedreservoirsusing different transfer functions are con ﬁ rmed and the results are veri ﬁ ed.


Introduction
Simulation of fluid flow in porous media in order to improve oil production from hydrocarbon reservoirs is vital. This includes proposing a Master Development Planning (MDP) for reservoir management (Siavashi et al., 2014).
Most countries that have a large oil and gas reserves are concerned with fractured reservoirs (Bourbiaux, 2010). Fractured porous media are characterized by substantial differences in the properties of the matrix and fracture regions (Noetinger et al., 2016). Recent research has shown that the ultimate recovery factor for a fractured reservoir is <10-60% or 70% (Bourbiaux, 2010). Simulation as a method of predicting fluid behavior using conventional simulators and the Eulerian approach is time-consuming and requires strong computing hardware (Biagi et al., 2016). Achieving an appropriate simulation method is essential. The streamline method is a powerful method with advantages and useful applications over conventional simulation methods based on finite volume and finite difference approaches. The streamline technique can be used to simulate hydrocarbon reservoirs at field scale with appropriate speed and accuracy (Batycky et al., 1997). In addition to its speed and precision, it has advantages which provide better flow visualization patterns and well intercommunication.
Fractured reservoirs can be divided into two connected media; the fracture network that provides the main path for fluid flow and matrix blocks, which are porous media and provide the bulk of reservoir volume and storage capacity. Three-dimensional images of fracture networks obtained by field fracturing data integration and stochastic modeling are highly complex and cannot be applied as direct input for a reservoir simulator (Bourbiaux et al., 1998). When fluid flow in reservoirs is simulated using mathematical models, the average properties of the porous media are used.
There are three general types of fractured reservoirs (Bourbiaux, 2010): Type A has a very tight matrix; thus, fluid flow and storage are not feasible in the matrix section. The permeability and porosity of the reservoir are determined by the fracture section (Single-Porosity Single-Permeability model (SPSP)). In type B, the matrix section has greater storage capacity for fluids and determines the main porosity of reservoir. The fluid flow rate is low in the matrix, and the matrix does not make a large contribution to the overall permeability of the reservoir (Dual-Porosity Single-Permeability (DPSP) model). In type C, the matrix section is porous and permeable and has an effective permeability similar to that of the fracture section (Dual-Porosity Dual-Permeability (DPDP) model).
The use of the SPSP, DPSP and DPDP models for simulation of fractured reservoirs are summarized in Figure 1 (Barenblatt et al., 1960;Di Donato et al., 2003;Datta-Gupta and King, 2007).
The DPSP model was used in this research because the permeability of matrix was considered to be small enough to ignore the Darcy flow (Jerbi et al., 2017). This model represents the matrix blocks as the source of the fracture network and fluid only flows in the fractures. This theory is reasonable because of the high-speed fluid flow in the fracture network with respect to the matrix blocks; thus, the fluid flow between the matrix blocks is neglected (Gilman, 2003).
A major challenge to simulating fractured reservoirs is the limited computing capacity of conventional simulators (Bourbiaux, 2010). The streamline method was derived from the streamtube method (LeBlanc and Caudle, 1971); however, the streamtube method has major problems and limitations. For example, changes in the rate and location of wells were not possible and, in most cases, the effect of gravity is ignored. The introduction and development of the streamlines concept solved many of these limitations (Baker, 2001).
Conventional simulation techniques based on the finite difference method have a two-step solution. In the first step, the pressure equation is solved and the pressure distribution in the domain is calculated. In the second step, the saturation equation is solved. In the streamline method, as in the IMplicit Pressure-Explicit Saturation (IMPES) reservoir simulation method, the pressure distribution is obtained implicitly over the Eulerian grid and then the streamlines are drawn using Pollock's semianalytical method. Afterwards, the transport equations are solved explicitly and linearly along the streamlines. Because the equations are one-dimensional along the streamlines, the number of calculations and, therefore, the amount of memory required for simulation is less than for conventional methods.
In conventional simulation methods, fluid flows between the simulation cells. In the streamline method, fluid flows in the direction of the pressure difference, which coincides with the physics of flow. The streamline method has two different timesteps. The total timestep (Dt total ) is used to solve the pressure equation implicitly. The local timestep (Dt sl ), or streamlines timestep, is used to solve the saturation equation explicitly along the streamlines. The implicit solution has no stability limitations in comparison with the explicit method when using large timesteps. Therefore, a larger total timestep size can be chosen to reduce the number of pressure solution updates and increase the simulation speed. These two timesteps are shown in Figure 2 (Doranehgard and Siavashi, 2018).
The increase in the speed of streamline method derives mainly from avoidance of the updating of the pressure field according to changes in the level of saturation. Another factor that increases the simulation speed is that every streamline is calculated separately, independent of other streamlines, with its own timestep size (Dt sl ).
The non-linear transport equations make the finite difference method very sensitive to changes in the size and orientation of the simulation cells and timestep (Baker, 2001). The streamline method is, therefore, more stable (Batycky et al., 1997). The streamlines never cross each other; thus, the amount of fluid produced is the sum of the streamlines . The merits of streamline method are as follows: high calculation speed and high accuracy; appropriate visualization of fluid flow and the relationship between the production and injection wells; better assessment of drained areas in the reservoir; high performance parallel processing option; help in history-matching.
The streamline method for field scale reservoirs is orders of magnitude faster than conventional simulators. In conventional simulators, increasing the number of simulation cells will almost exponentially increase the simulation run time. In the streamline method, this increase is nearlinear (Thiele, 2001); thus, this simulation method is more suitable for large scale models, as shown in Figure 3.

Governing equations for two-phase flow in porous media
The mass conservation equation for two-phase flow assuming a constant temperature is: where r a is density of a phase, S a is saturation of a phase, u ∝ is Darcy velocity of a phase (in vector form), q a is volumetric flow rate of a phase (for wells).
In fractured reservoirs fluid transfer between the matrix and fracture should be considered and is represented by the transfer function.
The Darcy velocity equation (momentum conservation equation) for the a phase in porous media  is: where m a is viscosity of a phase, P a is pressure of a phase, K ra is relative permeability of a phase, K is absolute permeability of rock, D is depth. The depth "D", which is equal to the distance between two geological layers in the z-direction and its unit for SI system is meter.
The pressure equation is derived by substituting Equation (2) into Equation (1) as: To find the pressure field in the streamline method, the saturation network should be considered constant and the pressure field should be calculated. In the DPSP model, pressure field calculation is only performed in the fracture part as calculation is not needed in the matrix part (Ahmadpour et al., 2016). The 2D pressure equation for two-phase flow of water and oil with the assumption of incompressible and immiscible fluids is as follows (Siavashi et al., 2014): where l t is the total mobility ratio, which is the summation of water and oil mobility ratios as: The water and oil mobility ratios (in the fracture) are expressed as: where K rwf and K rof are the water and oil relative permeability, respectively, in the fracture. When a fluid phase flows from matrix to fracture, the relative permeability of that phase depends on the saturation of the matrix cell (Lemonnier and Bourbiaux, 2010). In this research, the following equations were addressed (Brooks and Corey, 1964): where S wc is connate water saturation, S or is residual oil saturation.
For incompressible and two-phase flow (Ahmadpour et al., 2016): The equation of water saturation along a streamline for the matrix and fracture as a function of Time-Of-Flight (TOF), as shown in Equations (10) and (11), are as follows (Kazemi et al., 1976;Di Donato et al., 2003): Fig. 2. Schematic representation of total and local timesteps (Doranehgard and Siavashi, 2018).
TOF (t) (Batycky et al., 1997) is defined as the time required for a particle to travel distance S along a streamline traced from a source (e.g. an injection well) to a sink (e.g. a production well) and is written as : where ' is the porosity and d andũ are the location and velocity of the particle along the streamlines, respectively. This component separates the effects of geological heterogeneity from the transport equations (mass and momentum) and converts the 3D physical coordinates to multiple 1D TOF coordinates, allowing the transport equations to be solved simply. This is done by mapping the saturation equations from the 3D coordinates to the TOF coordinates. The multi-dimensional saturation equations are converted to a series of 1D equations for calculating variations in saturation along the streamlines. In addition, because only one streamline is loaded in the memory at any moment, the amount of memory required to perform the simulation is reduced significantly.
In Equation (10), f wf is the fractional flow of water in the fracture: The general steps of streamline simulation for a fractured reservoir using the DPSP model are as follows (Datta-Gupta and King, 2007;Hassane, 2013): 1) By using the numerical solution for the pressure equation and application of the Darcy law, the pressure distribution and fluid velocity for the fracture at one total timestep are achieved for the initial condition (for pressure and irreducible water saturation), porosity, permeability distribution and well and boundary conditions (Vitoonkijvanich et al., 2015).
2) The streamlines are drawn using Pollock's method (Pollock, 1988) considering the fluid velocity distribution in the fracture. It is assumed that the only connection between the injection and production wells is the fracture network because of the high fluid velocity of the fracture in comparison with the matrix.
3) When the streamlines are drawn, the TOF along them is calculated. The TOF contour locates and predicts the injected fluid front at the selected reservoir at different times and is a major advantage of the streamline method.
4) Fluid saturation and other fluid characteristics are transposed from solution grids to the streamlines. By solving the transport equations, fluid saturation can be calculated in matrix and fracture (Di Donato and Blunt, 2004). The transport equations are solved using a 1D technique along the streamlines in terms of TOF. In the streamline method, the explicit solution for the transport equations requires the use of a local timestep (Dt sl ) that is smaller than the total timestep (Dt total ) (Choobineh et al., 2015). 5) After solving the equations in a total timestep, the fluid saturation is transposed from the streamlines to the solution grids. A source of error in the streamline method occurs when transporting the solutions (fluid saturations) from the Eulerian grid to the streamlines and vice versa.
6) The operator-splitting technique is used for any effects running in different directions than the streamlines, such as the effect of gravity (Siavashi et al., 2014). 7) Repeat the previous steps until the simulation ends.

Matrix-fracture transfer functions
Transfer function defined as it specifically characterizes fractured reservoir flow behavior. From the dual porosity perspective, the fracture medium and the matrix medium act as flow regions in which the media interact. This interaction is defined mathematically by means of matrix-fracture transfer functions. The effect of the matrix-fracture mechanisms can differ and include variability of the flow properties in the matrix porous medium and the characteristics of the fracture network (Lemonnier and Bourbiaux, 2010). Several transfer functions for simulation of fractured reservoirs using dual porosity are discussed below. In these transfer functions, the effects of gravity and compressibility are neglected and fluid transfer between the matrix and fracture is controlled by imbibition process. The imbibition process plays a vital role in matrix-fracture fluid transfer (Sarma and Aziz, 2004). (Kazemi et al., 1976) The interaction (fluid transportation) between the matrix and fracture in the DPSP model is represented by the transfer functions (T i ). Kazemi et al. (1976) introduced the first multi-phase transfer function, known as the conventional transfer function, which is most commonly used in fractured reservoir simulation and many commercial simulator software packages (Gilman and Kazemi, 1983;Di Donato et al., 2003). They defined capillary pressure to introduce a transfer function that confirms the experimental results (Di Donato and Blunt, 2004).

Conventional transfer function
The simplified form of this transfer function that considers 2D Darcy flow (Gilman and Kazemi, 1983;Datta-Gupta and King, 2007) can be written as: where F s is the shape factor. These two equations are functions of the matrix and fracture phase saturation and pressure. The exchange rate between the matrix and fracture depends on the size of the matrix block. It considers the size of a single matrix block in the neighborhood of a single fracture for a dual-porosity model network (Jerbi et al., 2017).
The shape factor is a geometric coefficient used to calculate the matrix-fracture fluid transport coefficient as a function of fracture spacing (aperture). The shape factor is independent of time. Kazemi et al. used the following shape factor for the matrix-fracture transfer function for cubic blocks of a matrix (Gilman and Kazemi, 1983): where L is the length of the fracture in the x, y and z directions and is calculated experimentally and depends on the dimensions of the problem (Landereau et al., 2001). There is a pressure discontinuity at the interface of the two-phase flow of immiscible fluids, such as oil and water, that is caused by surface tension that leads to formation of curvature at the interface. The pressure of the wettingphase fluid is less than that of the non-wetting-phase fluid. This pressure difference is called capillary pressure that is denoted by P cm and P cf for the matrix and fracture, respectively (Chen et al., 2006).
Capillary pressure is a function of saturation. It also depends on the direction of change in saturation, meaning of drainage and imbibition. Drainage causes a decrease in wetting-phase saturation and imbibition causes an increase. The capillary pressure of the fracture is negligible for a low-permeable matrix block. The driving force for imbibition is the matrix capillary pressure. Imbibition continues until the difference in capillary pressure between the matrix and fracture falls to zero, at which water saturation will decrease to its minimal value and the water phase will become immobile. The capillary pressure in the matrix block is calculated as (Di Donato et al., 2003): where P max cm is the maximum capillary pressure that is depended on fluid and rock properties.
With the assumption of incompressible fluid: By ignoring the fracture capillary pressure, Equation (16) can be used to convert the conventional transfer function to (Al-Huthali and Datta-Gupta, 2004): Equations (10) and (11) can be discretized and rewritten in terms of local timestep (Dt sl ) as: In Equation (20), (S wf ) int denotes intermediate saturation and can be written as: This is a linear function of matrix saturation where b depends on the wettability and its value for a mixed-wet system is 100-10 000 times smaller than for strongly waterwet systems. Di Donato et al. (2003)

Proposed linear transfer function with constant coefficient
Generally, the transfer functions used in dual porosity simulators have a constant coefficient that is independent of time and another item, depending on forces such as fluid expansion, viscosity, capillarity, gravity and diffusion. Because the transfer function affects parameters such as the amount of oil recovered from the matrix block, the total recovery factor and amount of computations required, choosing an appropriate transfer function has a large effect on the performance of the selected reservoir simulator. The conventional transfer function is dependent on capillary pressure, which is a function of matrix water saturation. In this section, a new transfer function is introduced that is linearly proportional to the difference between the matrix and fracture water saturation values. This saturation difference causes water to be imbibed from the fracture to the matrix (capillary-controlled imbibition). This newly defined transfer function is simple and works with a constant coefficient using the saturation difference between the matrix and fracture: Parameter a is a constant coefficient which can be calculated from the experimental data and is independent of time. In the DPSP method, imbibition of water occurs from the fracture to the matrix and the drainage of oil is from the matrix to the fracture (counter-current imbibition). The reverse direction is not considered; therefore, the saturation difference between the matrix and fracture will change the transfer function.

Problem definition
Applying enhanced oil recovery methods, such as gas injection, water-alternating gas injection, chemical and thermal methods, are feasible for production from complex fractured reservoirs (Lemonnier and Bourbiaux, 2010). A 2D fractured reservoir sample with heterogeneous permeability is defined at this point to verify the streamline simulation model and compare the results of the transfer functions with those of commercial simulation software. In this sector model, a 5-spot pattern is considered with the injection well at the center and the production wells at the corners (Fig. 4). The injection rate is constant and the controlling parameter for the production wells is the bottomhole Dimensionless a = 1.18 Â 10 À8 1/day b = 1.5 Â 10 À8 1/day pressure (BHP). The sector model has 30 Â 40 Â 1 simulation cells and the total timestep size is 100 days. The fluid and reservoir properties are shown at Table 1. For permeability of the fracture, one layer of the SPE-10 model is used (Christie and Blunt, 2001). This heterogeneous permeability is shown in Figure 5. The simulation study was performed using the different types of transfer functions as introduced previously. The transfer function parameters were selected to allow comparison of the results. Equations (26) and (27) were used to evaluate the parameters: where S wm is average water saturation in matrix, S wðfÀmÞ is average water saturation difference between matrix and fracture. Table 2 shows the transfer functions parameters used in the simulation.   Overall, Figures 6-15 show the results of different transfer functions. It is clear that, in most cases, there are negligible discrepancies. The figures show that water saturation in the fracture is greater than in the matrix and can be attributed to the greater permeability of the fracture. The fluid front has same pattern for both the matrix and fracture.
The figures showing water saturation in the matrix reveal that the difference in saturation was high. In fact, one reason for the extreme saturation difference is the range selected for matrix saturation range (0.21-0.34). Figures 16-20 show the pressure distribution in the fracture.   It can be seen that the results of the commercial software were slightly different from those of the streamline method. With a decrease in the total timestep or size of the simulation grid cells, this difference becomes negligible. Figures 21-28 show the oil flow rate and cumulative volume of oil produced by production wells from the start of water injection for a period of 500 days. The curves of the different transfer function methods can be compared in these figures.
The results do not differ significantly. The greatest difference was for the results obtained from the Di Donato et al. transfer function and was dependent on the value of selected b.
The results of the proposed new transfer function are similar to that of the other solutions. The largest difference between the oil production rates obtained from the proposed transfer function and the commercial software response is 15% for well A at the first timestep. This difference gradually decreased in succeeding timesteps. Figure 29 shows the TOF contour and Figure 30 shows the streamlines after the first total timestep (100 days).
By investigation of the TOF grids (Fig. 29), the location of the fluid front at different times can be predicted. The streamlines grids (Fig. 30) represent the flow pattern and communication between the injection well and production wells. As mentioned, the direction of the streamlines is from the injection well to the production wells. The compaction of streamlines at one part of the sector-model shows the higher fluid velocity in that segment.
The TOF results shown in Figure 29 indicate that, after the injection of water, the breakthrough of water in well C occurred in less than 100 days and, after a short time, water breakthrough occurred in well A. The water front then reached wells B and D in 150 days. The streamlines in Figure 30 reveal that, the number of streamlines that end            in wells A and C are almost the same and are more than for wells B and D. This is consistent with the permeability of the fracture that is shown in Figure 5. Visualization of the streamlines also indicates that the flow path from the injection well to well A is more tortuous than the path from the injection well to well C and this is the reason for delayed breakthrough of the water front in well C. The results of the commercial software and streamline simulator based on the conventional Kazemi et al. transfer function should be the same because the applied transfer function in both simulators is the same. The oil flow rates reveal that the maximum differences between the commercial and streamline simulators for well A is 9% and that this difference decreases in the succeeding timesteps.

Conclusion
The current study simulated a heterogeneous fractured reservoir using several transfer functions for the streamline method. The results confirm the efficacy and capability of this simulation method. This method is equipped with streamlines and TOF and uses them to gain the following extra information over conventional simulation methods: the fluid flow pattern; the relationship between wells; analysis and interpretation of the flow path in the reservoir; location of injected fluid front (e.g. tracer injection) at different times.   The streamline simulation method is a complementary method for conventional simulators of hydrocarbon reservoirs. Comparison of the results of the presented transfer functions, which are all based on the imbibition mechanism, shows that they are similar. The differences observed relate to the amount of water saturation in the reservoir zones reached by the water front. The rate of change in water saturation in the linear transfer functions, especially for that of Di Donato et al., is greater than for that of Kazemi et al. This discrepancy does not create a large difference in the amount of produced oil. As shown in Figures 21-24, the greatest difference occurred in the first timestep. In succeeding timesteps, the results are closer and the graphs converge.
For single-layer heterogeneous fractured reservoirs, waterflooding was simulated using commercial software and an in-house streamline-based numerical code that uses the DPSP reservoir model. The Kazemi et al. List of symbols K permeability, m 2 q flow rate, m 3 S a a phase saturation, dimensionless p pressure, kpa f a fractional flow of a phase T transfer function, 1/s F s shape factor, 1/m 2 U Darcy velocity, m/s L a fracture distance in a direction P c capillary pressure, pa Greek letters r a density of a phase, kg/m 3 t time-of-flight Dt sl local timestep Dt t total timestep ' porosity l mobility ratio, 1/pa.s m viscosity, pa.s a, b transfer function constant coefficients, 1/s