Injection of water slug containing two polymers in porous media: Analytical solution for two-phase flow accounting for adsorption effects

Journal Pre-proof Injection of water slug containing two polymers in porous media: Analytical solution for two-phase flow accounting for adsorption effects Felipe de O. Apolinário, Anália S. de Paula, Adolfo P. Pires PII:

S0920-4105(20)30026-7

DOI:

https://doi.org/10.1016/j.petrol.2020.106927

Reference:

PETROL 106927

To appear in:

Journal of Petroleum Science and Engineering

Received Date: 20 August 2019 Revised Date:

2 January 2020

Accepted Date: 6 January 2020

Please cite this article as: de O. Apolinário, F., de Paula, Aná.S., Pires, A.P., Injection of water slug containing two polymers in porous media: Analytical solution for two-phase flow accounting for adsorption effects, Journal of Petroleum Science and Engineering (2020), doi: https://doi.org/10.1016/ j.petrol.2020.106927. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier B.V.

CRediT author statement Felipe de O. Apolinário: Methodology, Software, Formal Analysis, Investigation, Writing – Original Draft, Writing – Review & Editing Anália S. de Paula: Methodology, Formal Analysis, Investigation Adolfo P. Pires: Methodology, Formal Analysis, Investigation, Writing – Review & Editing, Supervision

Injection of water slug containing two polymers in porous media: analytical solution for two-phase flow accounting for adsorption effects Felipe de O. Apolinário1*, Anália S. de Paula2, Adolfo P. Pires1

5

1

Universidade Estadual do Norte Fluminense Darcy Ribeiro

2

Petróleo Brasileiro S.A.

*

Corresponding Author

Abstract 10

Polymer flooding is the most applied chemical method of enhanced oil recovery (EOR). Usually this process consists of injecting a slug containing dissolved chemicals displaced by water. This technique is modeled by a system of conservation laws with constant initial condition and not constant boundary conditions. In this paper we extend the two solute one-dimensional chromatography problem solution (Rhee et al., 2001)

15

for a two-phase environment. The exact solution for the injection of a slug containing two dissolved polymers driven by water in an oil reservoir is developed. It is considered that both polymers may adsorb in porous media following Langmuir’s adsorption isotherm. The solution is built splitting the original system into a purely chromatographic problem (auxiliary system), and a lifting equation that considers the

20

hydrodynamics properties of the system. Due to the chase water drive, interactions between waves arise along the space-time plane and change the path of the characteristics at the rear of the slug. The presented solution is composed by rarefaction and shock waves, and constant states. Moreover, a complete chromatography cycle develops in the porous media (complete separation of the chemical components).

25

Keywords: Enhanced Oil Recovery, Chemical Enhanced Oil Recovery, Polymer Flooding, Conservation Laws, Hyperbolic Systems of Partial Differential Equations 1

1. Introduction Injection of aqueous polymer solutions is the most applied chemical method of enhanced oil recovery (Sheng et al., 2015). Adding polymers to injection water increases the viscosity of the solution and reduces its mobility, avoiding viscous 5

fingering and early water breakthrough (Kargozarfard et al., 2018). The polymers commonly used in enhanced oil recovery (EOR) do not modify the interfacial tension between aqueous and oleic phase, and therefore do not change the residual oil saturation. In fact, the remarkable effect of polymers solution injection in a reservoir is the increase of oil production, i.e., a larger volume of oil is produced for a particular

10

volume of polymeric solution injected when compared to the same volume of pure water (Sorbie, 1991). This effect results in a higher oil recovery at the end of the production life of reservoirs subjected to polymer injection. The design of a polymer flooding process depends on the adsorption of the chemicals on the rock. The presence of divalent cations in the reservoir, either on the rock or in the

15

connate water, can increase the adsorption of the polymer by the rock, which reduces the viscosity of the polymeric solution and the efficiency of the process (Sorbie, 1991; Taber et al., 1997). The injection of a water pre-flush slug containing low concentration of divalent cations can minimize the adsorption of the polymers in sandstone reservoirs (Maitin & Volz, 1981; Davison & Mentzer, 1982; Algharaib et al., 2014). For carbonate

20

reservoirs, use of salinity resistant biopolymers or neutral pH water pre-flush with low concentration of
 and high concentration of  are the best options to avoid severe polymer adsorption (Ali & Barrufet, 1994; Lee et al., 2019).

Due to injectivity issues, to presence of divalent cations and to chemical components storage, polymer flooding is most used in onshore high permeable sandstone reservoirs 25

saturated with low/medium viscosity oil (Taber et al., 1997; Torrealba & Hoteit, 2019).

2

However, it has been reported successful implementation of this technique in different environments, including offshore reservoirs (Boardman et al., 1982). In Bohai Bay field, in China, the injection of polymers in an offshore reservoir containing heavy oil resulted in an increase of 3% in oil recovery, and a reduction of 5% in water cut (Zhou 5

et al., 2008). In a deep-water field in Angola, despite of the heavy oil, polymer injection resulted in 7% incremental oil production and 10% water cut reduction (Morel et al., 2012). Numerical simulation and analytical models are two of the most important tools to design polymer flooding projects. Numerical simulators can model complex scenarios

10

and reduce the uncertainty of polymer flooding projects design (Alsofi & Blunt, 2010; Alsofi & Blunt, 2014). On the other hand, analytical models can provide effective sensitive analysis on the performance of this technique under more restrictive assumptions (Bedrikovetsky, 1993). The mathematical model that governs the flow of water containing dissolved chemicals

15

is composed by the water volume conservation and mass conservation of each component. These equations result in a hyperbolic system of conservation laws. For a one-dimensional problem with constant initial and boundary conditions, the solution can be obtained by the method of characteristics and is composed by rarefaction and shock waves, and constant states (Bedrikovetsky, 1993).

20

For the case of water containing one polymer, the mathematical problem is composed

by a 2 × 2 hyperbolic system of conservation laws. The solution is obtained by an

extension of Buckley-Leverett (1942) theory, and its structure is formed by a saturation rarefaction wave followed by a concentration shock and the Buckley-Leverett shock (Patton et al., 1971). This structure can be generalized to include other enhanced oil

25

recovery methods (Pope, 1980), such as surfactant or alkali injection, and to different

3

adsorption isotherms (Langmuir, Henry or Freundlich adsorption isotherms) (Johansen & Winther, 1988).

The multicomponent polymer injection is modeled by a system of  + 1 ×  + 1

5

conservation laws, where is the number of polymers dissolved in water. The solution

to the continuous injection problem was presented in Johansen and Winther (1989), and it was developed from the associated chromatography problem (one-phase multicomponent problem), followed by an extension to the two-phase problem. It is important to point out that this solution is not applied for varying boundary conditions (slug injection).

10

Despite the efficiency of the continuous injection of chemicals to increase oil production, due to economic criteria and injectivity issues, usually a finite volume of polymeric solution is injected in the reservoir (polymer slug) followed by pure water (water drive) (Torrealba & Hoteit, 2019). In this case, a discontinuity arises in the boundary condition at the beginning of the water drive. At this point, interaction

15

between waves of different families occurs along the space-time plane. Therefore, in the case of slug injection of chemicals, the solution is not self-similar (Vicente et al., 2014). A theory to solve chemicals slug injection into oil reservoirs problems is presented in Pires et al. (2006). A potential function related to the conservation of water volume replaces the independent variable time in the system of hyperbolic equations.

20

Considering a slug containing chemical components dissolved, the potential function

splits the original  + 1 ×  + 1 system into an auxiliary system of equations and

a lifting equation. The auxiliary system depends on the thermodynamic model and

represents a pure chromatographic process (one-phase problem), whereas the lifting equation depends on the transport properties of the flow and on the solution of the 25

auxiliary system. For a multicomponent polymer slug injection, it can be noted that after

4

applying the splitting technique, the auxiliary system is similar to the one obtained in multicomponent chromatography (Rhee et al. 2001, Borazjani et al., 2015). The splitting technique has been applied to solve several problems involving polymer injection in porous media. Boa & Pires (2006) considered the case of one polymer 5

continuous injection where the polymer may adsorb according to Langmuir isotherm, and the amount adsorbed is affected by the water salt concentration. Silva et al. (2007) presented the solution for the continuous multicomponent polymer injection. Ribeiro & Pires (2008) considered different adsorption isotherms and fractional flow functions to solve the polymer slug injection problem.

10

Vicente et al. (2014) presented the solution to the problem of slug injection containing one polymer that adsorbs following Langmuir adsorption isotherm. The solution presented was built using the splitting technique and compared to the results obtained in numerical simulators. Borazjani et al. (2014) developed the exact solution to the case of polymer slug

15

injection considering salt effects, and the polymer adsorbs following Henry’s isotherm. De Paula & Pires (2015) considered Langmuir’s adsorption isotherm and salt effects in the adsorption curve to model polymer slug injection in porous media. Borazjani et al. (2016) applied the splitting technique to solve the two-phase problem of polymer slug injection with varying salinity. It was considered that the polymer follows

20

a linear adsorption isotherm and that the salt does not adsorb in porous media. The solution includes implicit formulae for saturation, polymer and salt concentrations and front trajectories of the components. Khorsandi et al. (2016) presented the exact solution to the problem of low salinity polymer slug injection considering cation exchange effects that lead to wettability

25

alteration. The results were compared to experimental data and numerical simulations.

5

De Paula et. al. (2019) presented the solution to the problem of slug injection containing

dissolved polymers that may adsorb following Henry´s isotherm. The solution was composed by a water saturation rarefaction, 2 concentration shock waves and the

Buckley-Leverett shock. The solution also included the effects of the interaction 5

between saturation and concentration waves. Due to the separation of the chemicals in the porous media, water banks appeared in the water saturation solution. It has been shown that it is also possible to build a solution to EOR problems considering advective transport, parabolic terms and relaxation non-equilibrium equations applying the splitting technique for the cases where the auxiliary system

10

allows the development of an analytical solution (Borazjani et al., 2015). In this paper we present the solution for the one-dimensional two-phase isothermal twocomponent polymer slug injection followed by water drive problem. It is considered that both polymers may be adsorbed by the porous media following Langmuir’s isotherm. It was also considered that the slug and water drive salinity are the same as

15

the connate water. The presented solution is a two-phase generalization of the onedimensional two-component chromatography solution presented by Rhee et al. (2001). This solution is general and may be applied to any problem with two dissolved chemical components, for instance: one polymer and one surfactant, one polymer and one salt, and so on. The restriction is that the adsorption must follow Langmuir adsorption

20

isotherm. Next section presents the mathematical description of the physical model followed by the detailed procedure of solution using the splitting technique (Pires et al., 2006) and the chromatography theory (Rhee et al., 2001). Then, the solution in auxiliary space is mapped onto the space-time plane. The paper ends with some discussions and

25

conclusions. 6

2. Mathematical Model We consider the injection of a water slug containing two dissolved polymers displaced by pure water. The following assumptions are adopted: 5

10

•

Two-phase, one-dimensional isothermal flow;

•

Incompressible system;

•

Homogeneous porous media;

•

Dispersion, gravity and capillarity are neglected;

•

Polymers are dissolved only in the aqueous phase;

•

Water density is not a function of polymer concentration;

•

Newtonian flow.

According to these hypotheses, the physical model is governed by the volume balance of water and by the mass conservation of each dissolved polymer in the slug. The resulting system of equations is given by:

15

  +   " +# 





, ,  

  , $

=0

+ 

" , , $

=0

    " +#  , $ +  " , , $ = 0    

(1)

where  is the water saturation, %& and % are the concentration of the components in the flowing phase, & and  the amount adsorbed by the rock,  is the water fractional

flow,  is the rock porosity and ' is the total flow velocity. We introduce the

following dimensionless variables in the system of equations (1): 20

( = ( =

)

*+ , 2

(2)

-3 ./ 010 * 4 + ,

(3)

7

where ( is the dimensionless position related to the length of the slug, ( represents

the number of slug volumes injected, Ω6 is the injected slug volume and 7 is the crosssectional area of the reservoir.

So, we can rewrite (1) in dimensionless form:

8:; + 8" 9# 89

5



8<9, ,  8);   , $

8:;  8" 9#     , $



8:;

=0

+ +

8" <9, , $ 8);

8" <9, , $ 8);

=0

(4)

=0

Water saturation is normalized using the relation: =

9); ,:; 9 = 9 > 9=

(5)

where  ? is the reservoir initial water saturation and  @ is the water saturation at the

injection point. 10

In this problem, the adsorption of both polymers is governed by Langmuir’s isotherm: A %B =

CD D

&∑FG CF F

(6)

where %B = %& , % .

When the slug injection begins (( = 0), there is no polymer in the reservoir (%& = ?

% = 0), and the water saturation is irreducible ( = 0). During the injection of the ?

15

slug, the water fractional flow is 1 at the inlet point (( = 0), and the injection

concentration is constant (%& , % ). After the injection of the slug (( = 0, ( = 1), the @

@

water drive begins (%& = % = 0). Thus, the initial and boundary conditions of the

problem are given by:

( , 0 = 0,

( = 0, H%& ( , 0 = %& , 0 < ( < ? % ( , 0 = % , ?

K

*+ ,

(7)

8

0, (  =  @  @  0,  %& ( = M%& , ( = 0, 0,  @ % 0, (  = M% ,  0,

( > 0 0 < ( < 1 ( > 1 0 < ( < 1 ( > 1

(8)

2.1. Splitting between thermodynamics and hydrodynamics

At this point we introduce the following potential function associated to water conservation: 5

NO = , %BN( − N(

(9)

in the system of conservation laws in dimensionless form (Equation 4), leading to: 8

R

9

S − 8) R<9,

8Q <9, ,  8#  ,  8Q

H8#  , 

10

8Q

8

&

+ 8) = 0 +

8

;

8

8);

S=0

 , 

;

(10)

= 0

(11)

This procedure splits the original 3 × 3 system of equations (Equation 4) into a 2 × 2

system (Equation 11), which is a function of thermodynamic properties only and is similar to the chromatography problem (Rhee et al., 2001), and a lifting equation (Equation 10), which is a function of hydrodynamic properties of the system and of the solution of the auxiliary system (Pires et al., 2006).

The initial and boundary conditions (Equations 7 and 8) in ( × O space become: 15

( , 0 = 0,

O = 0, H%& ( , 0 = %& , 0 < ( < % ( , 0 =

?

? % ,

0, O =  @ O = 1  @  0, %& O = M%& , ( = 0, 0,  @ % 0, O = M% ,  0,

K

*+ ,

O>0 01

01

(12)

(13)

9

Defining

as U, %& , %  and −

&

<9, , 

9

<9, , 

as VU, %& , % , the lifting equation

(Equation 11) becomes: 8 WX, ,  8Q

+

8X9, ,  8);

=0

(14)

The initial and boundary conditions (Equations 12 and 13) for the lifting equation 5

(Equation 14) are:

U → +∞ O = 0, Y V → −∞

(15)

U=1 ( = 0, Y V = −1

(16)

After both auxiliary system (Equation 10) and lifting equation (Equation 14) are solved

10

in ( × O plane, the inverse mapping to ( × ( plane is obtained from the following expression:

N( =

1Q

<"9); ,Q, ); ,Q, ); ,Q$

+

9

<"9); ,Q, ); ,Q, ); ,Q$

2.2. Solution of the auxiliary \ × \ problem

N(

(17)

The methodology to build the solution of the auxiliary problem follows the steps presented in Rhee et al. (2001) for the problem of two-component chromatography: 15

•

Apply the hodograph transformation to determine the concentration waves in the

•

Build the solution path in %& × % plane;

•

%& × % plane;

Map the solution from %& × % plane onto ( × O plane.

The derivation of the characteristic and shock waves in the hodograph plane for the 20

system of equations (11) is presented in Appendix A. The detailed theory for the twocomponent chromatography problem can be found in chapters 1 and 2 of Rhee et al. (2001). In the hodograph plane, the slopes of the characteristic curves are given by:

10

Γ : _ = R

1

& S = ` = & a&& −   + b&& −   + 4& & d &

1 



& Γ : e = R1 S = ` =  & a&& −   − b&& −   + 4& & d 1



&



(18) (19)

The characteristic velocities of the rarefaction waves are: f =

5

f =

1Q

1); 1Q

1);

= g hi 

= g h i

(20) (21)

where the parameters  and i are functions of _ and e, which are defined as:  = _ = jk

j&

(22)

i = ie = lk

l&

(23)

From Rankine-Hugoniot conditions, it is possible to determine the shock velocity: 10

m , i  , i   =

Q

);

= g hi  i 

m  ,  , i = ) = g hi  Q

;

(24) (25)

where the superscripts + and – represent the value of  or i before and after the shock,

respectively.

15

In equations 18-21, 24 and 25, the subscripts + and – denote the slow and fast wave family, respectively.

At this point we present a solution for an arbitrary initial and boundary condition. Figure 1 shows the solution path in the concentration plane. The injection condition is

represented by the point o, and the initial condition is represented by the point p. Points 7 and q are intermediary points along the solution path. We indicate a

20

rarefaction wave as “–“, and a shock wave as “→”.

11

Figure 1: Solution path in the concentration plane

During the slug injection (O < 1), the solution is self-similar and is composed by two 5

shock waves: o → 7 → p. The first shock velocity is denoted as m and is a i

characteristic parameter jump. On the other hand, m is the second shock velocity and is an  parameter jump in auxiliary space (Figure 2).

Across the shock o → 7 (i parameter jump), the concentration % jumps to the

initial condition, whereas the concentration %& increases to a value higher than the 10

injection condition (point 7). On the other hand, the component 1 concentration jumps to the initial condition through the shock 7 → p.

For O > 1, the injected fluid does not contain dissolved polymers, and the solution in

the concentration plane is composed by two rarefaction waves: p − q and q − o.

The first rarefaction wave is a Γ type family, whereas the second rarefaction is a Γ

characteristic family in concentration plane. 15

The characteristic curves
 , with slope f in ( × O plane, correspond to the segment

p − q in the hodograph plane; and the characteristic curves
 , with slope f in 12

( × O plane, correspond to the segment q − o in the hodograph plane. The slope

of the rarefaction waves, which is a function of polymers concentration, is given by (Figure 2):

5


 : f = g ries


 : f = g _i @

(26) (27)

where i @ denotes the value of i at injection condition. Along the rarefaction family
 , there is only component 2 in water (segment p − q in Figure 1, where %& =

%& = 0), and the concentration % is inversely proportional to the slope of the ?

10

rarefaction characteristic. Along the family
 , both chemicals are dissolved in water, and their concentrations are inversely proportional to the characteristic slope (Figure 2).

The water drive (used to displace the polymer slug) leads to interactions between the characteristic waves. Waves of different families (Γ and Γ ) are transmitted through

each other and their slopes change. Waves of same families (Γ and Γ or Γ and Γ )

adsorb each other, and a shock wave appears at the intersection point. The theory of 15

interaction between rarefaction and shock waves is detailed in Rhee et al. (2001).

Figure 2 shows the characteristic diagram in ( × O plane. The shock waves x and y are straight lines whose slope is given by:

m x" = @ , i  = 1, i  = i @ $ = g h@ i @ m y R  = k ,  =  @ , i = 1S = g h@

(29)

z = C

(30)

&

20

(28)

The coordinates of point x are: &

> > >  k# { &# 

Oz = &#> &

(31)

At point x the interaction between the shock m and rarefaction
 begins (Figure 2). At

this point the shock wave is transmitted through the rarefaction wave, and the 13

rarefaction wave through the shock wave. Each rarefaction characteristic
 carries a

constant  value, whereas i jumps from i @ to i ? through the shock wave. Therefore,

to calculate the shock path along the interaction, it is necessary to solve the Rankine-

5

Hugoniot conditions considering the value of  carried by each rarefaction wave. So, the new shock path is given by: 

(

 = z

O

z|  = }~g hi @ ( − ~

z|

"&#> $ &#



and O z| can be calculated by:

10

z|

(32)

&

#>



− 1 + 1

(33)

This shock path generated by its interaction with the rarefaction
 is part of a parabola

where the shock velocity decreases. At point
the interaction ends and the shock path

of m is a straight line with slope g i @ . Ahead of the shock m the component 2 no longer exists, i.e., % = 0. The coordinates of point
are (Figure 2):

| = C

k "&#> $

 k#

15

> { > k&

(34)

O| = #>&k

(35)

(

(36)

"&#> $


€ is: and the shock path  |‚

=C

QQƒ {

>

As the rarefaction
 intercepts shock m , the slope of the characteristics changes due to the jump of i. Thus, the new slope of rarefaction
 is given by:

20

 f  = g h

(37)

After the interaction with shock m , the rarefaction
 intercepts shock m . As both

waves belong to the same family, the shock adsorbs the rarefaction, but its slope changes. The first interaction between the shock m and the rarefaction
 occurs at

point y, whose coordinates are (Figure 2):

14

„ =



Q… C k"#> $ )… C #> "&#> $

(38)

O„ = g h"@ $ „ − z  

The shock path of m along the interaction with
 is:

„† (

5

= „

 ‡

R #> S  ‡

R #S

O „† = g (





(39)

(40)

„†

(41)

and when ( → +∞ the shock slope tends to the last rarefaction characteristic slope.

After the interaction with rarefaction family
 , shock m catches up rarefaction
 , and the shock adsorbs the rarefaction and its trajectory changes continuously (shock and rarefaction belong to the same family). The interaction begins at point € (Figure 2):

10

‚ =

"&Q; C {> )ƒ $ C { > "&{> $

O‚ = g i @ ‚ − |  + O|

After point € the shock path of m is given by:

‚† (

15

= ‚

"&{ > $ &{

O ‚† = g i  (



‚†

(42) (43)

(44) +1

(45)

Analogously to the case of the interaction between
 and m , when ( → +∞ the

shock slope tends to the last rarefaction characteristic slope.

15

Figure 2: Solution of the auxiliary system in plane ( × O

The solution of the auxiliary system is divided in 6 parts (%? , %?? , %??? , %?ˆ , %ˆ , and %ˆ? ),

5

which are bounded by the end of the polymer slug injection (O = 1), and by the first and the last points of the regions where waves interaction take place (points x, y,
,

and €). Thus, the overall solution is: %( , O

%? , % ,  ?? %??? , %?ˆ , %ˆ , %ˆ? ,

OJ1 1 J O J Oz Oz J O J O„ O„ J O J O| O| J O J O‚ O‚ J O

%&  %&  %  &

%& , %

(46)

where, %? ( , O is the self-similar part of the solution, given by: %? ( , O

10

@

%& , % ‹

%& , % ?

% , @

% , ?

% , ?

( J ˆ Q

ˆ‰ Šz Q ˆŒ Š„

Q

‰ Šz

J ( J

J (

Q

ˆŒ Š„

(47)

The other parts of the solution are:

16

%& %  & % &  %?? ( , O = %&  %&  %&  %& %??? ( , O = %& %  &  %&  %&  %&  %  &  %&



= =

‹ %& , % ? %& , %

=

=

 % ,

  = %&Ž ( , O, % = % ,

= %& , % = % , ?

?

= %& , % = % , ?

?

%& = 0, % = 0, % = 0, % = %  , O,  Ž (  &   %& = 0, % = % ,

Q&

Œ #= ,{ >  Q&

= %& , % = % , @

Œ #> ,{ >  Q

? % ,

ˆ‰ Šz Q

? % ,

   = %&Ž ( , O, % = %Ž ( , O, ‹

‰

= %&Ž ( , O, % = %Ž ( , O, @

ˆŒ Š„

( < 

‰

Q&

‰ #= ,{ >  Q& Œ #= ,{ > 

(

z|

< ( <

< ( <

Œ #, ,{ = {Œ  Q < ( ˆ Š„ Œ

( <  Q&

‰ #= ,{ =  Q&

‰ #= ,{ >  Q&

Œ #= ,{ >  Q&

Œ #> ,{>  Q

(48)

ˆ‰ Šz

Œ Š„

‰ #

= ,{ >  Q&

Œ #= ,{>  z| ( ( , O

< ( <

Q&

= = ‰ # ,{ 

Q

Q&

‰ #= ,{>  Q&

Q&

( , O < ( <

QQ …ƒ "#); ,Q$

<

< ( <

< (

< ( <  < ( <

<

< ( < ˆ

Q&

= = ‰ # ,{ 

#= ,{ =  Q&

Q&

‰ #= ,{ =  Q& < ( ‰ #= ,{ =  Q& < (  #= ,{ > 

= 0, % = % ,

= 0, % = %Ž ( , O, = 0, % =

( <

= 0, % = %Ž ( , O,

= 0, % = 0,

%?ˆ ( , O =

5

= 0, % = 0,

< ( < 

QQ…ƒ "#)

; ,Q$

(49)

Œ #, ,{ = { Œ  Q

ˆŒ Š„

Q&

= > ‰ # ,{  Q&

< ( < 

Π#

= ,{>  z| ( ( , O

% = %   , O, % = %   , O, < ( < &  &Ž ( Ž ( Œ #= ,{ >   z| „† % = %   , O, % = % ? , ( ( , O < ( < ( " ( , O$ &  &Ž (   „† ? ? ( " ( , O$ < ( %& = %& , % = % , (50)

17

( < %& = 0, % = 0, ‰ #= ,{ =  Q& Q& % = 0, % = %  , O, &  Ž ( = ,{ =  < ( <  #= ,{ >   #  ‰ ‰ Q& |‚ % = 0, % = %  , &  = > < ( < (  Q&

%ˆ ( , O =

%& = 0, % = 0,   %& = %&Ž ( , O, % = 0,  ? ? %& = %& , % = % ,

%& %  &  %ˆ? ( , O = %&  %&   %&

= 0, % = 0,

  = 0, % = %Ž ( , O,

= 0, % = 0,

 = %&Ž ( , O, % = 0,

=

? %& , %

=

? % ,

‰ # ,{  |‚ ( < ( QQƒ

"#= ,{ = $

<

QQƒ

Π"#

< ( < (

Œ „† ( " ( , O$

( <

= ,{ = $

Q&

„†

< (

(51)

" ( , O$

‰ #= ,{ =  Q& ‚† < ( < ( ‰ #= ,{=  QQ ‚† ( < ( <  "#=,{ƒ=$ (52) Œ QQƒ „† < ( < ( "  ( , O$  "#= ,{ = $ Œ

(

„†

" ( , O$ < (

In equations 47-52, %‘ ( , O denotes the rarefaction wave where component ’

5

concentration changes, subscripts “ + and “ − denote the rarefaction waves
 and
 ,

respectively, and the superscript in parenthesis denotes a constant state in the phase

plane. Figures 3-8 show the concentration profiles for each part of the solution (Equation 46) along the ( × O plane.

18

Figure 3: Concentration profile for solution %? ( , O

Figure 4: Concentration profile for solution %?? ( , O

19

Figure 5: Concentration profile for solution %??? ( , O

Figure 6: Concentration profile for solution %?ˆ ( , O

20

Figure 7: Concentration profile for solution %ˆ ( , O

Figure 8: Concentration profile for solution %ˆ? ( , O

21

In figures 3-8 it is possible to see the development of the chromatographic cycle. Note that component 1 front travels ahead of component 2 due to its lower adsorption rate. In Figure 3 there is a region where the concentration of both components is the injection @

@

condition concentration (%& , % ). This region is followed by a shock wave where the 5

concentration of both components change. After the shock wave, the concentration of component 2 jumps to its initial condition, whereas the concentration of component 1 increases to a higher value than the injection condition (see Figure 1). In figure 4, the water drive has begun, and a rarefaction wave appears at the rear of the polymer slug. Note that at this part of the solution the chromatographic cycle begins.

10

Figures 5 and 6 stress, besides the separation of the components, the development of a sharp edge on the concentration profile of both components. In figure 7, components 1 and 2 are completely separated and travel in porous media as two single component slugs separated by a pure water region. In figure 8, it is possible to see the spreading of the concentration at the rear part of each component slug. This

15

effect is dependent on the adsorption isotherm.

When O → +∞, the rarefaction waves are completely absorbed by the shock waves,

thus components 1 and 2 concentration along the porous media is equal to their water drive concentration (Rhee et al., 2001). 2.3. Solution of the lifting equation

20

In this subsection we present the solution of the lifting equation, which depends on the solution of the auxiliary system (previous section) and the transport properties (relative permeability and viscosity of each phase), using the method of characteristics. Applying the chain rule in equation 10 we find: 8W 8X

8X 8Q

25

+ 8) = − 8 8X

;

8W 8  8Q

− 8

8W 8 

8Q



In constant concentration regions of the solution, we have

(53) 8 8Q

=

8 8Q

= 0. Thus,

22

8W 8X

8X 8Q

+

8X

8);

= 0

(54)

In these regions each characteristic curve carries a constant value of U, and its velocity is given by: 1Q

1);

5

=

8WX, ,  8X

(55)

In regions where %& and/or % vary, U is no longer constant along the characteristic. In

this case, U can be found through the following equation: 1X

1);

= − 8

8W 8 

8Q

− 8

8W 8 

8Q



(56)

The viscosity of the polymer solution is calculated by:

10

– 1 ”• %& , %  = ”• + —& %& + — % 

(57)

– where ”• is the viscosity of pure water, and the coefficient —A is an empiric parameter

that represents the effect of the polymer concentration in the solution viscosity. In this work, we assume —& > — .

Figure 9 presents the solution of the lifting equation in the ( × O plane, which is divided in 9 regions:

15

-

Region (I): constant state: U = U ? , %& = % = 0 (( axis);

Region (1): constant state: U = U & , %& = % = 0; Region (2): constant state: U = U  , %& = %&

‹

and % = 0 (Figure 1);

Region (3): U˜ rarefaction, %& = %& and % = % ; @

@

-

Region (4): interaction between U rarefaction and
 , concentrations vary from

-

Region (5): U™ rarefaction, %& = 0 and % = % ;

20

-

injection condition o to the intermediate state q (Figure 1); 

Region (6): interaction between Uš rarefaction and
 , %& = 0 and % varies

from %



to % = 0 (Figure 1); ?

23

-

Region (7): U› rarefaction, %& = % = 0;

Region (4-): U rarefaction, %& varies from %&

Region (5-): constant state: U

U ™ , %&

‹

%

to %& and % ?

0.

%

?

0;

Figure 9: Lifting equation solution in auxiliary plane

5

The solution of U( , O is presented in six parts. When O J 1, the solution is self-

similar. The remaining regions are bounded by O interaction (points x, y,
, and €):

U( , O

10

U? , U ,  ?? U??? , U?ˆ , Uˆ , Uˆ? ,

OJ1 1 J O J Oz Oz J O J O„ O„ J O J O| O| J O J O‚ O‚ J O

1 (part II), and by the waves

(66)

The self-similar part is given by:

24

œ > >

U? ( , O

&

8WRX , , S U @ , ( J O }  8X   & œ > > U  , O, O }8WRX , , S J  J ˜ ( (

U  ,  U & ,  U ? ,

Q

ˆ‰ Q ˆŒ

8X

J ( J

Q

ˆŒ

Q

ˆ‰

(67)

J ( J ∞

( → ∞

The structural formula for this part is p → 1 → 2 → 3 P o (Figure 10).

Figure 10: Solution U? in V U plane

5

The solution U?? ( , O is:

U?? ( , O

&

  U @ , ( J O P 1  Ÿ 8X  & Q& U  , O, O P 1 8W"X ž, –, –$Ÿ J  J ( 8X ‰ #= ,{ =   › ( Q& Q& U  , O, š ( = ,{ =  J ( J  #= ,{>   # ‰ ‰  Q& Q& U™ ( , O, J ( J  #= ,{>  #= ,{ > 

8W"X ž , –, –$

‰

Œ

U ( , O,  #=,{> J ( J  #> ,{> Œ Œ  Q& Q U˜ ( , O, Œ#> ,{> J ( J ˆ‰   Q Q J ( J ˆ U , ˆ‰ Œ Q U & , J ( J ∞ ˆŒ  U ? , ( → ∞ Q&

Q&

(68)

and its structural formula is p → 1 → 2 → 3   P 4   P 5   P 6   P 7   P o. The superscript ′′ denotes the first rarefaction wave in the region.

25

For U??? ( , O we have

U??? ( , O

5

Figure 11: Solution U?? in V × U plane

U @ ,  U  , O,  › ( U  , O, š (  U™ ( , O,

( J O P 1  O P 1  Q&

8X

8W"X ž , –, –$

‰ #= ,{ =  Q&

‰ #= ,{ >  Q&

8X

J ( J

&

Ÿ

Q&

‰ #= ,{ >  Q&

J ( J 

Π#

= ,{ > 

U ( , O,  #=,{> J ( J ( Œ   QQ… z| U ( , O, ( J ( J Œ #> ,{=   QQ… Q J ( J > = U , Œ # ,{  ˆŒ Q U & , J ( J ∞ ˆŒ  U ? , ( → ∞ z|

&

Ÿ

8W"X ž , –, –$

J ( J 

Q&

‰ #

= ,{ = 

(69)

The structural formula is p → 1 → 2 P 4   P 4  → 4   P 5′′ P 6   P 7   P o

(Figure 12). The superscript ′ denotes the last rarefaction wave in the region. We also

add superscripts “ ” and “P“ in the rarefaction notation (either ¥′ or ¥′′) to indicate that the rarefaction is a left or a right state of a shock wave, respectively. In figure 12,

¥ ¦ ( , O denotes the concentration rarefaction waves in region ¥ connected through a

10

shock. This structure appears in concentration profiles that cross the shock waves O z| ,

O „† and O ‚† (regions where concentration rarefactions waves interact with shock waves).

26

Figure 12: Solution U??? in V × U plane

The solution U?ˆ ( , O is given by

5

U?ˆ ( , O

U @ ,  U  , O,  › ( U  , O,  š ( U™ ( , O,  U ( , O,  U ( , O,  & U , U ? ,

( J O P 1  O P 1  Q&

8W"X ž , –, –$ 8X

8W"X ž , –, –$

‰ #= ,{ =  Q&

‰ #= ,{ >  Q&

8X

J ( J 

Q&

‰ #

J ( J 

&

Ÿ

= ,{ > 

&

Ÿ

J ( J 

Q&

Π#

Q&

‰ #

= ,{> 

z| J ( J ( Œ #= ,{ >  z| „† ( J ( J ( „† ( J ( J ∞

= ,{ = 

(70)

( → ∞

and its path in V U plane is shown in figure 13. The structural formula for this part of

the solution is p → 1 → 4   P 4  → 4   P 5   P 6′′ P 7   P o.

27

Figure 13: Solution U?ˆ in V × U plane

For the region of solution Uˆ ( , O, we have:

5

Uˆ ( , O

U @ ,  U  , O, › (  Uš ( , O,  U™ ( , O,  ™ , U   U ( , O,  & U , U ? ,

( J O P 1 

8W"X ž , –, –$

O P 1  Q&

8X

8W"X ž , –, –$

‰ #= ,{ =  Q&

8X

J ( J 

J ( J (

‰ #= ,{ >  |‚ ( J (

QQƒ Œ = = Œ # ,{ 

J

J

Q&

‰ #

= ,{> 

|‚

QQƒ

&

Ÿ

&

Ÿ

J ( J 

Q&

‰ #

= ,{ = 

(71)

Œ

Œ #= ,{ =  „† ( J (

( J ( J ∞ ( → ∞ „†

Figure 14 presents the solution path of Uˆ in V U plane, whose structural formula is p → 1 → 4   P 5  → 5   P 6   P 7′′ P o.

28

Figure 14: Solution Uˆ ( , O in V U plane

Finally, the solution Uˆ? ( , O is given by

5

Uˆ? ( , O

U @ ,   U› ( , O,  Uš ( , O,

( J O P 1  O P 1  Q&

8W"X ž , –, –$ 8X

8W"X ž , –, –$ 8X

J ( J (

‰ #= ,{ =  ‚† ( J (

‚†

J Œ =ƒ = U ™ , Œ # ,{   QQ „† U ( , O, Œ =ƒ J ( J ( = Œ # ,{    & „† ( J ( J ∞ U , U ? , ( → ∞ QQ

&

Ÿ

&

Ÿ

J ( J 

Q&

‰ #

= ,{= 

(72)

Figure 15 shows the solution path of Uˆ? in V U plane. Its structural formula is p → 1 → 4   P 5  → 6   P 7′′ P o.

29

Figure 15: Solution Uˆ? in V × U plane

5

From the definition of the variables U and V, the water fractional flow function

( , O, %& , %  and saturation ( , O can be easily determined through the following expressions: 

&

X); ,Q

( , O

(73) P

WX); ,Q, ,  X); ,Q

(74)

2.4. Inverse mapping to time domain

10

At this point we have already calculated the solution of , %& and % in the ( O plane. The next step is the inverse mapping of these solutions onto the space-time plane through the relation: N(

1Q <"9); ,Q, ); ,Q, ); ,Q$

9 N( <"9); ,Q, ); ,Q, ); ,Q$

(75)

30

Rarefaction waves and constant states are calculated directly from expression (75), whereas the shock trajectories are determined by (Pires et al., 2006): §A =

<¦

9 ¦ ˆD

,¨

, P

(76)

where mA is the respective shock velocity in the auxiliary plane (Equations 24-25).

5

The complete mathematical derivation and the exact equations for the inverse mapping are presented in Appendix B. Figure 16 presents the characteristic diagram of the solution in ( ( plane.

Figure 16: Characteristic diagram in ( ( plane

10

The solution is also divided in six regions. For ( J 1 the solution is self-similar. Beyond this point, the solution is bounded by the start and end of the waves interactions regions (points x, y,
, and €): ( , O

? ,  ,  ?? ??? , ?ˆ , ˆ , ˆ? ,

( J 1 1 J ( J z z J ( J „ „ J ( J | | J ( J ‚ ‚ J (

(77)

31

The self-similar part of the solution (? ) is given by:  ,

? ( , ( 

>

>

   @ , ( J ( 9>  > >  , , S  ˜ ( , ( . ( J ( J ( § 9 >

  ,  & ,   ? ,

,

S

( § J ( J ( § ( § J ( J §zK ( §zK ( J (

(78)

Figure 17 presents the self-similar part of the solution for three different times ((& J

( J (˜ J 1).

5

Figure 17: Self-similar part of the solution

The path of the solution ? in the   plane (phase plane) is shown in Figure 18

(structural formula: o P 3 → 2 → 1 → p). This part of the solution is equivalent to the solution of continuous water injection containing two dissolved 10

polymers.

32

Figure 18: Solution path for ? in phase plane

The solution ?? is

?? ( , ( 

5

= =

,

S



(

   @ , ( P 1 ( J 9 >  > = = ªª = =   ,  ,  9žªª C k#= «{ = ¬  = = = ®   ¬  = ® = ®
™

 ,

(

(

9­ªª C k#= «{ > ¬

(



9¯ªª C k«#= ¬ {>

(

= ® > >   9 œ C k«#> ¬ { >  > >  ¬ { >    , ( § J ( J ( §  & ( § J ( J §zK (  , ?  , §zK ( J (

(79)

The profile of the solution ?? is presented in figure 19 and its path in the phase plane is

shown in figure 20. The structural formula is o P 7   P 6   P 5   P 4   P 3 → 2 →

1 → p. The same notation previously used for the V U solution is followed for

  plane.

33

Figure 19: Solution profile for ?? ( , (  for ( (1 J ( J z )

Figure 20: Solution path for ?? in phase plane

5

At the beginning of water drive, for 1 J ( J z (Figure 19), two concentration rarefaction waves appear at the rear of the slug. These rarefactions interact with the saturation rarefaction waves from region 3 (Figure 16). Note that the separation of the components begins (chromatographic cycle), where component 1 travels ahead of component 2 in porous media. 34

The solution ??? is given by: ??? ( , ( 

 @ ,    ,  ,  › ( (  š ( , ( ,   ™ ( , ( ,

( J

= =

 , , S 9 >

= =

 , , S 9 >

= =

( P 1 J ( J




9žªª C k#= «{ = ¬ = ®


9­ªª C k#= «{

( P 1

9žªª C k#= «{ = ¬

( P 1 J ( J

( P 1 J ( J >  ¬

= =




= ®

( P 1




9­ªª C k#= «{ > ¬ = ®


9¯ªª C k«#= ¬ { >

( P 1

( P 1

 ( P 1 J ( J R ªª = 1:; z| (  ¬ 9¯ C k«# {  =
= ®

(80)

 is the value of  on a rarefaction wave
 after its interaction with shock § . where Ž

5

Figure 21 shows the profile of the solution ??? , and in figure 22 we present its path in

the phase plane. The structural formula for this region is: o P 7   P 6   P 5   P 4   → 4  P 2 → 1 → p.

Figure 21: Solution profile for ??? ( , (  for (™ (z J (™ J „ ) 35

Figure 22: Solution path for ??? in phase plane

In this part of the solution, a lower viscosity water region between two regions of higher viscosity appears (Figure 21). This effect is caused by the separation of the dissolved 5

polymers. In the low viscosity region, a peak in water saturation can be observed in the

solution profile (see beginning of  ( , (  in figure 21).

For the solution ?ˆ we have ?ˆ ( , ( 

 @ ,    ,  ,  › ( (  š ( , ( ,   ™ ( , ( ,

( J

= =

 , , S 9 >

= =

 , , S 9 >

= =

( P 1 J ( J




9žªª C k#= «{ = ¬ = ®


9­ªª C k#= «{

( P 1

9žªª C k#= «{ = ¬

( P 1 J ( J

( P 1 J ( J >  ¬

= =




= ®

( P 1




9­ªª C k#= «{ > ¬ = ®


9¯ªª C k«#= ¬ { >

( P 1

( P 1

 = ®  ( P 1 J ( J R ªª = 1:; z| ( ¬ 9¯ C k«# {  Œ , = S  
(81)

36

Figure 23 presents the saturation profile for the solution ?ˆ and figure 24 shows its path in  ×  plane. For this case the structural expression is o − 7   − 6   − 5   − 4   →

4  P 4   → 1 → p.

5

Figure 23: Solution profile for ?ˆ ( , (  for (š („ J (š J | )

Figure 24: Solution path for ?ˆ in phase plane At this point the polymers are almost separated, and there is only a small region where

both components coexist ( ( , (  in figure 23). The peak of water saturation 10

increases and the constant state (2) no longer appears in the solution. The solution ˆ is

37

ˆ ( , ( 

 @ ,    ,  , › ( (   š ( , ( ,   ™ ( , ( ,

( J

= =

 , , S 9 >

= =

 , , S 9 >

= =

( P 1 J ( J




9žªª C k#= «{ = ¬ = ®


9­ªª C k#= «{

( P 1

9žªª C k#= «{ = ¬

( P 1 J ( J

( P 1 J ( J >  ¬

= =




= ®

( P 1




9­ªª C k#= «{ > ¬ = ®



( P 1

( P |  |

 = ® = =  9 C k«#= ¬ { =  = = Œ , = S
Figure 25 shows the profile for ˆ for (› , where | J (› J ‚ . The path of ˆ in the 5

phase plane is presented in figure 26. Its structural representation is o P 7   P 6   P

5   → 5  P 4   → 1 → p.

Figure 25: Solution profile for ˆ ( , (  for (› | J (› J ‚  38

Figure 26: Solution path for ˆ in  ×  plane

Note that the components are completely separated in porous media (full chromatographic cycle) (Figure 25). The chemical components are traveling in two 5

different slugs separated by a pure water bank (constant state  ™ ).

Finally, ˆ? is given by: ˆ? ( , ( 

 @ ,    ,  ,  › ( (  š ( , ( , 

( J

= =

 , , S 9 >

= =

 , , S 9 >

= =

( P 1 J ( J




9žªª C k#= «{ = ¬


( P 1

= ‰ ,±‰ S

= =


9žªª C k#= «{ = ¬

( P 1 J ( J

=



( P 1

‰ S
( P ‚  ‚


= =

,

S

   ™ ( P ‚  ‚ J ( J ¯Œ ( P |  | ,   ‰ 9­ªª‰ C {±‰ 9 C k«#= ¬ { =  = = Œ , = S
 J  J  ( | | ( Œ 9°ªªŒ C k#±Œ   ( ( 9 ¯Œ C k«#= ¬ { =  Œ , = S
(83)

10

 where iŽ is the value of i on a rarefaction wave
 before its interaction with the

shock § . The saturation profile of the solution ˆ? is shown in figure 27.

39

Figure 27: Solution profile for ˆ? ( , (  for (² ((² L ‚ )

In figure 28 we present the solution path of ˆ? ( , (  in   plane. The structural

representation of the solution in this region is o P 7   P 6   → 5  P 4   → 1 → 5

p.

Figure 28: Solution path for ˆ? in phase plane In this part of the solution, the slugs containing the chemical components travel with different velocities in the porous media, and the high water saturation region (water 40

bank) becomes larger. Moreover, the constant state (5) no longer exists, which implies in no constant concentration regions in porous media.

For ( → +∞, both concentration rarefaction waves are completely adsorbed by the

respective shock waves of the same family. Thus, the chemical components 5

concentrations are %& = % = 0 and the water saturation is  =  @ (boundary condition for ( > 1) (Rhee et al., 2001).

3. Conclusions

In this paper we derive the analytical solution to the one-dimensional two-phase water slug injection containing two dissolved polymers problem. The solution was built using 10

the splitting technique, which consists in applying a potential function that splits the original system of equations into an auxiliary system and a lifting equation. The solution of the problem was obtained in the auxiliary plane and subsequently mapped to the space-time plane. Effects of interactions between waves of different families and between waves of same family were considered in the construction of the solution.

15

The solution is divided in six regions bounded by the crossing points of the waves in the space-time plane. Analytical expressions and profiles for water saturation and polymer concentrations were presented for each region. It was shown the development of a full chromatographic cycle, and how the separation of the chemicals influences the saturation profile. It is highlighted that in regions where

20

a lower viscosity aqueous solution appears between two regions of higher viscosity aqueous solution, water banks without polymers appear. We also showed that the water bank region increases with time due to the different velocities of the polymer slugs. The presented solution can be used to validate reservoir simulators and to evaluate the most important parameters to polymer flooding projects design. It is also possible to

25

forecast fluids production, and water and polymer breakthrough.

41

The solution methodology shown in this paper can be applied to solve similar slug injection problems, such as surfactant-polymer systems or low-salinity polymer systems. The only restriction is that the adsorption of the chemicals must follow Langmuir adsorption isotherm. 5

4. Acknowledgements This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. 5. References Algharaib, M., Alajmi, A., & Gharbi, R. (2014). Improving Polymer Flood Performance

10

in High Salinity Reservoirs. Journal of Petroleum Science and Engineering, 115, pp. 17-23. doi: http://dx.doi.org/10.1016/j.petrol.2014.02.003 Ali, L., & Barrufet, M.A. (1984). Profile Modification Due to Polymer Adsorption in Reservoir

Rocks.

Energy

&

Fuels,

8(6),

pp.

1217-1222.

doi:

http://dx.doi.org/10.1021/ef00048a008 15

Alsofi, A.M., & Blunt, M.J. (2010). Streamline-Based Simulation of Non-Newtonian Polymer Flooding. SPE Journal, 15(04), 901-911. doi: https://doi.org/10.2118/123971PA Alsofi, A.M., & Blunt, M.J. (2014). Polymer Flooding Design and Optimization Under Economic Uncertainty. Journal of Petroleum Science and Engineering, 124, 46-59. doi:

20

https://doi.org/10.1016/j.petrol.2014.10.014 Bedrikovetsky, P. G. (1993). Mathematical Theory of Oil and Gas Recovery. London: Kluwer Academic Publishers. Boa, P. M. F., & Pires, A. P. (2006). Salt Effects on Polymer Adsorption in Chemical Flooding of Oil Reservoirs. In 11th Brazilian Congress of Thermal Sciences and

25

Engineering, Curitiba, Brazil.

42

Boardman, R. S., Moore, L. J., Julian, M. H., Bilbrey, D.G., & Moore, J. S. (1982). Design and Implementation of Four Enhanced Oil Recovery in Bay Fields of South Louisiana. In SPE Enhanced Oil Recovery Symposium. Tulsa, Oklahoma, USA. Borazjani, S., Bedrikovetsky, P. G., & Farajzadeh, R. (2014). Exact Solution for non5

Self-Similar Wave-Interaction Problem During Two-Phase Four-Component Flow in Porous

Media.

Abstract

and

Applied

Analysis,

2014.

doi:

http://dx.doi.org/10.1155/2014/731567 Borazjani, S., Bedrikovetsky, P. G., & Farajzadeh, R. (2016). Analytical Solutions of Oil Displacement by a Polymer Slug with Varying Salinity. Journal of Petroleum 10

Science and Engineering, 140, 28-40. doi: https://doi.org/10.1016/j.petrol.2016.01.001 Borazjani, S., Roberts, A.J., & Bedrikovetsky, P.G. (2016). Splitting in Systems of PDEs for Two-Phase Multicomponent Flow in Porous Media. Applied Mathematics Letters, 53, 25-32. doi: https://doi.org/10.1016/j.aml.2015.09.014 Buckley, S.E. & Leverett, M.C. (1942). Mechanism of Fluid Displacement in Sands.

15

Transactions of the AIME, 146(01), pp. 107-116. doi: https://doi.org/10.2118/942107-G Davison, P & Mentzer, E. (1982). Polymer Flooding in North Sea Reservoirs. SPE Journal, 22(03), pp. 353-362. doi: https://doi.org/10.2118/9300-PA Johansen, T., & Winther, R. (1988). The Solution of the Riemann Problem for a Hyperbolic System of Conservation Laws Modeling Polymer Flooding. SIAM Journal

20

of Mathematical Analysis, 19(3), 541-566. doi: https://doi.org/10.1137/0519039 Johansen, T., & Winther, R. (1989). The Riemann Problem for Multicomponent Polymer Flooding. SIAM Journal of Mathematical Analysis, 20(4), 908-929. doi: https://doi.org/10.1137/0520061

43

Kargozarfard, Z., Ruai, M., & Ayatollahi, S. (2018). Viscous Fingering and Its Effect on Areal Sweep Efficiency During Waterflooding: an Experimental Study. Petroleum Science, 16(7), 105-116. doi: https://doi.org/10.1007/s12182-018-0258-6 Khorsandi, S., Changhe, Q., & Johns, R. T. (2016). Displacement Efficiency for Low 5

Salinity Polymer Flooding Including Wettability Alteration. In SPE Improved Oil Recovery Conference, Tulsa, USA. Lee, Y., Lee, W., Jang, Y., & Sung, W. (2019). Oil Recovery by Low-Salinity Polymer Flooding in Carbonate Oil Reservoirs. Journal of Petroleum Science and Engineering, 181, 1-9. doi: https://doi.org/10.1016/j.petrol.2019.106211

10

Mandal, A. (2015). Chemical Flood Enhanced Oil Recovery: A Review. International Journal

of

Oil,

Gas

and

Coal

Technology,

9(3),

241-264.

doi:

http://doi.org/10.1504/IJOGCT.2015.069001 Maitin, B.K. & Volz, H. (1981). Performance of Deutsche Texaco Ag’s Oerrel and Hankensbuettel Polymer Floods. In SPE/DOE Enhanced Oil Recovery Symposium, 15

Tulsa, USA. doi: http://doi.org/10.2118/9794-MS Morel, D., Vert, M., Jouenne, S., Gauchet, R., & Bouger, Y. (2012). First Polymer Injection in Deep Offshore Field Angola: Recent Advances in the Dalia/Carmelia Field Case. SPE Oil and Gas Facilities, 1(2), 43-52. doi: https://doi.org/10.2118/135735-PA Patton, J. T., Coats, K. H., & Colegrove, G. T. (1971). Prediction of Polymer Flood

20

Performance. SPE Journal, 11(1), 72-84. doi: https://doi.org/10.2118/2546-PA de Paula, A. S, & Pires, A. P. (2015). Analytical Solution for Oil Displacement by Polymer Slugs Containing Salt in Porous Media. Journal of Petroleum Science and Engineering, 135, 323-335. doi: https://doi.org/10.1016/j.petrol.2015.09.001 Pires, A. P., Bedrikovetsky, P. G., & Shapiro, A. A., 2006. A Splitting Technique for

25

Analytical Modeling of Two-Phase Multicomponent Flow in Porous Media. Journal of 44

Petroleum

Science

and

Engineering,

51(1-2),

54-67.

doi:

https://doi.org/10.1016/j.petrol.2005.11.009 Pope, G. A. (1980). The Application of Fractional Flow Theory to Enhanced Oil Recovery. SPE Journal, 20(3), 191-205. doi: https://doi.org/10.2118/7660-PA 5

Rhee, H. -K, Aris, R., & Amundson, N. R. (2001). First-Order Partial Differential Equations: Theory and Application of Hyperbolic Systems of Quasilinear Equations, vol. 2. New Jersey, USA: Prentience-Hall, Englewood Cliffs. Ribeiro, P. M., & Pires, A. P. (2008). The Displacement of Oil by Polymer Slugs Considering Adsorption Effects. In SPE Annual Technical Conference and Exhibition,

10

Denver, USA. Sheng, J. J., Leonhardt, B., & Azri, N. (2015). Status of Polymer-Flooding Technology. Journal

of

Canadian

Petroleum

Technology,

54(2),

116-126.

doi:

https://doi.org/10.2118/174541-PA Silva, R. C. A., Cardoso, C. B., & Pires, A. P., 2007. The Role of Adsorption Isotherms 15

on Chemical Flooding Oil Recovery. In SPE Annual Technical Conference and Exhibition, Anaheim, USA. Sorbie, K. S. (1991). Polymer-Improved Oil Recovery. New York, USA: Springer Science and Business Media. de Paula, A.S., Apolinário, F.O., & Pires, A.P. (2019). Water Slug Injection Containing

20

n

Polymers

in

Porous

Media.

AIChE

Journal,

65(11).

doi:

https://doi.org/10.1002/aic.16735 Taber, J. J., Martin, F. D, & Seright, R. S. (1997). EOR Screening Criteria Revisited – Part 1: Introduction to Screening Criteria and Enhanced Oil Recovery Field Projects. SPE Reservoir Engineering, 12(3), 189-198. doi: https://doi.org/10.2118/35385-PA

45

Torrealba, V.A. & Hoteit, H. (2019). Improved Polymer Flooding Injectivity and Displacement by Considering Compositionally-Tuned Slugs. Journal of Petroleum Science and Engineering, 178, 14-26. doi: https://doi.org/10.1016/j.petrol.2019.03.019 Vicente, B. J., Viatcheslav, I. P., & Pires, A. P. (2014). Semi-Analytical Solution for a 5

Hyperbolic System Modeling 1D Polymer Slug Flow in Porous Media. Journal of Petroleum

Science

and

Engineering,

115,

102-109.

doi:

https://doi.org/10.1016/j.petrol.2014.02.005 Zhou, W., Zhang, J., Feng, G., Jiang, W., Sun, F., Zhou, S., & Liu, Y. (2008). Key Technologies of Polymer Flooding in Offshore Oilfield of Bohai Bay. In SPE Asia 10

Pacific Oil and Gas Conference and Exhibition. Perth, Australia.

46

Appendix A In this appendix we present the detailed derivation of the solution of the auxiliary problem in the hodograph space, following the procedure presented in chapters 1 and 2 of Rhee et al. (2001). For the sake of clarity, we adopted the same notation as the one of 5

Rhee et al. (2001). Applying the chain rule in the auxiliary system (Equation 11), we find: H

&& 8Q + & 8Q + 8) = 0 8

8

8

;

& 8Q +  8Q + 8) = 0 ; 8

8

8

(A.1)

where,

A‘ = 8 D , ¨ = 1,2; ’ = 1,2 8#

F

10

(A.2)

It is possible to write ( and O as a function of %& and % applying the hodograph

transformation. This procedure can only be applied if the Jacobian matrix does not vanish nor approaches infinite for any concentration pair: o=

15

8



8

8Q ´ 8  8Q

8); ´ 8 8);

=

8 8

8Q 8);

− 8)

8 8 ;

8Q

≠ 0, o ∈ ℝ

(A.3)

If equation A.3 is satisfied, we can write ( = ( %& , %  and O = O%& , % , and apply

the chain rule to find the relation between the partial derivatives of %& and % , and the

hodograph variables: 8

=−

8

@

8); 8); 8Q

8

20

8Q

8

=

=

¸¹ ¸º

¸¹ ¸º

@

(A.5)

¸¹ ¸»;

@

=−

(A.4)

(A.6)

¸¹ ¸»;

@

(A.7) 47

Replacing equations A.4-A.7 in equation A.1, we find the system of equations in the hodograph space:

− + && ; − & ; = 0 8 8 8 ¼ 8Q 8); 8); + & −  =0 8Q

8)

8

8

8)

8

(A.8)

The characteristic velocity of system of equations (A.8) are the roots of: 5

& `  − && −  ` − & = 0

where ` = 1 (see section 1.3 of Rhee et al. (2001) for details).

(A.9)

1



Solving equation A.9 for `, we find

_ & & ` = Ye ½ =  & a&& −   ¦ b&& −   + 4& & d

(A.10)

& Γ : _ = R1 S = ` =  & a&& −   − b&& −   + 4& & d

(A.11)

The characteristic curves Γ and Γ are denoted as: 10

1









&

& Γ : e = R1 S = ` =  & a&& −   + b&& −   + 4& & d 1

&

(A.12)

In the phase plane (plane %& × % ), the characteristic parameter _ varies along Γ and it is always negative, whereas e varies along Γ and it is positive.

From Langmuir’s adsorption isotherm, we calculate the partial derivatives of the

15

adsorption isotherm with respect to the dissolved polymer concentrations: && = &C

C &C      C  

(A.13)

& = &C

C C 

   C  

(A.14)

& = &C

C C 

   C  

(A.15)

 = &C

C &C      C  

20

(A.16)

Substituting equations A.13-A.16 in equation A.9 it is possible to obtain: 

% R1 S − ¾ − % + %&  1 − %& = 0 1



1



(A.17) 48

where, ¾=

C C C C

(A.18)

Differentiating equation A.17 with respect to % leads to: 1  1

5

a%

1

1

− ¾ − % + %& d = 0

(A.19)

Solving the ordinary differential equation A.19 we find:

%& = `% −

¿À

À&

(A.20)

¾ + %& − %  + 4%& % = 0

(A.21)

Equation A.20 represents the characteristic curves (straight lines) in the phase plane.

10

The equations for the two characteristic waves are found applying the parameters _ and

e in equation A.20: Γ : %& = e% − l& ¿l

(A.22)

Γ : %& = _% − j& ¿j

(A.23)

As the parameter _ is negative and e is positive, the family Γ is composed by straight

15

lines with positive slopes, and the family Γ by straight lines with negative slopes.

For any pair _, e, the constant concentration state can be determined by the following relations (derived from equations A.22 and A.23): %& = − j&l& ¿jl

(A.24)

% = j&l& ¿

20

(A.25)

We now establish a relation between the phase plane and the ( × O plane. Consider the linear combination of the auxiliary system equations (Equation A.1):

f& 8) + f& && + f &  8Q + f 8) + f& & + f   8Q = 0 8

;

8

8

;

8

(A.26)

49

where f& and f are the eigenvalues of equation A.1. If Á represents a characteristic parameter (either _ or e), the derivatives

8 8Â

and

8 8Â

will be in the same direction if the

following system of equations (Rhee et al., 2001):

f& 5

f f&

8Q

8Â

8Q

8Â

8 8Â

− f& && + f & 

− f& & + f   + f

8 8Â

=0

f& && + f & 

8); 8Â

8); 8Â

=0

(A.27)

=0

(A.28) (A.29)

8 8Â

+ f& & + f   8Â = 0 8

(A.30)

is satisfied. From equations A.27 and A.29, one can find that: 10

8Q

8Â

− && − `& 

8); 8Â

=0

(A.31)

Moreover, we consider that the independent variables ( and O can be written as a

function of the parameter Á. The eigenvalues of the auxiliary system are given by f = 8) = 8Q

;

Thus,

15

−f 8Â 8Q

8); 8Â

¸º ¸Ã ¸»; ¸Ã

(A.32)

=0

(A.33)

Defining e as the characteristic parameter for the family
 (Γ family on the hodograph plane), and _ for the family
 (family Γ on the hodograph plane), we have

8Q 8l

− f

− f 8j

8Q

20

8); 8l

8); 8j

=0 =0

(A.34) (A.35)

Comparing equations A.34 and A.35 with equation A.31, we can establish the relation

between the hodograph plane and the ( × O plane: f = 1) = && − e& 1Q

;

(A.36) 50

f =

1Q

1);

= && − _&

(A.37)

Replacing equations A13.-A.16 in equations A.36 and A.37, we find l& 

f = g h RjkS RlkS j&

f = g h R 5

where h =

j&  l&

jk

C C

 = _ = i = ie =

S R

lk

(A.38)

S

(A.39)

. For convenience, we denote

j&

jk

(A.40)

l&

lk

(A.41)

The functions _ and ie represent a pair of characteristic parameters and their

values lie between 0 ≤ _ ≤ k ≤ ie ≤ 1. The plane _ × ie is spanned by &

10

orthogonal straight lines, where the  waves are vertical and i waves are horizontal. The characteristic velocities are now recast as: f = g hi 

(A.42)

f = g h i

(A.43)

After deriving the expressions for the characteristic waves of the auxiliary problem, the

15

shock velocities are found from the Rankine-Hugoniot conditions:

m , i  , i   =

Q

);

= g hi  i 

m  ,  , i = ) = g hi  Q

;

(A.44) (A.45)

where the superscripts + and – represent the value of  or i before and after the shock

wave, respectively. 20

Now we analyze how f and f change along the phase plane. We denote § the

derivative of f with respect to the concentration along the family Γ and § the

51

derivative of f with respect to the concentration along the family Γ . We will consider g > g& , and therefore, h > 1. Thus, we have: (Œ ‰ (

(‰ Œ (

5

¸Å‰ ¸Æ ¸¹ ¸Æ

=

= −2g h

¸Å‰ ¸Ç ¸¹ ¸Ç

=

= −2g h

j& l&œ jklkœ

j&œ l& jkœ lk

= −2g h #&  i ˜ < 0 &k

= −2g h {& ˜ i  < 0 &k

(A.46)

(A.47)

Since both derivatives (Equations A.46 and A.47) are negative, we conclude that f decreases along Γ and Γ as % increases.

Appendix B

10

In this appendix the inverse mapping of the waves from the ( × O plane to the ( × (

plane is presented. We also derive the exact coordinates of the crossing points of these waves related to the beginning and the ending points of the interactions between the waves.

The relation between the shock waves in the ( × O plane and ( × ( plane is given by (Pires et al., 2006):

15

§A = 9¦ˆ , ¨ = +, − <¦

D

(B.1)

where §A is the shock velocity in the ( × ( plane and mA is the shock velocity in the ( × O plane.

From equations (75) and (A.38), the velocity of the rarefaction family
 in ( × (

plane is: 20

1);

=C

<

1);

=C

<

1:;

 k#{

 9



(B.2)



(B.3)

and the velocity of the rarefaction family
 (Equations 75 and A.39): 1:;

 k#

 {9

 is given by: From equation (B.1), the shock path x

52

( = ( C



<‰

k#> { > 9 ‰

Following the same procedure, the shock path  y is: ( = (

C k#> 9‰

R1) S

z|

= g hi

z|

=

<‰

(B.4)

(B.5)

From equations (A.43) and (A.44), we find the shock velocity x
in ( × O plane: 5

1Q

;

@

−



~ > & È

…ƒ

C k{ > ~C k{ > );

(B.6)

Then, through equation (B.1), the velocity of shock x
in ( × ( plane is: R 1: ;S 1)

;

<‰

ɺ Ÿ É»; …ƒ

9 ‰ 

(B.7)

Therefore, we can write the shock path x
as:

( = R 1: ;S 1)

;

10

z|

( − z  + z

(B.8)

 intercepts the first rarefaction where z , z  is the point where the shock x

characteristic of
 , given by: z = C

C k#> { > 9‰

 k#

z = C



> { > "&#> $

<‰

k#> { > 9‰

z

(B.9) (B.10)

The rarefaction
 is transmitted through the shock x
, its new slope is:

15

1); 1:;

=C



<

k# { = 9



(B.11)

The shock x
finishes at point
, which is the point where the last characteristic of the rarefaction
 intercepts the shock x
. The coordinates of point
are:

| =

| =

ʉ ʉ   > ‰ :… )… Ë ‡aÈ= d Ì> ‰+‰ Ë Ì ‰+ ʉ ʉ   > ‰+‰ Ë Ì = > ‰  Ë ‡aÈ d Ì ‰+

C

<‰

 k«#= ¬ { > 9 ‰

| − 1

(B.12)

(B.13)

53

After the interaction with the rarefaction
 , the shock path continues from point
as a straight line up to point €, which is the point where the interaction between shock m

 is: with the rarefaction
 begins. The shock path
€ ( =

5

<‰

C { > 9 ‰

( − |  + |

(B.14)

and point € coordinates are: ‚ =

‚ =

ʉ ʉ :ƒ  )ƒ  > ‰ Ë Ì ‰+ Ë ‡È= aÌ> d ‰+‰ ʉ ʉ  Ë Ì> ‰+‰ Ë ‡È= aÌ> d ‰+‰ 

<‰



C k#= «{ > ¬ 9‰

r‚ − 1s

(B.15)

(B.16)

After point €, the rarefaction
 is absorbed and the path of shock €∞ is:

( = C

<‰

 {9

10

‰

( − ‚  + ‚

(B.17)

y . The last interaction occurs between the transmitted rarefaction
 and the shock 

The rarefaction is absorbed by the shock wave, and the new shock path is defined as: ( = C

<‰

 k#9

‰

( − „  + „

(B.18)

where the point „ , „  is given by:

„ = 15

ʉ :… )…  Ë ‡aÈ> d Ì= ‰+‰ ʉ ʉ   > ‰+‰ Ë ‡È > = ‰  Ë ‡aÈ d Ì ‰+

„ = C

<‰

 k#

> 9 ‰

„

(B.19)

(B.20)

54

Highlights •

The problem of two-component slug injection in an oil reservoir is solved;

•

The problem is decoupled into a thermodynamic system and a transport equation;

•

A full chromatographic cycle is developed in the porous media;

•

Pure water banks appear due to the separation of the chemicals in the reservoir.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: