\DeclareBoldMathCommand\balpha

α \DeclareBoldMathCommand\btauτ \DeclareBoldMathCommand\bgammaγ \DeclareBoldMathCommand\bbetaβ \DeclareBoldMathCommand\bgrad∇ \DeclareBoldMathCommand\bXX

Modelling shear thinning polymer flooding using a dynamic viscosity model

Prabir

{\rm{Daripa}}^{\rm a}

and Rohit

{\rm{Mishra}}^{\rm b}

^aDepartment of Mathematics
^bMechanical Engineering
Texas A&M University
College Station, TX 77843-3368 Author for correspondence (email:[email protected])email:[email protected]

Abstract

Two distinct effects that polymers exhibit are shear thinning and viscoelasticity. The shear thinning effect is important as the polymers used in chemical enhanced oil recovery usually have this property. We propose a novel approach to incorporate this shear thinning effect through an effective dynamic viscosity of the shear thinning polysolution. The procedure of viscosity calculation of the polysolution, although based on a very basic power law model, is based on empiric coefficients which depends on a spatio-temporally evolving variable namely concentration of polymer. Since viscosity calculation is done pointwise, the accuracy of the model is higher than what exists in literature. This method has been integrated with an existing method for a Newtonian physics based model of porous media flows. The solver uses a hybrid numerical method developed by Daripa & Dutta [1, 2, 3]. The above method solves a system of coupled elliptic and transport equations modelling Darcy’s law based polymer flooding process using a discontinuous finite element method and a modified method of characteristics. Simulations show (i) competing effects of shear thinning and mobility ratio; (ii) injection conditions such as injection rate and injected polymer concentration influence the choice of polymers to optimise cumulative oil recovery; (iii) permeability affects the choice of polymer; (iii) dynamically evolving travelling viscosity waves; and (v) shallow mixing regions of small scale viscous fingers in homogeneous porous media. This work shows an effective yet easy approach to make design choices of polymers in any given flooding condition.

Article Highlights

•

A dynamic viscosity framework for modeling enhanced oil recovery by shear thinning polymer flooding is proposed.
•

This framework is applied to an in-house surfactant-polymer flooding code.
•

A new software that solves shear-thinning polymer flooding problem with or without surfactant has been developed.

Keywords Shear thinning polymer flooding, Enhanced oil recovery, Power law fluids, Data driven modelling, Immiscible two-phase flow

1 Introduction

Multiphase multicomponent porous media flows are ubiquitous. They occur in environment, climate, subsurface, biomedical sciences, etc. just to name a few. They also occur in chemical enhanced oil recovery (CEOR) by polymer and/or polymer-surfactant flooding, a topic related to the subject of this paper. In these chemical floods, a polymer-thickened aqueous phase which we will henceforth call polysolution is injected to displace oil. Assuming polysolution as a Newtonian fluid and neglecting capillary pressure, these flooding processes are modeled by a nonlinear system of elliptic equation for pressure and hyperbolic equations for pressure and saturation (see Daripa et al. [4]). An increase in viscosity, due to polymer, of the displacing fluid inhibits fingering instability and increases the saturation behind the front sweeping the oil. In that study, it has been established through nonlinear wave analysis and numerical simulation using a front tracking method [5] that both of these effects of polymer improve oil recovery. This method is not applicable as it is when capillary pressure is taken into account (in particular when surfactant is also used).

In Daripa & Dutta [2], authors extended the polymer model of Daripa et. al [4] to include the effects of capillary pressure and also of another component, namely surfactant. This new model which we call SP-model (‘S’ stands for surfactant and ‘P’ stands for polymer), is a coupled nonlinear system of elliptic and transport equations. This SP model and the numerical method are very briefly reviewed in §2.1 and §3 which is adequate for the development of the remaining part of this paper. Convergence of this numerical method has been proven in Daripa & Dutta [3]. A Matlab code [1] available on Github has been developed implementing this method. Using this code, numerical simulations have been performed to validate this method qualitatively and quantitatively (see Daripa & Dutta [3]), and to evaluate the relative performance of several floods in homogeneous and heterogeneous reservoirs. As we will see below, this model and the numerical method are the building blocks for the data driven model and the method that incorporates the shear thinning effect of polymer. In this paper, we are interested in including non-Newtonian effects of polymer in SP-flooding model, develop a numerical method to numerically solve this new model, and then perform numerical simulations to validate this method. In particular, we perform simulations with polymer-flooding exclusively by setting the concentration of surfactant to zero in the SP model. This provides an understanding of data driven non-Newtonian effect of polymers on viscous fingering and oil recovery properties of polymer floods. In turn, as we will see, this is helpful in the selection of a polymer for maximum oil recovery at any given flooding condition.

Polymers that are commonly used in chemical enhanced oil recovery are usually stiff polymers of shear thinning type. Thus, stability of such a displacement process where the displacing fluid is shear thinning plays an important role. It is known from theoretical studies [6, 7] that shear thinning behavior of non-Newtonian displacing fluids suppresses fingering instability. In these and many other studies, shear thinning fluids are generally modeled as power-law fluids which depend on two indices, shear rate and density of polysolution. These are taken as constants in linear stability studies. However, in practice such stability results are of little use in CEOR since power law indices are functions of the type of polymer and its concentration which change in space and time in the flow domain (see §2 and §3 below). Stability of such a displacement process with time and space dependent parameters in the power law has not been studied to-date. Even if it were, such stability results are likely to be of little value for practical applications to CEOR since values of these parameters evolve in space-time and are not known a priori. For this and many applications related reasons, there is a need to better understand and model the displacement processes with non-Newtonian fluids in general. It appears that simulation is the only viable alternative for accurate prediction of flow features, flow instabilities and other oil recovery performance measures.

Towards this end, it is worth mentioning some other works from a huge body of literature on CEOR. CEOR methods increase oil recovery by altering fluid-fluid and/or fluid-rock interaction in the reservoir by reducing interfacial tension (IFT) and/or altering the viscosity of the injected fluid for mobility control. Another way the injected chemicals increase recovery is by altering the wettability of the rock to increase oil permeability [8][9]. Pope et al. [10] reported that polymer flooding can increase oil recovery up to 12-30 % original oil initially in-place (OOIP). The idea behind polymer flooding is to introduce an intermediate layer of higher viscosity fluid which increases the cumulative oil recovered by reducing mobility ratio and thereby reducing the finger formation of less viscous fluid (water) displacing the high viscous fluid (oil). Although this method has been used widely in industry, it still lacks a fundamental understanding of the role of polymers in CEOR. Sheng et al. [11] in their review paper on the status of polymer flooding listed different challenges associated with polymer flooding. The major challenge is associated with the basic underlying physics of polymer and its behavior at variable rheological conditions such as temperature, salinity, permeability among others. Due to this highly variable nature of polymers it becomes a challenge in predicting oil recovery by simulation or by experiments.

A polymer working in a specific set of conditions might not be viable in a different oil field completely. The complexity around usage of polymer is so profound that there is no general consensus on a very basic injection parameter, namely, the amount of polymer injected. Between 1970s to 1980s the amount of polymer injected is about 100-200 ppm PV (pore volume) in Chinese projects which changed to 500 to 600 ppm PV in early 1990s. This changed again in 2000s to 400-500ppm PV [12]. Niu et al. [13] reports that when the amount of polymer injected is larger than 400ppm PV, incremental oil recovery becomes less sensitive to the amount of polymer injected. Levitt et al. [14] finds a similar trend, namely increasing polymer (HPAM) viscosity from 3 cP to 60 cP does not significantly change recovery from two pore volumes (PV) of tertiary polymer injection. Apart from the above mentioned uncertainty, another reason polymer flooding may not be the ideal flooding technique is the higher injection rate required to inject polymers. This is overcome by yet another advanced CEOR method termed as the Surfactant-Polymer (SP) flooding. SP flooding is a chemical oil recovery method that introduces a layer of surfactant before injecting the polymer which mobilizes the oil by reducing the capillary pressure. However, to better understand the non-Newtonian effects of polymers we are only focusing on shear thinning polymer flooding in this paper.

Prior research conducted in the field of Saffman-Taylor instabilities suggests that instabilities arise when a less viscous fluid displaces a fluid with higher viscosity. It has been seen from theoretical studies [7, 6] that shear thinning behavior of non-Newtonian displacing fluids may minimize the instability effects (finger formation) on the interface separating the displaced and displacing fluids. The polymers that are currently used in field operations are considered inelastic shear thinning fluids. There have been numerical studies analyzing the inelastic behavior of polymers. The study by Durst et al. [15] shows that elastic behavior is observed only above a critical Deborah (De) number. This number is very small when we consider velocity of shear thinning fluid in porous medium. This has been shown experimentally by Marshall & Metzner [16]. Therefore, the focus of efforts in the past has been to analyze the shear thinning behavior of non-Newtonian fluids as the elastic behavior has less significance in this field of application.

After establishing the importance of shear thinning effects of polymer, the next big question is which polymer is the best for oil recovery. While there is no one good answer to this question multiple studies indicate a couple of good choices. The choice of polymer plays a crucial role in CEOR. The benefits of polymer flood is contingent on the extent of polymer retention in the oil field with minimum retention favorable for high recovery [17]. The three major mechanisms identified by Willhite et al. [18] for polymer flow in porous media are polymer adsorption, mechanical entrapment and hydrodynamic retention.

Apart from the performance of polymer in oil recovery, the choice must also be made based on its environmental impacts. Natural polymers such as Xanthane and Schizophyllan are less detrimental to the environment when compared to hydrolyzed polyacrylamide (HPAM), a widely used industrial polymer which can cause environmental problems (increases the difficulty in oil-water separation, degrades naturally to produce toxic acrylamide and endanger local ecosystem). Agi et al. [19] reviewed different natural polymers in reference to its application in enhanced oil recovery. Many previous studies show non-Newtonian behavior of natural polymer. For example, Liu et al. [20] in their rheology study of Cassava starch finds that the polymer exhibits shear thinning behavior i.e. viscosity decreases with increasing strain rate. Most of the other natural polymers show similar behavior. These natural polymers in an eco-friendly way enhances oil recovery with higher sweep efficiency. This sweep efficiency can be improved if these are used as nanofluids [21]. Nanoparticles have advantages of being tolerant to high salinity, high temperature and retention in highly permeable reservoir. However, there are concerns on the cost of full-scale field implementation and toxicity of nanofluids. Agi et al [22] identified Cissus populnea(CP) as an important biopolymer and showed an efficient process to synthesize Cissus populnea nanofluid (CPNF). They showed at same concentrations the viscosity of CPNF is higher than CP which in turn is higher than commonly used natural polymer – Xanthane. While Xanthane shows a decrease in viscosity with increasing temperatures, CP and CPNF show an evident increase in viscosity. All these recent studies show promise on the possibility of improving enhanced oil recovery if the properties of the polymers are understood and applied carefully. Introduction of biopolymer nanoparticles is another emerging field where focus must be given to analyze the fluid properties not just experimentally but through polymer flood simulations.

Extensive numerical and computational studies for tertiary oil recovery is available in the literature. Daripa et al. [4] provids a comprehensive analytical study of instability control in tertiary oil recovery. Afsharpoor et al. [23] uses upper convected Maxwell equation to model strong extensional flow effect to relate the shear stress with the shear rate. Clemens et al. [24] performs a pore-scale evaluation of polymers displacing oil using experiments and CFD simulation. The experiments show shear thinning behavior of the polymer. CFD simulations show that viscosity is lower at pore throats than that in the pore. This shear thinning behavior affects the displacement efficiency of the polymer flooding. However, the power law model with fixed parameters $(n,\epsilon)$ is usually considered in CFD simulations when in relality these should be functions of polymer concentration. This assumption of constant parameters becomes even more detrimental in predicting recovery in real flooding simulations as the concentration of the polymer changes in space and time. Nandwani et al. [25] models surfactant flooding process using ANSYS FLUENT (a commercial CFD code) and shows that ultralow IFT, minimal fingering and low diffusion rate of the surfactant in oil are responsible for higher oil recovery. For Newtonian fluids, Daripa [26] studies a system of three fluids where the driving fluid drives the fluid in the middle which simultaneously drives the third fluid. He shows that there exists a critical viscosity for the middle fluid which leads to the most stable setup i.e. least finger growth. Recently, Manzoor et. al. [27] tested the effects of pressure variations on heavy oil recovery process in a homogeneous glass beads-packed physical model. The study compares the simulation results with the experiments. The main conclusions are that maximum pressure and periodic pressure variations have the potential to enhance heavy-oil recovery. Moreover, they report that the numerical simulations are highly sensitive to uncertainty in permeability, porosity, diffusion coefficient of polymer, number of grid points and heavy-oil density.

Xie et al. [28] perform experiments on a microchip with heterogeneous porous structures where oil is displaced by dispersed polymer. They show dispersion effect even when the polymer particles are smaller than the pore throat. So the plugging effect is not the major mechanism for preferential flow control by dispersed polymers. Simulations [28][29] for this type of flow show that dispersed polymers smartly controls the preferential flow by inducing pressure fluctuations and thereby increasing efficiency. Xie et al. [30] uses Herschel-Bulkley model to find viscosity of the polymer in a two phase (polymer aqueous solution and oil) Lattice-Boltzman simulation. They show that apart from modelling shear rheological behavior of polymer it is also important to model the three phases (polymer, water and oil) individually. Although the Herschel-Bulkley model has been widely used in the literature to model pseudo-plastics, we use the power law model. The major difference between the two is that the power law model has one less empirical parameter. The reason for this choice is the widespread availability of the power-law coefficients for different polymers used in this study. The different rheological models important for polymer flooding simulations are shown in figure 1.

Refer to caption — Figure 1: Rheological models

Surfactant-polymer flooding simulations helps us in understanding the complex phenomenon of multi-component fluids moving across a porous region. The simulations get even more challenging as one of the components (polymer) behaves as a non-Newtonian fluid i.e. shear rate is not a linear function of strain rate. To account for this, we propose and numerically study a data driven approach to incorporate shear thinning effect of polymer (shear thickening case is similar) in a hybrid numerical method developed by Daripa & Dutta [2, 3]. The data driven approach proposed here is based on the recognition that (i) shear thinning property can be modeled by an effective viscosity of polysolution based on a strain rate using a power law model, (ii) values of the power law parameters in the power law model that fit for the given shear rate and concentration of polymer vary and can be estimated from curve fit based on experimental results, (iii) effective viscosity can be used in an otherwise Darcy’s law based Newtonian model of SP floods without changing the fundamental equations yet implementing an accurate physical model which reduces to the Newtonian rheology with a constant viscosity in regions where the polymer concentration is zero, and (iv) shear thinning effect of any displacing fluid is easily incorporated in any existing Darcy’s law and Newtonian fluid mechanics based code. This paper is laid out as follows. Data driven mathematical model of SP flooding, which consists of the governing equations for SP flooding based on Newtonian fluid dynamics and Ostwald de Waele power law model for shear thinning polysolution, is described in §2. Numerical method is described in §3. Numerical results in rectilinear and quarter five-spot geometries are presented in §4. Finally we conclude in §5. The algorithm and the flow chart for the method are given in Appendix A.1. The method for calculation of finger width is given in Apprndix A.2.

2 Non-Newtonian Mathematical model

The non-Newtonian mathematical model is built from combining Newtonian model of porous media flow for multi-phase multi-component porous media flow involving polymer as a component among others with power law model for the shear thinning fluid. This is discussed in this section. Next section discuses the numerical method.

2.1 Newtonian model of SP flooding

The Newtonian model of Surfactant-Polymer (SP) flooding in porous media is used here as a prototype model. In this section, we recall from Daripa & Dutta [2] some relevant facts about this model. It is given by the following coupled nonlinear system of elliptic equation $\eqref{eq:presssure}_{2}$ for global pressure $p$ and transport equations (2), (3) and (4) for saturation $s$ , concentration $c$ of polymer, and concentration $\Gamma$ of surfactant respectively.

{\bf v}=-{\bf K}({\bf x})\lambda(s,c,\Gamma){\bf\nabla}p,\quad-{\bf\nabla}\cdot\left({\bf K}({\bf x})\lambda(s,c,\Gamma){\bf\nabla}p\right)=q_{a}+q_{o},

(1)

\phi\frac{\partial s}{\partial t}+\frac{\partial f_{a}}{\partial s}{{\bf v}}\cdot{\bf\nabla}s+{\bf\nabla}\cdot\left(D\frac{\partial p_{c}}{\partial s}{\bf\nabla}s\right)\\ =g_{s}-\frac{\partial f_{a}}{\partial c}{{\bf v}}\cdot{\bf\nabla}c-\frac{\partial f_{a}}{\partial\Gamma}{{\bf v}}\cdot{\bf\nabla}\Gamma-{\bf\nabla}\cdot\left(D\frac{\partial p_{c}}{\partial\Gamma}{\bf\nabla}\Gamma\right),

(2)

\phi\frac{\partial c}{\partial t}+\left(\frac{f_{a}}{s}{{\bf v}}+\frac{D}{s}\frac{\partial p_{c}}{\partial s}{\bf\nabla}s+\frac{D}{s}\frac{\partial p_{c}}{\partial\Gamma}{\bf\nabla}\Gamma\right)\cdot{\bf\nabla}c=g_{c},

(3)

\phi\frac{\partial\Gamma}{\partial t}+\left(\frac{f_{a}}{s}{{\bf v}}+\frac{D}{s}\frac{\partial p_{c}}{\partial s}{\bf\nabla}s+\frac{D}{s}\frac{\partial p_{c}}{\partial\Gamma}{\bf\nabla}\Gamma\right)\cdot{\bf\nabla}\Gamma=g_{\Gamma}.

(4)

The terms in the above equations are defined in Table-1 and after Table-1.

Term	Description
${\bf K}({\bf x})$	Absolute Permeability Tensor
$\lambda(s,c)$	Total Mobility
$q_{a},q_{o}$	Source/sink terms for global pressure equation
$s,c,\Gamma$	Saturation (volume fraction), Concentration of polymer, Concentration of surfactant
$\phi$	Porosity of the medium
$p_{c}$	Capillary pressure
$k_{rj}$	(j=o) Relative permeability of oil; (j=a) Relative permeability of aqueous phase
$\mu_{j}$	(j=o) Viscosity of oil (constant); (j=a) Viscosity of aqueous phase
$f_{j}$	(j=o) fractional flow function of oil; (j=a) fractional flow function of the aqueous phase
$g_{s},g_{c},g_{\Gamma}$	Source terms for saturation, concentration and surfactant transport

Table 1: Nomenclature

The elliptic equation $\eqref{eq:presssure}_{2}$ solves for a new function called global pressure which is defined in terms of a thermodynamic pressure and capillary phase pressures so that we have a canonical elliptic equation (1) for the global pressure $p$ . The following relation between this global pressure and phase pressures is taken from Daripa & Dutta [2].

$\displaystyle p=\frac{1}{2}(p_{o}+p_{a})$	$\displaystyle+\frac{1}{2}\int^{s}_{s_{c}}(f_{o}-f_{a})(\zeta,c,\Gamma)\frac{dp_{c}}{d\zeta}(\zeta,\Gamma)d\zeta+\frac{1}{2}\int^{\Gamma}_{\Gamma_{c}}(f_{o}-f_{a})(s,c,\xi)\frac{dp_{c}}{d\xi}(s,\xi)d\xi$
	$\displaystyle-\frac{1}{2}\int\eta^{c}(s,c,\Gamma)\left(\frac{\partial c}{\partial x}dx+\frac{\partial c}{\partial y}dy\right)-\frac{1}{2}\int\eta^{s}(s,c,\Gamma)\left(\frac{\partial s}{\partial x}dx+\frac{\partial s}{\partial y}dy\right)$
	$\displaystyle-\frac{1}{2}\int\eta^{\Gamma}(s,c,\Gamma)\left(\frac{\partial\Gamma}{\partial x}dx+\frac{\partial\Gamma}{\partial y}dy\right)+C,$	(5)

where $C$ is a reference pressure which takes the place of an integration constant in the calculations. In the above,


$\displaystyle\eta^{c}(s,c,\Gamma)$	$\displaystyle=\int^{s}_{s_{c}}\frac{\partial}{\partial c}\left(f_{o}-f_{a}\right)(\zeta,c,\Gamma)\frac{dp_{c}}{d\zeta}(\zeta,\Gamma)d\zeta+\int^{\Gamma}_{\Gamma_{c}}\frac{\partial}{\partial c}\left(f_{o}-f_{a}\right)(s,c,\xi)\frac{dp_{c}}{d\xi}(s,\xi)d\xi,$	(6a)
$\displaystyle\eta^{s}(s,c,\Gamma)$	$\displaystyle=\int^{\Gamma}_{\Gamma_{c}}\frac{\partial}{\partial s}\left(f_{o}-f_{a}\right)(s,c,\xi)\frac{dp_{c}}{d\xi}(s,\xi)d\xi,$	(6b)
$\displaystyle\eta^{\Gamma}(s,c,\Gamma)$	$\displaystyle=\int^{s}_{s_{c}}\frac{\partial}{\partial\Gamma}\left(f_{o}-f_{a}\right)(\zeta,c,\Gamma)\frac{dp_{c}}{d\zeta}(\zeta,\Gamma)d\zeta,$	(6c)

where $s_{c}$ is the value of the aqueous phase saturation at which $p_{c}(s_{c},\Gamma)=0$ and $\Gamma_{c}$ is the surfactant concentration value defined similarly by $p_{c}(s,\Gamma_{c})=0$ .

The aqueous phase viscosity $\mu_{a}$ in this Newtonian model depends linearly on $c$ : $\mu_{a}=\mu_{w}+\beta c$ where $\mu_{w}$ is the constant viscosity of water and $\beta$ is a polymer specific constant. Thus, this system does not include non-Newtonian effect of polymer. The phase relative permeability $k_{rj}$ depends on saturation $s$ and capillary pressure dependent residual saturation $s_{r}$ (see [2]). Since capillary pressure depends on concentration $\Gamma$ of surfactant, the phase relative permeability $k_{rj}$ depends on both $s$ and $\Gamma$ . The total mobility $\lambda=\lambda_{a}+\lambda_{o}$ where $\lambda_{j}=k_{rj}/\mu_{j}$ is the mobility of phase $j$ . Therefore $\lambda_{a}$ depends on the triplet ( $s,c,\Gamma$ ) but $\lambda_{o}$ depends only on ( $s,\Gamma$ ). However, both of these mobility we write generically as $\lambda_{j}(s,c,\Gamma)$ for convenience. The function $f_{j}(s,c,\Gamma)=\lambda_{j}(s,c,\Gamma)/\lambda(s,c,\Gamma)$ is the fractional flow function of the phase $j$ and the coefficient $D(s,c,\Gamma)={\bf K}(\mathbf{x})\lambda_{o}(s,\Gamma)f_{a}(s,c,\Gamma)$ .

The source/sink terms in the elliptic equation $\eqref{eq:presssure}_{2}$ are given by (see Daripa & Dutta [2])

\displaystyle q_{a}=\begin{cases}Q\\ -(\lambda_{a}/\lambda)\,Q\,;\\ 0\end{cases}q_{o}=\begin{cases}0\\ -(\lambda_{o}/\lambda)\,Q\\ 0\end{cases}\text{ at}\,\begin{cases}{\bf x}^{i}=(0,0)\qquad\qquad\qquad\text{ (Source)}\\ {\bf x}^{p}=(1,1)\qquad\qquad\qquad\text{ (Sink)}\\ {\bf x}\in\Omega\setminus\{(0,0)\cup(1,1)\}\,\text{ (Elsewhere)}\end{cases},

(7)

and three source terms in the transport equations are defined as (see Daripa & Dutta [2])

\displaystyle g_{s}=\begin{cases}(1-f_{a})Q\\ 0\end{cases}g_{c}=\begin{cases}(c^{i}-c)Q/s\\ 0\end{cases}g_{\Gamma}=\begin{cases}(\Gamma^{i}-\Gamma)Q/s\\ 0\end{cases}\text{ at }{\bf x}=\begin{cases}{\bf x}^{i}\qquad\quad\text{ (source)}\\ \Omega\setminus\{{\bf x}^{i}\}\text{ (elsewhere)}\end{cases},

where $Q$ is the volumetric injection/production rate. The following initial and boundary conditions are prescribed.


$\displaystyle\quad s({\bf x})=s_{0}({\bf x}),\quad c({\bf x})=c_{0}({\bf x}),\quad\Gamma({\bf x})=\Gamma_{0}({\bf x});\quad$	$\displaystyle{\bf x}\in\Omega,\quad t=0,$	(8a)
$\displaystyle\quad\bgrad s\cdot\hat{{\bf n}}=0,\quad\bgrad c\cdot\hat{{\bf n}}=0,\quad\bgrad\Gamma\cdot\hat{{\bf n}}=0\quad\&\quad{\bf v}_{j}\cdot\hat{{\bf n}}=0;\qquad$	$\displaystyle{\bf x}\in\partial\Omega,\quad t>0,\quad(j=a,o),$	(8b)

where $\hat{{\bf n}}$ is a unit vector normal to $\partial\Omega$ .

This SP-flooding model reduces to polymer flooding model if we set surfactant concentration $\Gamma=0$ in the above equations. This eliminates the transport equation (4) for surfactant, the source term $g_{\Gamma}$ , the terms involving $\Gamma$ from equations (2) and (3), and dependency of $D(s,c,\Gamma)$ , $\lambda(s,c,\Gamma)$ and fractional flow functions $f_{j}(s,c,\Gamma)$ on $\Gamma$ . This completes the description of the Newtonian model of Surfactant-Polymer (SP) flooding. Finer details can be found in Daripa & Dutta [2].

2.2 Model of shear thinning

To include shear thinning effect of polymers in the above model, following Ostwald de Waele power law model (9) is used for the calculation of the effective viscosity $\mu_{a}$ of the aqueous phase.

\mu_{a}=\rho\varepsilon\dot{\gamma}^{n-1},

(9)

where shear rate $\dot{\gamma}$ is given by

\dot{\gamma}=2\sqrt{\mid\Pi_{D}\mid},

(10)

in terms of the second invariant $\Pi_{D}$ of the strain rate tensor given by [31]

\Pi_{D}=-\frac{1}{4}\Bigg{[}\Bigg{(}\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\Bigg{)}\Bigg{]}^{2}+\frac{\partial u}{\partial x}\frac{\partial v}{\partial y}.

(11)

The power law (9) depends on density $\rho$ , strain rate $\dot{\gamma}$ , and two parameters $(\varepsilon,n)$ . The density $\rho(\bf{x},t)$ is the weighted sum of water and polymer densities based on the local in time and space value of concentration, $c(\bf{x},t)$ , of polymer which evolves according to its transport equation (3). The strain rate $\dot{\gamma}$ relates to the second invariant of the strain rate tensor, $\Pi_{D}$ , according to (10) [31].

The strain rate $\dot{\gamma}$ is calculated from the velocity field $(u,v)={\bf v}$ which is obtained as part of the solution of model (1) through (4) described in the next section §3. The parameters $(\varepsilon,n)$ depend on the type and concentration of polymers. This is exemplified in Fig. 2 in which experimental values of $\varepsilon$ and $n$ versus concentration for two shear thinning polymers, namely Xanthane and Schizophyllan, are shown from Hatscher [32] and Lindner et al. [7]. Values of these two parameters in the flow domain at any specific time are found from curve fitting local values of $c$ with experimental data similar to the ones shown in Fig. 2 for the specific polymers in use. Local values of density $\rho$ of polysolution is also found from local values of $c$ , concentration of polymer. Similarly, values of shear rate $\dot{\gamma}$ in the flow domain at any specific time are found from velocity field $(u,v)={\bf v}$ which is obtained as part of the solution of the model (1) through (4) described in §3. All these values are then used in the power law (9) to find the local effective viscosity $\mu_{a}({\bf x},t)$ . This affects the model equations (1)-(4) through the functions $\lambda_{a},\lambda,f_{a}$ and $D$ appearing in these equations which depend on $\mu_{a}({\bf x},t)$ . This procedure is a significant improvement and provides accurate values of interest than conventional practices in which usually power law parameter values are taken at fixed values of polymer concentration at the injection point. In §3, we show this by integrating shear thinning model in the numerical method.

This nonlinear coupling of effective viscosity $\mu_{a}({\bf x},t)$ with the system (1)-(4) and influence of experimental data make it increasingly hard to predict associated response of viscosity variation in the flow domain, flow features, fingering instability and oil recovery performance measures to shear thinning. Therefore, accurate physical interpretation of numerical solution of this data driven model problem is not going to be easy. Regardless, various numerical results obtained using this data driven numerical method discussed in the next section are analyzed and interpreted later in terms of physics.

3 Numerical Method

3.1 A Brief Review of the Numerical Method for the Newtonian Model

We first briefly review from Daripa & Dutta [2] the hybrid numerical method for solving the mathematical model of §2. In this method, global pressure equation (1) is solved using a Discontinuous Finite Element Method (DFEM) [33]. The system of transport equations (2), (3), & (4) is solved by a time-implicit finite difference method based on the Modified Method Of Characteristics (MMOC) [34, 35]. Two different types of grid are used in the numerical method: the finite difference grid shown in Fig. 3(a) for the transport equations and the finite element grid shown in Fig. 3(b) for the pressure equation. Thus, transfer of data between these two grid system during numerical solution process is required in this numerical method because of the coupling of elliptic and transport equations. Details about the computational grid, the numerical method, and the entire computational algorithm for the above Newtonian SP flooding model are given in Daripa & Dutta [daripa2017modelling]. The convergence of the method has been proven in one-dimension in Daripa & Dutta [3].

(a) Uniform FD grid for the transport equations

(b) Uniform FE mesh for the pressure equation

Figure 3: Discretization of the computational domain for (a) the Transport Equations and (b) the Pressure Equation.

3.2 Numerical Method for the non-Newtonian Model

The numerical method discussed above in §3.1 is adapted for this shear thinning model by first calculating the aqueous phase (polysolution) viscosity $\mu_{a}({\bf x},t)$ using data driven power law model (9) as discussed in §2.2. The effective viscosity and appropriate associated quantities which depend on this viscosity are calculated at all grid points of both types of grid. For example, total mobility $\lambda$ which the elliptic equation $\eqref{eq:presssure}$ depends on is calculated at all finite element grid points at every time level since effective viscosity $\mu_{a}$ changes in time which $\lambda$ depends on as explained in §2.2. Similarly, $\lambda_{a},\lambda,f_{a}$ and the function $D$ which the transport equations (2), (3), & (4) depend on are calculated at all finite difference grid points at every time level. All this adds to the computational cost to study the shear thinning effect of polymer. In appendix A.1, the entire algorithm that implements this numerical method and the flow chart are given.

4 Results

Three sets with two simulations per set (for two polymers) are carried out for homogeneous and two different heterogeneous permeability fields. Extensive comparisons have been made to validate this approach. Results are summarised below.

4.1 Shear thinning induced nonuniform viscosity and travelling waves

As discussed in §2.2, variable viscosity profile in a flow domain due to shear thinning effect is very hard to predict a priori due to the influence of experimental data and permeability field of the porous media. It is even harder to predict its development in time a priori as the aqueous phase viscosity $\mu_{a}$ , calculated in space and time from the shear thinning power law model (9), is nonlinearly coupled with the elliptic equation for pressure and the transport equations for saturation and polymer as discussed earlier. However, a posteriori one can analyze such profiles to get an insight into the development of any phenomena and/or distinct patterns discovered in the process. Here we exemplify this by analyzing the evolution of such profiles in a quarter five spot simulation in heterogeneous porous media.

Simulations were carried out in a quarter five spot geometry with heterogeneity shown in Fig. 4. Fig. 5(a) and Fig. 5(b) show plots of viscosity profile at the horizontal mid-section at three different times for the two polymers. We see in these figures several localized travelling waves (peaks) in viscosity at each time level for both the polymers. Development of these patterns appear to be a consequence of shear thinning property of the displacing fluids since we notice in these two figures that these patterns and their speeds are different for these two polymers even though simulations for both the polymers were carried with the same injected polymer concentration (IPC) and same injection rate (IR). These patterns and trends also change as we change the IPC and IR which will be discussed in the future when we have mathematically analyzed this complex model and have a better quantitative understanding of the nonlinear complex dynamics hidden in this model.

4.2 Competing effects: viscosity ratio vs shear thinning

Shear thinning polymers during its passage through the pores undergo shear stresses which lead to a reduction in viscosity. This viscosity reduction is proportional to the local velocity gradients which relates to the injection rates. So, if the polymer were to behave as a Newtonian fluid, high injection rates would mean high oil recovery. But due to the shear effect, a high injection rate may not necessarily correspond to high oil recovery. Therefore, modelling polymer as a non-Newtonian fluid gives a very practical benefit of assessing different injection rates and testing an ideal injection rate which helps in maintaining the viscosity ratio while avoiding the viscosity reduction due to shear thinning. This competing effect can be seen in Fig. 6 where the injection rates were varied and cumulative oil recovery (COR) recorded for time step t=100.

It can be seen in the figure that at lower injection rates Xanthane has slightly higher COR as compared to Schizophyllan at high IPC. But at these low injection rates, Schizophyllan has a higher COR at low IPC. As we increase the injection rates, both the polymers are subject to higher shear rates. Due to this, the advantage they provide by maintaining the viscosity ratio will start to diminish depending on the power law parameters of the polymer. We see that reduction in viscosity is higher for Schizophyllan because of diminishing returns in COR with increasing injection rate. For Xanthane, on the other hand, reduction in viscosity is not at the same rate as for Schizophyllan. Another interesting observation is flipping trend in COR for Schizophyllan. At low IR, higher IPC results in higher COR but at high IR the trend is flipped and lower IPC results in higher COR. This non-linearity is a clear indication of the two competing effects of mobility ratio and viscosity reduction due to shear thinning.

This competitive behavior is also observed in Fig. 7 for simulations in rectilinear geometry with heterogeneity generated using the following permeability tensor ${\bf K}({\bf x})$ taken from [36] and shown in Fig. 10.

\displaystyle{\bf K}({\bf x})=50\left[0.5(1-10^{-7})(\sin{(6\pi\cos{x})}\cos{(4\pi\sin{(3y)})}-1)+1\right]{\mathit{1}},

(12)

where ${\mathit{1}}$ is the second order identity tensor. The shear rate is higher in regions of high permeability. This directly relates to the ease with which the flow moves in these regions. Due to high velocity gradient, the flow experiences higher shear which eventually affects the viscosity of polysolution. Even when the concentration is same, viscosity varies for the four cases due to differences in the shear. Therefore, it can be safely concluded that while high concentration leads to higher viscosity, it is the permeability field that affects the flow movement through shear. Therefore, a shear thinning fluid will perform differently in different permeability fields.

4.3 Simulation in rectilinear geometry

4.3.1 homogeneous porous media

First, simulations in rectilinear geometry are carried out for a homogeneous permeability field with polymers Xanthane and Schizophyllan in separate experiments. Fig. 8 shows the saturation fields obtained with Xanthane and Schizophyllan at four time levels. This shows comparative development of the level sets of saturation due to the effect of two polymers that respond differently to the shear in the field at flooding conditions IR=120,000 and IPC=1500 wppm. It shows that the flow with polymer Xanthane moves faster than polymer Schizophyllan in this case as IR and IPC are set to their respective highest values in the range of values tested. At such low concentration of polymer Xanthane, power law index $n$ is closer to 1 and hence the viscosity has very little shear rate dependency (see (9)). It behaves almost like a non-uniform Newtonian fluid, non-uniform embodies the fact that aqueous phase viscosity will still vary in space and time due to the data and time dependency of $\rho$ and $\varepsilon$ which the power law (9) depends on. On the other hand, for Schizophyllan viscosity depends on shear rate with the index $n$ around 0.7 and affects the viscosity profoundly in comparison to Xanthane. This difference justifies the difference in the speed with which the saturation fronts move in these two cases. Perhaps polymer Xanthane undergoes significant shear thinning in comparison to Schizophyllan and hence moves faster due to reduction in viscosity.

Fig. 8 also shows narrow mixing regions of mild viscous fingers for both the polymers, though the mild fingers in the case of Schizophyllan are somewhat more pronounced. These fingers are mild in comparison to classical ones for displacements involving Newtonian flows because mobility ratio across the interfaces in this figure is position dependent and can change from favorable to unfavorable constantly along an interface because of data dependent shear thinning property of the fluids as discussed earlier. Finger widths in these cases have been computed using the procedure described in appendix A.2. Fig. 9(a) shows finger width growth for polymers Xanthane and Schizophyllan at IR=10,000 and IPC=300 wppm. Fig. 9(b) shows the effect of finger width growth on cumulative oil recovery (COR). It is evident that Schizophyllan with higher finger width shows higher recovery as compared to Xanthane. However, this trend is reversed for higher values of IPC and IR. Plots in Fig. 9(c),(d) show Xanthane to have higher finger width and COR at IR=120,000 and IPC=1500 wppm. This clearly shows the importance of choosing a polymer based on the IR and IPC conditions in a given polymer flood.

4.3.2 heterogeneous porous media

Next, heterogeneous permeability field shown in Fig. 10 is used for simulation in core flood. The saturation fields and COR are shown in Fig. 11 and Fig. 12 respectively. Although most of the flooding is now affected by this heterogeneous permeability field, but for the same IR and IPC clear differences can be noticed on the predicted COR for Schizophyllan and Xanthane. This is a clear indicator that choice of polymer also depends on the permeability field. For the given case, Xanthane seems to outperform Schizophyllan in terms of COR. But given a different set of IR, IPC and permeability field, the choice of polymer may change. In practice, therefore permeability fields also influence the choice of polymer and other outcomes from shear thinning polymer flooding. The fingering nature of the flow pattern seen in Fig. 11 is due to channeling effect as there are regions of high permeability.

4.4 Quarter five spot simulation

Finally, a quarter five spot geometry with a heterogeneous logarithmic permeability field shown in Fig. 4 is used for simulation. Fig. 13 and Fig. 14 show the saturation field and COR for two different polymers. Here again the influence of the permeability field is evident. Also the comparative simulations show significant difference as the variable permeability field leads to higher shear and subsequent variation for different polymers used in the flooding. In simulations with heterogeneous permeability fields, mean finger width does not really help in understanding the underlying physics as the flow becomes very complex due to the interaction of heterogeneous permeability, time dependent variable viscosity field and travelling viscosity waves generated by the shear thinning effect which has been alluded to earlier. The only metric that can help understand the flow then becomes the COR shown in Fig. 14 for both the polymers.

4.5 Effect of injection rate (IR) and injected polymer concentration (IPC) on Cumulative Oil Recovered (COR)

The effect of IR and IPC is polymer specific and can be quantified on a common parameter namely COR. Figures 15 and 16 show the COR for varying IR and IPC conditions for homogeneous rectilinear polymer flooding simulation. It can be seen that both the polymers predict COR with a strong dependence on IR. However, IPC seems to be a less important parameter and the degree to which it influences the COR differs for the two polymers. Xanthane has a stronger dependence on IPC for predicting the COR as compared to Schizophyllan. This information enables correct estimation of polymer required for any flooding process and also indicates the importance of IR that will be required for optimum recovery.

4.6 Sensitivity of the power law parameters ( $\varepsilon$ and $\eta$ )

Significance and implications of the way the power law viscosity model (9) is implemented here are exemplified in this section. This model is not an ordinary power law model with fixed parameters, rather values of the parameters $n$ and $\varepsilon$ are based on the type and concentration $c(x,t)$ of polymer. Since the concentration of polymer is evolving in space and time according to its own transport equation and the power law viscosity model coefficients at every point in space and time are derived from experimental data, effectively this means the classical shear thinning fluid, if thought to have constant values of the power law parameters, itself is changing in space and time. In this way we treat the flow in the most physical way. To the best of the authors’ knowledge, to-date this type of space and time dependent viscosity model has not been implemented in a full scale shear thinning polymer flood simulations.

Fig. 17 shows plots of cumulative oil recovery (COR) versus time for constant values of $n$ and $\varepsilon$ (average values over the entire range chosen) and also for experimental data driven values of these two parameters. Interestingly, difference is observed right from the beginning and is very significant at breakthrough. The graph has three distinct regions. The first region from t=0 to t=300. The second one from t=300 to t=1600 and the third one from t=1600 to t=2500. We see in this figure that simulation with experimental data dependent parameters results in significant variation in COR rate over all three regions as compared to a more or less constant COR rate for constant parameter simulation. In the first region, COR rate is lower for the data dependent parameters than that for constant values of parameters indicating that higher effective polymer viscosity makes the flow slow. This is also seen in Fig. 18 where level sets of saturation at four different time levels with constant power law coefficients (top) and with variable power law coefficients (bottom) in rectilinear homogeneous polymer flooding are compared. The parameters chosen are for high polymer concentration relating to high viscosity. In the second region, due to now mobilized flow which corresponds to higher shear rates, the viscosity starts to decrease even if the concentration stays the same leading to higher COR rate. Finally, in the third region the COR rate decreases due to reduction in sweep efficiency as some aqueous phase reaches the production well. In this region, both constant and data dependent parameter simulations show similar trend. The maximum variation at water breakthrough is around 33 $\%$ . This is a clear indication of the importance of data dependent shear thinning polymer flood simulations. It is worth observing in Fig. 18 that saturation in the data dependent parameter case is more diffused than in the constant parameter case. This can be attributed to the viscosity which is not only varying in space but is also dependent on the local concentration of polymer and therefore smearing the so-called interface.

4.7 Sensitivity to ${\rm K}_{max}$

Any given polymer flooding simulation is highly sensitive to the permeability field. In this section, effect of varying the maximum value, ${\rm K}_{max}$ , of permeability has been discussed. Multiple simulations with the two polymers at different injection rates and ${\rm K}_{max}$ values were run. Fig. 19 shows the cumulative oil recovered at water breakthrough for different conditions. A clear trend is not seen which is expected given the highly non-linear and data driven nature of the problem. However there are some trends that are seen across the two polymers. For lower injection rates, increasing the ${\rm K}_{max}$ value is seen to increase the COR. At some particular cutoff this trend flips which is expected to be lower than IR=60000 for Xanthane but higher for Schizophyllan. As seen from the IR=60000 the overall effect of increasing ${\rm K}_{max}$ is different for Xanthane and Schizophyllan, however the trend is decreasing for both the polymers at the maximum injection rate (IR=120000). So to generalize, increasing the value of ${\rm K}_{max}$ increases (decreases) COR for lower (higher) injection rates. This can be explained from the relation of shear rate with viscosity for complex fluids. At lower injection rates the flow experiences lower shear rate resulting in higher viscosity of the shear thinning fluid. Effect of this higher viscosity is opposed by the lower resistance to the flow from the medium at higher ${\rm K}_{max}$ value. It appears that this second effect wins over the first one with the net effect of increasing overall COR at water breakthrough.

A similar trend is seen for the homogeneous case. Fig. 18 shows the cumulative oil recovered for different polymers and injection rate across ${\rm K}_{max}$ values. However, in this case for the chosen range of ${\rm K}_{max}$ values, the cutoff value is not seen for Xanthane but is seen to be between IR=60,0000 and IR=120,000 for Schizophyllan.

5 Conclusions

A dynamic empiric coefficient based shear thinning model of polymer flooding has been implemented in an in-house code for modelling multi-component multi-phase fluid flow in porous media. Simulations with this data driven model using this code have been performed for two shear thinning polymers, Xanthane and Schizophyllan, at various values of the parameters of the problem. Numerical results which are qualitatively consistent with physics show the merits of this easy-to-implement inexpensive and fast method. Time dependent viscosity profiles in the flow field show a train of travelling viscosity waves which are not yet fully understood theoretically and opens the door for future research for a better theoretical understanding of the complex PDE model for the shear tinning polymer flooding. The simulations also show (i) competing effects of shear thinning and mobility ratio; (ii) injection conditions such as injection rate and injected polymer concentration influence the choice of polymers to optimise cumulative oil recovery (these effects are relatively much more significant in the case of Schizophyllan compared to Xanthane); (iii) permeability field also affects the choice of polymer as polymers show varying movement for different shear rates that are caused by heterogeneity; and (iv) dynamically evolving travelling wave patterns of viscosity profiles which interact with the underlying flow and the heterogeneity field to generate narrow region of fingers which need to be further probed in the future. The significance of this work is that it shows a simple way, yet accurate enough, to incorporate shear thinning effect of polymer in an otherwise Newtonian model of multiphase multicomponent porous media flow model and provides an effective yet easy approach to make design choices of polymers for CEOR in any given flooding condition. We summarize below some of the specific results.

1.

Schizophyllan polymer viscosity decreases faster with increasing shear rate as compared to Xanthane.
2.

Xanthane performs better in terms of COR for higher IPC and IR conditions while Schizophyllan shows higher COR for lower IPC and IR. This fact is reflected from the MFW which in turn shows that the mean finger width is higher for the polymer producing higher COR for a given condition.
3.

Competing effects of viscosity ratio and shear thinning behavior is seen in Schizophyllan polymer for the range of IR and IPC tested. COR is higher for higher IPC at lower injection rates which reduces to a COR less than the case of lower IPC at higher injection rates. This indicates that lower IPC might be preferable for some polymers at higher injection rates.
4.

Permeability field directly affects the shear rates and therefore the viscosity of the aqueous solution. Polymers might perform differently in different permeability fields when shear thinning behavior is accounted for.
5.

The percentage difference for COR is at $33\%$ for shear thinning model when comparing varying parameters and constant power law parameters. This shows the importance of these models to predict the overall recovery in a given polymer flood simulation. These percentages may vary depending on the IPC, IR and permeability.

Appendix A Appendix

A.1 Algorithm

Here we outline the algorithm for the SP-flood simulation that includes shear thinning effect of polymer. The algorithm for the polymer flood is essentially a special case of the same with zero concentration for the surfactant. The step-by-step algorithm built on the one described in Daripa & Dutta [2] is given below. To accomplish some of the steps below, one needs to refer Daripa & Dutta [2].

1.

Define the Cartesian grid in the domain using equal, uniform grid sizes in both the axes. Generate the finite element mesh.
2.

Generate a heterogeneity field on this grid.
3.

Choose an initial interface separating the injected fluid from the resident fluid.
4.

Set the model parameters: $\mu_{o}$ , $\mu_{w}$ , $s_{ro}^{\sigma 0}$ , $s_{ra}^{\sigma 0}$ .

Initialize the state variables $s$ , $c$ and $\Gamma$ as

\displaystyle s_{0}=\begin{cases}1\hfill&x\in\Omega^{+}\\ s_{0}^{\sigma 0}\hfill&x\in\Omega^{-}\\ \end{cases},\qquad c_{0}=\begin{cases}0.01\hfill&x\in\Omega^{+}\\ 0\hfill&x\in\Omega^{-}\\ \end{cases},\qquad\Gamma_{0}=\begin{cases}0.005\hfill&x\in\Omega^{+}\\ 0\hfill&x\in\Omega^{-}\\ \end{cases}.

6.

Calculate $\sigma(s^{n},\Gamma^{n})$ , $\mu_{a}(c^{n},\nabla{\bf v}^{n-1})$ , $s_{ra}(s^{n},\Gamma^{n})$ , $s_{ro}(s^{n},\Gamma^{n})$ , $\lambda_{a}(s^{n},c^{n},\Gamma^{n})$ , $\lambda_{o}(s^{n},c^{n},\Gamma^{n})$ , $\lambda(s^{n},c^{n},\Gamma^{n})$ using $s^{n}$ , $c^{n}$ , $\Gamma^{n}$ which are values of $s$ , $c$ and $\Gamma$ respectively at the $n^{th}$ time level.
7.

Solve the global pressure equation to get $p^{n}$ and subsequently compute ${\bf v}^{n}$ .
8.

Use ${\bf v}^{n}$ , $s^{n}$ , $c^{n}$ , $\Gamma^{n}$ and the quantities calculated in Step 6, to solve for $s^{n+1}$ , $c^{n+1}$ and $\Gamma^{n+1}$ , thus completing a full time step.
9.

If breakthrough is achieved: then stop; else update $n=n+1$ and repeat from Step 6 above.

In order to reduce computational cost, a few iterations of Steps 6 and 8 are done before updating the pressure in Step 7. The pseudocode (see Algorithm 1) and flow-chart (see Fig. 21) for the procedure are given here.

/* Set up Cartesian grid, FE Mesh, permeability field and model parameters */

Set

i,j=0,\ldots,M;\;{\bf x}_{ij}=\left(\dfrac{i}{M},\dfrac{j}{M}\right);

/* (

M\times M

is the grid size) */

Set

\Sigma=

Initial interface;

\Sigma=\partial\Omega^{+}\cap\partial\Omega^{-}

1Set

{\bf K}({\bf x})

Set

\mu_{o},\,\mu_{w},\,s_{ro}^{\sigma 0},\,s_{ra}^{\sigma 0},\,\tilde{q}=\textit{values from \ref{table:1}}

;

/* Initialization */

Set

t=0;\;\Delta t=\dfrac{1}{N};\;Tstop=N\Delta t

;

N

chosen for accuracy */

2 for i = 0, …, M do

3 for j = 0, …, M do

4 Set

(s,\,c,\,\Gamma)({\bf x}_{ij},0)=\begin{cases}(1,0,0)\hfill&{\bf x}_{ij}\in\Omega^{+}\\ (s_{0}^{\sigma 0},c_{0},\Gamma_{0})\hfill&{\bf x}_{ij}\in\Omega^{-}\\ \end{cases}

;

5 end for

7 end for

/* Computation loop */

8 while $\left(s({\bf x}_{M,M},t)\leq 1-s_{0}^{\sigma 0}\;\;\&\&\;\;t<Tstop\right)$ do

9 Compute

\{\sigma,\,\mu_{a},\,s_{ra},\,s_{ro},\,\lambda_{a},\,\lambda_{o},\,\lambda,\,p_{c}\}\;\text{using}\;(s^{n},c^{n},\Gamma^{n},{\bf v}^{n-1})

;

10 Use data driven power law model to calculate

\mu_{a}

;

11 Update

\lambda_{a},\lambda,f_{a},D(s,c,\Gamma)

using

\mu_{a}

;

13 Solve global pressure equation for

p^{n},{\bf v}^{n}

;

14 Recompute

\{s_{ra},\,s_{ro},\,\lambda_{a},\,\lambda_{o},\,\lambda\}\;\text{using}\;(s^{n},c^{n},\Gamma^{n},{\bf v}^{n})

;

15 Solve transport equations for

s^{n+1}

c^{n+1}

and

\Gamma^{n+1}

;

16 Set

t=t+\Delta t

;

18 end while

Algorithm 1 SP flooding simulation

Figure 21: Flow-chart for SP flooding simulation

A.2 Post processing for mean finger width calculation

This appendix has bearing for section 4.3.1. Following procedure is used to calculate the mean finger width (flowchart shown in Fig. 22) which is based on identifying the mixing layer, fingers and numerical interface shown in Fig. 23.

1.

Identify the mixing layer: The 2D domain is scanned at each time step to locate a reduction in saturation. For each $y$ level, let $x_{0}$ denote the start of mixing layer and $x_{1}$ denote the end of mixing layer.
2.

Find the mean saturation for the mixing layer: The mean saturation is found from the mixing region by taking a mean of saturation at each $y$ level. All the means are then averaged to find a level set based on that cutoff concentration. As the mean concentration goes below this set value, the $x$ value corresponding to it is defined as $x_{\rm MFW}$ (shown in green in Fig. 23).
3.

Find the number of fingers: After establishing the location of the interface, number of fingers across that interface is counted. This is done by assigning binary coefficients to each cell based on whether or not they are higher than the cutoff concentration. A counter keeps track of the number of cells that are present before the cutoff concentration limit is not met at which point the counter is re-initialized to zero to count the next finger. This way we can find the count of fingers in the interface from the number of time this flip happens.
4.

Find mean finger width: The mean finger width is found by adding the counter as mentioned in step-3 and dividing that value by the number of flips (number of fingers).

Author Contributions PD lead in the overall planning and formulation of the study outline. PD also developed the theory, algorithms and implemented the algorithms in an in-house code. RM helped in performing the simulations using advanced computing resources provided by Texas A&M High Performance Research Computing during the semester he took a course on Hydrodynamic Stability with the first author PD.

Funding Financial support from the U.S. National Science Foundation through grant DMS-1522782 and from TAMU internal grant T3 from the office of Vice President of Research to the first author PD is gratefully acknowledged. Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.

Data Availability Statement The datasets generated during the current study are available from the corresponding author on reasonable request.

Software Availability The software generated during the current study are available from the corresponding author on reasonable request. The software will be available on GitHub upon publication of this paper.

Conficts of interest The authors declare no confict of interest.

References

[1] P. Daripa and S. Dutta. DFEM-MMOC based EOR Code in MATLAB. https://github.com/daripa8371/EOR/, 2020.
[2] P. Daripa and S. Dutta. Modeling and simulation of surfactant–polymer flooding using a new hybrid method. Journal of Computational Physics, 335:249–282, 2017.
[3] P. Daripa and S. Dutta. On the convergence analysis of a hybrid numerical method for multicomponent transport in porous media. Applied Numerical Mathematics, 146:199–220, 2019.
[4] P. Daripa, J. Glimm, B. Lindquist, and O. McBryan. Polymer floods: a case study of nonlinear wave analysis and of instability control in tertiary oil recovery. SIAM Journal on Applied Mathematics, 48(2):353–373, 1988.
[5] P. Daripa, J. Glimm, J. Grove, B. Lindquist, and O. McBryan. Reservoir simulation by the method of front tracking. In Proc. of the IFE/SSI seminar on Reservoir Description and Simulation with Emphasis on EOR, pages 1–18, 1986.
[6] D. Bonn and J. Meunier. Viscoelastic free-boundary problems: non-newtonian viscosity vs normal stress effects. Physical Review Letters, 79(14):2662, 1997.
[7] A. Lindner, D. Bonn, and J. Meunier. Viscous fingering in a shear-thinning fluid. Physics of Fluids, 12(2):256–261, 2000.
[8] M. Zhou, J. Bu, Y. Ma, J. Zou, H. Fu, and F. Yang. Synthesis of new sulfobetaine gemini surfactants with hydroxyls and their effects on surface-active properties. Journal of Surfactants and Detergents, 21(6):867–877, 2018.
[9] M. Zhou, L. Xia, Y. He, L. Zhang, X. Qiao, and X. Zhong. Synthesis of new salt-resistant sulfonate gemini surfactants with hydroxyl groups. Journal of Surfactants and Detergents, 18(2):303–308, 2015.
[10] G. A. Pope et al. Recent developments and remaining challenges of enhanced oil recovery. Journal of Petroleum Technology, 63(07):65–68, 2011.
[11] J. J. Sheng, B. Leonhardt, and N. Azri. Status of polymer-flooding technology. Journal of Canadian petroleum technology, 54(02):116–126, 2015.
[12] J. J. Sheng. Modern chemical enhanced oil recovery: theory and practice. Gulf Professional Publishing, 2010.
[13] J. G. Niu, P. Chen, Z.B. Shao, D.M. Wang, G. Sun, and Y. Li. Research and development of polymer enhanced oil recovery. Research and development of enhanced oil recovery in Daqing. Petroleum Industry Press, Beijing, pages 227–325, 2006.
[14] D. Levitt, S. Jouenne, I. Bondino, E. Santanach-Carreras, M. Bourrel, et al. Polymer flooding of heavy oil under adverse mobility conditions. In SPE enhanced oil recovery conference. Society of Petroleum Engineers, 2013.
[15] F. Durst, R. Haas, and W. Interthal. Laminar and turbulent flows of dilute polymer solutions: a physical model. In Progress and Trends in Rheology, pages 218–223. Springer, 1982.
[16] R. J. Marshall and A. B. Metzner. Flow of viscoelastic fluids through porous media. Industrial & Engineering Chemistry Fundamentals, 6(3):393–400, 1967.
[17] K. S. Sorbie. Introduction to polymer flooding. In Polymer-Improved Oil Recovery, pages 1–5. Springer, 1991.
[18] G. P. Willhite and J. G. Dominguez. Mechanisms of polymer retention in porous media. In Improved oil recovery by surfactant and polymer flooding, pages 511–554. Elsevier, 1977.
[19] A. Agi, R. Junin, J. Gbonhinbor, and M. Onyekonwu. Natural polymer flow behaviour in porous media for enhanced oil recovery applications: a review. Journal of Petroleum Exploration and Production Technology, 8(4):1349–1362, 2018.
[20] J. Liu, L. Guo, L. Yang, Z. Liu, and C. He. Study on the rheological property of cassava starch adhesives. Advance Journal of Food Science and Technology, 6(3):374–377, 2014.
[21] A. Agi, R. Junin, and A. Gbadamosi. Mechanism governing nanoparticle flow behaviour in porous media: insight for enhanced oil recovery applications. International Nano Letters, 8(2):49–77, 2018.
[22] A. Agi, R. Junin, M. O. Abdullah, M. Z. Jaafar, A. Arsad, W. R. W. Sulaiman, M. Norddin, M. Abdurrahman, A. Abbas, A. Gbadamosi, et al. Application of polymeric nanofluid in enhancing oil recovery at reservoir condition. Journal of Petroleum Science and Engineering, 194:107476, 2020.
[23] A. Afsharpoor, K. Ma, K. Mateen, A. Duboin, S. Jouenne, and P. Cordelier. Micro-scale experiment and cfd modeling of viscoelastic polymer; trapped oil displacement and deformation at the dead-end. In SPE improved oil recovery symposium. OnePetro, 2014.
[24] T. Clemens, K. Tsikouris, M. Buchgraber, L. Castanier, and A. Kovscek. Pore-scale evaluation of polymers displacing viscous oil computational fluid dynamics simulation of micromodel experiments. SPE Reservoir Evaluation & Engineering, 16(02):144–154, 2013.
[25] S. K. Nandwani, M. Chakraborty, and S. Gupta. Chemical flooding with ionic liquid and nonionic surfactant mixture in artificially prepared carbonate cores: a diffusion controlled cfd simulation. Journal of Petroleum Science and Engineering, 173:835–843, 2019.
[26] P. Daripa. Studies on stability in three-layer hele-shaw flows. Physics of Fluids, 20(11):112101, 2008.
[27] A. A. Manzoor. Modeling and simulation of polymer flooding with time-varying injection pressure. ACS omega, 5(10):5258–5269, 2020.
[28] C. Xie, W. Lei, M. T. Balhoff, M. Wang, and S. Chen. Self-adaptive preferential flow control using displacing fluid with dispersed polymers in heterogeneous porous media. Journal of Fluid Mechanics, 906:A10, 2021.
[29] C. Xie, W. Lei, and M. Wang. Lattice boltzmann model for three-phase viscoelastic fluid flow. Physical Review E, 97(2):023312, 2018.
[30] C. Xie, J. Zhang, V. Bertola, and M. Wang. Lattice boltzmann modeling for multiphase viscoplastic fluid flow. Journal of Non-Newtonian Fluid Mechanics, 234:118–128, 2016.
[31] M. T. Schobeiri. Applied fluid mechanics for engineers. McGraw-Hill Education, 2014.
[32] S. Hatscher. Schizophyllan as a biopolymer for eor lab and field results. winterfall. International Energy Agency (IEA), Germany, 2016.
[33] S. Hou, W. Wang, and L. Wang. Numerical method for solving matrix coefficient elliptic equation with sharp-edged interfaces. J. Comput. Phys., 229:7162–7179, 2010.
[34] J. Douglas Jr and T. F. Russell. Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures. SIAM J. Numer. Anal., 19(5):871–885, 1982.
[35] J. Douglas Jr. Finite difference methods for two-phase incompressible flow in porous media. SIAM J. Numer. Anal., 20(4):681–696, 1983.
[36] J. A. Ferreira, G. Pena, and G. Romanazzi. Anomalous diffusion in porous media. Appl. Math. Model., 40:1850–1862, 2015.


$\mathbf{(a)}$	$\mathbf{(b)}$

$\mathbf{(c)}$	$\mathbf{(d)}$


$\mathbf{(a)}$	$\mathbf{(b)}$

$\mathbf{(c)}$	$\mathbf{(d)}$