On the Modelling of Immiscible Viscous Fingering in Two-Phase Flow in Porous Media

Viscous fingering in porous media is an instability which occurs when a low-viscosity injected fluid displaces a much more viscous resident fluid, under miscible or immiscible conditions. Immiscible viscous fingering is more complex and has been found to be difficult to simulate numerically and is the main focus of this paper. Many researchers have identified the source of the problem of simulating realistic immiscible fingering as being in the numerics of the process, and a large number of studies have appeared applying high-order numerical schemes to the problem with some limited success. We believe that this view is incorrect and that the solution to the problem of modelling immiscible viscous fingering lies in the physics and related mathematical formulation of the problem. At the heart of our approach is what we describe as the resolution of the “M-paradox”, where M is the mobility ratio, as explained below. In this paper, we present a new 4-stage approach to the modelling of realistic two-phase immiscible viscous fingering by (1) formulating the problem based on the experimentally observed fractional flows in the fingers, which we denote as fw∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ f_{\rm w}^{*} $$\end{document}, and which is the chosen simulation input; (2) from the infinite choice of relative permeability (RP) functions, krw∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ k_{\rm rw}^{*} $$\end{document} and kro∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ k_{\rm ro}^{*} $$\end{document}, which yield the same fw∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ f_{\rm w}^{*} $$\end{document}, we choose the set which maximises the total mobility function, λT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \lambda_{\text{T}}^{{}} $$\end{document} (where λT=λo+λw\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \lambda_{\text{T}}^{{}} = \lambda_{\text{o}}^{{}} + \lambda_{\text{w}}^{{}} $$\end{document}), i.e. minimises the pressure drop across the fingering system; (3) the permeability structure of the heterogeneous domain (the porous medium) is then chosen based on a random correlated field (RCF) in this case; and finally, (4) using a sufficiently fine numerical grid, but with simple transport numerics. Using our approach, realistic immiscible fingering can be simulated using elementary numerical methods (e.g. single-point upstreaming) for the solution of the two-phase fluid transport equations. The method is illustrated by simulating the type of immiscible viscous fingering observed in many experiments in 2D slabs of rock where water displaces very viscous oil where the oil/water viscosity ratio is (μo/μw)=1600\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ (\mu_{\text{o}} /\mu_{\text{w}} ) = 1600 $$\end{document}. Simulations are presented for two example cases, for different levels of water saturation in the main viscous finger (i.e. for 2 different underlying fw∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ f_{\rm w}^{*} $$\end{document} functions) produce very realistic fingering patterns which are qualitatively similar to observations in several respects, as discussed. Additional simulations of tertiary polymer flooding are also presented for which good experimental data are available for displacements in 2D rock slabs (Skauge et al., in: Presented at SPE Improved Oil Recovery Symposium, 14–18 April, Tulsa, Oklahoma, USA, SPE-154292-MS, 2012. 10.2118/154292-MS, EAGE 17th European Symposium on Improved Oil Recovery, St. Petersburg, Russia, 2013; Vik et al., in: Presented at SPE Europec featured at 80th EAGE Conference and Exhibition, Copenhagen, Denmark, SPE-190866-MS, 2018. 10.2118/190866-MS). The finger patterns for the polymer displacements and the magnitude and timing of the oil displacement response show excellent qualitative agreement with experiment, and indeed, they fully explain the observations in terms of an enhanced viscous crossflow mechanism (Sorbie and Skauge, in: Proceedings of the EAGE 20th Symposium on IOR, Pau, France, 2019). As a sensitivity, we also present some example results where the adjusted fractional flow (fw∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ f_{\rm w}^{*} $$\end{document}) can give a chosen frontal shock saturation, Swf∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ S_{\rm wf}^{*} $$\end{document}, but at different frontal mobility ratios, M(Swf∗)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ M(S_{\rm wf}^{*} ) $$\end{document}. Finally, two tests on the robustness of the method are presented on the effect of both rescaling the permeability field and on grid coarsening. It is demonstrated that our approach is very robust to both permeability field rescaling, i.e. where the (kmax/kmin) ratio in the RCF goes from 100 to 3, and also under numerical grid coarsening.


Introduction
Viscous instability in fluid displacements in porous media is a classical physical problem with a long and rich history. Here, we focus on viscous instability where the source of the effect is that the displacing phase viscosity (μ 1 ) is a much lower than that of the displaced phase (μ 2 ), i.e. (μ 2 /μ 1 ) >> 1. Both miscible and immiscible viscous unstable displacements have been extensively described in the literature in both Hele-Shaw cells and also in porous media (see review by Homsy 1987).
The continuum transport equations describing both miscible and immiscible displacement in porous media are quite well established. The general derivation of these equations can be carried out by applying material balance and then, in a porous media, invoking a form of the Darcy Law, either in its single-phase or two-phase forms to obtain (in 2D or 3D) a pressure equation. Solution of this pressure equation yields the velocity field, which may then be used in a transport equation to describe either the single-phase or two-phase transport. In the former miscible case, the two fluids are viewed to be miscible in all proportions, and the viscosity of any mixture of the 2 fluids (the less and the more viscous) is described by a mixing rule for the viscosity, μ mix = μ mix (f, μ 1 , μ 2 , a, b, …) where f is the fraction of one of the fluids, a, b,… are some model parameters and the viscosity tends correctly to the limits (μ 1 or μ 2 ) where f = 0 or 1. Experimental and numerical modelling results of miscible viscous fingering in 2D (rectangular and five-spot) models show excellent agreement (Sorbie et al. 1995;Zhang et al. 1997).
The generalised two-phase equations are more complex in that they contain both viscous and capillary pressure terms, and these are described in several books and papers (e.g. Peaceman 1977;Aziz and Settari 1979;Pinder and Gray 2008). In the context of immiscible viscous fingering, these equations have been laid out very clearly in the recent study by Bakharev et al. (2020). We take the immiscible displacing phase as water ( 1 = w ) and the displaced phase as a viscous oil ( 2 = o ). Henceforth, we will consider only water and oil as the immiscible displacing and displaced phases, respectively. We assume an incompressible fluid-porous medium system, i.e. fluid and rock densities are independent of pressure. Then, in the absence of gravity, the coupled two-phase equations for pressure (say P o ) and saturation (say S o ) in 3D are as follows: where the quantities, o and w are the oil and water phase mobilities, o = (k ro ∕ o ) and w = (k rw ∕ w ) , k ro and k rw are the relative permeability (RP) functions, o and w are the fluid viscosities, where the subscripts o and w refer to oil and water, respectively. The absolute permeability is represented by the tensor, k , although in this work we will take this to be diagonal and isotropic; that is, in 2D: where k xx = k yy . A full two-phase tensor formulation is presented elsewhere (Pickup and Sorbie 1996), but all of the findings of this paper can be carried over into such a generalised formulation. The capillary pressure, P c , is introduced to close the equations by relating the pressures in the phases (which are different) through a nonlinear capillary pressure term, P c (S w ) = P o − P w , which is a function of one of the phase saturations, say S w , but by definition, S w + S o = 1 , so either phase saturation can be used. In the viscous dominated limit, then P c (S w ) = 0 , which assumes that there is a single pressure, P = P o = P w , and the above equations reduce accordingly.
The vast majority of two-phase immiscible simulations in porous media, in applications such as aquifer remediation, packed bed flow, in CO 2 storage or in oil reservoirs, are carried out including only viscous forces and usually gravity. Capillarity is neglected-that is it is assumed that, P c (S w ) = 0-as it is viewed to be a "small scale" (mm-cm) or even a "pore scale" effect. However, the two-phase transport, described by Eq. (2), strictly must have a lower capillary length scale, L cap , which is related to capillary diffusion. This is discussed further below, but the issue should be kept in mind for the development of our approach to immiscible viscous fingering below.

Immiscible Viscous Fingering: Experimental Review
There is an extensive literature on viscous fingering from the initial work of Engelberts and Klinkenberg (1951) and Hill (1952) and the early classical work of Taylor and Saffman (1959). Indeed, we find that the earliest explicit mention of the term "viscous fingering" is in the paper by Engleberts and Klinkenberg (1951), and this is for immiscible (oil/water) fingering. Focusing only on immiscible fingering, this earlier work up to the late 1980s has been reviewed by Homsy (1987) and Kueper and Frind (1988). However, there is vastly more work on the theory (stability analysis, numerics, etc.) than on experiments (e.g. Chuoke et al. 1959; Alemán and Slattery 1988;Chikhliwala and Yortsos 1988). One of the earliest studies showing immiscible fingering was in the work of van Meurs and van der Poel (1958) reproduced in Fig. 1. In this work, refractive index matched oil and beads were used to construct high permeability bead packs and thus when water was injected into the oil (viscosity ratio, ( o ∕ w ) = 80 ) the resulting fingers could be seen and photographed directly. The viscous fingers in their experiment are quite beautiful, and they also give a very early but unheeded clue on how to model them; in their experiments, the authors could measure oil recovery but they could not directly measure the water saturation, S w , in the fingers. However, it is evident visually that S w in the main fingering region is quite high, certainly in the range S w ≈ 0.3-0.8, and this is important, as we explain below. More recent work on immiscible viscous fingering in 2D rock slabs has been reported by researchers in the group of Skauge and co-workers (Skauge et al. 2012(Skauge et al. , 2013(Skauge et al. , 2020. Two examples of immiscible fingering from this school are shown in Fig. 2 from the work of (a) Vik et al. (2018) where ( o ∕ w ) ≈ 500 and (b) Skauge et al. (2012) where ( o ∕ w ) = 2000 , although both of these groups looked at a wide range of viscosity ratios. In these experiments, X-ray imaging and X-ray saturation measurements were used both to visualise the viscous fingers and also to measure the local water (and oil) saturations in the 2D rock slabs, as shown in Fig. 2c. It is this direct calibrated measurement of water saturation in the viscous fingers that is crucial to the modelling developments in this paper. Again, this work confirms directly that the water saturation in the fingers, S w , is "high"; direct evidence of this is given in several of the papers from this group (Skauge et al. 2012(Skauge et al. , 2013Vik et al. 2018). For example, point A in Fig. 2c shows that the water saturation in the finger can be S w ~ 0.58.

1D Buckley-Leverett Analysis and the "M-Paradox" for Unstable Displacement
Linear two-phase immiscible displacement is frequently analysed in 1D using the Buckley-Leverett (BL) equation which applies in the viscous dominated limit of the flow equations, i.e. when capillary pressure, P c = 0 (Buckley and Leverett 1942). In 1D incompressible flow, there is no requirement for a pressure equation (only pressure differences can be calculated from the flow solutions), and the BL equation for the development of the S w (x, t) profile is given by: where v t is the total velocity, v t = Q∕ (A. ) , Q is the volumetric flow rate, A the crosssectional area and ϕ the porosity. The function f w S w is the fractional flow of the water (displacing) phase within the porous medium which is given by the following well-known variant forms of the same equation (Dake 1978):  (1958). N p is the number of PV injection of the water phase An important quantity governing viscous instability in an immiscible system is the mobility ratio, M S w = w o = k rw . o k ro . w (Dake 1978). Note that, if the magnitudes of the relative permeabilities were similar (i.e. k rw ≈ k ro ), then M would be Vik et al, 2018 (SPE 190866) Immiscible viscous fingering in a high permeability Bentheimer rock slab where the images are captured by X-ray imaging and X-ray local saturation measurements and, by suitable calibration, water saturations can be measured; a is from Vik et al. (2018), b is reconstructed from Skauge et al. (2012) and c is reconstructed from Skauge et al. (2013) (almost) independent of S w and M ≈ o w >> 1 for unstable immiscible displacement of oil by water. Indeed, there is one common convention to take the relative permeability terms in M to be at their respective relative permeability (RP) endpoints where k rw and k ro are evaluated at S w = 1 − S or and S w = S wi , respectively, that is: where S wi and S or are the initial water saturation and the residual oil, respectively. Defining M S w in this way would give a value of M close to M ≈ o w >> 1 ; the authors would describe this as "conventional wisdom".
It is well known that the BL equation leads to a shock front solution with a water saturation shock front height, S w = S wf , which can be found from the fractional flow, as shown in Fig. 3 where the tangent construction gives the S wf as shown (Dake 1978). In this figure, "normal" relative permeabilities based on the widely used Brooks-Corey (BC) form are given, i.e. k rw S w = w . S wn 1∕ w and k ro S w = o (1 − S wn ) 1∕ o , where the normalised water saturation, S wn , is given by S wn = S w − S wi 1 − S or − S wi Corey 1964, 1966).
To illustrate what is predicted by 1D BL theory, Fig. 3 shows "conventional" (Brooks-Corey, BC) oil/water relative permeabilities, k rw S w and k ro S w , as functions of S w , and the corresponding fractional flow curves f w S w using these relative permeability functions for oil/water viscosity ratios, ( o ∕ w ) = 1600 and 5. We would expect the higher viscosity ratio case to lead to viscous instability and the lower case to be much more stable or at least should show less pronounced fingering. The S w (x, t 1 ) profiles along the system at dimensionless time t = t 1 = 0.05 predicted by BL theory are shown in Fig. 1b for each viscosity ratio, and the calculated values of shock front height,S wf , and the value of the mobility ratio at this shock front, i.e. M(S wf ) , are shown in this figure. The S w (x, t) profiles in Fig. 1b show that shock fronts form for each viscosity ratio and, as is well known, a higher S wf is observed for the ( o ∕ w ) = 5 case than for the 1600 ratio case. However, as a result of the calculated values of S wf , the rather counter-intuitive result is that the value of S wf for the high viscosity contrast case yields a value of M(S wf ) at the shock front which is quite low; M(S wf ) = 1.22 in this case, which is "stable". On the other hand, the more "stable" case for the viscosity ratio of 5 gives M(S wf ) = 6.60. We refer to this behaviour of the two ( o ∕ w ) cases here, as the "M-paradox". Fig. 3 a Conventional relative permeabilities,k rw and k ro , and the associated fractional flow curves, f w,5 and f w,1600 for o ∕ w = 5 and 1600, respectively, and b the corresponding 1D water saturation profiles, S w (x, t = t 1 ), where t 1 = 0.1; the corresponding Buckley-Leverett shock front saturations, S wf , and the mobility ratio at this front, M S wf , are also shown The name "M-paradox" becomes even more well-deserved when we consider the results of a very fine grid 2D numerical simulation of "viscous fingering" using the conventional relative permeabilities in Fig. 3 for the high viscosity contrast case, ( o ∕ w ) = 1600. An example 2 D simulation for this two-phase case of a water → oil displacement with this high viscosity contrast is shown in Fig. 4. In this figure, the initial water saturation is S wi = 0.17 (where k rw S wi = 0 ) and in the location in the middle of the finger (shown by the small arrow), then S w = 0.205. Thus, in the middle of the finger, the water saturation is barely above the initial value ( ΔS w = S w − S wi ≈ 0.035 ). The simulation results in very low water saturations in the fingers yielding "wispy" or "ghost" fingers which do not at all resemble experiments such as those shown in Fig. 1 or 2 which show much higher S w values in the main fingers. Later in this paper, we will offer a resolution of the M-paradox.
Thus, the M-paradox is that, for very adverse viscosity ratio, ( o ∕ w ) >> 1, a low BL shock is formed in 1D two-phase flow which has a low value of M(S wf ) ≈ 1, which is in turn "stable". The corollary is that as we attempt to make the displacement more stable by greatly reducing ( o ∕ w ) , for example to ~ 5, the M(S wf ) becomes greater, i.e. about ~ 6.6 in this case, which is hence apparently less stable.
The BL solution is, as discussed, a 1D solution and viscous fingering is an inherently 2D (or 3D) phenomenon. Thus, the BL solution misses any 2D flow mechanisms but "compensates" for the fingering in an "averaged" manner by giving the low saturation "pseudo-BL finger" shown in Fig. 3b. In fact, we conclude that 1D BL theory is at best somewhat misleading, and at worst simply incorrect, when viscous fingering occurs since this is a 2D phenomenon, and this has also been suggested recently by others (Luo et al. 2017;Sorbie and Skauge 2019). A "middle road" interpretation is that BL theory is adequate, but some care must be taken in one's choice of relative permeability functions to ensure that physically sensible results are obtained.

Fig. 4
Simulation of "viscous fingering" with conventional relative permeabilities showing the "wispy" or "ghost" fingers where the S w in the water finger (shown as red) is just above ( ΔS w ~ 0.02-0.03) the initial (immobile) water saturation, S wi = 0.17 1 3

Method for Modelling Immiscible Viscous Fingering
A number of studies modelling immiscible fingering using direct numerical simulation have appeared in the literature, which were nicely summarised recently in the work on immiscible fingering by Bakharev et al. (2020). Although this particular study was both careful and thorough, the forms of the experimental relative permeabilities do not lead easily to fingering of the sort observed in the work quoted in Figs. 1 and 2 of this paper, especially since a homogeneous randomly perturbed permeability field was used in their work (Bakharev et al. 2020).
To develop an improved approach for the direct numerical simulation of immiscible viscous fingering, we start from the conclusion that 1D BL theory is inappropriate, or at least should be used with caution, for any process or mechanism which is fundamentally 2D (or 3D) in nature, as viscous fingering clearly is. We also note that, proceeding with "normal" relative permeabilities leads to the incorrect fractional flows observed in the various regions in viscous fingering, i.e. in the fingers and in the bypassed zones. Physically (see Figs. 1 and 2), we observe that the water saturation S w is quite high in the fingers and is very low in the bypassed regions of high S o (where S w = S wi ). The 4-part approach we propose to correct these shortcomings, which arise from the M-paradox, is summarised as follows: (i) We start from a chosen fractional flow function, denoted f * w , where we specify a higher "shock front" saturation,S * wf , than is found from the conventional relative permeabilities. However, from f * w alone, only the ratio of relative permeabilities can be found, and it is straightforward to show that However, if the form of one of the RPs is assumed, for example, k * rw , then the other RP, k * ro in this case, is given by Alternatively, we can choose generalised forms of both k * rw and k * ro such that we obtain the f * w (and thus S * wf ) we require. But this still gives an infinite set of potential RP pairs ( k * rw and k * ro ) which give the same f * w , and this is addressed next. (ii) The new fractional flow function, f * w , can be chosen to agree with experiment in terms of the water saturation in the finger, i.e. to avoid the "ghost" fingering situation illustrated in Fig. 4. However, we note in (i) above that there is an infinite choice of RP functions which yield the same f * w . Below, we describe a very general parameterisation of the RP functions used to represent k * rw and k * ro , which enables us to generate a wide range of such function (for the same f * w ). The clue to which is the "correct" set of RP functions comes from two connected sources, (a) the conjecture that the "correct" set of RPs is that which maximises the total mobility function, T (where T = o + w ), i.e. minimises the pressure drop across the fingering system, and (b) from linear stability analysis of the two-phase unstable system. This is discussed in more detail below.
(iii) We then choose a heterogeneous flow domain, which we define here as a random correlated field (RCF) with a given level of permeability (k) heterogeneity (i.e. the variance or range of k, which can be described by a Dykstra-Parsons coefficient or simply as a (k max /k min ) ratio) and a correlation structure defined by its dimensionless correlation length, D = ( ∕L) , where and and L are the correlation length of the permeability field and the system length, respectively. Other choices are also possible to match the actual permeability fields and correlation structures for specific cases. However, the RCF incorporates both heterogeneity and structure in a simple and quantifiable manner. [*NB The Dykstra-Parsons coefficient, V DP , is a measure of the degree of permeability heterogeneity. If permeability is plotted on a log probability plot then, V DP = (k 50 − k 84.1 )∕k 50 , where k 50 is the 50% probability (median) permeability value (mD) and k 84.1 is the permeability (in mD) with probability of 84.1% (at one standard deviation).] (iv) We then choose a sufficiently fine spatial grid (in 2D Δx and Δy ) for the numerical simulations where (Δx∕L) ~ (Δy∕L) << D (where ~ implies "of order"). If a conventional numerical reservoir simulator is used using very simple single-point upstream weighting for the spatial derivatives in the two-phase transport equation, then the level of dispersivity in the calculation is of order, ≈ (Δx∕2) (Peaceman 1977).
The above "recipe" for modelling immiscible viscous fingering may at first appear somewhat arbitrary, but it has a very strong physical rationale, and each aspect of this approach is discussed in turn below. The first point on the choice of fractional flow, f * w , is to "correct" the "error" in the 1D BL equation, if conventional RP functions are used. This is similar in some respects to the approach reported in Luo et al. (2017) who also adjust the fractional flow to give a more appropriate curve to describe the averaged behaviour of the system in the presence of viscous fingering. Luo et al. (2016) then proceeded to match this adjusted fractional flow curve to production data for unstable immiscible water flooding and also polymer flooding. They did not take the following steps to actually simulate the fingering directly, as we do in this work.
Secondly, from the infinite set of possible RP functions, k * rw and k * ro , consistent with the chosen f * w , we choose the set which maximises the T function, this ensures that the resulting finger pattern has the lowest pressure drop possible. The pressure drop across the system has the formal structure and so it is strictly the functional ∫ dx T which should be minimised. However, this is not really a 1D integral and the actual calculation required for any 2D (or 3D) spatial saturation distribution, S w (x, y), involves the pressure drop calculation, i.e. which is the pressure equation presented in Sect. 1 (with Pc = 0). For the ΔP to be a minimum across the system (i.e. that the path of least resistance is followed by the fingering system), then T must be a maximum. We regard this as a type of "variational principle" for the fluid flow and present it here as a conjecture. It is consistent with the very sharp pressure drops that are observed in all immiscible fingering experiments. Further evidence that the maximum mobility assumption is correct is provided by the results from linear stability analysis for immiscible fingering by Chikhliwala and Yortsos (1988); they show that the criterion for unstable finger growth are equivalent to conditions on the relative permeabilities which are met by the maximum mobility principle. This issue of maximum T will also be supported by further evidence in a forthcoming publication. It also turns out that for the analytic forms for the RPs that we choose, it appears we can always find this optimal T function. Chikhliwala and Yortsos (1988) derived some conditions, based on linear stability theory, on the form of the RP curves on the onset of viscous instability in two-phase flow. We have found that their conditions indicate the same choice of RP as would lead to this same condition on T . Intuitively, it might be thought that this condition would be best met by having the largest mobility (i.e. the largest k * rw ) for the water phase and this is indeed partly true, but the fact that the 2 RP functions are constrained (by f * w ), means that there is in fact a maximum physical T that is achievable. This issue will be developed in much more detail with supporting calculations in a forthcoming paper.
Thirdly, in point (iii) above, we recognise that the fingering process is quite random, but in a 2D domain we can match the correlation structure to the possible numbers of fingers that initially emerge. Thus, if we take a smaller D , then we expect to see more initial fingers and we would expect the actual number of fingers in the early time, N f to be of order, N f ≈ 1∕ D . If this latter conjecture is approximately true, then taking (Δx∕L) = (Δy∕L) << D , ensures that the grid Peclet number, N Pe,grid , fully resolves the emerging fingers; this quantity is given as follows: where D num is the level of numerical dispersion ( D num = .v ), and since ≈ Δx∕2 , then N Pe,grid ≈ 2N x where N x is the number of grid blocks (in the flow direction), Thus, by definition, we are simulating the system at a high mesh Peclet number which is more than adequate to resolve the emergent fingering, since we have ~ D L∕Δx grid blocks per average finger. A second way of viewing this matter follows the discussion introduced at the end of Sect. 1 of this paper, and concerns the capillary length, L cap . This is the length scale below which capillarity (from P c ) would "wash out" very small viscous fingers, and hence is equivalent to a certain level of capillary dispersion, D cap . It is shown elsewhere (Stephen et al. 2001) that the form of this nonlinear capillary dispersivity is: where we note it is the slope of the capillary pressure that governs the capillary dispersion. In our approach, we essentially identify D num ≈ D cap . In fact, D cap is a nonlinear function which is zero at both S w = S wi and at S w = (1 − S or ) (since the mobilities λ w and λ o are, respectively, zero at these limits), and hence, it has a maximum value of D max . The quantity D max can be calculated from the relative permeabilities and the slope of the Pc, and hence, the corresponding L cap (dispersivity) can be found. The fine grid defines the level below which capillary effects resolve the immiscible fingering; thus, by definition, the simulation of viscous dominated fingering may be carried out, as shown below.
For a given permeability field, we find that this is perfectly adequate down to dimensionless correlations lengths of D ≈ 1∕100 since we can obtain very similar behaviour (number of emergent macroscopic fingers, breakthrough times and recovery profiles) under grid coarsening (as demonstrated later in this paper). In a later paper, we will describe the effects of capillary pressure explicitly for both water-wet and oil-wet systems, since rather different behaviour is observed for each of these cases (Stokes et al. 1986).
There is an additional (unforeseen) effect which also helps in the numerical resolution of the emerging immiscible fingers. From point (i) above, when we choose the fractional flow function, f * w , and then derive the underlying RP functions, it turns out that these functions are "forced" to have a rather restrictive form which also helps to better define the fingers and to suppress most of the numerical dispersion arising in the system. This is not obvious until some of the numerical results are presented below.

and Grid Resolution
In the immiscible water → oil displacement simulations presented here, we take the very adverse viscosity ratio, ( o ∕ w ) = 1600 . The general analytical forms of the relative permeability curves used to construct the fractional flow function, f * w , used in the immiscible fingering simulations are as follows: where S wn and S on are the normalised water and oil saturation, respectively, where S on = 1 − S wn and S wn = ; the constants w , o , w and o and the normal Brooks-Corey parameters and w , o , w and o are additional constants which allow us to have some flexibility in the forms of the RPs, and hence on the f * w function. We have found that these empirical forms of the RP functions give us sufficient flexibility to generate any chosen modified fraction flow function, f * w , which we wish. In addition, the above RP analytical equations reduce to the Brooks and Corey forms when the simulation in Fig. 4). For the base case simulation of fingering (Case 1), we construct the fractional flow function, f * w , shown in Fig. 5c which is based on the RP functions which are also shown in Fig. 5a and b on linear and a semi-log scales, respectively; as noted above, ( o ∕ w ) = 1600 . The various parameters used to generate the RPs are also given in this figure. In this case, we chose the RPs by trial and error in order to obtain a flood front saturation height of S wf = 0.275 which is sufficiently well above the initial water saturation, S wi = 0.17 to generate a higher water saturation in the fingers. Note also here that the M S w value at the 1D shock front, where S w = S wf , is M S wf = 19 , which is sufficiently high to be unstable, but is far below the "conventional" value of mobility ratio M ≈ o ∕ w ≈ 1600 . An additional Case 2 was also chosen with the higher shock front saturation 0.34S wf = 0.34 , all functions and parameters for this case are shown in Fig. 6. We note that these calculated RPs may now be thought of as "pseudo-relative permeabilities" which arise from choosing the fw* function and then maximising the λ T . Features of the RPs such as crossover point and end point levels are "lost" since these characteristics refer to "normal" relative permeabilities measured under viscous stable conditions. An additional fractional flow plot is shown in both Fig. 5c (Case 1) and Fig. 6c (Case 2) where a viscous "polymer" with w replaced by p = 25 w is shown. The shock front tangent construction is slightly different in this case, as discussed in Pope (1980) and Sorbie (1991), and for Case 1, this results in a slightly higher shock front S w = S wf = 0.306 and a mobility ratio of M S wf = 14.625 . This is a little more stable than the waterflood case (where M S wf = 19 ) but it does not appear that injecting polymer would make a very significant difference in this case to increasing the displacement of oil. We will return to this issue in the results section of this paper. For Case 2, the corresponding polymer values are S w = S wf = 0.395 and a mobility ratio of M S wf = 20.96 . It may be surprising to find that the frontal mobility value for the polymer in Case 2 is actually higher than that of water, but this is partly to do with the manner of the polymer construction (see Pope (1980) and Sorbie (1991)) and the frontal S wf value for the polymer is still above that of the water. The chosen Case 2 RP functions, k * rw and k * ro on a a linear and b a log-linear scale along with c the generated modified "fingering" fractional flow function, f * w , giving a higher shock front saturation S w = S wf = 0.34, in this case; initial water saturation is S wi = 0.17 and residual oil (where k ro = 0) is S or = 0.2 For the permeability field, we choose a square domain where L = L x = L y , where the dimensionless length, L = 1. Three random correlated fields are chosen with the same base case level of heterogeneity in terms of the Dykstra-Parson coefficient ( V DP = 0.66) with dimensionless correlation lengths, D = 1/30, 1/15 and 1/10 as shown in Fig. 7. The permeability range is shown in inset of Fig. 7 and ranges from k min = 10 mD to k max 1000 mD, and thus, the base case k max /k min = 100; this is reduced in some of the simulations presented below. Specifying the grid structure with N x = 2000 in the direction of flow and N y = 1000 in the transverse direction fully defines the system. A standard numerical reservoir simulation model was used for these calculations which employed a very simple numerical scheme (single-point upstreaming) for the two-phase transport equations. Sufficiently, short time steps were taken to ensure the time truncation error was negligible.

Water → Oil Viscous Fingering Simulations
The system for viscous dominated two-phase immiscible flow is now fully defined; the relative permeabilities (RPs) are given for Case 1 in Fig. 5 and for Case 2 in Fig. 6. All RPs were derived from the desired f * w curves also shown in Figs. 5c ( S wf = 0.275) and 6c ( S wf = 0.34), the viscosity ratio is ( o ∕ w ) = 1600 and the 3 heterogeneous permeability RCFs, with λ D = (1/10), (1/15) and (1/30), are given in Fig. 7. The 2D system is square with dimensionless side lengths L = L x = L y = 1 and the grid size is N x = 2000 in the direction of flow and N y = 1000 (perpendicular to it). Water injection starts with system at initial water saturation, S wi = 0.17 (where k rw = 0).

Case 1 Simulations
The results of the simulation of these unstable water → oil displacements for Case 1 are shown in Fig. 8, at 5 times (PV injected). We recall that the water shock front saturation ( S wf ) from the fractional flow curve for Case 1 in Fig. 5 is S wf = 0.275 , where the mobility ratio, M S wf = 19 . In Fig. 8, the water saturation colour scheme shows red for Fig. 7 The 3 heterogeneous permeability random correlated fields (RCF) used in the viscous fingering simulations, where the dimensionless correlation lengths are, from left to right, D = (1/10), (1/15) and (1/30). The permeability range is shown in inset and ranges from k min = 10 to k max 1000 S w = 0.3-0.4 and black is the initial water saturation, S wi = 0.17 (where k rw = 0 ). It is clear that the viscous fingers are very clearly defined and sharply resolved in these simulations, and several observations can be made. As conjectured earlier in this paper, the early time results in Fig. 8a at T = 0.005PV show the emergence of many small fingers with clearly more of these in the shortest correlation length case where λ D = (1/30) and fewest in the largest λ D = (1/10) case. As the displacement progresses, then for the longest correlation length case (λ D = (1/10)), fewer dominant fingers emerge more rapidly and break through earlier, as shown in Fig. 8b-d (from 0.01 to 0.045PV). The corollary is that, for the shortest λ D = (1/30) case, a much larger number of smaller fingers persist for longer and, as a result, breakthrough is later. At the latest time shown here in Fig. 8e at time = 0.08PV, it is clear that there are larger regions of defending fluid (oil) which have been bypassed by the water fingers for the longer correlations length cases with λ D = (1/10). A consequence of all this behaviour is observed in the oil recovery (%) versus time (PV) results and the water production profiles versus time (PV), both shown in Fig. 9. These results show that the longer the correlation length, the less efficient the oil recovery is, and earlier the breakthrough and rise of water fraction in the produced fluids from the system. Figure 9 also shows that, after 0.5PV of injection, the oil recovery is quite low at ~ 12-14% of the original oil in the system; after 2PV of injection, this only rises to ~ 14.5-16% (not shown).
The corresponding dimensionless pressure drop (ΔP D ) across the whole system is shown for the 3 simulation cases in Fig. 10, where ΔP D is defined as, ΔP D = ΔP/ΔP 0 , where ΔP is the calculated pressure drop as a function of PV injected and ΔP 0 is the pressure drop at the same injection rate initially when only (viscous) oil is flowing. Again, this agrees well qualitatively with experiment in form (Skauge et al. 2012(Skauge et al. , 2013Vik et al. 2018), but in our calculation we find, as expected, the longer the dimensionless correlation length, the earlier the breakthrough (Fig. 8) and the sharper the drop in ΔP D .

Case 2 Simulations
A second case is presented for the details described in Fig. 6 (Case 2) for an f * w with higher frontal saturations of S wf = 0.34 . The purpose of this second case is first to show an example with a higher water saturation in the fingers than in Case 1. In addition, the Case 1 set of RPs in Fig. 5 was "close" to the lowest pressure (maximum mobility, T ) criterion discussed in Sect. 4, but the Case 2 set of RPs (shown in Fig. 6) has genuinely maximised T . This is shown in Fig. 11 where the results for 4 sets of RPs (numbered 1-4) constrained to have the same fractional flow function, f * w (Fig. 11b), are used to calculate the total mobility functions shown in Fig. 11a. The Case 2 RPs in Fig. 6 give a T function which coincide with the dashed line, which is the maximum possible mobility that is achievable for this case. To find this maximum mobility, the end points of the RPs were fixed (at S wi and S or ) and all other parameters in the RPs were allowed to vary, until a max λ T was found consistent with the given fw*. For a given experimental case, the fw* would be chosen to match the saturations in the fingers, the RPs were then found by maximising mobility (λ T ) and with these RP functions, all other quantities  Fig. 7 with λ D = (1/10), (1/15) and (1/30). The saturation patterns are shown at the dimensionless times (in PV) indicated in (a-e). The water shock front saturation ( S wf ) from the fractional flow curve in Fig. 5c is S wf = 0.275, and the water saturation colour scheme shows red for S w = 0.3-0.4 and black is the initial water saturation, S wi = 0.17 (where k rw = 0) emerge from the normal immiscible simulations, i.e. the fingering patterns, oil recoveries, watercuts and pressure drops. This matter will be discussed in more detail and its relation with linear stability analysis will be explained further in a forthcoming paper.
For the Case 2 RPs, and with all other data as in Case 1, the water → oil displacement patterns are shown at 4 times in Fig. 12 for the 3 permeability models in Fig. 7. Note that the red colour now corresponds to a higher saturation than in the Case 1 displacements in Fig. 12, but that the same qualitative features are observed for Case 2 as for Case 1, e.g. more earlier  Fig. 9, for each of the 3 permeability fields in Fig. 7 fingers for the shorter correlation length permeability field, etc. Since the water saturation is higher in the fingers in Case 2, then a more efficient oil displacement by water and a later breakthrough time is expected and this is confirmed by the results in Fig. 13. The oil recoveries now range from ~18 to 22% after 0.5 PV of injection which is again quite typical for such a viscous oil. Also, as a result of the higher Case 2 water saturations, the pressure profiles, plotted as dimensionless pressure drops (defined above) in Fig. 14, are higher than observed for Case 1, although they have the same qualitative order for the different λ D fields.

Simulation of Polymer Flooding of the Unstable Water → Oil Displacement
The injection of low-viscosity water into a very viscous oil results in unstable displacement, as in the two examples presented above (Cases 1 and 2) where the ratio ( o ∕ w ) = 1600 . Polymeric solutions of a few 100 s to a few 1000 s of ppm (mg/L in Fig. 11 a The total mobility functions, f * w , and b the fixed fractional flow, f * w , for the various example RP curves that give this f * w ; the Case 2 RP curves, k * rw and k * ro , (see Fig. 6) give the maximum f * w function possible (minimum pressure drop) as shown in (a) where the blue line and dashed line coincide  Fig. 5c is S wf = 0.34, and the water saturation colour scheme shows red for S w = 0.35-0.5 and black is the initial water saturation, S wi = 0.17 (where k rw = 0) water) can greatly viscosify the aqueous phase and we may ask if increasing the viscosity of the aqueous phase by a factor of × 25 or so would make any difference. In Figs. 5c and 6c, the construction to the Case 1 and Case 2 fractional flows is shown, respectively, for a viscous "polymer" with w replaced by p = 25 w (Pope 1980;Sorbie 1991) . For Case 1, this results in a slightly higher shock front S w = S wf = 0.306 and a mobility ratio at the front of M S wf = 14.625 ; cf. the waterflood M S wf = 19 . This suggests that injecting polymer would not make much difference to increasing the displacement of oil, in this case. Indeed, when applying extended (1D) BL theory to test this out, it  Fig. 12, for each of the 3 permeability fields in Fig. 7 appears that a modest improvement in displacement efficiency (~15 to 25%) is achievable as shown in Sorbie and Skauge (2019). However, in a series of experiments performed in sandstone slabs, Skauge et al. (2012) showed that the level of improvement by injecting polymer was really quite remarkable (+ 50 to 70%). The X-ray analysis made it possible to visualise changes in local saturation, and accelerated oil production was explained, and visualised as crossflow of oil that moved into the established water channels (i.e. into the path of least resistance). Salmo andSkauge (2015, 2017) derived relative permeabilities by history matching experiments dominated by crossflow. Seright et al. (2018) also noted that a polymer with p = 25 w greatly improved the oil displacement by ~ 60%+ for an oil with o = 1600 w . Our own analysis of Seright et al.'s result showed that about 1/3 of this increase can be explained by extended BL theory (Sorbie and Skauge 2019). Another important observation from the work of Vik et al. (2018) was that, when polymer was injected, fingering still occurred (i.e. it was not completely suppressed, although it was a little less than in the waterflood). In addition, this experimental oil recovery increase was observed very rapidly (in PV) after tertiary* polymer injection (*when the polymer solution was injected after the original water fingers had broken through at the system outlet). The viscous crossflow phenomenon is discussed in more detail in Sorbie and Skauge (2019).
We now test if our approach can actually predict such a response by taking one of our viscous fingering simulations for the Case 1 RPs, the D = (1/10) case, since this gives the least efficient recovery, and injecting a viscosified polymer of viscosity, p = 25 w , thus reducing the ratio ( p ∕ w ) = 64 . However, the fingers have already broken through by 0.045PV of water injection (see Fig. 7d) and, by 0.1PV of injection, the watercut was already ~ 0.85 (i.e. 85% of the produced fluid is water) as shown in Fig. 8. In the polymer fluid injection at 0.1PV, it is stressed that nothing else is changed in this simulation except the viscosity of the injected aqueous (polymer) phase.
The initial pattern of viscous fingers for the water displacing oil up to T = 0.1PV is exactly as already shown for Case 1 in Fig. 8. Results of the polymer injection (at 0.1PV) are compared with those for the continued injection of water in Fig. 15. The distribution of the (red) water fingers is shown after 0.1PV just before polymer injection in Fig. 15a; the subsequent development of the fingering/water saturation distribution is shown in Fig. 15b-f at the PV of water/polymer solution injection indicated. Figure 16 shows the oil recovery versus PV and the watercut versus PV injected for both the water injection only and the 0.1PV of water followed by polymer injection.
The results showing the water saturations in Fig. 15 show that even after polymer injection, the fingering pattern persists (as observed experimentally), but the water saturation distributions are certainly modified in 4 important ways, as follows: (i) Firstly, the saturations along the finger "backbones" increase, as evidence by the higher S w "tendrils" along the main fingers clearly seen in Fig. 15d at 0.5PV to 15f at 1PV. This is due to direct polymer displacement along the highly conducting fingers. (ii) The more viscous polymer moving along these fingers causes sweep of the bypassed oil into the finger ahead of the polymer, from whence it is produced very rapidly. This is very clear evident of the viscous crossflow mechanism (see Sorbie and Skauge 2019). Direct evidence of viscous crossflow can be seen visually by carefully comparing the S w values at the edges of the waterflooding and polymer flooding fingers in Fig. 15. In particular, observe the edges of the "black" (high So) zone where this effect is clearest. This viscous crossflow effect is even more clearly evident in the animations of the dynamic displacements.
(iii) The bypassed region shrinks and is more invaded by water flowing out of the main flowing channels. Again, this is the other part of the viscous crossflow mechanism. (iv) Close to the inlet of the system (Fig. 15d-f), S w increases significantly (white region) due to direct displacement by the viscous polymer. This region effectively "feeds" the downstream region with fluid.
As a result of the above fluid flow process-i.e. the superposition of the viscous crossflow mechanism and the viscous fingering mechanism-the incremental oil recovery is very significant as seen in Fig. 16. This figure shows that the % improvement in oil displacement relative to the waterflood alone is 61.5% at 1PV and this goes to a 73.8% increase at 2PV. Both the extremely large response in terms of oil recovery and the speed of this response are very close to the values of oil recovery and the response times observed experimentally (Skauge et al. 2012(Skauge et al. , 2013Seright et al. 2018).
We repeat that no further adjustment of the RP functions ( k * rw and k * ro ) was required to achieve this result, i.e. no special "polymer RP functions" were required. Simply changing the water viscosity was sufficient to obtain the calculated "polymer" responses observed.

Simulations Adjusting the Underlying Fractional Flows, f * w
In the sensitivity calculations presented here, we take an example with an even higher water saturation within the finger (i.e. a higher S wf ) than the 2 examples presented above (where S wf = 0.275 or 0.34, for Cases 1 and 2, respectively). In this case, we choose S wf = 0.46 but this can be done using an infinite range of fractional flow functions as indicated in Fig. 17, where 4 possible f * w versus S w are given (Fig. 17a) along with the corresponding values of the frontal mobility ratio, M(S wf ) , shown in the semilog plot of the universal curve of f * w versus M (Fig. 17b). The underlying relative permeabilities from the f * w functions in Fig. 17a are not presented here. The 4 examples in Fig. 17 all have S wf = 0.46 ( S wf =0.45 for no. 3) but they have different frontal M values of, M(S wf ) = 1.38, 3.66, 6.60 and 11.09 for nos. 1-4 in Fig. 17, respectively.
Simulations were carried out in the CRF permeability model with λ D = (1/15) ( Fig. 7), again with the viscosity ratio ( o ∕ w ) = 1600 and the same grid as in the base case for these 4 examples, and results are presented in Fig. 18. Full simulations of these 4 cases have been carried out but here only a snapshot of the phase displacement fingering patterns is shown at 0.120PV for each numbered example in Fig. 18. As expected, the no. 1 simulation with M = 1.38 is stable and has a high ( S w > 0.55) and "compact" water saturation behind the front, and hence, the observed retardation of this front relative to the other cases in Fig. 18 at the same time. The variation of the front in no. 1 is simply due to permeability heterogeneity of the CRF, with a "frontal roughness" of order ~ D . Some degree of actual viscous instability is seen in the other 3 examples, although no. 2 with M = 3.66 still shows relatively high water saturation levels behind the front. Numbers 3 and 4, with M = 6.60 and 11.09, respectively, in Fig. 18 show gradually increasing instability and fingering with well-defined fingers forming in each case.
The main objective of these calculations is to demonstrate that, by manipulating the input fractional flow function, f * w , then we can effectively chose the S wf which we believe gives the correct saturations within the fingers and also the frontal mobility ratio, M S wf , which characterises the degree of immiscible instability.

Viscous Fingering Simulation Results Under Permeability Field Rescaling and Grid Coarsening
In order to test the robustness of the viscous fingering methodology proposed in this paper, and to support the arguments presented in Sect. 4, we apply two numerical tests to our method, using Case 1 as an example, as follows: (i) Firstly, we take the D = 1/15 RCF which, in the Case 1 model, has a permeability minimum of k min,1 = 10 and maximum of k max1 = 1000, and we perform an order preserving rescaling of the permeability field maintaining its exact structure, first to k min,2 = 100 and k max,2 = 1000, and then to k min,3 = 333 and k max,3 = 1000. To achieve this, the simple rescaling to the k i2 distribution, for example, is given by k i2 = k min,2 + k i1 −k min,1 k max,1 −k min,1 k max,2 − k min,2 ; and this is applied, where the old value of k i,1 in distribution 1 is replaced by k i2 in the new distribution 2. This is to ensure that it is not the permeability contrast that is dominating the fingers, but the field structure ( D ), and some degree of heterogeneity. (ii) Secondly, we again take Case 1 with D = 1/15, and gradually 2 × 2 coarsen it from the base case simulation with N x = 2000, N y = 1000, to N x = 1000/N y = 500, and to N x = 500/N y = 250. The coarsening scheme takes the arithmetic average of the cells as the coarsened block permeability, and for transmissibility calculations, it takes the harmonic average across the adjacent coarse block interfaces.
For permeability field (CRF) rescaling, the CRFs look identical (not shown). The simulations of the actual fingering patterns look very similar in each case but are not identical, as shown by the example fingering patterns. The (k max /k min ) = 100 and 10 cases in Fig. 19a and b compared at times T = 0.025 and 0.08PV, respectively, look very similar; but the (k max /k min ) = 3 case in Fig. 19a and b look somewhat different with more "linear" fingers. At the earlier time before breakthrough (0.025PV), the number of fingers is very similar in all 3 cases and the major fingers are in the same places, but the fine details differ slightly. After breakthrough at 0.08PV of water injection, the pattern and quantity of bypassed oil are again approximately the same. However, the oil recoveries and watercuts shown in Fig. 19c are very close, with recoveries at 0.5PV being within ~ 2% of each other, despite the fact that the k max ∕k min has been reduced from 100 to 3. Likewise, the profile of normalised pressure drop, ΔP∕ΔP 0 versus PV injected for all (k max /k min ) = 100, 10 and 3 cases are almost identical as shown in Fig. 19d. Figure 20 shows the effect of grid coarsening on the finger patterns and on the oil recovery and watercut profiles for the D = 1/15 model (Fig. 7). The finger patterns are shown at times T = 0.025 and 0.08 PV in Fig. 20a and b, respectively, for the 3 numerical grids N x /N y = 2000/1000, 1000/500 and 500/250. These patterns are again very similar but not absolutely identical in detail. However, the oil recoveries and watercut profiles versus PV shown in Fig. 20c are almost identical. The conclusion from these results are that as long as we preserve the conditions in our 4-part approach, viz. (i) the form of the f * w (derived from the k * rw and k * rw ), (ii) the maximum mobility ( T ) RP functions honouring f * w , (iii) the permeability field as a CRF with a level of heterogeneity and structure ( D ) and (iv) a sufficiently fine grid, then the procedure is robust. Furthermore, it captures the main physics of the viscous fingering process.

Summary, Discussion and Conclusions
Generating realistic immiscible viscous fingering by direct numerical simulation discretising the PDEs (partial differential equations) describing the process has long been recognised as being a difficult task. This has been viewed by many as a problem of "the numerics", and many studies have appeared using higher-order numerical schemes of various types, with some limited success. However, the present authors do not view this as being the central issue in modelling immiscible fingering. We view it as being a problem of establishing the correct physics and related mathematics. Specifically, the problem arises due to a shortcoming in 1D fractional flow theory in describing fingering systems with "conventional" RPs, where the mechanism is inherently 2D (or 3D). We have described the consequences of this as the "M-paradox" which in turn when using "conventional" relative permeability (RP) functions reproduces the erroneous (or at best "averaged") Buckley-Leverett 1D result.
In this paper, we propose a novel 4-step methodology for simulating immiscible twophase viscous fingering in porous media. This is based on recognising the M-paradox and the experimental observation that the water saturation in displacing fingers is relatively high. Therefore, the 4 stages of this approach, involve (i) choosing the form of the f * w to give sufficiently high S wf in the fingering zone, (ii) recognising that this still only constrains Fig. 19 Water displacing oil viscous fingering simulations for the base case data in the D = (1/15) CRF for the original permeability range (k min1 = 10 and k max1 = 1000) and the rescaled ranges (k min2 = 100 and k max2 = 1000) and (k min3 = 333 and k max3 = 1000) at a T = 0.020PV, b T = 0.075PV along with c the corresponding oil recovery versus PV and watercut versus PV and d the normalised pressure drops, (ΔP/ΔP o ), for each case the ratio k * ro ∕k * rw which is satisfied by an infinite number of possible sets of RP, then we choose the set which shows the maximum total mobility function, T (iii) establishing the appropriate permeability field as a RCF (in this case) with a level of heterogeneity ( k variation) and structure ( D ) and (iv) simulating the process with a sufficiently fine grid, using very simple conventional numerical techniques (e.g. single-point upstreaming). It is shown that following this scheme we can generate very realistic immiscible fingering which shows the detailed qualitative features of immiscible fingering patterns observed in 2D experiments.
Considering the choice of f * w for a given simulation of immiscible fingering in experiments, we believe that it is probably best to match f * w to the observed saturation levels in the fingers through a "history matching" procedure. The chosen S wf in f * w can be successively increased in order to match the (relatively high) water saturations in the experimentally observed viscous fingers. This approach is currently being pursued to model the experiments of Skauge et al. Results are illustrated by performing 2D ( L x = L y ) water → oil viscous fingering calculations for 2 example cases (Case 1 and Case 2) with ( o ∕ w ) = 1600 , in 3 permeability fields of the same heterogeneity (range of permeabilities) but different correlation lengths, D = 1/10, 1/15 and 1/30. The Case 1 f * w has a shock saturation of S wf = 0.275 and Case 2 has the higher shock saturation of S wf = 0.34 ; for both cases, S wi = 0.17 (where k rw = 0 ). Case 2 also has RPs (see Figs. 6,11) where the mobility function T S w (pressure drop) has be accurately maximised (minimised); the Case 1 example is very "close" to this mobility maximum. As expected, many more proto-fingers are observed in the short correlation case ( D = 1∕30 ) in both Cases 1 and 2 and more major fingers persist for longer until later in the displacement sequence. In the long correlation length case ( D = 1/10), fewer proto-fingers develop into fewer major fingers more quickly. The oil recover/watercut versus PV injected results reflect this showing higher recoveries and later breakthrough for the low correlation case and vice versa for the long correlation case (with the middle correlation length case being in the middle). The pressure drop behaviour across the system also agrees well qualitatively with the experimental observations. Experimental observations of the displacement by water of very viscous oils performed by Skauge et al. (2012Skauge et al. ( , 2013, Vik et al. (2018) and Seright et al. (2018) showed the very surprising result that oil recovery was greatly improved by viscosifying the injected water using polymer such that the viscosity of the polymer was, p = 25 w . The direct experimental observation was that the polymer did not stop fingering from occurring, despite the fact that the polymer application greatly increased oil recovery. The actual mechanism turned out to be viscous crossflow which is analysed and explained in Sorbie and Skauge (2019). The results in this paper fully confirm this but the full details are not presented here. Using the Case 1 simulation for the D = 1/10 case (which showed the least efficient oil recovery), a simulation of polymer flooding after 0.1PV of water displacement was performed, making no further assumptions other than the viscosity of the "polymer" solution ( p = 25 w ). The prediction of the simulation was that oil recovery by polymer injection did indeed improve very significantly relative to the water injection by 61.5% at 1PV and 73.8% at 2PV; these are very significant increases which agree very well with the experimentally observed improvements referred to above. All other simulated quantities also agreed well with the experiment, viz. the timing of the incremental oil recovery was very "fast", the fingering was still present, clear evidence for viscous crossflow was observed and the ΔP versus PV behaviour was similar (not presented here).
A further example was presented where the modified fractional flow, f * w , could be specified to give a fixed even higher S wf (0.46), but at various potential frontal mobility ratio values, M S wf . These are precisely the parameters which are required to model both experimental and field-scale systems where viscous fingering occurs. This sensitivity calculation showed that we can go from quite (frontal) stable to more unstable displacements.
Finally, two quite stringent numerical tests were applied to the basic method of simulating immiscible fingering in water → oil displacements, using Case 1 and D = 1/15 to illustrate our findings. The first used a permeability field where the base case range of k min to k max = 10-1000, i.e. k max ∕k min = 100 , was rescaled to values of k max ∕k min = 10 and 3 . This was shown to make quite a small difference to the overall statistic of the fingering which were very similar (not absolutely identical) and the recoveries/watercuts were very close. Secondly, one of the base case models was coarsened from the original grid with N x /N y = 2000/1000 successively to N x /N y = 1000/500 and then N x /N y = 500/250. Again, the fingering pattern was very similar (not identical) and the recoveries and watercut profiles were almost identical. This paper is far from being the last word on the numerical modelling of immiscible viscous fingering, but we hope that it will stimulate some interesting discussion and further promising research approaches.