Numerical Simulation of Two-Phase Flows in Heterogeneous Porous Media

In this work, we present three-dimensional numerical simulations of water-oil flow in porous media in order to analyze the influence of the heterogeneities in the porosity and permeability fields and, mainly, their relationships upon the phenomenon known in the literature as viscous fingering. For this, typical scenarios of heterogeneous reservoirs submitted to water injection (secondary recovery method) are considered. The results show that the porosity heterogeneities have a markable influence in the flow behavior when the permeability is closely related with porosity, for example, by the Kozeny-Carman (KC) relation. This kind of positive relation leads to a larger oil recovery, as the areas of high permeability (higher flow velocities) are associated with areas of high porosity (higher volume of pores), causing a delay in the breakthrough time. On the other hand, when both fields (porosity and permeability) are heterogeneous but independent of each other the influence of the porosity heterogeneities is smaller and may be negligible.


INTRODUCTION
Every oil field has a life cycle that goes from discovery to abandonment. Therefore, predicting when to leave the field or when to apply a secondary or tertiary recovery method to increase the volume of the recovered oil is crucial for the oil industry. This prediction is usually made using reservoir simulations, which consists of the elaboration of numerical models that take into account the physical phenomenology typical of flows in porous media, in the search for approximate solutions as close as possible to reality [38]. From the mathematical model of the problemdescribed by a system of partial differential equations of the hyperbolic-parabolic type -computational models are used to analyze the behavior of the reservoir by means of predictions, which use experimental data to determine the parameters employed in the model [19].
Oil recovery methods are classified into primary, secondary and tertiary. The primary consists of the volume of oil produced through the natural energy given by the decompression of the fluids [43]. When reservoirs retain a large amount of oil, after primary recovery the secondary or tertiary methods are used to obtain additional production. Secondary recovery is the amount of hydrocarbons obtained by supplementing the primary energy with the injection of fluids into strategically selected wells. On the other hand, the tertiary recovery uses physical, chemical or biological products (polymers, steam, solvents, gases, bacteria, etc.) to increase the recovery factor to around 60%. It is used in reservoirs where the secondary recovery method fails or would fails if employed; however, are expensive methods [41,43]. One of the most employed secondary method is the injection of water, which generally produces good results in oil reservoirs, increasing the recovery factor to about 30%. However, when we inject water (less viscous) into an oil-filled (more viscous) porous medium, an unstable interface (morphologically) occurs between the two fluids, creating a finger profile. In addition, the heterogeneities present in the rocks, mainly in the porosity and permeability, affect the appearance and behavior of this finger growth [5]. This phenomenon is known in the literature as viscous fingering [27]. When it happens, some of the injected water advances at velocities higher than the mid-wave front, not by sweeping the oil and arriving early to the production wells [25].
Many authors have been studied the impact of heterogeneities in the permeability or porosity fields on the porous media flows ( [3,4,6,8,12,13,16,17,20,21,24,26,28,29,30,35,36,39,42,44]). Typically, the permeability in natural reservoirs may vary in space by 3 to 4 orders of magnitude whereas porosity vary by one order [20]. In [35] the permeability values, in square centimeters, varies from 10 −16 to 10 −3 , a range that covers 13 orders of magnitude, while the variation in porosity values is generally no more than 0.2. As a consequence, in studies about the influence of the heterogeneities in porous media hydraulic properties on the flow, it is usually assumed heterogeneity acting only in the permeability field and, the variations in the porosity field are often neglected ( [3,6,8,12,16,17,20,28,29,35,39,44]). On the other hand, some authors have considered both porosity and permeability as heterogeneous fields, but independent of each other (unrelated). For example, Amaziane et. al. [4] have studied numerical methods for two-phase flow in heterogeneous porous media assuming that porosity and permeability are periodic oscillatory functions in space and, Gundersen et. al. [21] have solved a 2D problem considering a log-normal permeability field and a Gaussian porosity field. However, Sperl et al. [42] have experimentally observed, that the permeability is determined by the micro-structure of the material, that is, it depends on the porosity, mainly, on the arrangement and size of the pores. Therefore, it is important to study the relationships between the porosity and permeability fields on the water-oil flow in heterogeneous reservoirs.
Hassan et al. [22,23] have concluded that direct relations between porosity and permeability fields have pronounced influence on the macrodispersion of solutes in heterogeneous porous media. More recent, Correa et al. [13] have analized the relationship between the porosity and permeability fields on the two-phase flow (water-oil) in a 2D toy problem using the Kozeny-Carman relation. It was observed a delay in the breakthrough time when the heterogeneous porosity and Tend. Mat. Apl. Comput., 21, N. 2 (2020) CARNEIRO, BORGES and MALTA 341 permeability fields are related in comparison with the case where the same permeability fields were used but the porosity field considered a constant (equal to the mean of heterogeneous one). Finally, we comment that Henderson et al. [24] have used a generalized three-parameter Kozeny-Carman model (TPKCG) to obtain the permeability field from the Gaussian distributed fractal porosity field to study the viscous fingering phenomena during a water injection process in oil reservoir simulations. They have found that this generalization of the Kozeny-Carman equation can be used in numerical simulations of oil recovery processes susceptible to hydrodynamic instability phenomena.
Motivated by the previous discussion, the objective of this work is to study the influence of the heterogeneities present in the porosity and permeability fields and their relations on the behavior of the water front, taking into account the appearance and evolution of the viscous fingering, in reservoirs submitted to the secondary recovery process. For this end, we performed and analized numerical simulations in tri-dimensional scenarios based on the SPE benchmark model 2 (SPE2) [11].
The outline of this work is as follows. Section 2 shows the mathematical model of the problem, as well as the KT method to approximate the saturation equation. Numerical experiments are reported in Section 3. Finally, conclusions are given in Section 4.

THE MODEL PROBLEM
The mathematical model that describes an immiscible two-phase flow in a rigid saturated porous medium, with incompressible fluids, and in the absence of capillary forces, consists of a system of equations formed by the water saturation transport equation (2.1), Darcy's Law (2.2) and the condition of incompressibility (2.3), that is, and where S w is the water saturation, φ is the porosity, k is the permeability, f w is the flux function of water, v t = v w + v o is the total velocity, ρ o is the density of crude oil, ρ w is the density of water, g is the acceleration of gravity, z is the depth, λ o = k ro /µ o , and λ w = k rw /µ w are the mobilities of oil and water, respectively, with λ t = λ o + λ w . Here k rw and k ro are the relative permeabilities of water and oil and µ w , µ o the respective viscosities.
For the problem defined by equations (2.1), (2.2) and (2.3) to be complete, it is necessary to establish initial and boundary conditions. To do this, first we exhibit in Figure 1 the 1/4 fivespot geometry, that will be used in the numerical examples. Considering the symmetry of the problem, our study is restricted to the hatched area where we have one injector well (in the lower left corner) and one producer well (in the upper right corner). The domain to be studied (1/4 model) consists of a parallelepiped Ω ⊂ R 3 with injector (Γ in , black ) and producer (Γ out , gray) wells located in opposite corners, as shown in Figure 2. The null-flow condition is imposed across the ∂ Ω (boundary of Ω), except in Γ in ∪ Γ out where we consider v t · n = q in in Γ in and pressure P = 0 in Γ out . For the initial time condition, t = 0, we assume the reservoir is completely filled with oil [13].
In order to increase computational efficiency and allow that different numerical methods may be used, we decouple the system (2.1)-(2.3), and the velocity-pressure and saturation problems are solved separately. Therefore, we take macro time steps, in which the velocity-pressure problem is solved using a mixed method with lower-order Raviart-Thomas elements [2,18]. Between two macro-steps, with the velocity frozen, transport micro-steps are performed, respecting the CFL condition, to update the saturation [37]. In this step we apply an extension of the Kurganov and Tadmor (KT) method [32] for heterogeneous porosity fields [13], that is briefly described in the next section.

The Kurganov and Tadmor (KT) method
For the sake of simplicity, we present the KT method for a unidimensional transport problem, considering the porosity φ = 1. To this end, let the conservation law equation given by subject to the following initial and boundary conditions It is important to note the relation between equations (2.1) and (2.4). To define the KT method, we integrate the hyperbolic differential equation (2.4) at each time step over the intervals , as illustrated in Figure 3. Next, we use the REA (Reconstruct-Evolve-Average) algorithm, in the form Figure 3: KT scheme adapted from [32].

Reconstruction (R):
The values of S n j are reconstructed by piecewise linear polynomial functions, that is We observe that to ensure mass conservation, linear reconstructions must satisfy In addition, to avoid the appearance of spurious oscillations, due to proposed linear reconstruction by parts, the slopes must be chosen conveniently, that is, Thus, the KT method has the TVD (Total Variation Diminishing) property, described below [33].
First of all, we consider the S function defined in a mesh with M cells, a Total Variation TV is set as One method is called TVD if, for any S n , its values for S n+1 satisfy the following inequality TV(S n+1 ) TV(S n ).
(2.10) Therefore, the KT method uses a slope limiter to determine (S x ) n j in order to satisfy the TVD relation (2.10). To do this, a MC limiter (monotonized central-difference limiter) is employed. It is given by a family of discrete derivatives, parameterized by θ , with 1 ≤ θ ≤ 2 , that is, where the MinMod is defined by Therefore, the semi-discrete KT formulation is given by (2.14) In (2.14), the coefficients a n j±1/2 are the local propagation velocities a n j±1/2 := max u∈C(s − with ρ the spectral ratio. If f is a convex function we have a n j±1/2 : where s + j−1/2 and s − j−1/2 are the values of s(x,t) in the point (x j−1/2 ,t n ) on the linear piecewise reconstructions S n j (x) and S n j−1 (x), respectively, given by and Similary, we have The ordinary differential equation (2.13) defines the KT method for a hyperbolic differential equation of the type (2.4). This equation will be solved using the explicit Euler method, that is, denoting the right hand side of equation ( where s min and s max are the maximum and minimum values of a scalar quantity for a given t n .

The Kozeny-Carman (KC) relation
As we have already mentioned, the main objective of this work is to study the effects of heterogeneities in the permeability and porosity fields as well as possible relationships between them on the viscous fingering process. Several laboratory studies have shown that the permeability mainly depends on the pore size and shape of the pore [7,10,15,40,42]. One of the widely used models associating the porosity (φ ) and the permeability k(φ ) is the Kozeny-Carman correlation [9,31] given by where c is the Kozeny constant representing the particle shape factor and S is the specific surface. In the following experiments, equation (2.22) will be take into account as an example, for the construction of related permeability and porosity fields. It should be noted that any other type of model that relates these two fields could be used. We are not presenting any defence of the KC model.

NUMERICAL RESULTS
All numerical simulations are based on the SPE benchmark 2 (SPE2) [11] adapted to simulate the secondary recovery process (water-oil) in a three-dimensional domain with 304.80 × 15.24 × 304.80 m 3 . Figures 4 and 5 illustrate the permeability (k) and porosity (φ ) fields of the original SPE2 geological model. In particular, we consider the Upper Ness formation cut from the Brent Group field, since this region has well-defined permeable zones (channels) that accentuate the viscous fingering process. Specifically, in this work we choose sections of the original geological models that are displayed in Figures 6 and 7 Figure 5).
In the following simulations we analyze five scenarios: • Case 0: k and φ homogeneous; • Case 1: k heterogeneous and φ homogeneous; • Case 2: k homogeneous and φ heterogeneous; • Case 3: k and φ heterogeneous, and; • Case 4: the inverse of the relation KC (2.22) with c = 5 and S = 4.47 × 10 4 are used to generate the φ from k ( Figure 6).
The parameters c and S are chosen such that the average of the porosity field is equal to the original one. The fluid properties were adapted from SPE2 benchmark and are presented in Table 1. The parameters associated with the five-spot problem described at section 2 are given in Table 2. In addition, we present the relative permeability functions where S rw and S ro are the residual saturations of water and oil, respectively. Figure 8 exhibits the water saturation results at t = 180, 360, 540, 720 and 900 days for Case 0 (k and φ homogeneous) which are used only for comparison with heterogeneous cases. These results are consistent with those found in [41].
Figures 9 and 10 display the water front behavior at t = 180, 360, 540, 720 and 900 days for Cases 1 and 2, respectively. We notice that the viscous fingering process is more pronounced in Case 1, where the permeability is heterogeneous, which provides an early breakthrough when compared to Case 2 (homogeneous permeability and heterogeneous porosity) as it can be seen in Figure  11 where it is shown the oil production curves for both cases. These results explain, in part, the reason why many studies on stochastic porous media flows consider heterogeneities only in the permeability fields (see [8,14,16,17,20]). Porosity (φ ) 0.125 Table 2: Five-spot parameters.

Parameter Value
Water flow in injector well q in 198.7341 (m 3 /d) Producer well pressure (P) 0.0(Pa) Saturation in injection (S) 1.0 The water saturation at t = 180, 360, 540, 720 and 900 days for Case 3 (k and φ heterogeneous, without any relation) is exbihited in Figure 12. These results are very similar to those obtained for Case 1 (Figure 9), where the porosity field is homogeneous. The porosity heterogeneities have a small impact on the viscous fingering process when the porosity and permeability fields are not related. This fact is confirmed by the similarity of the production curves from Cases 1 and 3 ( Figure 13). Once again, these behaviors support the hypothesis assumed in the most of stochastic porous media flows studies, that only the permeability field is usually considered heterogeneous.
In Figure 14 the water saturation at t = 180, 360, 540, 720 and 900 days is showed to Case 4 (KC relation), where the porosity and permeability fields strongly related. Unlike the previous cases, we can observe a strong influence of the porosity heterogeneities on the viscous fingering process.
(e) t = 900d To better compare Cases 1, 2, 3 and 4 the water saturation plots at 720 days are displayed in Figure 15. It is possible to observe that the variation (variance) in the porosity fields is typically smaller than that in the permeability fields. Therefore, when there is no relationship between the permeability and porosity fields, the porosity heterogeneities yield little influence on the viscous fingering process, if compared to the effects produced by the heterogeneities in the permeability field. These observations justify the hypothesis that the porosity field is homogeneous, taken into account by many works in the literature. On the other hand, if there is a strong relationship between permeability and porosity fields (Figure 15-c) the porosity heterogeneities play an important role in the formation and development of the viscous fingering process. In this case, high permeability zones with higher local velocities are associated to high porosity zones, i.e., zones with larger capacitance. This scenario promotes a desirable delay in the water front and, consequently, in the breakthrough time. After the breakthrough time, part of the injected water displaces the oil while the other part starts to be produced, generating a decline in the production of the reservoir. [41]. To illustrate this issue, Figure 16 displays the oil production as a function of time (days) for all simulated scenarios. The constant part of the graph shows the period before (e) t = 900d Figure 9: Water saturation at different times for Case 1.
the breakthrough time, that is, a period in which all the injected water remains in the porous medium producing the same amount of oil. Although the same permeability field is considered in Cases 1, 3 and 4, when the KC relation (Case 4) is applied, it is observed a delay of about 395 days in the breakthrough time if compared with homogeneous and uncorrelated porosity fields (Cases 1 and 3, respectively).

CONCLUSION
The Numerical simulations of the secondary oil recovery process from a highly heterogeneous reservoir presented in this work showed that not only the existence of heterogeneities in the permeability and porosity fields are important in the flow pattener, but also the possible relations between these fields can play a relevant role in the development of the viscous fingering phenomenon with pronounced impact on the oil production. For example, when the well known Kozeny-Carman relation was applied in the experiments, a delay in the advance of the water front, causing a larger oil recovery, was observed compared to the case in which the porosity (e) t = 900d    field is considered constant or it is not related to permeability. Finally, we comment that these observations suggest that the simplifying hypothesis, widely used in stochastic studies, that only the permeability field is heterogeneous should be applied only in specific cases.

ACKNOWLEDGMENT
The first author thanks the National Council for Scientific and Technological Development (CNPq) for the scholarship that supported the execution of this work.