Regular Article
Streamline simulation of wateroil displacement in a heterogeneous fractured reservoir using different transfer functions
^{1}
PhD Candidate of Petroleum Engineering, School of Chemical Engineering, Institute of Petroleum Engineering, University of Tehran,
Tehran, Iran
^{2}
Professor of Chemical Engineering, Institute of LNG, College of Engineering, University of Tehran,
Tehran, Iran
^{3}
Assistant Professor of Mechanical Engineering, Applied MultiPhase Fluid Dynamics Laboratory, School of Mechanical Engineering, Iran University of Science and Technology,
Tehran, Iran
^{*} Corresponding author: avatani@ut.ac.ir
Received:
17
November
2017
Accepted:
1
February
2018
Main parts of oil and gas reserves are stored in fractured reservoirs. Simulation of multiphase flow in fractured reservoirs requires a large amount of calculations due to the complexity, reservoir scale and heterogeneity of the rock properties. The accuracy and speed of the streamline method for simulating hydrocarbon reservoirs at field scale make it more applicable than conventional Eulerian simulators using finite difference and finite element techniques. Conventional simulators for fractured reservoirs consume a great deal of time and expense and require powerful CPUs like supercomputers. This makes the development of a fast, powerful and precise simulation method of great importance. The present study was undertaken to develop a computational code as a streamline simulator to study waterflooding in a twodimensional fractured reservoir with heterogeneous permeability using the Dual PorositySingle Permeability (DPSP) model. In this simulator, the pressure equation is solved implicitly over an Eulerian grid and then the streamlines are generated using Pollock's semianalytical method and are traced. At this point, the TimeOfFlight (TOF) is developed and the saturation equations are mapped and solved explicitly along the streamlines. Next, the results are transferred back onto the Eulerian grid and the calculations are repeated until the simulation end time. In fractured reservoirs, the interaction between the matrix and fracture is included in the transfer functions. Transfer functions model fluid flow and production mechanisms between the matrix and fracture. They introduce source/sink equations between the matrix and fracture and they are distributed throughout the media. In the current study, a problem is simulated using streamline method and several important transfer functions. A new linear transfer function with a constant coefficient 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 with those of the commercial software and showed appropriate accuracy for the newly introduced transfer function. The accuracy and efficiency of the streamline simulator for simulation of twophase flow in fractured reservoirs using different transfer functions are confirmed and the results are verified.
© M. Mesbah et al., published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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 timeconsuming 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. Threedimensional 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 (SinglePorosity SinglePermeability 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 (DualPorosity SinglePermeability (DPSP) model). In type C, the matrix section is porous and permeable and has an effective permeability similar to that of the fracture section (DualPorosity DualPermeability (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; DattaGupta 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 highspeed 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 twostep 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 PressureExplicit 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 onedimensional 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 (Δt_{total}) is used to solve the pressure equation implicitly. The local timestep (Δt_{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 (Δt_{sl}).
The nonlinear 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 (Siavashi et al., 2016). 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 historymatching.
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.
Fig. 1 Schematics of different fracture − models for simulating fractured reservoirs (DattaGupta and King, 2007). 
Fig. 2 Schematic representation of total and local timesteps (Doranehgard and Siavashi, 2018). 
Fig. 3 Simulation run time ratio for various number of grids for streamline and conventional (IMPES) method. 
2 Governing equations for twophase flow in porous media
The mass conservation equation for twophase flow assuming a constant temperature is: $$\frac{\partial}{\partial t}({\rho}_{\alpha}{S}_{\alpha}\mathrm{\varnothing})+\nabla \cdot ({\rho}_{\alpha}{\displaystyle \stackrel{\rightharpoonup}{{u}_{\alpha}}})={\rho}_{\alpha}{q}_{\alpha}\text{,}$$(1) where ρ_{α} is density of α phase, S_{α} is saturation of α phase, ${{\displaystyle \overrightarrow{u}}}_{\propto}$ is Darcy velocity of α phase (in vector form), q_{α} is volumetric flow rate of α 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 α phase in porous media (Siavashi et al., 2016) is: $$\stackrel{\rightharpoonup}{{u}_{\alpha}}}=\frac{K{K}_{r\alpha}}{{\mu}_{\alpha}}(\nabla {P}_{\alpha}{\rho}_{\alpha}g\nabla D)\text{,$$(2) where μ_{α} is viscosity of α phase, P_{α} is pressure of α phase, K_{rα} is relative permeability of α 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 zdirection and its unit for SI system is meter.
The pressure equation is derived by substituting Equation (2) into Equation (1) as: $$\frac{\partial}{\partial t}({\rho}_{\alpha}{S}_{\alpha}\mathrm{\varnothing})+\nabla \cdot \left({\rho}_{\alpha}\frac{K{K}_{r\alpha}}{{\mu}_{\alpha}}(\nabla {P}_{\alpha}{\rho}_{\alpha}g\nabla D)\right)={\rho}_{\alpha}{q}_{\alpha}\text{.}$$(3)
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 twophase flow of water and oil with the assumption of incompressible and immiscible fluids is as follows (Siavashi et al., 2014): $$\left[\frac{\partial}{\partial x}\left({K}_{f}{\lambda}_{t}\frac{\partial {P}_{f}}{\partial x}\right)+\frac{\partial}{\partial y}\left({K}_{f}{\lambda}_{t}\frac{\partial {P}_{f}}{\partial y}\right)\right]={q}_{w}+{q}_{o}\text{,}$$(4) where λ_{t} is the total mobility ratio, which is the summation of water and oil mobility ratios as: $${\lambda}_{t}={\lambda}_{wf}+{\lambda}_{of}\text{.}$$(5)The water and oil mobility ratios (in the fracture) are expressed as: $${\lambda}_{wf}=\frac{{K}_{rwf}}{{\mu}_{w}},\text{\hspace{0.22em}}{\lambda}_{of}=\frac{{K}_{rof}}{{\mu}_{o}}\text{,}$$(6)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): $${K}_{rw}={\left(\frac{{S}_{w}{S}_{wc}}{1{S}_{wc}{S}_{or}}\right)}^{2}\text{,}$$(7) $${K}_{ro}={\left(\frac{{S}_{o}{S}_{or}}{1{S}_{wc}{S}_{or}}\right)}^{2}\text{,}$$(8) where S_{wc} is connate water saturation, S_{or} is residual oil saturation.
For incompressible and twophase flow (Ahmadpour et al., 2016): $$\begin{array}{c}\hfill {S}_{of}+{S}_{wf}=1\hfill \\ \hfill {S}_{om}+{S}_{wm}=1\hfill \end{array}\text{.}$$(9)
The equation of water saturation along a streamline for the matrix and fracture as a function of TimeOfFlight (TOF), as shown in Equations (10) and (11), are as follows (Kazemi et al., 1976; Di Donato et al., 2003): $$\frac{\partial {S}_{wf}}{\partial t}+\frac{\partial {f}_{wf}}{\partial \tau}=\frac{T}{{\phi}_{f}}\text{,}$$(10) $$\frac{\partial {S}_{wm}}{\partial t}=\frac{T}{{\phi}_{m}}\text{.}$$(11)
TOF (τ) (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 (Siavashi et al., 2016): $$\tau ={\displaystyle \underset{0}{\overset{s}{{\displaystyle \int}}}}\frac{\phi}{\left{\displaystyle \overrightarrow{u}}(\delta )\right}d\delta \text{,}$$(12) where φ is the porosity and δ and $\overrightarrow{u}$ 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 multidimensional 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 as follows: $${f}_{wf}=\frac{{\lambda}_{wf}}{{\lambda}_{wf}+{\lambda}_{of}}\text{.}$$(13)
The general steps of streamline simulation for a fractured reservoir using the DPSP model are as follows (DattaGupta 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 (Δt_{sl}) that is smaller than the total timestep (Δt_{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 operatorsplitting 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.
3 Matrixfracture 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 matrixfracture transfer functions. The effect of the matrixfracture 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 matrixfracture fluid transfer (Sarma and Aziz, 2004).
3.1 Conventional transfer function (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 multiphase 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).
The simplified form of this transfer function that considers 2D Darcy flow (Gilman and Kazemi, 1983; DattaGupta and King, 2007) can be written as: $$\begin{array}{c}\hfill {T}_{w}={F}_{s}{K}_{m}{\lambda}_{wf}({P}_{wf}{P}_{wm})\hfill \\ \hfill {T}_{o}={F}_{s}{K}_{m}{\lambda}_{of}({P}_{of}{P}_{om})\hfill \end{array}\text{,}$$(14) 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 dualporosity model network (Jerbi et al., 2017).
The shape factor is a geometric coefficient used to calculate the matrixfracture 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 matrixfracture transfer function for cubic blocks of a matrix (Gilman and Kazemi, 1983): $$\overrightarrow{{F}_{s}}}=4\left(\frac{1}{{L}_{x}^{2}}{\displaystyle \overrightarrow{i}}+\frac{1}{{L}_{y}^{2}}{\displaystyle \overrightarrow{j}}+\frac{1}{{L}_{z}^{2}}{\displaystyle \overrightarrow{k}}\right)\text{,$$(15) 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 twophase 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 nonwettingphase 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). $$\begin{array}{c}\hfill {P}_{cm}={P}_{om}{P}_{wm}\hfill \\ \hfill {P}_{cf}={P}_{of}{P}_{wf}\hfill \end{array}\text{.}$$(16)
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 wettingphase saturation and imbibition causes an increase. The capillary pressure of the fracture is negligible for a lowpermeable 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): $${P}_{cm}={P}_{cm}^{\mathrm{max}}\frac{{S}_{wm}^{1}{(1{S}_{orm})}^{1}}{{S}_{wcm}^{1}{(1{S}_{orm})}^{1}}\text{,}$$(17) where ${P}_{cm}^{\mathrm{max}}$ is the maximum capillary pressure that is depended on fluid and rock properties.
With the assumption of incompressible fluid: $${T}_{o}+{T}_{w}=0\text{.}$$(18)
By ignoring the fracture capillary pressure, Equation (16) can be used to convert the conventional transfer function to (AlHuthali and DattaGupta, 2004): $$T=F{K}_{m}\frac{{\lambda}_{wf}{\lambda}_{om}}{{\lambda}_{wf}+{\lambda}_{om}}{P}_{cm}\text{.}$$(19) Equations (10) and (11) can be discretized and rewritten in terms of local timestep (Δt_{sl}) as: $${({S}_{wf})}_{i}^{\mathit{int}}={({S}_{wf})}_{i}^{\mathit{int}}\frac{{T}_{i}^{n}}{{\varphi}_{f}}\mathrm{\Delta}{t}_{sl}\text{,}$$(20) $${({S}_{wm})}_{i}^{n+1}={({S}_{wm})}_{i}^{n}+\frac{{T}_{i}^{n}}{{\phi}_{m}}\mathrm{\Delta}{t}_{sl}\text{.}$$(21)In Equation (20), (S_{wf})^{int} denotes intermediate saturation and can be written as: $${({\mathit{S}}_{\mathit{w}\mathit{f}})}_{\mathit{i}}^{\mathit{int}}={\left({\mathit{S}}_{\mathit{w}\mathit{f}}\right)}_{\mathit{i}}^{\mathit{n}\frac{(\mathrm{\Delta}\mathit{t}\mathit{s}\mathit{l})}{\mathrm{\Delta}\mathit{T}}}({f}_{wf}({S}_{wf,\mathit{i}}^{n}){f}_{wf}({S}_{wf,\mathit{i}1}^{n})))$$(22)
3.2 Linear transfer function (Di Donato et al.)
Di Donato et al. (2003) proposed a linear transfer function for matrix blocks surrounded by water as: $$\begin{array}{cc}\hfill T=\beta {\phi}_{m}(1{S}_{orm}{S}_{wm})\hfill & \hfill {S}_{wf}>0\hfill \\ \hfill T=0\hfill & \hfill {S}_{wf}=0\hfill \end{array}\text{.}$$(23) This is a linear function of matrix saturation where β depends on the wettability and its value for a mixedwet system is 100–10 000 times smaller than for strongly waterwet systems. Di Donato et al. (2003) described the method for obtaining this transfer function. (Di Donato et al., 2003).
3.3 Linear transfer function (Kazemi et al.)
Many researchers have proposed other transfer functions which are consistent with the imbibition experiments and saturation curves of fractured reservoirs. de Swaan introduced a function similar to Equation (10) for S_{wf} = 1 and utilized a convolution integral to calculate matrix saturation for S_{wf} < 1 (de Swaan, 1978). Kazemi et al. extracted a linear transfer function from the convolution integral of de Swaan (Di Donato et al., 2003): $$T=\beta {\mathrm{\varnothing}}_{m}({S}_{wf}(1{S}_{orm}{S}_{wm})({S}_{wm}{S}_{wim}))\text{.}$$(24)
3.4 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 (capillarycontrolled imbibition). This newly defined transfer function is simple and works with a constant coefficient using the saturation difference between the matrix and fracture: $$\begin{array}{cc}\hfill T=0\hfill & \hfill {S}_{wf}<{S}_{wm}\hfill \\ \hfill T=\alpha ({S}_{wf}{S}_{wm})\hfill & \hfill {S}_{wf}\ge {S}_{wm}\hfill \end{array}\text{,}$$(25) Parameter α 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 (countercurrent imbibition). The reverse direction is not considered; therefore, the saturation difference between the matrix and fracture will change the transfer function.
4 Problem definition
Applying enhanced oil recovery methods, such as gas injection, wateralternating 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 5spot 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 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 SPE10 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: $$\beta {\phi}_{m}(1{S}_{orm}{{\displaystyle \overline{S}}}_{wm})=F{K}_{m}{P}_{cm}({{\displaystyle \overline{S}}}_{wm})\frac{{k}_{rom}({{\displaystyle \overline{S}}}_{wm})}{{\mu}_{o}}\text{,}$$(26) $$\alpha \mathrm{\Delta}{{\displaystyle \overline{S}}}_{w(fm)}=\beta {\phi}_{m}(1{S}_{orm}{{\displaystyle \overline{S}}}_{wm})\text{,}$$(27) where ${{\displaystyle \stackrel{\u203e}{S}}}_{wm}$ is average water saturation in matrix, ${{\displaystyle \overline{S}}}_{w(fm)}$ is average water saturation difference between matrix and fracture.
Table 2 shows the transfer functions parameters used in the simulation.
Fig. 4 5Spot pattern sector model. 
Fluid and reservoir properties.
Fig. 5 Permeability of fracture (heterogeneous fractured reservoir). 
Transfer functions parameters.
5 Results and discussion
The results of the simulation studies and comparison of different transfer functions are summarized in this section. To verify the results, IMEXCMG a fully implicit simulation software was applied. This simulation software uses the conventional transfer function (ComputerModelingGroup, 2007). The water saturation and pressure distribution in the fracture and matrix after 300 days of water injection are shown in the sectormodel grids. The oil flow rate and cumulative volume of oil produced by the production wells using different transfer functions are represented. The TOF and streamlines are also shown.
Figures 6–10 represent water saturation in the fracture.
Figures 11–15 represent water saturation in the matrix.
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 β.
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 sectormodel 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.
Fig. 6 Water saturation in fracture by commercial software. 
Fig. 7 Water saturation in fracture by conventional transfer function. 
Fig. 8 Water saturation in fracture by Di Donato et al. linear transfer function. 
Fig. 9 Water saturation in fracture by Kazemi et al. linear transfer function. 
Fig. 10 Water saturation in fracture by new linear transfer function with constant coefficient. 
Fig. 11 Water saturation in matrix by commercial software. 
Fig. 12 Water saturation in matrix by conventional transfer function. 
Fig. 13 Water saturation in matrix by Di Donato et al. linear transfer function. 
Fig. 14 Water saturation in matrix by Kazemi et al. linear transfer function. 
Fig. 15 Water saturation in matrix by new linear transfer function with constant coefficient. 
Fig. 16 Pressure distributions in fracture by commercial software. 
Fig. 17 Pressure distributions in fracture by conventional transfer function. 
Fig. 18 Pressure distributions in fracture by Di Donato et al. linear transfer function. 
Fig. 19 Pressure distributions in fracture by Kazemi et al. linear transfer function. 
Fig. 20 Pressure distributions in fracture by new linear transfer function with constant coefficient. 
Fig. 21 Oil flow rate at production well A. 
Fig. 22 Oil flow rate at production well B. 
Fig. 23 Oil flow rate at production well C. 
Fig. 24 Oil flow rate at production well D. 
Fig. 25 Cumulative volume of produced oil at production well A. 
Fig. 26 Cumulative volume of produced oil at production well B. 
Fig. 27 Cumulative volume of produced oil at production well C. 
Fig. 28 Cumulative volume of produced oil at production well D. 
Fig. 29 Timeofflight contour after the first total timestep (100 days). 
Fig. 30 Streamlines after the first total timestep (100 days). 
6 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 singlelayer heterogeneous fractured reservoirs, waterflooding was simulated using commercial software and an inhouse streamlinebased numerical code that uses the DPSP reservoir model. The Kazemi et al. conventional transfer function was used in the commercial software. In the streamline simulator, in addition to the conventional transfer function, three linear transfer functions were reviewed and employed. The commercial software verified the results of the streamline simulator, especially for the Kazemi et al. conventional transfer function. By adjusting the transfer function parameters, the results became more similar and all transfer functions were fairly acceptable.
The transfer functions discussed in this research are the Kazemi et al. conventional transfer function, the Di Donato et al. linear transfer function, the Kazemi et al. linear transfer function and the proposed linear transfer function with constant coefficients. The proposed linear transfer function is simple to implement, has less calculation complexity, is appropriate for the streamline method and offers acceptable accuracy when compared with other transfer functions.
List of symbols
S_{α}: α phase saturation, dimensionless
f_{α}: fractional flow of α phase
L_{α}: fracture distance in α direction
Greek letters
ρ_{α}: density of α phase, kg/m^{3}
α, β: transfer function constant coefficients, 1/s
References
 Ahmadpour M., Siavashi M, Doranehgard M.H. (2016) Numerical simulation of twophase flow in fractured porous media using streamline simulation and IMPES methods and comparing results with a commercial software, J. Cent. South Univ. 23, 10, 2630–2637. [CrossRef] [Google Scholar]
 AlHuthali A., DattaGupta A. (2004) Streamline simulation of countercurrent imbibition in naturally fractured reservoirs, J. Pet. Sci. Eng. 43, 3–4, 271–300. [CrossRef] [Google Scholar]
 Baker R. (2001) Streamline technology: reservoir history matching and forecasting its success, limitations, and future, J. Can. Pet. Technol. 40, 4, 23–27. [CrossRef] [Google Scholar]
 Barenblatt G.I., Zheltov I.P., Kochina I.N. (1960) Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata], J. Appl. Math. Mech. 24, 5, 1286–1303. [CrossRef] [Google Scholar]
 Batycky R.P, Blunt M.J., Thiele M.R. (1997) A 3D fieldscale streamlinebased reservoir simulator, SPE Reserv. Eng. 12, 04, 246–254. [CrossRef] [Google Scholar]
 Biagi J., Agarwal R., Zhang Z. (2016) Simulation and optimization of enhanced gas recovery utilizing CO_{2}, Energy 94, 7, 78–86. [CrossRef] [Google Scholar]
 Bourbiaux B. (2010) Fractured reservoir simulation: a challenging and rewarding issue, Oil Gas Sci. Technol. − Rev. IFP 65, 2, 227–238. [Google Scholar]
 Bourbiaux B., Cacas M.C., Sarda S., Sabathier J.C. (1998) A rapid and efficient methodology to convert fractured reservoir images into a dualporosity model, Rev. IFP 53, 6, 785–799. [CrossRef] [Google Scholar]
 Brooks R., Corey T. (1964) HYDRAU uc properties of porous media, Hydrology Papers, Colorado State University 24. [Google Scholar]
 Chen Z., Huan G., Ma Y. (2006) Computational methods for multiphase flows in porous media, SIAM, Philadelphia, USA. [Google Scholar]
 Choobineh M.J., Siavashi M., Nakhaee A. (2015) Optimization of oil production in water injection process using ABC and SQP algorithms employing streamline simulation technique, Modares Mech. Eng. 15, 5, 227–238. [Google Scholar]
 Christie M.A., Blunt M.J. (2001) Tenth SPE comparative solution project: a comparison of upscaling techniques, in: SPE Reservoir Simulation Symposium, Society of Petroleum Engineers. [Google Scholar]
 ComputerModelingGroup (2007) Launcher user's guide, computer modeling group ltd, calgary, advanced oil/gas reservoir simulator. [Google Scholar]
 DattaGupta A., King M.J. (2007). Streamline simulation: theory and practice, Society of Petroleum Engineers Richardson, Richardson, TX, USA. [Google Scholar]
 de Swaan A. (1978) Theory of waterflooding in fractured reservoirs, Soc. Pet. Eng. J. 18, 02, 117–122. [CrossRef] [Google Scholar]
 Di Donato G., Blunt M.J. (2004) Streamlinebased dualporosity simulation of reactive transport and flow in fractured reservoirs, Water Resour. Res. 40, 4, W04203. [CrossRef] [Google Scholar]
 Di Donato G., Huang W., Blunt M. (2003). Streamlinebased dual porosity simulation of fractured reservoirs, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
 Doranehgard M.H., Siavashi M. (2018) The effect of temperature dependent relative permeability on heavy oil recovery during hot water injection process using streamlinebased simulation, Appl. Therm. Eng. 129, 11, 106–116. [CrossRef] [Google Scholar]
 Gilman J.R. (2003) Practical aspects of simulation of fractured reservoirs, in: International Forum on Reservoir Simulation, Buhl, BadenBaden, Germany. [Google Scholar]
 Gilman J.R., Kazemi H. (1983) Improvements in simulation of naturally fractured reservoirs, SPE J. 23, 4, 695–707. [Google Scholar]
 Hassane T.F.R. (2013) The application of streamline reservoir simulation calculations to the management of oilfield scale, HeriotWatt University, Edinburgh, Scotland, UK. [Google Scholar]
 Jerbi C., Fourno A., Noetinger B., Delay F. (2017) A new estimation of equivalent matrix block sizes in fractured media with twophase flow applications in dual porosity models, J. Hydrol. 548, 40, 508–523. [CrossRef] [Google Scholar]
 Kazemi H., Merrill L.S. Jr., Porterfield K.L., Zeman P. (1976) Numerical simulation of wateroil flow in naturally fractured reservoirs, Soc. Pet. Eng. J. 16, 06, 317–326. [CrossRef] [Google Scholar]
 Landereau P., Noetinger B., Quintard M. (2001) Quasisteady twoequation models for diffusive transport in fractured porous media: largescale properties for densely fractured systems, Adv. Water Resour. 24, 8, 863–876. [CrossRef] [Google Scholar]
 LeBlanc J.L., Caudle B.H. (1971) A streamline model for secondary recovery, Soc. Pet. Eng. J. 11, 01, 7–12. [CrossRef] [Google Scholar]
 Lemonnier P., Bourbiaux B. (2010) Simulation of naturally fractured reservoirs. State of the art, Part 1Physical mechanisms and simulator formulation, Oil Gas Sci. Technol. − Rev. IFP 65, 2, 239–262. [Google Scholar]
 Noetinger B., Roubinet D., Russian A., Le Borgne T., Delay F., Dentz M., de Dreuzy J.R., Gouze P. (2016) Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale, Transp. Porous Media 115, 2, 345–385. [Google Scholar]
 Pollock D.W. (1988) Semianalytical computation of path lines for finitedifference models, Groundwater 26, 6, 743–750. [CrossRef] [Google Scholar]
 Sarma P., Aziz K. (2004). New transfer functions for simulation of naturally fractured reservoirs with dual porosity models, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
 Siavashi M., Blunt M.J., Raisee M., Pourafshary P. (2014) Threedimensional streamlinebased simulation of nonisothermal twophase flow in heterogeneous porous media, Comput. Fluids 103, 10, 116–131. [CrossRef] [Google Scholar]
 Siavashi M., Tehrani M.R., Nakhaee A. (2016) Efficient particle swarm optimization of well placement to enhance oil recovery using a novel streamlinebased objective function, J. Energy Resour. Technol. 138, 5, 052903. [CrossRef] [Google Scholar]
 Thiele M.R. (2001) Streamline simulation, in: Proc. Sixth International Forum on Reservoir Simulation, citeseer. [Google Scholar]
 Vitoonkijvanich S., AlSofi A.M., Blunt M.J. (2015) Design of foamassisted carbon dioxide storage in a North Sea aquifer using streamlinebased simulation, Int. J. Greenh. Gas Control 33, 12, 113–121. [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Schematics of different fracture − models for simulating fractured reservoirs (DattaGupta and King, 2007). 

In the text 
Fig. 2 Schematic representation of total and local timesteps (Doranehgard and Siavashi, 2018). 

In the text 
Fig. 3 Simulation run time ratio for various number of grids for streamline and conventional (IMPES) method. 

In the text 
Fig. 4 5Spot pattern sector model. 

In the text 
Fig. 5 Permeability of fracture (heterogeneous fractured reservoir). 

In the text 
Fig. 6 Water saturation in fracture by commercial software. 

In the text 
Fig. 7 Water saturation in fracture by conventional transfer function. 

In the text 
Fig. 8 Water saturation in fracture by Di Donato et al. linear transfer function. 

In the text 
Fig. 9 Water saturation in fracture by Kazemi et al. linear transfer function. 

In the text 
Fig. 10 Water saturation in fracture by new linear transfer function with constant coefficient. 

In the text 
Fig. 11 Water saturation in matrix by commercial software. 

In the text 
Fig. 12 Water saturation in matrix by conventional transfer function. 

In the text 
Fig. 13 Water saturation in matrix by Di Donato et al. linear transfer function. 

In the text 
Fig. 14 Water saturation in matrix by Kazemi et al. linear transfer function. 

In the text 
Fig. 15 Water saturation in matrix by new linear transfer function with constant coefficient. 

In the text 
Fig. 16 Pressure distributions in fracture by commercial software. 

In the text 
Fig. 17 Pressure distributions in fracture by conventional transfer function. 

In the text 
Fig. 18 Pressure distributions in fracture by Di Donato et al. linear transfer function. 

In the text 
Fig. 19 Pressure distributions in fracture by Kazemi et al. linear transfer function. 

In the text 
Fig. 20 Pressure distributions in fracture by new linear transfer function with constant coefficient. 

In the text 
Fig. 21 Oil flow rate at production well A. 

In the text 
Fig. 22 Oil flow rate at production well B. 

In the text 
Fig. 23 Oil flow rate at production well C. 

In the text 
Fig. 24 Oil flow rate at production well D. 

In the text 
Fig. 25 Cumulative volume of produced oil at production well A. 

In the text 
Fig. 26 Cumulative volume of produced oil at production well B. 

In the text 
Fig. 27 Cumulative volume of produced oil at production well C. 

In the text 
Fig. 28 Cumulative volume of produced oil at production well D. 

In the text 
Fig. 29 Timeofflight contour after the first total timestep (100 days). 

In the text 
Fig. 30 Streamlines after the first total timestep (100 days). 

In the text 