=Paper= {{Paper |id=Vol-1623/paperapp1 |storemode=property |title=Numerical Simulation of Chemical Enhanced Oil Recovery Processes |pdfUrl=https://ceur-ws.org/Vol-1623/paperapp1.pdf |volume=Vol-1623 |authors=Bakhbergen Bekbauov, Abdumauvlen Berdyshev, Zharasbek Baishemirov |dblpUrl=https://dblp.org/rec/conf/door/BekbauovBB16 }} ==Numerical Simulation of Chemical Enhanced Oil Recovery Processes== https://ceur-ws.org/Vol-1623/paperapp1.pdf
     Numerical Simulation of Chemical Enhanced Oil
                  Recovery Processes

    Bakhbergen Bekbauov1⋆ , Abdumauvlen Berdyshev2 , and Zharasbek Baishemirov2
                 1
                 Al-farabi Kazakh National University, Almaty, Kazakhstan
                             bakhbergen.bekbauov@kaznu.kz
            2
              Abai Kazakh National Pedagogical University, Almaty, Kazakhstan
                        berdyshev@mail.ru, zbai.kz@gmail.com



        Abstract. In this paper we develop a new mathematical formulation for chem-
        ical compositional reservoir simulation, and provide a comparison of its results
        on alkaline-surfactant-polymer flooding with those of UTCHEM simulator. Our
        research has found that the existing chemical compositional model estimates the
        adsorption effect on the transport of a component reasonably well but it does
        not satisfy the principle of mass conservation. Since the total mass conserva-
        tion equation follows from summing the species-conservation equations over all
        components, the obtained equation violates the principle of total mass conser-
        vation as well. With these partial differential equations as governing equations,
        several simulators have been developed. In this work, we propose an approach
        to model the change in pore volume due to adsorption that satisfies the mass
        conservation law, and allows applying a sequential solution approach.

        Keywords: chemical compositional model,surfactant, porosity, adsorption


1     Introduction

Chemical flooding is one of the most promising and broadly applied enhanced oil re-
covery (EOR) processes. Chemical flooding can be further subdivided into alkaline
flooding, surfactant flooding, polymer flooding, and alkaline-surfactant-polymer (ASP)
flooding. Alkali reduces adsorption of the surfactant on the rock surfaces and reacts
with acids in the oil to create natural surfactant. Surfactants are chemicals that used
to reduce the interfacial tension between the involved fluids, increasing oil mobility.
ASP flooding is a form of chemical EOR method that can allow operators to extend
reservoir pool life and extract incremental reserves currently inaccessible by conven-
tional methods. While ASP flooding has a high efficiency it is technical, costly, and
risky. Model studies can assist in this evaluation.
    Most multiphase compositional models reported in the literature [1], [2], [3], [4] and
[5] are limited in their applicability in one way or another (single species, equilibrium
mass transfer, and lack of miscibility modeling etc). The mathematical formulation
⋆
    Corresponding author
Copyright ⃝
          c by the paper’s authors. Copying permitted for private and academic purposes.
In: A. Kononov et al. (eds.): DOOR 2016, Vladivostok, Russia, published at http://ceur-ws.org
               Numerical Simulation of Chemical Enhanced Oil Recovery Processes        665

developed in this work is extended from the UTCHEM model formulation for use in
chemical flooding studies that does not have these common limitations.
    Sequential schemes are very suitable for chemical compositional flow problems which
include a large number of chemical components. Only the implicit pressure and explicit
composition (IMPEC) formulation was used for chemical compositional reservoir sim-
ulation so far, but there is no obvious reason why the sequential formulation can’t
be used as well. Because of explicit solution of compositions, the size of time steps is
limited to stabilize the general procedure.
    Chen et al. [6] presented a numerical approach that solves both pressure and com-
positions implicitly. Though the approach was claimed to be sequential and extended
from the IMPEC approach used in UTCHEM model [7], the mathematical formulations
for the governing equations did not undergo any change in their model.
    The basic equations used in UTCHEM model that describe multiphase, multi-
component flow in permeable media are the species-conservation, pressure (an over-
all mass-continuity), and energy conservation equations. Accumulation terms in the
species-conservation equations used in UTCHEM model account for the reduction in
pore volume caused by adsorption.
    During the process of this research, it was revealed that this commonly used ap-
proach estimates the adsorption effect on the transport of a component reasonably
well but it does not satisfy the species-conservation equation. Since the total mass
conservation equation follows from summing the species-conservation equations over
all components, the obtained equation violates the principle of total mass conservation
as well. In recent years with use of these governing equations several simulators were
developed for simulation of the chemical flooding processes [6], [7], [8] and [9].
    In this work we introduce a new approach to model the reduction in pore volume
due to adsorption that satisfies the conservation equations. In certain situations, such as
significant change in the effective pore size due to adsorption, these enhancements are
essential to properly model the physical phenomena occurring in petroleum reservoirs.
In addition, this new approach for modeling the adsorption effect on the transport
of a component makes it possible to develop a new mathematical formulation for the
sequential chemical compositional reservoir simulation.



2   Mathematical model

Consider a bulk volume Vb at some point within a porous medium domain. Let us
assume that this representative elementary volume (REV) is made up of np + 1 phases
(np fluid phases and a solid phase consisting of rock grains or soil) with nc chemical
species. Conceivably, at least, each species can exist in any phase and can transfer
between phases via evaporation, condensation, dissolution, adsorption and so forth.
    In our model formalism, each pair (i, α), with i chosen from the species indices and
α chosen from the phases, is a constituent. Each constituent (i, α) has its own intrinsic
mass density ρiα , measured as mass of i per unit volume of α, and its own average
velocity −
         →u iα . Each phase α has its own volume fraction ϕα . The volume fraction of
phase α, ϕα , is the volume of phase α divided by the bulk volume Vb .
666     B. Bekbauov et al.

  If the index nc represents the species and the index np + 1 represents the phases
making up the solid phase, then in terms of the above defined mechanical variables the
mass balance for each constituent (i, α) is
                                                             {                     }
       ∂                        →
                                −                              i = 1, ..., nc ,
         (ϕα ρiα ) + ∇ · (ϕα ρiα u iα ) = Riα + rmiα + qiα ,                         (1)
      ∂t                                                       α = 1, ..., np + 1.

From left to right in Eq. (1), the terms are now the accumulation, transport, and source
terms, the last consisting of three types. The mass fraction of component i in phase α
in Vb is defined to be ωiα . The parameter∑nc ωiα is the mass of component i in phase α
divided by mass of the phase. Hence, i=1         ωiα = 1. With that definition,
                                      {                                      }
                  ρiα = ρα ωiα ,        i = 1, ..., nc , α = 1, ..., np + 1,         (2)

where ρα is intrinsic mass density of phase α.
   The source term Riα accounts for the rate of mass generation (Riα > 0) and
consumption (Riα < 0) of component i in phase α, either through chemical or biological
reactions. There is no general function for Riα . An example of a first-order reaction
rate for radioactive decay or biodegradation is
                                       {                                      }
              Riα = −ki ϕα ρα ωiα ,      i = 1, ..., nc , α = 1, ..., np + 1,      (3)

where ki is the decay constant or reaction rate coefficient in units of inverse time.
   The second source term rmiα expresses the rate of mass transfer of component i
from or into the phase α owing to vaporization or condensation. Adsorption is described
through isotherms.
   The last source term in Eq. (1) qiα represents physical sources (wells).
   Substitution of Eqs.(2) and (3) into Eq. (1) gives:

               (ϕα ρα ωiα ) + ∇ · (ϕα ρα ωiα −
                                             →
            ∂
                                             u iα ) = −ki ϕα ρα ωiα + rmiα + qiα ,
            ∂t                                                                             (4)
                                              {i = 1, ..., nc , α = 1, ..., np + 1.}
                                                                  ∑np +1
From definition, volume fractions must obey the constraint α=1            ϕα = 1. It is well
known that the porosity ϕ is defined as the fraction of the bulk permeable medium
that is pore space, that is, the pore volume Vp divided by the bulk volume   ∑Vnpb . The fact
that the all fluid phases jointly fill the voids (pores) implies the relation α=1    ϕα = ϕ.
    The phase saturation Sα is defined as the fraction of the pore volume occupied by
phase α, that is, volume of phase α Vα divided by the pore volume Vp . The saturation
of fluid phase α can also be defined as Sα = ϕα /ϕ. For fluid phases such as liquids and
vapors, ϕα = ϕSα , α = 1, ..., np , where ϕSα also called the fluid content. For the solid
(s) phase ϕs = 1 − ϕ, which is the grain volume divided by the bulk volume Vb . We can
rewrite equation (4) in the following form by noting that the porosity is ϕ = 1 − ϕs
and defining the fluid saturations Sα = ϕα /ϕ, α = 1, ..., np :

             (ϕSα ρα ωiα ) + ∇ · (ϕSα ρα ωiα →
                                             −
          ∂
                                             u iα ) = −ki ϕSα ρα ωiα + rmiα + qiα ,
          ∂t                                                                               (5)
                                                    {i = 1, ..., nc , α = 1, ..., np , }
                Numerical Simulation of Chemical Enhanced Oil Recovery Processes              667

for the fluids, and if we fix a coordinate system in which −
                                                           →
                                                           u is = 0, and note that qis = 0,
then the momentum balance for the solid phase reduces to
             ∂
                ((1 − ϕ)ρs ωis ) = −ki (1 − ϕ)ρs ωis + rmis ,       {i = 1, ..., nc .}         (6)
             ∂t
The statistical average apparent velocity of constituent (i, α) owing to both convection
and dispersion is the sum of the barycentric velocity of phase α and the diffusion
velocity of species i in phase α:
                                −
                                →′
                  −
                  →
                  u iα = −→
                          u α + u iα ,   {i = 1, ..., nc , α = 1, ..., np .}         (7)

Since phase velocities are typically more accessible to measurement than species veloc-
ities, it is convenient to rewrite the constituent mass balance equation (5) as
                                                  −
                                                  →
           (ϕSα ρiα ) + ∇ · (ϕSα ρiα →
                                     −
        ∂
                                     u iα ) + ∇ · J Diα = −ki ϕSα ρiα + rmiα + qiα ,
        ∂t                                                                                     (8)
                                                      {i = 1, ..., nc , α = 1, ..., np , }
       −
       →                −
                        →′
where J Diα = ϕSα ρiα u iα stands for the diffusive flux of constituent (i, α).
    So far, the mathematical formulation of the mass conservation equations developed
above is essentially the same as the standard formulation described in [7]; where it
differs is in the treatment of average velocity in the governing equations. Here we start
to deviate from the standard formulation.
    The fluxes of component i in phase α with respect to volume-averaged velocity
−
→                ¯ · ∇(ρ ω ) and mass-averaged velocity −        →               ¯ · ∇(ω )
J Diα = −ϕSα K̄   iα     α iα                                    J Diα = −ϕSα K̄    iα  iα
owing to hydrodynamic dispersion alone were presented in [10]. The flux with respect
to bulk volume-averaged velocity is proposed in this work:
            −
            →          ¯ · ∇(ϕS ρ ω ),
            J Diα = −K̄ iα        α α iα        {i = 1, ..., nc , α = 1, ..., np , }    (9)
                   ¯ for a homogeneous, isotropic permeable medium [11] are
Two components of K̄iα

                      Diα    αlα u2xα + αtα (u2yα + u2zα )         {                      }
          (Kxx )iα =       +
                                        |−
                                         →                 ,
                       τ                 u α|                          i = 1, ..., nc ,
                                                                                              (10)
                      (αlα − αtα )uxα uyα                              α = 1, ..., np ,
           (Kxy )iα =
                             |→
                              −           ,
                              u α|
where the subscript l refers to the spatial coordinate in the direction parallel, or longi-
tudinal, to bulk flow, and t is any direction perpendicular, or transverse, to l. Diα is the
effective binary diffusion coefficient of component i in phase α [12], αlα and αtα are the
longitudinal and transverse dispersivities, and τ is the permeable medium tortuosity.
    A general set of partial differential equations (11) for the conservation of component
i in fluid phase α is obtained upon substitution of the definition for flux (Eq. (9)) into
Eq. (8):

              (ϕSα ρα ωiα ) + ∇ · (ϕSα ρα ωiα −
                                              →
           ∂                                            ¯ · ∇(ϕS ρ ω ) =
                                              u iα ) − K̄iα         α α iα
           ∂t                                                                                 (11)
           = −ki ϕSα ρα ωiα + rmiα + qiα ,        {i = 1, ..., nc , α = 1, ..., np , }
668      B. Bekbauov et al.

The term rmiα is difficult to calculate without detailed analysis of the transport oc-
curring within the phases. One typically simplifies the equations by using overall com-
positional balance equations. Overall compositional balance equations can be obtained
by summing Eqs. (6) and (11) over the solid and np fluid phases:
             [ np                      ]
           ∂ ∑                                ∑np
                   ϕSα ρiα + (1 − ϕ)ρis + ∇ ·     (ϕSα ρiα −
                                                           →
                                                           u α) −
          ∂t α=1                              α=1
                                            [ np                      ]
                 ∑np
                                              ∑                                     (12)
           −∇ ·        ¯ · ∇(ϕS ρ )] = −k
                     [K̄                           ϕS ρ + (1 − ϕ)ρ      +Q ,
                          iα          α iα            i             α iα                is         i
                 α=1                                      α=1
                                                                               {i = 1, ..., nc , }
            ∑np
where Qi = α=1    q is the injection/production rate for component i per bulk volume.
         ∑np +1 iα
We have α=1 rmiα = 0, a relation following from the inability to accumulate mass
at a volumeless phase interface.
    In equation (12), expressing the content of component i in phase α in terms of
volume fraction
                                                   }
                  ρα ωiα = ρi ciα , α = 1, ..., np
                                                      {i = 1, ..., nc , }         (13)
                  (1 − ϕ)ρs ωis = ϕρi ĉi ,

we get
                [ ( np                  )]
                       ∑                             ∑np
                                                         (Sα ciα −
                                                                 →
             ∂
                 ϕρi       Sα ciα + ĉi    + ∇ · ϕρi             u α )−
             ∂t        α=1                           α=1
                     np [
                                                         ( np           )
                   ∑                        ]             ∑                                            (14)
              −∇ ·        ¯
                         K̄ ∇(ϕρ S c ) = −k ϕρ                 S c + ĉ + Q ,
                               iα·     i α iα              i    i          α iα     i        i
                    α=1                                             α=1
                                                                             {i = 1, ..., nc , }

where ρi is the component mass density in units of mass of i per unit volume of i, ciα
is the component concentration in units of volume of i in phase α per unit volume of α,
and ĉi is the adsorbed component concentration, measured as volume of i in phase α
per unit pore volume. The linear, Freundlich and Langmuir adsorption isotherm models
are applied to calculate the adsorbed concentrations ĉi . In our definition
                                  ( np             )
                               ∑
                               nc   ∑
                                       Sα ciα + cˆi = 1,                          (15)
                                     i=1    α=1

but
                                           ∑
                                           nc ∑
                                              np
                                                     Sα ciα ̸= 1,                                      (16)
                                           i=1 α=1

unlike existing model. To∑account for the reduction in pore volume caused by adsorp-
                           ncv
tion, the coefficient (1 − i=1   ĉi ) is introduced into the overall compositional balance
               Numerical Simulation of Chemical Enhanced Oil Recovery Processes                669

equation (14) in UTCHEM model. The coefficient represents reduction in pore volume
due to adsorption, ĉi is the adsorbed concentration of species i, and ncv is the total
number of volume-occupying components. During the process of this research, it was re-
vealed that even though this approach estimates the adsorption effect on the transport
of a component reasonably well, it does not satisfy the species-conservation equation
since the coefficient is multiplied only to the first summand of the accumulation term
in Eq. (14). It is well known that an equation remains balanced when both sides of an
equation are multiplied by the same nonzero quantity.
    In the present work we introduce a new approach to model the reduction in pore
volume due to adsorption that satisfies the continuity equation. Let us denote the
modified volume fraction of phase α due to adsorption by ϕ̂α . Porosity ϕ̂ is defined as the
fraction of the bulk permeable medium that is pore space remaining after adsorption.
This porosity is related to the original porosity ϕ as follows:
                                         (            )
                                               ∑
                                               ncv
                                   ϕ̂ = ϕ 1 −      ĉi .                                (17)
                                                 i=1

The saturation of fluid phase α is defined as Sα = ϕ̂α /ϕ̂. Now using the same derivation
procedure as carried out above, we obtain
             [ ( np                  )]
                    ∑                              ∑np
                                                       (Sα ciα →
                                                               −
          ∂
              ϕ̂ρi      Sα ciα + ĉi    + ∇ · ϕ̂ρi             u α) −
          ∂t        α=1                            α=1
                    np [
                                                         ( np           )
                   ∑                       ]                ∑                        (18)
             −∇ ·        ¯
                        K̄ · ∇(ϕ̂ρ S c ) = −k ϕ̂ρ               S c + ĉ + Q ,
                          iα        i α iα             i   i          α iα       i       i
                   α=1                                          α=1
                                                                         {i = 1, ..., nc , }
   The phase flux from Darcy’s law is

                  −
                  →       k̄¯krα
                  uα =−          (∇pα − γα ∇z),                {α = 1, ..., np , }             (19)
                        ϕ̂Sα µα

where k̄¯ is the permeability tensor, krα is the relative permeability of fluid phase α, µα
is the dynamic viscosity of fluid phase α, pα is the pressure in fluid phase α, γα is the
specific weight for fluid phase α, and z represents depth.
    Variation of pore volume with pore pressure p can be taken into account by the
pressure dependence of porosity. The porosity depends on pressure due to rock com-
pressibility, which is often assumed to be constant and can be defined as

                                 ϕ = ϕR [1 + cf (p1 − ps )],                                   (20)

where ϕR is the porosity at a specific pressure ps , p1 is the water phase pressure, and
cf is the pore compressibility at ps .
    A slightly compressible fluid has a small but constant compressibility. For a slightly
compressible fluid, the component density ρi can be written as:

                     ρi = ρiR [1 + c0i (p1 − pR )],        {i = 1, ..., nc , }                 (21)
670     B. Bekbauov et al.

where ρiR is the density of component i at the standard pressure pR , a constant value,
c0i is the compressibility of component i.
     Now since reference density ρiR is constant for each component we can divide
through both sides of Eq. (18) by ρiR . In terms of the dimensionless density ρ̄i = ρi /ρiR
Eq. (18) can be written as:
              [ ( np                   )]
                      ∑                                 ∑np
                                                            (Sα ciα −
                                                                    →
            ∂
               ϕ̂ρ̄i      Sα ciα + ĉi      + ∇ · ϕ̂ρ̄i             u α) −
           ∂t         α=1                               α=1
                      np [
                                                               ( np              )
                     ∑                            ]              ∑                  Qi    (22)
              −∇ ·          ¯
                           K̄iα · ∇(ϕ̂ρ̄i Sα ciα ) = −ki ϕ̂ρ̄i       Sα ciα + ĉi +     ,
                     α=1                                         α=1
                                                                                    ρiR
                                                                                              {i = 1, ..., nc , }

   We sum the mass balance equations above over the nc components to obtain the
equation of continuity, or conservation of total mass. The equation of continuity is
                                             { np                        }
                                                 ∑          ∑nc             ∑nc
          ϕFt (c̃i , ĉi ) + ϕ̂R ct
                                    ∂p1
                                        + ∇ · ϕ̂       −
                                                       →
                                                    (Sα u α     ρ̄i ciα ) =
                                                                                Qi
                                                                                   , (23)
                                     ∂t          α=1        i=1
                                                                                ρ
                                                                            i=1 iR

where we used that
                            ∑
                            nc
                                      −
                                      →
                                  ∇ · J Diα = 0,              {α = 1, ..., np , }                                   (24)
                            i=1

(net dispersive flux in a phase is zero), and according to the total reaction definition
                                   ( np            )
                               ∑nc   ∑
                                          Riα + Ris = 0.                             (25)
                                    i=1       α=1

The total compressibility, ct , is
                (               ){                                           }
                        ∑
                        ncv            [                        ]∑
                                                                 nc
            ct = 1 −        ĉi    cf + 1 + cf (2p1 − ps − pR )       0
                                                                    (ci c̃i ) ,                                     (26)
                           i=1                                                               i=1

and                                           [(                    ) n                ]       (n      )
                                         ∂              ∑
                                                        ncv          ∑ c
                                                                                             ∂ ∑ cv

            Ft (c̃i , ĉi ) = (p1 − pR )           1−         ĉi            c0i c̃i       −        ĉi .           (27)
                                         ∂t             i=1            i=1
                                                                                             ∂t i=1
We define the overall concentration c̃i as
                                   ∑
                                   np
                           c̃i =         Sα ciα + ĉi ,             {i = 1, ..., nc , }                             (28)
                                   α=1

and by definition
                                                ∑
                                                nc
                                                      c̃i = 1.                                                      (29)
                                                i=1
                Numerical Simulation of Chemical Enhanced Oil Recovery Processes                        671

     The pressure equation is developed by substituting Darcy’s law (Eq. (19)) for the
phase flux term of Eq. (23), using the definition of capillary pressure pcα1 = pα −p1 , α =
2, ..., np . The pressure equation in terms of the reference phase (phase 1) pressure is

                        ∂p1      (            )     ( ∑  np           )
                 ϕ̂R ct            ¯                 ¯
                            − ∇ · k̄ λrT c ∇p1 = ∇ · k̄     λrαc ∇pcα1 −
                         ∂t                             α=2
                                                                                                        (30)
                                ( ∑   np           )                     ∑nc
                                                                             Qi
                            −∇ · k̄¯     λrαc γα ∇z − ϕFt (c̃i , ĉi ) +        ,
                                     α=1
                                                                             ρ
                                                                         i=1 iR

where
                                        ∑
                                        nc
                           λrαc = λrα         ρ¯i ciα ,         {α = 1, ..., np , }                     (31)
                                        i=1
and total relative mobility is
                                                     ∑
                                                     np
                                        λrT c =               λrαc .                                    (32)
                                                     α=1
The extension of the LET correlations is used to represent the relative permeability
and capillary pressure curves [13], [14].
   Applying the mean value estimate for character sums in Eq. (22), we can write
                     ∑
                     np
                                        → ∑
                                        −
                                               np
                           Sα ciα −
                                  →
                                  u α = ũ i Sα ciα ,                  {i = 1, ..., nc , }              (33)
                     α=1                      α=1

and
          ∑                               ¯ ∑ (                 )
          np                                   np
                ¯
               K̄iα · ∇(ϕ̂ρ̄i Sα ciα ) = K̃i ·    ∇ ϕ̂ρ̄i Sα ciα ,                {i = 1, ..., nc , }   (34)
         α=1                                    α=1
      −
      →         ¯
where ũ i and K̃i can be defined as some averages. Since differentiation and summation
are interchangeable operations in this system, the sum of the gradients can be calculated
as the gradient of the sum
                                          (              )
           ¯ ∑ (               )                ∑
               np                               np
                                     ¯
          K̃ ·
           i       ∇ ϕ̂ρ̄ S c
                           i α iα= K̃ · ∇ ϕ̂ρ̄
                                          i         S c   i,     {i = 1, ..., n .}
                                                                    α iα             (35)         c
               α=1                                            α=1

   Equation (22) can be written using Eqs. (33), (34), and (35) as below:
        [ ( np                  )]        (                    )
     ∂         ∑                              → ∑
                                              −
                                                     np
         ϕ̂ρ̄i     Sα ciα + ĉi    + ∇ · ϕ̂ρ̄i ũ i      Sα ciα −
     ∂t        α=1                                  α=1
               [                          ]             ( np              )
                 ¯     (     ∑np        )                 ∑                  Qi                         (36)
         −∇ · K̃i · ∇ ϕ̂ρ̄i      Sα ciα     = −ki ϕ̂ρ̄i       Sα ciα + ĉi +     ,
                             α=1                          α=1
                                                                             ρiR

                                                                                  {i = 1, ..., nc , }
    This new mathematical formulation of species conservation equations makes it pos-
sible to apply a sequential solution approach to solve these equations implicitly for the
672     B. Bekbauov et al.

                   ∑np
total concentration α=1   Sα ciα of each component. A flash calculation is then per-
formed to obtain the phase saturations and the concentrations of components in each
                           −
                           →         ¯
phase. Numerical values of ũ i and K̃i can be most simply calculated as the weighted
averages

                ∑
                np                         ∑
                                           np
                     Sα ciα −
                            →
                            uα                 ¯ · ∇(ϕ̂ρ S c )
                                              K̄iα      i α iα
        −
        →                             ¯
        ũ i = α=1np             ,   K̃i = α=1 (              )   ,   {i = 1, ..., nc , }   (37)
                 ∑                                  ∑
                                                    np
                       Sα ciα                ∇ ϕ̂ρi    Sα ciα
                 α=1                                α=1

obtained from previous time-step values.
    The sequential solution procedure is carried out in the following order: (a) solution
of the pressure equation (30) implicitly, (b) solution of the transport system of equations
(36) implicitly for the total concentration of each component.
    Details about the derivation of the mathematical model formulation are provided
in our previous publication [15].




3     Numerical Results

ASP flooding is the most promising EOR solution for one of the greatest challenges
facing the oil industry worldwide: after conventional water flooding the residual oil
(drops trapped by capillary forces) in reservoirs around the world is likely to be around
70% of the original oil in place. The mathematical formulation is evaluated in the
modeling of a field scale ASP EOR process.




               Fig. 1. Computational domain and well pattern illustration
               Numerical Simulation of Chemical Enhanced Oil Recovery Processes        673

    As illustrated in Fig. 1, the ASP flooding pilot has 4 injection wells and 9 production
wells in an inverted five-spot well pattern. The ASP process was conducted in a 4-slug
sequence: pre-flush polymer flood, alkaline/surfactant slug, alkaline/surfactant/polymer
slug, and a polymer drive. Total simulation time is 551 days. Reservoir properties in-
clude heterogeneous permeability and initial water saturation fields. The reservoir is
at a depth of 4150 ft., has an average initial pressure of 1770 psi, and the porosity is
assumed to be constant throughout the reservoir and equal to 0.3. Grid dimensions are
19 × 19 × 3. The OOIP is 395,427 bbls, the crude oil viscosity is 40 cp, the initial brine
salinity is 0.0583 meq/ml and the initial brine divalent cation concentration is 0.0025
meq/ml.
    We use S3GRAF software, developed and licensed by Sciencesoft Ltd., for post-
processing the output data.
    Three flowing phases and eleven components are considered in the numerical simu-
lations. The phases are water, oil and microemulsion, while the components are water,
oil, surfactant, polymer, chloride anions, divalent cations (Ca++, Mg++), carbonate,
sodium, hydrogen ion, and oil acid. The ASP interactions are modeled using the re-
actions: in situ generated surfactant, precipitation and dissolution of minerals, cation
exchange with clay and micelle, and chemical adsorption. Note the detailed chemical
reaction modeling, and the heterogeneous and multiphase petroleum reservoir under
consideration.
    A comparison with UTCHEM has also been performed. The matches between old
and new formulations’ numerical results for the matched variables are shown in Figs. 2-
5 for the injected pore volume in the range 0 - 1.0 PV. Comparative studies show that




                  Fig. 2. Average pressure vs. total injected pore volume



the results obtained from IMPEC implementation of the newly proposed formulation
are in a good agreement with that of UTCHEM simulator. In the scope of this research
work, through its application to the above-mentioned numerical experiment and com-
parisons with UTCHEM model results, the newly developed formulation has proven to
be reliable, practical, and accurate. The mathematical model and numerical simulation
developed in this work can also be used to study the transport of contaminants and
674    B. Bekbauov et al.




            Fig. 3. Oil and water saturations vs. total injected pore volume




               Fig. 4. Microemulsion saturation vs. injected pore volume




 Fig. 5. Adsorbed surfactant ratio (ML per ML of pore volume) vs. injected pore volume
               Numerical Simulation of Chemical Enhanced Oil Recovery Processes         675

remediation of contaminated aquifers surfactants.



4     Conclusion

In the scope of this research work, a new mathematical model formulation for multi-
component, multiphase flow in porous media has been developed. During the process
of this research, it was revealed that commonly used approach estimates the adsorption
effect on the transport of a component reasonably well but it does not satisfy the mass
conservation or continuity equation. In the present work we introduce a new approach
to model the reduction in pore volume due to adsorption that satisfies the continuity
equation. The mathematical formulation developed in the scope of this work is ex-
tended from the UTCHEM model formulation for use in chemical flooding studies. A
comparison with UTCHEM has also been performed. Comparative studies show that
the results obtained from IMPEC implementation of the newly proposed formulation
are in a good agreement with that of UTCHEM simulator. The implementation of a
sequential solution approach for chemical compositional reservoir simulation based on
the formulation described in this paper is scheduled for the future.



Acknowledgments. This paper was supported by the Ministry of Science and Ed-
ucation of the Republic of Kazakhstan under grants No. 1735/GF4 and No. 0128/GF4.



References

1. Baehr, A.L. and Corapcioglu, M.Y.: Ground water Contamination by Petroleum Products:
   2. Numerical Solution. Water Resour. Res., 23(10), 201 (1987)
2. Mayer, A.S. and Miller, C.T.: A Compositional Model for Simulating Multiphase Flow,
   Transport and Mass Transfer in Groundwater Systems. Paper presented at the Eighth
   International Conference on Computational Methods in Water Resources, Venice, Italy,
   June 11-15. (1990)
3. Sleep, B.E. and Sykes, J.F.: Compositional Simulation of Groundwater Contamination by
   Organic Compounds: 1. Model Development and Verification. Water Resour. Res., 29(6),
   1697-1708, June. (1993)
4. Abriola, L.M. and Pinder, G.F.: A Multiphase Approach to the Modeling of Porous Media
   Contamination by Organic Compounds: 2. Numerical Simulation. Water Resources Res.,
   21 (1), (1985a)
5. Kalurachchi, J.J. and Parker, J.C.: Modeling Multicomponent Organic Chemical Transport
   in Three-Phase Porous Media. J. of Contaminant Transport, 5, 349. (1990)
6. Chen, Z., Ma, Y., and Chen, G.: A sequential numerical chemical compositional simulator.
   Transport in Porous Media 68, 389-411 (2007)
7. Delshad, M., Pope, G.A., Sepehrnoori, K.: UTCHEM Version-9.0, Technical Documenta-
   tion, Center for Petroleum and Geosystems Engineering. The University of Texas at Austin,
   Texas. (2000)
676     B. Bekbauov et al.

8. Luo, H., Al-Shalabi, E.W., Delshad, M., Panthi, K., and Sepehrnoori, K.: A Robust
   Geochemical Simulator to Model Improved Oil Recovery Methods. SPEJ (Preprint,
   SPE173211-PA). (2015)
9. Luo, H., Delshad, M., Li, Z., and Shahmoradi, A.: Numerical simulation of the impact of
   polymer rheology on polymer injectivity using a multilevel local grid refinement method.
   Petroleum Science, 13(1): 110-125. (2016)
10. Lake, L.W.: Enhanced Oil Recovery” (Englewood Cliffs, NJ: Prentice Hall Inc.), 550 pp.
   (1989)
11. Bear, J.: Dynamics of Fluids in Porous Media, Dover, New York. (1972)
12. Bird, R.B., Stewart, W.E., and Lightfoot, E.N.: Transport Phenomena, 2nd edition, John
   Wiley and Sons, New York. (2002)
13. Bekbauov, B. E., Kaltayev, A., Wojtanowicz, A. K., Panfilov, M.: Numerical Modeling
   of the Effects of Disproportionate Permeability Reduction Water-Shutoff Treatments on
   Water Coning. J. Energy Resour. Technol. 135(1), 011101. Paper No: JERT-12-1134; doi:
   10.1115/1.4007913 (2013)
14. Bekbauov, B. E., Kaltayev, A., and Nagy, S.: Three-Dimensional Thermal Petroleum
   Filtration Study of Water Coning. Arch. Min. Sci., 55(1), pp. 201215. (2010)
15. Bekbauov, B. E., Kaltayev, A., Berdyshev, A.: A New Mathematical Formulation of
   the Governing Equations for the Chemical Compositional Simulation // arXiv:1512.08170
   [physics.flu-dyn] (2015)