[1]A.C. Alvarez These authors contributed equally to this work.

These authors contributed equally to this work.

These authors contributed equally to this work.

[1]Instituto de Computação, Universidade Federal do Rio de Janeiro, Avenida
Athos da Silveira Ramos, 274, Rio de Janeiro, 21941-590, Rio de Janeiro,
Brasil

2]Lab. Fluid Dynamics, IMPA, Estrada Dona Castorina, 110, Rio de Janeiro,
22460-320, Rio de Janeiro, Brasil

3] Civil Engineering and Geosciences, TU Delft, Stevinweg 1, Delft, 2628 CE,
Delft, The Netherlands

# Modeling low saline carbonated water flooding
including surface complexes

<EMAIL_ADDRESS>J. Bruining<EMAIL_ADDRESS>D. Marchesin
<EMAIL_ADDRESS>*  [  [

###### Abstract

Carbonated water flooding (CWI) increases oil production due to favorable
dissolution effects and viscosity reduction. Accurate modeling of CWI
performance requires a simulator with the ability to capture the true physics
of such process. In this study, compositional modeling coupled with surface
complexation modeling (SCM) are done, allowing a unified study of the
influence in oil recovery of reduction of salt concentration in water. The
compositional model consists of the conservation equations of total carbon,
hydrogen, oxygen, chloride and decane. The coefficients of such equations are
obtained from the equilibrium partition of chemical species that are soluble
both in oleic and the aqueous phases. SCM is done by using the PHREEQC
program, which determines concentration of the master species. Estimation of
the wettability as a function of the Total Bound Product (TBP) that takes into
account the concentration of the complexes in the aqueous, oleic phases and in
the rock walls is performed. We solve analytically and numerically these
equations in $\displaystyle 1-$D in order to elucidate the effects of the
injection of low salinity carbonated water into a reservoir containing oil
equilibrated with high salinity carbonated water.

###### keywords:

surface complexation modeling,wettability, carbonate water flooding,
conservation laws

##  1 Introduction

Carbonated Water Injection (CWI) is an oil recovery technique that increases
the oil production due to favorable dissolution effects causing viscosity
reduction [19, 9]. When this method is combined with low salinity brine
injection, the brine behaves like a natural solvent that enhances oil recovery
[17, 54]. In this way CWI takes advantage of brine flooding to improve oil
recovery from 5 to 20% of the oil initially in place (OIIP) [27, 40, 63]. This
happens because geochemical reactions between the injected carbonated brine
and rock can alter the petrophysical properties of the reservoir.

While there are several studies showing the relevance of CWI coreflood for
enhanced oil recovery, understanding such process is still a challenge. Among
such processes are the effects of wettability, the flow of mineral salts in
water and flow of CO2 in the aqueous and oleic phases. In this paper, we use
multiphase compositional modeling to quantify the oil and water saturations as
well as the CO2 transfer between the aqueous and oleic phases. As novelty we
take into account the formation of surface complexes, which participate in the
mechanism of wettability modification due to charge transfer [35, 7]. With
this new effect we can quantify the intrinsic nonlinear relation leading to
optimal recovery conditions. We also evaluate the modifications of
$\displaystyle pH$, of the velocity of the saline front and of the magnitude
of the saturation shock. These changes are the main factors for the increase
in oil productivity.

Modeling enhanced oil recovery (EOR) by conventional
fractional$\displaystyle-$flow theory [16, 45] has had numerous extensions to
take into account different mechanisms that contribute to injectivity. Some
model extensions have included CO2 [66, 18, 59] and molecular diffusion [30].
Thus compositional modeling has been developed to include master chemical
species in the aqueous and oleic phases [64, 26]. More recently some authors
have included in the model the transport of ions [1, 6, 61]. Some authors
emphasize that surface complex formation is a mechanism responsible for
increased oil recovery [3]. It has often been stated that the injection of
water with low salinity and CO2 produces chemical reactions that change
wettability of rock walls favorable to enhanced oil recovery [36].

Some research underlines that wettability alteration changes the oil recovery
when carbonated water at low salinity is injected, see, e.g., [70, 69].
Understanding this recovery method is one of the purposes of this work, for
conditions typical of the Brazilian pre-salt basin.

In [33] it is highlighted that the formation of surface complexes contributes
to changes in wettability. This happens because charges on the calcite and
oleic surfaces have opposite signs, which contributes to squeezing out the
water film; this effect forms an oil-wet surface [31]. This effect can be
explained in terms of changes in $\displaystyle pH$ values and of larger shock
amplitudes in the saturation profile. Understanding this phenomenon is an
expected result of this work.

To incorporate this process, we use multiphase compositional modeling with
master species. We consider the concentrations of ions generated by all
equilibrium reactions present in the carbonated brine and CO2 in the aqueous
and oleic phases. To do so, we use a one dimensional incompressible flow model
that describes two-phase flow with geochemical modeling ([26]). These
processes are studied by means of a system of balance laws for the transported
quantities. Using this model we study the flow of oil, water and dissolved
carbon dioxide in a sandstone rock.

The geochemical data are obtained by utilizing the modeling capabilities of
the program PHREEQC (acronym of pH-REdox-Equilibrium C-program). This is a
computer program written in C++ designed to perform a wide variety of aqueous
geochemical calculations (see details about its implementation in [58, 57, 8,
7]). Furthermore, we use the program for calculation of the surface
complexation in calcite rock [65]. The program determines the concentration of
ions and molecules dissolved in water inclusive carbon dioxide (CO2), which is
the only compound that exists in both phases. By assuming chemical
equilibrium, we define the behavior of all dissolved compounds by means of
four transport equations of master species, which are chloride, oil, twice
oxygen minus hydrogen and inorganic carbon.

We also consider the ions ($\displaystyle\mbox{H}^{+}$,
$\displaystyle\mbox{OH}^{-}$, $\displaystyle\mbox{CO}_{3}^{2-}$,
$\displaystyle\mbox{HCO}_{3}^{-}$, $\displaystyle\mbox{Cl}^{-}$), water (H2O)
and the sorbed species $\displaystyle oil_{s,NH^{+}}$, $\displaystyle
oil_{w,COOH}$, $\displaystyle Cal_{s,OH}$, $\displaystyle Cal_{w,CO_{3}H}$,
$\displaystyle oil_{s,N}$, $\displaystyle oil_{w,COO^{-}}$, $\displaystyle
oil_{w,COOCa^{+}}$, $\displaystyle Cal_{s,OH_{2}^{+}}$, $\displaystyle
Cal_{s,CO_{3}^{-}}$, $\displaystyle Cal_{w,CO_{3}^{-}}$, $\displaystyle
Cal_{w,CO_{3}Ca^{+}}$ [12]. The species with complexes lead to what we call
the surface complexes$\displaystyle-$chloride ionic carbon
dioxide$\displaystyle-$oil$\displaystyle-$water SC-CLICDOW model. This model
summarizes the most relevant concepts described in [14, 15, 32, 41, 42]. To
illustrate these effects we inject water with the same $\displaystyle pH$ but
with lower salinity than the initial salinity in the core.

In this study, we use SCM through the geochemistry PHREEQC program. This
solver is utilized to calculate the concentration of the complexes and to
predict the wettability alteration of minerals through the Total Bound Product
(TBP). Details of the procedure employed here can found in [51]. SCM is a
chemical equilibrium technique used to model the interactions of water, oil,
brine and rock [49, 23, 23]. Such a method has been used in several studies
for characterizing the surface adsorption phenomenon [50, 60] and determining
the wettability of minerals at reservoir conditions [28, 11, 51].

In these works the model is based on adsorption of aqueous solute into surface
functional groups characterized by a set of chemical reactions. Moreover, some
experiments indicate that since carbonates are salt-type minerals, their
surface reactivity is different from sandstone and clays. Another conclusion
is that dissolution and precipitation interfere on ion adsorption dynamics at
functional sites. The main effects that influence the oil-brine-rock
interaction are the brine chemistry and the oil composition, i.e., acid and
base number [20].

Several works describe how to quantity the relationship between salt
concentration and water flooding process, e.g., [35, 56, 2]. In this paper we
use a similar idea, which is proposed by [35] utilizing the relative
permeability functions based on Corey’s correlation.

Based on Gibbs rule and assuming that the sodium and chloride concentrations
are approximately equal, we consider four balance laws, i.e., total carbon,
hydrogen$\displaystyle-$oxygen, chloride and decane. Each one consists of
three terms, i.e., accumulation, convection and the combination of molecular
and capillary diffusion. In this way, we obtain analytically the Riemann
solution, which consists basically in applying the method of characteristics
(MOC) through the wave curve method. In this work we use this method to seek
analytical solutions of the SC$\displaystyle-$CLICDOW model similar to the
model treated in [18, 66, 16, 45, 59, 21, 37, 1]. The Riemann solution
consists of a concatenation of spreading and shock waves, implementing certain
admissibility conditions ([44, 55, 29, 46, 47]).

A Riemann solver for the proposed geochemical model is developed to quantify
the geochemical processes of water injection with CO2 and low salinity in a
carbonated reservoir. We also take into account the surface complexes
formation as a mechanism to change wettability.

A Riemann solver is developed to automate the construction of solution paths.
To do so, we take into account the bifurcation structures, which are not part
of the classical fractional flow method used by [59]. We also provide
comparisons with numerical solutions obtained by means of a commercial program
(COMSOL). This procedure is extremely useful because it allows to include in a
unified manner the geochemistry, the equilibrium reactions and the charge
balances. Also, the method serves to study different situations and to include
more chemical species in the system. Using the extended Gibbs rule (see Eq.
(1)), we reduce the mathematical complexity associated with considering the
large number of physical constraints and parameters that are included in the
geochemistry program PHREEQC.

We use a similar methodology as developed in [6, 1] where the extended Gibbs
phase rule was used to focus on the principal chemical species and to
incorporate the geochemistry of the oil recovery. This circumvents
simultaneous solution of the transport equations, of the equilibrium relations
and of the effect of low salt concentrations.

The aims on this paper are (1) to quantify how the presence of surface
complexes affects wettability (2) to understand how changes in the salt
concentration affects relative permeability and as a consequence of the
injectivity (3) to analyze the wave structure of the solutions (e.g., the
occurrence of a pH wave embedded in a constant pH flood), front salt formation
and a jump in water saturation. The Riemann solution confirms that pH
variations occur with various numerical schemes, inclusive discontinuous
Galerkin, which is expected when surface complexes are formed.

The additional advantage of the Riemann solution is that it can be used to
perform a bifurcation analysis, and to make an inventory of the possible
qualitatively different solutions. The bifurcations occur at coincidence and
inflection loci. Clearly, bifurcations are essential to build the analytical
solution, as well as to determine the location where qualitative changes of
the behavior of the solution are expected. This determination is a useful tool
for mathematical modeling in oil recovery.

By means of numerical and analytical methods we aim in this paper at
quantifying the recovery improvement when carbonated water at low salinity is
injected in a reservoir that contains carbonated brine in equilibrium with an
oleic phase and carbon dioxide.

For the fulfillment of the objectives our paper is organized as follows.
Section II gives the physical model and the equilibrium equations considered
in this study. In Section III the derivation of the mass balance equations for
master species are presented. Moreover, Corey parameters of the flow functions
depending on salt concentration is presented. In Section IV a summary of the
surface complexation model is presented together with values of the parameters
considered in this study. Besides, regression formulas for parameters of the
system of conservation laws depending on pH and chloride concentration are
presented. Furthermore, a procedure for the calculation of wettability
depending on TBP is described. Section V describes the Riemann solver and the
strategy to obtain the Riemann solution. Section VI gives the results in terms
of the $\displaystyle pH$, the chloride concentration, the water saturation
and the total velocity. The calculation suggests that a low salinity
carbonated water flood improves the recovery because it admits a high
dissolved concentration of carbon dioxide.

##  2 Physical model

We consider the injection of low salinity brine (0.5 $\displaystyle mol/liter$
NaCl, saturated with CO2 at a $\displaystyle pH=4.0$) into an inert rock
filled with an oleic phase. Injection and initial fluids contain carbon
dioxide and other unrelated ions such as sodium chloride. We assume chemical
equilibrium in both the aqueous phase and the oleic phase [1]. The ions and
water are only present in the aqueous phase, decane is only present in the
oleic phase. Indeed, dissolution of oil in the aqueous phase is disregarded.
The solubility of carbon dioxide decreases dramatically at high salt
concentration. We assume that the flow is governed by Darcy’s law. The
temperature is chosen to be 39oC because literature data is available and the
pressure is chosen to be well above the pressure at which a gaseous phase can
exist.

We apply the extended Gibbs phase rule to determine the number of degrees of
freedom $\displaystyle n_{f}$. This rule states (see, e.g., [34, 53]) that
this number is given by

$n_{f}=N_{s}+n_{s}-N_{r}-n_{c}+2-p,$ (1)

where $\displaystyle N_{s}$ is the number of dissolved chemical species,
$\displaystyle n_{s}$ is the number of surface species, $\displaystyle N_{r}$
is the number of chemical reactions and $\displaystyle n_{c}$ is the number of
constraints, e.g., the charge balance. The number $\displaystyle 2$ represents
the temperature and pressure and $\displaystyle p$ the number of phases.

We follow Appelo and Parkhurst [58] and [7] and obtain with the geochemistry
program PHREEQC, when we add water, $\displaystyle CaCO_{3}\left(solid\right)$
and $\displaystyle NaCl$, that there are fifteen different chemical species
$\displaystyle\left(N_{s}=15\right)$ with molar concentrations in the aqueous
phase: $\displaystyle c_{a,CO_{2}},$ $\displaystyle c_{a,CO_{3}^{2-}},$
$\displaystyle c_{a,HCO_{3}^{-}},$ $\displaystyle c_{a,CaHCO_{3}^{+}},$
$\displaystyle c_{a,CaCO_{3}},$ $\displaystyle c_{a,NaCO_{3}^{-}},$
$\displaystyle c_{a,NaHCO_{3}},c_{a,H_{2}O},$ $\displaystyle c_{a,H^{+}},$
$\displaystyle c_{a,OH^{-}},$ $\displaystyle
c_{a,CaOH^{+}},c_{a,Ca^{2+}},c_{a,Cl^{-}},c_{a,Na^{+}}$ and the alkane
$\displaystyle\left(A\right)$ concentration in the oleic phase $\displaystyle
c_{o,A}$ . The alkane only occurs in the oleic phase, whereas all the other
components occur only in the aqueous phase. $\displaystyle CaCO_{3}$ occurs
both in the solid phase with concentration $\displaystyle c_{r,CaCO_{3}}$ and
in the aqueous phase with concentration $\displaystyle c_{a,CaCO_{3}}$. In
addition we have $\displaystyle n_{s}=11$ sorbed species $\displaystyle
oil_{s,NH^{+}}$, $\displaystyle oil_{w,COOH}$, $\displaystyle Cal_{s,OH}$,
$\displaystyle Cal_{w,CO_{3}H}$, $\displaystyle oil_{s,N}$, $\displaystyle
oil_{w,COO^{-}}$, $\displaystyle oil_{w,COOCa^{+}}$, $\displaystyle
Cal_{s,OH_{2}^{+}}$, $\displaystyle Cal_{s,CO_{3}^{-}}$, $\displaystyle
Cal_{w,CO_{3}^{-}}$, $\displaystyle Cal_{w,CO_{3}Ca^{+}}.$ There are seven
surface reactions taking into account complexes. These surface reactions are
[12]

$\displaystyle\displaystyle oil_{s,}{}_{{}^{NH^{+}}}\rightleftharpoons
oil_{s,N}+H^{+}$ $\displaystyle\displaystyle\quad\log K_{-11}=-6.0$
$\displaystyle\displaystyle oil_{w,COOH}\rightleftharpoons
oil_{w,COO^{-}}+H^{+}$ $\displaystyle\displaystyle\quad\log K_{-21}=-5.0$
$\displaystyle\displaystyle oil_{w,COOH}+Ca^{2+}\rightleftharpoons
oil_{w,COOCa^{+}}+H^{+}$ $\displaystyle\displaystyle\quad\log K_{22}=-3.8$
$\displaystyle\displaystyle Cal_{s,OH}+H^{+}\rightleftharpoons
Cal_{s,OH_{2}^{+}}$ $\displaystyle\displaystyle\quad\log K_{13}=11.8$
$\displaystyle\displaystyle Cal_{s,OH}+HCO_{3}^{-}\rightleftharpoons
Cal_{s,CO_{3}^{-}}+H_{2}O$ $\displaystyle\displaystyle\quad\log K_{33}=5.8$
$\displaystyle\displaystyle Cal_{w,CO_{3}H}\rightleftharpoons
Cal_{w,CO_{3}^{-}}+H^{+}$ $\displaystyle\displaystyle\quad\log K_{-14}=-5.1$
$\displaystyle\displaystyle Cal_{w,CO_{3}H}+Ca^{2+}\rightleftharpoons
Cal_{w,CO_{3}Ca^{+}}+H^{+}$ $\displaystyle\displaystyle\quad\log K_{24}=-2.6$
(2)

We use four different sorption sites. Each of the sorption sites can receive
sorption molecules of two types, called $\displaystyle weak$ and
$\displaystyle strong$. For instance the strong adsorption sites in oil can
receive both $\displaystyle oil_{s,NH^{+}}$ and $\displaystyle oil_{s,N}$, so
that the sum of the sorbed concentration at each sorption site is a given
constant. A similar situation occurs for the weak adsorption sites on oil
$\displaystyle\left(oil_{s}\right),$ but now with three surface species. The
concentration of $\displaystyle CaCO_{3}$ in the solid phase is known and
constant. Thus taking into account the charge balance equation we have
$\displaystyle n_{c}=5$ constraints.

We notice that the species can be derived from seven master species, i.e.,
$\displaystyle C\left(4\right),$ $\displaystyle H\left(1\right),$
$\displaystyle O\left(-2\right),$ $\displaystyle Ca(2),$ $\displaystyle
Na(1),$ $\displaystyle Cl\left(-1\right)$, $\displaystyle C\left(-4\right)$,
where the value in parentheses denotes the constant valence of the species. We
note that $\displaystyle C\left(-4\right)$ is used to denote organic carbon as
opposed to $\displaystyle C\left(4\right),$ which denotes inorganic carbon.

We consider the following eight equilibrium reactions in the aqueous phase
(thus we have $\displaystyle N_{r}=15$)

$\displaystyle\left(CO_{2}\right)_{aq}+H_{2}O\rightleftharpoons
HCO_{3}^{-}+H^{+}$ $\displaystyle HCO_{3}^{-}\rightleftharpoons
CO_{3}^{2-}+H^{+}$ $\displaystyle H_{2}O\rightleftharpoons OH^{-}+H^{+}$
$\displaystyle\left(CaCO_{3}\right)_{aq}\rightleftharpoons
Ca^{2+}+CO_{3}^{2-}$ $\displaystyle Ca^{2+}+H_{2}O\rightleftharpoons
CaOH^{+}+H^{+}$ $\displaystyle CO_{3}^{2-}+Ca^{2+}+H^{+}\rightleftharpoons
CaHCO_{3}^{+}$ $\displaystyle Na^{+}+CO_{3}^{2-}\rightleftharpoons
NaCO_{3}^{-}$ $\displaystyle Na^{+}+HCO_{3}^{-}=NaHCO_{3}$

We dropped the subscript $\displaystyle\left(aq\right)$ on all compounds
except $\displaystyle CaCO_{3}$ and $\displaystyle CO_{2}$ as we assume that
they only occur in the aqueous phase. All possible equilibrium equations can
be derived from these eight equilibrium reactions.

Thermodynamic equilibrium requires that the chemical potential of
$\displaystyle\left(CaCO_{3}\right)_{r}$ in the solid phase is equal to the
chemical potential of $\displaystyle\left(CaCO_{3}\right)_{aq}$ in the aqueous
phase. This can be represented by

$\left(CaCO_{3}\right)_{r}\rightleftharpoons\left(CaCO_{3}\right)_{aq}.$ (3)

In the same way the chemical potential of carbon dioxide in the aqueous phase
is equal to the chemical potential in the oleic phase. This can be represented
as

$\left(CO_{2}\right)_{o}\rightleftharpoons\left(CO_{2}\right)_{aq}.$ (4)

As we consider solid, aqueous and oleic phases, the number of phases
$\displaystyle p$ is $\displaystyle 3$. We have the charge balance equation
given by

$\left(\begin{array}[c]{c}2c_{a,CO_{3}}+oil_{w,COO^{-}}+c_{a,HCO_{3}}+c_{a,OH}+c_{a,NaCO3}+c_{a,Cl}+Cal_{s,CO_{3}^{-}}+Cal_{w,CO_{3}^{-}}\\\
=2c_{a,Ca}+oil_{s,}{}_{{}^{NH^{+}}}+c_{a,H}+c_{a,Ca(HCO_{3})}\\\
+c_{a,Na}+c_{a,CaOH}+oil_{w,COOCa^{+}}+Cal_{s,OH_{2}^{+}}+Cal_{w,CO_{3}Ca^{+}}\end{array}\right).$

This charge balance equation can also be derived from the mass balance
equations. Thus, we adopt the charge balance equation allowing one mass
balance equation to be removed.

Following Gibbs rule described in (1) the number $\displaystyle n_{f}$ of
degrees of freedom is

$n_{f}=N_{s}+n_{s}-N_{r}-n_{c}-p+2=15+11-15-5-3+2=5.$ (5)

In Eq. 5 we use that $\displaystyle N_{s}+n_{s}-N_{r}-n_{c}=6$. Given the
temperature and pressure, we only need three concentrations to specify the
composition of the three phase system.

For the independent concentrations we choose the hydrogen ion $\displaystyle
c_{a,H^{+}}c,$ concentration the chloride ion $\displaystyle c_{a,Cl^{-}}$
concentration and the sodium ion concentration $\displaystyle c_{a,Na^{+}}$.
Another assumption consists in taking the chloride and sodium ion
concentration as equal. This assumption is to simplify our model.

##  3 Mass balance equations for porous medium flow with surface species

All compounds that are used in the model are built with the atoms
$\displaystyle C\left(4\right),$ $\displaystyle H\left(1\right),$
$\displaystyle O\left(-2\right),$ $\displaystyle Ca\left(2\right)$,
$\displaystyle Na\left(1\right)$, $\displaystyle Cl\left(-1\right)$ and
$\displaystyle C\left(-4\right)$. For instance carbon dioxide consists of one
atom of carbon with valence four and two oxygen atoms with valence minus two.
The uncompensated valences result in the charge of the ion, e.g.,
$\displaystyle CO_{3}^{-2}$ consists of one atom $\displaystyle
C\left(4\right)$ and three oxygen atoms $\displaystyle O\left(-2\right).$
Hence $\displaystyle CO_{3}^{2-}$ has a charge of minus 2. Organic carbon is
denoted by $\displaystyle C\left(-4\right)$, leading to for instance
$\displaystyle CH_{4},$ where the four valent carbon is compensated by four
monovalent hydrogens. Two carbons may combine to $\displaystyle
C_{2}\left(-6\right)$, and $\displaystyle n$ carbons to $\displaystyle
C_{n}\left(-2n-2\right),$ but we denote all organic carbon by $\displaystyle
C\left(-4\right).$

The concentrations of all species are expressed in term of activities and
activity coefficients, where the activity coefficients $\displaystyle\gamma$
depend only on the ionic strength $\displaystyle\mu.$

###  3.1 Master species in terms of activities in the aqueous phase

For the derivation of the mass balance equation we first derive the equations
for the dissolved master species, i.e., $\displaystyle C\left(4\right)$,
$\displaystyle H\left(1\right)$, $\displaystyle O\left(-2\right),$
$\displaystyle Ca\left(2\right),$ $\displaystyle Cl\left(-1\right),$
$\displaystyle Na\left(1\right)$ concentrations in the aqueous solution, which
do not include the surface master species.

For convenience we first define the total aqueous inorganic carbon (i.e.,
without hydrocarbon) concentration $\displaystyle C_{a,C\left(4\right)}$ by

$\displaystyle\displaystyle C_{a,C\left(4\right)}$
$\displaystyle\displaystyle:=\frac{a_{a,CO_{2}}}{\gamma_{a,CO_{2}}\left(\mu\right)}+\frac{a_{a,CO_{3}^{2-}}}{\gamma_{a,CO_{3}^{2-}}\left(\mu\right)}+\frac{a_{a,HCO_{3}^{-}}}{\gamma_{a,HCO_{3}^{-}}\left(\mu\right)}+\frac{a_{a,CaHCO_{3}^{+}}}{\gamma_{a,CaHCO_{3}^{+}}\left(\mu\right)}$
$\displaystyle\displaystyle+\frac{a_{a,CaCO_{3}}}{\gamma_{a,CaCO_{3}}\left(\mu\right)}+\frac{a_{a,NaCO_{3}^{-}}}{\gamma_{a,NaCO_{3}^{-}}\left(\mu\right)}+\frac{a_{a,NaHCO_{3}}}{\gamma_{a,NaHCO_{3}}\left(\mu\right)}.$
(6)

In the same way we define the total aqueous hydrogen concentration by

$C_{a,H\left(1\right)}:=\left(\begin{array}[c]{c}2\frac{a_{a,H_{2}O}}{\gamma_{a,H_{2}O}}+\frac{a_{a,H^{+}}}{\gamma_{a,H^{+}}\left(\mu\right)}+\frac{a_{a,OH^{-}}}{\gamma_{a,OH^{-}}\left(\mu\right)}\\\
+\frac{a_{a,HCO3^{-}}}{\gamma_{a,OH^{-}}\left(\mu\right)}+\frac{a_{a,CaHCO_{3}^{+}}}{\gamma_{a,CaHCO_{3}^{+}}\left(\mu\right)}+\frac{a_{a,CaOH^{+}}}{\gamma_{a,CaHCO_{3}^{+}}\left(\mu\right)}+\\\
\frac{a_{a,NaHCO_{3}}}{\gamma_{a,CaHCO_{3}^{+}}\left(\mu\right)}\end{array}\right):=2\frac{a_{a,H_{2}O}}{\gamma_{a,H_{2}O}}+\delta
C_{a,H\left(1\right)}.$ (7)

where the concentration of water $\displaystyle
c_{a,H_{2}O}=a_{a,H_{2}O}/\gamma_{a,H_{2}O}$ is much larger than the other
concentrations.

The total aqueous oxygen concentration can be written as

$C_{a,O\left(-2\right)}:=\left(\begin{array}[c]{c}\frac{a_{a,H_{2}O}}{\gamma_{a,H_{2}O}}+\frac{3a_{a,CO_{3}^{2-}}}{\gamma_{a,CO_{3}^{2-}}\left(\mu\right)}+\frac{a_{a,OH^{-}}}{\gamma_{a,OH^{-}}\left(\mu\right)}+\frac{2a_{a,CO2}}{\gamma_{a,CO2}\left(\mu\right)}\\\
+\frac{3a_{a,HCO3^{-}}}{\gamma_{a,HCO3^{-}}\left(\mu\right)}+\frac{3a_{a,CaHCO_{3}^{+}}}{\gamma_{a,CaHCO_{3}^{+}}\left(\mu\right)}+\frac{a_{a,CaOH^{+}}}{\gamma_{a,CaOH^{+}}\left(\mu\right)}+\\\
+\frac{3a_{a,CaCO_{3}}}{\gamma_{a,CaCO_{3}}\left(\mu\right)}+\frac{3a_{a,NaCO_{3}^{-}}}{\gamma_{a,NaCO_{3}^{-}}\left(\mu\right)}+\frac{3a_{a,NaHCO_{3}}}{\gamma_{a,NaHCO_{3}}\left(\mu\right)}\end{array}\right):=\frac{a_{a,H_{2}O}}{\gamma_{a,H_{2}O}}+\delta
C_{a,O\left(-2\right)}.$ (8)

The total aqueous calcium concentration is

$C_{a,Ca\left(2\right)}:=\frac{a_{a,Ca^{2+}}}{\gamma_{a,Ca^{2+}}\left(\mu\right)}+\frac{a_{a,CaHCO_{3}^{+}}}{\gamma_{a,CaHCO_{3}^{+}}\left(\mu\right)}+\frac{a_{a,CaOH^{+}}}{\gamma_{a,CaOH^{+}}\left(\mu\right)}+\frac{a_{a,CaCO_{3}}}{\gamma_{a,CaCO_{3}}\left(\mu\right)}+\beta_{Ca}/2.$
(9)

where $\displaystyle\beta_{Ca}$ is the equivalent fraction of calcium (see
Definition in [7]). The total aqueous sodium concentration is

$C_{a,Na\left(1\right)}:=\frac{a_{a,Na^{+}}}{\gamma_{a,Na^{+}}\left(\mu\right)}+\frac{a_{a,NaCO_{3}^{-}}}{\gamma_{a,NaCO_{3}^{-}}\left(\mu\right)}+\frac{a_{a,NaHCO_{3}}}{\gamma_{a,NaHCO_{3}}\left(\mu\right)}+\beta_{Na}.$
(10)

where $\displaystyle\beta_{Na}$ is the equivalent fraction of sodium.

The total aqueous chloride concentration is

$C_{a,Cl\left(-1\right)}:=\frac{a_{a,Cl^{-}}}{\gamma_{a,Cl^{-}}\left(\mu\right)}.$
(11)

###  3.2 Multiphase mass-balance equations with surface complexes

#### Carbon balance

We can write the mass balance equation for carbon as

$\displaystyle\displaystyle\partial_{t}\left(\varphi
C_{a,C\left(4\right)}S_{w}+\varphi
c_{o,CO_{2}}S_{o}+\left(1-\varphi\right)\left(c_{r,CaCO_{3}}+Cal_{s,CO_{3}^{-}}\right)\right)$
$\displaystyle\displaystyle+\partial_{x}u\left(C_{a,C\left(4\right)}f_{w}+c_{o,CO_{2}}f_{o}\right)$
$\displaystyle\displaystyle=\partial_{x}\left(\varphi\left(D_{w}S_{w}\partial_{x}C_{a,C\left(4\right)}+D_{o}S_{o}\partial_{x}c_{o,CO_{2}}\right)\right)+$
$\displaystyle\displaystyle\partial_{x}\left(\mathcal{D}C_{a,C\left(4\right)}\partial_{x}S_{w}\right)+\partial_{x}\left(\mathcal{D}c_{o,CO_{2}}\partial_{x}S_{o}\right),$
(12)

where $\displaystyle c_{o,CO_{2}}$ is the concentration of CO2 in the oil
phase. The surface master species are arbitrarily chosen to be $\displaystyle
oil_{s,N},$ $\displaystyle oil_{w,COO},Cal_{s,O},$ $\displaystyle
Cal_{w,CO_{3}}$. This choice uses that the sum of the species derived from the
master species is equal to the number of active sites, which is considered to
be constant for two sets of adsorbed sites on the oil surface and two sets
adsorbed onto the calcite surface. Note that $\displaystyle S_{o}=1-S_{w}.$
The equilibrium conditions of the sorbed species are taken from [12].

#### Hydrogen balance

In the same way we find for the total hydrogen balance

$\displaystyle\displaystyle\partial_{t}\left(\varphi
C_{a,H\left(1\right)}S_{w}\right)+\partial_{t}\varphi\left(oil_{s,}{}_{{}^{N-H^{+}}}+oil_{w,COO-H}\right)S_{o}$
$\displaystyle\displaystyle+\partial_{t}\left(1-\varphi\right)\left(2Cal_{s,O-2H^{+}}+Cal_{s,O-H}+Cal_{w,CO_{3}^{-}-H}\right)$
$\displaystyle\displaystyle+\partial_{x}\left(uC_{a,H\left(1\right)}f_{w}\right)+\partial_{x}u\left(oil_{s,}{}_{{}^{N-H^{+}}}+oil_{w,COO-H}\right)f_{o}=$
$\displaystyle\displaystyle\partial_{x}\left(\mathcal{D}C_{a,H\left(1\right)}\partial_{x}S_{w}\right)+\partial_{x}\left(\mathcal{D}\left(oil_{s,}{}_{{}^{N-H^{+}}}+oil_{w,COO-H}\right)\partial_{x}S_{o}\right)$
$\displaystyle\displaystyle+\partial_{x}\left(\varphi
D_{w}S_{w}\partial_{x}C_{a,H\left(1\right)}\right)+\partial_{x}\left(\varphi
D_{o}S_{o}\partial_{x}\left(oil_{s,}{}_{{}^{N-H^{+}}}+oil_{w,COO-H}\right)\right),$
(13)

where the surface master species are chosen to be $\displaystyle oil_{s,N-},$
$\displaystyle oil_{w,COO-},Cal_{s,O-},$ $\displaystyle Cal_{w,CO_{3}^{-}-}.$

#### Oxygen balance

In the same way we find for the total oxygen balance

$\displaystyle\displaystyle\partial_{t}\left(\varphi\left(C_{a,O\left(-2\right)}S_{w}+2c_{o,CO_{2}}S_{o}\right)+\left(1-\varphi\right)\left(3c_{r,CaCO_{3}}+2\left(Cal_{s,CO_{3}^{-}}\right)\right)\right)$
$\displaystyle\displaystyle+\partial_{x}u\left(C_{a,O\left(-2\right)}f_{w}+2c_{o,CO_{2}}f_{o}\right)=\partial_{x}\left(\mathcal{D}C_{a,O\left(-2\right)}\partial_{x}S_{w}\right)+2\partial_{x}\left(\mathcal{D}c_{o,CO_{2}}S_{o}\partial_{x}S_{o}\right)$
$\displaystyle\displaystyle+\partial_{x}\varphi\left(D_{w}S_{w}\partial_{x}C_{a,O\left(-2\right)}+2D_{o}S_{o}\partial_{x}c_{o,CO_{2}}\right).$
(14)

#### Calcium balance

For the total calcium we find

$\scriptstyle\displaystyle\mathtt{\partial_{t}\left(\varphi
S_{w}C_{a,Ca\left(2\right)}+\varphi S_{o}\left(oil_{w,COO-
Ca^{+}}\right)+\left(1-\varphi\right)\left(c_{r,CaCO_{3}}+Cal_{w,CO_{3}-Ca^{+}}\right)\right)+\partial_{t}\beta}_{Ca}$
$\scriptstyle\displaystyle\partial_{x}u\left(C_{a,Ca\left(2\right)}f_{w}\right)+\partial_{x}u\left(oil_{w,COO-
Ca^{+}}\right)f_{o}=\partial_{x}\left(\mathcal{D}C_{a,Ca\left(2\right)}\partial_{x}S_{w}\right)$
$\scriptstyle\displaystyle+\partial_{x}\left(\mathcal{D}oil_{w,COO-
Ca^{+}}\partial_{x}S_{o}\right)+\partial_{x}\left(\varphi
D_{w}S_{w}\partial_{x}C_{a,Ca\left(2\right)}\right)+\partial_{x}\left(\varphi
D_{o}S_{o}\partial_{x}oil_{w,COO-Ca^{+}}\right).$ (15)

#### Sodium, Chlorine and total oil equation

For the sodium equation we obtain

$\partial_{t}\left(\varphi
S_{w}C_{a,Na\left(1\right)}\right)+\partial_{t}\beta_{Na}+\partial_{x}\left(uC_{a,Na\left(1\right)}f_{w}\right)=\partial_{x}\left(\mathcal{D}C_{a,Na\left(1\right)}\partial_{x}S_{w}\right)+\partial_{x}\left(\varphi
D_{w}S_{w}\partial_{x}C_{a,Na\left(1\right)}\right).$ (16)

In the same way we obtain for the first order terms of the chlorine equation

$\partial_{t}\left(\varphi
S_{w}C_{a,Cl\left(-1\right)}\right)+\partial_{x}\left(uC_{a,Cl\left(-1\right)}f_{w}\right)=\partial_{x}\left(\mathcal{D}C_{a,Cl\left(-1\right)}\partial_{x}S_{w}\right)+\partial_{x}\left(\varphi
D_{w}S_{w}\partial_{x}C_{a,Cl\left(-1\right)}\right).$ (17)

For the total oil, e.g., heptane we retain the first order terms

$\partial_{t}\left(\varphi
S_{o}c_{o,C\left(-4\right)}\right)+\partial_{x}\left(uc_{o,C\left(-4\right)}f_{o}\right)=\partial_{x}\left(\mathcal{D}c_{o,C\left(-4\right)}\partial_{x}S_{o}\right)+\partial_{x}\left(\varphi
D_{o}S_{o}\partial_{x}c_{o,C\left(-4\right)}\right).$ (18)

The oil concentration can be obtained from the dissolved carbon dioxide
concentration $\displaystyle c_{o,CO_{2}}$ (proportional to $\displaystyle
c_{a,CO_{2}}$) with the EOS

$\frac{c_{o,C\left(-4\right)}}{c_{O,C\left(-4\right)}}+\frac{c_{o,CO_{2}}}{c_{O,CO2}}=1.$
(19)

###  3.3 System of conservation laws

Based on Gibbs rule and assuming that sodium and chloride concentration are
similar we consider four conservation laws, namely for total carbon,
hydrogen$\displaystyle-$oxygen, chloride and decane. Each one consists of four
terms, i.e. accumulation, convection, molecular diffusion and capillary
diffusion. We neglect diffusion and capillarity effects.

We can write the mass balance equation for carbon as

$\displaystyle\displaystyle\partial_{t}\left(\varphi
C_{a,C\left(4\right)}S_{w}+\varphi
c_{o,CO_{2}}S_{o}+\left(1-\varphi\right)\left(Cal_{s,CO_{3}^{-}}\right)\right)$
$\displaystyle\displaystyle+\partial_{x}\left(u\left(C_{a,C\left(4\right)}f_{w}+c_{o,CO_{2}}f_{o}\right)\right)=0.$
(20)

For the total oil, e.g., heptane we obtain

$\partial_{t}\left(\varphi
S_{o}c_{o,C\left(-4\right)}\right)+\partial_{x}\left(uc_{o,C\left(-4\right)}f_{o}\right)=0.$
(21)

In the same way we obtain the chlorine equation

$\partial_{t}\varphi
S_{w}C_{a,Cl\left(-1\right)}+\partial_{x}\left(uC_{a,Cl\left(-1\right)}f_{w}\right)=0.$
(22)

To remove one balance equation, we combine the hydrogen and oxygen balance
equations to a single equation, in such a way that the water concentration is
eliminated. We do so because the water concentration is much higher than the
other concentrations.

We subtract twice the oxygen equation from the hydrogen equation and obtain
after substitution of Eq. (13) and Eq. (14)

$\displaystyle\displaystyle\partial_{t}\left(2\delta
C_{a,O\left(-2\right)}-\delta
C_{a,H\left(1\right)}\right)S_{w}-\partial_{t}\varphi\left(oil_{s,}{}_{{}^{N-H^{+}}}+oil_{w,COO-H}\right)S_{o}$
$\displaystyle\displaystyle-\partial_{t}\left(1-\varphi\right)\left(2Cal_{s,O-2H^{+}}+Cal_{s,O-H}+Cal_{w,CO_{3}^{-}-H}\right)$
$\displaystyle\displaystyle
2\partial_{t}\left(\varphi\left(2c_{o,CO_{2}}S_{o}\right)+\left(1-\varphi\right)\left(3c_{r,CaCO_{3}}+2\left(Cal_{s,CO_{3}^{-}}\right)\right)\right)$
$\displaystyle\displaystyle+\partial_{x}u\left(4c_{o,CO_{2}}-oil_{s,}{}_{{}^{N-H^{+}}}-oil_{w,COO-H}\right)f_{o}+\partial_{x}u\left(\left(2\delta
C_{a,O\left(-2\right)}-\delta C_{a,H\left(1\right)}\right)f_{w}\right)=0.$
(23)

In more compact way the system of conservation laws (20) and (23) can be
written

$\displaystyle\displaystyle\partial_{t}\left(\varphi\rho_{w1}S_{w}+\varphi\rho_{o1}S_{o}+\rho_{r1}\right)+\partial_{x}\left(u\left(\rho_{w1}f_{w}+\rho_{o1}f_{o}\right)\right)=0,$
(24)
$\displaystyle\displaystyle\partial_{t}\left(\varphi\rho_{o2}S_{o}\right)+\partial_{x}\left(u\left(\rho_{o2}f_{o}\right)\right)=0,$
(25)
$\displaystyle\displaystyle\partial_{t}\left(\varphi\rho_{w3}S_{w}\right)+\partial_{x}\left(u\left(\rho_{w3}f_{w}\right)\right)=0,$
(26)
$\displaystyle\displaystyle\partial_{t}\left(\varphi\rho_{w4}S_{w}+\varphi\rho_{o4}S_{o}+\rho_{r4}\right)+\partial_{x}\left(u\left(\rho_{w4}f_{w}+\rho_{o4}f_{o}\right)\right)=0,$
(27)

where the coefficient functions are defined by

$\displaystyle\displaystyle\rho_{w1}=C_{a,C\left(4\right)},~{}\rho_{o1}=c_{o,CO_{2}},~{}\rho_{r1}=(1-\varphi)\left(Cal_{s,CO_{3}^{-}}\right),$
(28)
$\displaystyle\displaystyle\rho_{o2}=C_{a,C\left(-4\right)},~{}~{}\rho_{w3}=C_{a,Cl},$
(29) $\displaystyle\displaystyle\rho_{w4}=\left(2\delta
C_{a,O\left(-2\right)}-\delta
C_{a,H\left(1\right)}\right),~{}\rho_{o4}=\left(4c_{o,CO_{2}}-\left(oil_{s,}{}_{{}^{N-H^{+}}}+oil_{w,COO-H}\right)\right),$
(30)
$\displaystyle\displaystyle\rho_{r4}=(1-\varphi)\left(4Cal_{s,CO_{3}^{-}}-\left(2Cal_{s,O-2H^{+}}+Cal_{s,O-H}+Cal_{w,CO_{3}^{-}-H}\right)\right).$
(31)

where $\displaystyle\rho_{wi}$, $\displaystyle\rho_{oi}$ and
$\displaystyle\rho_{ri}$ depend on $\displaystyle pH$ and $\displaystyle[Cl]$.
It is possible to verify that when the chemical surface complexes are not
present then system (24)-(27) is reduced to the system studied in [6].

###  3.4 Fractional flow

The fractional flows for water and oil are saturation-dependent functions
defined as follows. We denote $\displaystyle
S_{e}(S_{w})=(S_{w}-S_{wr})/(1-S_{wr}-S_{or}),~{}\text{for}~{}S_{w}\geq
S_{wr}~{}\text{and}~{}S_{e}=0~{}\text{for}$ $\displaystyle
S_{w}<S_{wr};k_{rw}(S_{w})=k_{w}(S_{e}(S_{w}))^{n_{w}}~{}~{}\text{ and
}~{}~{}k_{ro}(S_{w})=k_{o}(1-S_{e}(S_{w}))^{n_{o}}$ ($\displaystyle
S_{o}=1-S_{w}$) ([25]). Here $\displaystyle n_{w}$ and $\displaystyle n_{o}$
are the Corey exponents. The parameters $\displaystyle k_{w}$ and
$\displaystyle k_{o}$ are end point relative permeabilities of water and oil
phases, respectively.

The water viscosity is taken as $\displaystyle\mu_{w}=0.001$ and the oil
viscosity as $\displaystyle\mu_{o}=0.002$ when they are constant; then we have
the fractional flow functions for the aqueous and oleic phases

$f_{w}(S_{w})=\frac{k_{rw}(S_{w})/\mu_{w}}{(k_{rw}(S_{w})/\mu_{w}+k_{ro}(1-S_{w})/\mu_{o})},\quad\text{and}\quad
f_{o}(S_{w})=1-f_{w}(S_{w}),$ (32)

where the water and oil permeabilities $\displaystyle k_{rw}(S_{w})$ and
$\displaystyle k_{ro}(S_{o})$ are expressed in terms of their saturations;
$\displaystyle\mu_{w}$ and $\displaystyle\mu_{o}$ are the viscosities of the
aqueous and oleic phases. We disregard capillarity and diffusive effects.

In this work we consider the relation between relative permeability and low
salinity from the ideas described in [35]. Such a relation is based on the
definition of the weight parameter $\displaystyle\theta$ to modify rock and
fluid properties regarding salinity level. This model considers salt as an
aqueous component that can be transported and traced. In this relation the
relative permeability functions depend on water saturation $\displaystyle
S_{w}$ and salinity. This relationship is defined by the upper ($\displaystyle
HL$) and lower ($\displaystyle LS$) limits for salinity called high and low
salinity limits, respectively. To do that Corey’s coefficients are adjusted
for low and high salt concentration.

The weighting function $\displaystyle\theta$ is introduced by a linear
relationship among two values for residual oil saturation corresponding to
high and low salinity limits

$\theta=\frac{S_{or}-S_{orw}^{LS}}{S_{orw}^{HL}-S_{orw}^{LS}},$ (33)

where the parameter $\displaystyle\theta$ is used for interpolating between
preset high and low salinity curves for relative permeability and capillary
pressure

$K_{wi}(S_{w})=\theta K_{ri}^{HL}(S_{w})+(1-\theta)K_{ri}^{LS}(S_{w}),$ (34)

where $\displaystyle i=w,o$ for water and oil permeability, respectively. For
high salinity regime we take $\displaystyle k_{w}=0.25$, $\displaystyle
k_{o}=0.5$, $\displaystyle n_{w}=3$ and $\displaystyle n_{o}=2$, while for low
salinity the values $\displaystyle k_{w}=0.41$, $\displaystyle k_{o}=0.5$,
$\displaystyle n_{w}=3$ and $\displaystyle n_{o}=2$ are used. Furthermore, in
this study, the upper limit used for salinity corresponds to maximum residual
oil saturation ($\displaystyle S_{or}^{H}=0.35$). The lower limit corresponds
to minimum residual oil saturation ($\displaystyle S_{or}^{L}=0.15$).

##  4 Surface complexation modeling

In this study, we use SCM to predict the wettability alteration of minerals
through TBP. Details of the procedure employed here can be found in [28, 51].
The geochemistry solver PHREEQC was utilized to calculate the concentrations
of the complexes. This program uses as input the properties of formation water
and oil components of the crude oil. The properties of the minerals for
several formation waters and for two polar oil components in the crude oil
acid (e.g., COOH and NH+) are used here. We take the input data described in
[24]. Input properties of Oil are depicted in Table 1 and 2.

Oil | Density (g/cm³) at 20oC | TAN (mg KOH/g oil) | TBN (mg KOH/g oil)
---|---|---|---
1 | 0.86 | 0.1 | 1.9
2 | 0.9 | 0.38 | 2.3
Table 1: Properties of the minerals and the polar oil components.

The total acidic number (TAN) and the total base number (TBN) are used to
calculate the oil site density using the formulas (see [10])

$N_{S,COOH}=0.602\times 10^{6}\frac{TAN}{1000~{}a_{oil}~{}MW_{KOH}},$ (35)

$N_{S,NH^{+}}=0.602\times 10^{6}\frac{TBN}{1000~{}a_{oil}~{}MW_{KOH}},$ (36)

where $\displaystyle MW_{KOH}=56.1g/mol$ molecular weight of potassium
hydroxide. Here $\displaystyle a_{oil}$ [$\displaystyle m^{2}/g$] denotes the
specific area of oil, which we assume to be the same as its respective
carbonate minerals in aqueous solutions taken from [67].

Surface | Site Density (site/$\displaystyle nm^{2}$) | Surface Area ($\displaystyle m^{2}$/g) | Mass (g)
---|---|---|---
Calcite | 4.9 | 2.0 | 0.2
Oil 1=COOH | 0.54 | 2.0 | 2.59
Oil 1=NH+ | 10.2 | 2.0 | 2.59
Oil 2=COOH | 2.04 | 2.0 | 2.71
Oil 2=NH+ | 12.34 | 2.0 | 2.71
Table 2: Surface density and area of the two types of oil.

Input data sets corresponding to salt are chosen for several formation waters
where the ion concentrations of $\displaystyle Na^{2+}$ and $\displaystyle
Cl^{-}$ varies. The data are taken to yield the coefficients of system
(24)-(27), which depend on the ion concentration chosen in this work to
represent the variability of state space, i.e., $\displaystyle pH$ and
$\displaystyle Cl^{-}$. Since we assume that the ion concentrations of
$\displaystyle Na^{2+}$ and $\displaystyle Cl^{-}$ are similar in the
simulations, these ion concentrations vary from 30 to 1290 mmol/kgw. In turn,
the ion concentrations of $\displaystyle Ca^{2+}$, $\displaystyle Mg^{2+}$ and
$\displaystyle SO_{4}^{2-}$ are taken as constant. The $\displaystyle pH$ is
taken to vary from 2.7 to 8. Carbon is also taken as a constant in this study.
The data are summarized in Table 3.

Ion | Synthetic injected water( mmol/kgw)
---|---
$\displaystyle Na^{+}$ | 30-1290
$\displaystyle Mg^{2+}$ | 40
$\displaystyle Ca^{2+}$ | 20
$\displaystyle Cl^{-}$ | 30-1290
$\displaystyle SO_{4}^{2-}$ | 20
Table 3: Formation waters where chloride and sodium ion concentrations vary.
Other minerals are fixed.

###  4.1 Regression formulas Matlab, Eureqa and PHREEQC program

The expressions for the coefficients $\displaystyle\rho_{wi},~{}\rho_{oi}$ and
$\displaystyle\rho_{ri}$ ($\displaystyle i=1,\ldots,4$) in system (24)-(27)
are obtained using the ion concentrations of the complexes with formulas given
in 28-31. These coefficients can also be determined by allowing the other
minerals as $\displaystyle Ca^{2+}$ and $\displaystyle Mg^{2+}$ to vary as
well. In such a case, and in accordance with the generalized Gibbs rule, it is
necessary to add other equations to the model and to increase the number of
degrees of freedom. This useful procedure can be used to study enhanced oil
recovery by Smart Water injection, which we will study in future work with the
methodology explained here.

Processing of the PHREEQC output data is done with MATLAB program, by using
the curve fitting tool.

Figure 1: Coefficients $\displaystyle\rho_{w1}$(left) and
$\displaystyle\rho_{o1}$(right).

Figure 2: Coefficients $\displaystyle\rho_{r1}$(left) and
$\displaystyle\rho_{o2}$(right).

Figure 3: Coefficients $\displaystyle\rho_{w3}$(left) and
$\displaystyle\rho_{w4}$(right).

Figure 4: Coefficients $\displaystyle\rho_{o4}$(left) and
$\displaystyle\rho_{r4}$(right).

Denoting the ion concentration of hydrogen and chloride by $\displaystyle
x=pH$ and $\displaystyle y=[Cl]$ the formula of the coefficients for the case
of Oil $\displaystyle 1$ are

$\displaystyle\displaystyle\rho_{w1}=-8.247\cdot 10^{-7}\,x^{5}+8.841\cdot
10^{-9}\,x^{4}\,y+2.472\cdot 10^{-5}\,x^{4}+9.621\cdot 10^{-9}\,x^{3}\,y^{2}$
$\displaystyle\displaystyle+3.387\cdot 10^{-9}\,x^{3}\,y-2.803\cdot
10^{-4}\,x^{3}+9.781\cdot 10^{-6}\,x^{2}\,y^{3}-2.523\cdot
10^{-5}\,x^{2}\,y^{2}-$ $\displaystyle\displaystyle 2.843\cdot
10^{-5}\,x^{2}\,y+0.0015\,x^{2}-5.102\cdot 10^{-5}\,x\,y^{4}+4.901\cdot
10^{-5}\,x\,y^{3}$ $\displaystyle\displaystyle+8.39\cdot
10^{-5}\,x\,y^{2}+9.35\cdot 10^{-5}\,x\,y-0.0036\,x+1.49\cdot
10^{-4}\,y^{5}-2.383\cdot 10^{-4}\,y^{4}$
$\displaystyle\displaystyle+1.341\cdot 10^{-4}\,y^{3}-1.758\cdot
10^{-4}\,y^{2}-1.601\cdot 10^{-4}\,y+0.007,$ (37)
$\displaystyle\displaystyle\rho_{o1}=-1.575\cdot 10^{-5}x^{4}y+1.01\cdot
10^{-4}x^{4}-1.83\cdot 10^{-5}x^{3}*y^{2}+3.95\cdot 10^{-4}x^{3}y-$
$\displaystyle\displaystyle 0.002x^{3}+3.9010\cdot
10^{-6}x^{2}y^{3}+2.644\cdot 10^{-4}x^{2}y^{2}-0.0034x^{2}y+0.0132x^{2},$
$\displaystyle\displaystyle+1.747\cdot 10^{-5}xy^{4}-1.0440\cdot
10^{-4}xy^{3}-0.0011xy^{2}+0.0115xy-0.0369x+$ $\displaystyle\displaystyle
4.917\cdot 10^{-4}y^{5}-0.0014y^{4}+0.0015y^{3}+7.16\cdot
10^{-4}y^{2}-0.0134y+0.046,$ (38)
$\displaystyle\displaystyle\rho_{r1}=-1.056\cdot 10^{-10}x^{5}+1.526\cdot
10^{-10}x^{4}y+3.667\cdot 10^{-9}x^{4}+3.284\cdot 10^{-10}x^{3}y^{2},$
$\displaystyle\displaystyle-3.606\cdot 10^{-9}x^{3}y-4.844\cdot
10^{-8}x^{3}-6.403\cdot 10^{-10}x^{2}y^{3}-3.359\cdot 10^{-9}x^{2}y^{2},$
$\displaystyle\displaystyle+2.797\cdot 10^{-8}x^{2}y+2.980\cdot
10^{-7}x^{2}-3.184\cdot 10^{-9}xy^{4}+1.668\cdot 10^{-8}xy^{3},$
$\displaystyle\displaystyle-6.75\cdot 10^{-9}xy^{2}-7.32\cdot
10^{-8}xy-8.22\cdot 10^{-7}x+1.03\cdot 10^{-8}y^{4}-4.66\cdot 10^{-8}y^{3}$
$\displaystyle\displaystyle+4.384\cdot 10^{-8}y^{2}+5.35\cdot
10^{-8}y+8.19\cdot 10^{-7},$ (39)
$\displaystyle\displaystyle\rho_{o2}=3.742\cdot 10^{-6}x^{5}+4.911\cdot
10^{-6}x^{4}*y-1.278\cdot 10^{-4}x^{4}-4.653\cdot 10^{-6}x^{3}y^{2},$
$\displaystyle\displaystyle-1.09\cdot 10^{-4}x^{3}y+0.0016x^{3}-3.25\cdot
10^{-5}x^{2}y^{3}+1.2\cdot 10^{-4}x^{2}y^{2}+8.35\cdot 10^{-4}x^{2}y,$
$\displaystyle\displaystyle-0.0087*x^{2}+9.718\cdot 10^{-5}xy^{4}+1.862\cdot
10^{-4}xy^{3}-7.897\cdot 10^{-4}xy^{2},$
$\displaystyle\displaystyle-0.0025xy+0.022x-6.073\cdot 10^{-4}y^{4}+4.18\cdot
10^{-4}y^{3}+9.108\cdot 10^{-4}y^{2}$
$\displaystyle\displaystyle+0.0027*y+1.812,$ (40) $\rho_{w3}=0.37y,$ (41)
$\displaystyle\displaystyle\rho_{w4}=2.61\cdot 10^{-4}x^{5}+1.67\cdot
10^{-4}x^{4}y-0.0079x^{4}+2.78\cdot 10^{-5}x^{3}y^{2},$
$\displaystyle\displaystyle-0.003x^{3}y+0.0914x^{3}+8.04\cdot
10^{-4}x^{2}y^{3}-0.0025x^{2}y^{2}+0.025x^{2}y-0.51x^{2},$
$\displaystyle\displaystyle-0.0038xy^{3}+0.0113xy^{2}-0.0744xy+1.3520x+0.0021y^{3}-0.0072y^{2}$
$\displaystyle\displaystyle+0.0616y-2.5890,$ (42)
$\displaystyle\displaystyle\rho_{o4}=-5.701\cdot 10^{-5}x^{5}-6.681\cdot
10^{-5}x^{4}*y+0.0019*x^{4}-7.327\cdot 10^{-5}x^{3}y^{2}$
$\displaystyle\displaystyle+0.0017x^{3}y-0.0231x^{3}+9.061\cdot
10^{-6}x^{2}y^{3}+0.0011x^{2}y^{2}-0.0142x^{2}y$
$\displaystyle\displaystyle+0.1273x^{2}+2.365\cdot 10^{-5}xy^{4}-1.965\cdot
10^{-4}xy^{3}-0.0046xy^{2}+0.0483xy$
$\displaystyle\displaystyle-0.322x+9.37\cdot
10^{-4}y^{5}-0.0023y^{4}+0.002y^{3}+0.006y^{2}-0.057y+0.34,$ (43)
$\displaystyle\displaystyle\rho_{r4}=5.86\cdot 10^{-9}x^{4}-3.762\cdot
10^{-10}x^{3}y-1.33\cdot 10^{-7}x^{3}+1.228\cdot 10^{-10}x^{2}y^{2}$
$\displaystyle\displaystyle+3.162\cdot 10^{-9}x^{2}y+1.055\cdot
10^{-6}x^{2}+8.124\cdot 10^{-9}xy^{3}-2.369\cdot 10^{-8}xy^{2}$
$\displaystyle\displaystyle+1.87\cdot 10^{-8}xy-3.282\cdot 10^{-6}x-4.032\cdot
10^{-8}y^{3}+1.093\cdot 10^{-7}y^{2}$ $\displaystyle\displaystyle-1.056\cdot
10^{-7}y+3.072\cdot 10^{-6}.$ (44)

The form of the function described for the above coefficients is the same for
the two types of oil used in this study. These formulas show how much the pH
affects the coefficients (see Figures 1-4). In all coefficients the
significant changes happen around $\displaystyle pH=5$.

The coefficients in (4.1)-(4.1) are used as input to solve the system of
conservation laws (24)-(27). First we use the numerical solution from the
solver COMSOL, which is based on the finite element method. Second, we obtain
the solution by the in house Riemann solver RPNfilho. To obtain a numerical
solution we provide the coefficients and explore a set of parameters which
permit to obtain stable solutions.

###  4.2 Wettability estimation

As suggested in [24], we assume that wettability is proportional to the Total
Bound product (TBP). TBP estimation is based on the quantification with SCM of
the attractive electrostatic forces existing between the rock-brine and the
oil-brine interfaces. These forces are represented by Bound Product (BP),
which is the product of the mole fractions of the oil ($\displaystyle O_{i}$)
and mineral ($\displaystyle m_{i}$) sites with unlike charges [24] and it is
given by the Equation:

$BP=O_{i}*m_{i}$ (45)

For a given mineral or rock, the sum of all the BP is the Total Bond Product
(TBP) [24]. The SCM using PHREEQC solver predicts the oil adhesion onto a
mineral or rock surface using the electrostatic pair linkages existing between
the mineral–brine and oil–brine interfaces. It has been reported in the
literature that the higher TBP corresponds to the higher tendency for oil to
be adsorbed onto the surface and vice-versa. The TBP can be written as:

$\displaystyle\displaystyle
TBP_{1}=(oil_{w,COO-})*(Cal_{s,CaOH2+}+Cal_{w,CO3Ca+}+Ca_{w,CO3Mg+}),$
$\displaystyle\displaystyle
TBP_{2}=(oil_{s,NH+})*(Cal_{w,CO3-}+Ca,_{s,CaO-}+Cal_{s,CaCO3-}+Cal_{s,CaSO4-}),$
$\displaystyle\displaystyle
TBP_{3}=oil_{w,COOCa+}*(Cal_{w,CO3-}+Cal_{s,CaO-}+Cal_{s,CaCO3-}+Cal_{s,CaSO4-}),$
$\displaystyle\displaystyle
TBP_{4}=oil_{w,COOMg+}*(Cal_{w,CO3-}+Cal_{s,CaO-}+Cal_{s,CaCO3-}+Cal_{s,CaSO4-}),$
$\displaystyle\displaystyle TBP=TBP_{1}+TBP_{2}+TBP_{3}+TBP_{4}.$

The complexes of oil and rock are calculated with the PHREEQC solver with
input values given in Tables 1-3.

Figure 5: TBP vs pH and Cl for temperature of T=39oC.  Figure 6: TBP vs pH and
Cl for temperature of T=39oC.

Figures 5 and 6 shows TBP for oil type 1 and 2 respectively (see Table 1). In
both cases, the TBP takes higher values for higher chloride concentration,
which corresponds to what was reported in [62]. Their results indicate that a
wettability shift from oil-wet to water-wet occurs when salinity is reduced.

The relation between TBP and $\displaystyle pH$ is slightly different with
respect to salinity ion concentration. Figure 6 shows that TBP has an
oscillation with respect to $\displaystyle pH$, i.e., for $\displaystyle
pH=2.7$ a peak is observed in TBP, to later decrease and obtain a minimum at
$\displaystyle pH=3$. Subsequently, TBP increases until it reaches a peak
above $\displaystyle pH=5$, and then decreases for higher pH values. This
behavior is maintained for temperatures of 95oC but the peaks are more
accentuated, see Figure 7. These observed patterns of the TBP relationship
with respect to pH are somewhat different from those reported by [62], where a
wettability shift from water-wet to oil-wet can be achieved with increasing
brine pH.

Our numerical experiment and the experimental result from [62] suggest that
the relationship between TBP and wettability is not proportional to the
quantity $\displaystyle pH$.

Figure 7: TBP vs pH and Cl for temperature of T=95oC

$\displaystyle\displaystyle
TBP(x,y)=0.0054x^{4}+0.0091x^{3}y-0.1177x^{3}-0.0041x^{2}y^{2}-0.1433x^{2}y$
$\displaystyle\displaystyle+0.9304x^{2}-0.0469xy^{3}+0.1905xy^{2}+0.5342xy-3.269x-0.099y^{4}$
$\displaystyle\displaystyle+0.615y^{3}-1.422y^{2}+0.2537y+4.72.$ (46)

Figure 8: Residual Oil Saturation as a function of the USBM (United States
Bureau of Mines) wetting index. [22]

Once wettability is estimated we use the procedure proposed in [13] to
estimate the residual oil $\displaystyle S_{or}$. Given the wettability index
we obtain $\displaystyle S_{or}$ from Figure 8. Next, we obtain the initial
water $\displaystyle S_{wi}$ using the empirical relationship ([48])

$S_{or}=2.0698S_{wi}^{3}-4.3857S_{wi}^{2}+2.17419S_{wi}+0.148.$ (47)

For estimating the residual oil $\displaystyle S_{or}$, we use the linear
relationship with salt concentration given in [35].

Here, we propose a new formula for the parameter $\displaystyle\theta$ in
equation (33) that takes into account the relationship between brine and pH as
well as the wettability. The formula provides a direct method for calculating
the parameter $\displaystyle\theta$, allowing to consider the concentration of
complexes in the fractional flow. The use of TBP in this context is a novel
approach, as previous works have used total ionic strength instead (see e.g.,
[38]).

Formula for the parameter $\displaystyle\theta$ in (33) in term of TBP is
given by

$\theta=\frac{TBP(pH,Cl)-TBP_{l}}{TBP_{h}-TBP_{l}},$ (48)

where TBP is given in (46), $\displaystyle TBP_{l}$ and $\displaystyle
TBP_{h}$ represent the TBP for the lower and higher brine concentration for a
fixed $\displaystyle pH$. For example for Oil 1 in Table 1 and $\displaystyle
pH=4$ we have $\displaystyle TBP_{l}=0.3841$ and $\displaystyle
TBP_{h}=0.534$. This formulation takes advantage of the established
relationship between TBP and wettability. Therefore, once TBP is scaled, it
serves to ensure the effect of salinity on the fractional flow function in an
appropriate way.

##  5 Solving the Riemann-Goursat problem by MOC

Once we have determined the coefficients of the system of conservation laws
(24)-(27) and fixed the fractional flow function $\displaystyle f_{w}$ in (32)
with $\displaystyle\theta$ given in (48), we can find the solution of such a
system. Notice that the fractional flow has been determined taking into
account the wettability values (through TBP), which in turn are estimated from
the complexes. This solution of the system will allow us to determine the
influence of the reduction of salt concentration at the point of injection.

The basis for MOC for the Riemann-Goursat problem is to assume that the
independent variables $\displaystyle(s_{w},pH,Cl,u)$ are functions of the
variable $\displaystyle\xi=x/t$, which is possible because the Riemann
solution are scale invariant under the map
$\displaystyle(t,x)\rightarrow(\alpha t,\alpha x)$; then we can take
$\displaystyle
U=(s_{w},pH,Cl,u)(x,t)=(s_{w},pH,Cl,u)(x(\xi),t(\xi))=(s_{w},y,u)(\xi)$ with
characteristic lines $\displaystyle\xi=x/t$. If we assume that the functions
$\displaystyle U(\xi)=(s_{w}(\xi),y(\xi),u(\xi))$ are differentiable along the
characteristic lines, the following generalized eigenvalue equation for
characteristics values holds

$Ar=\lambda Br,\quad\text{where}\quad A={\partial F}/{\partial U},\quad
B={\partial G}/{\partial U}.$ (49)

The eigenvector $\displaystyle r$ is parallel to $\displaystyle d{U}/d\xi$, so
the rarefaction curves are tangent to the characteristic field given by the
normalized eigenvector $\displaystyle r$.

The physical model consists of four unknowns states variables $\displaystyle
s_{w},pH,[Cl],u$ satisfying system (24)-(27). The system has four equations
with the unknown variables $\displaystyle(s,y,u)$ with $\displaystyle
y=(y_{1},y_{2})=(pH,[Cl])$ given by

$\displaystyle\displaystyle\frac{\partial}{\partial
t}(\varphi\rho_{wj}(y)s_{w}+\varphi\rho_{oj}(y)s_{o}+(1-\varphi)\rho_{rj}(y))+\frac{\partial}{\partial
x}((u)(\rho_{wj}(y)f_{w}+\rho_{oj}(y)f_{o})=0.$ (50)

We take $\displaystyle s=s_{w}$, $\displaystyle f=f_{w}$, $\displaystyle
s_{o}=1-s$ and $\displaystyle f_{o}=1-f_{w}$. The accumulation $\displaystyle
G$ and flux functions $\displaystyle\hat{F}=(u/\phi)F$ are given by

$\displaystyle\displaystyle
G_{j}=\varphi\rho_{wj}(y)s_{w}+\varphi\rho_{oj}(y)s_{o}+(1-\varphi)\rho_{rj}(y)$
(51) $\displaystyle\displaystyle F_{j}=\rho_{wj}(y)f_{w}+\rho_{oj}(y)f_{o}$
(52)

The index $\displaystyle w$ (water) is often replace by the index
$\displaystyle a$ (aqueous phase), and the index $\displaystyle j=1,\cdots,4$
is used to denote chemical species.

The main features in the $\displaystyle 1$-D Riemann solutions of hyperbolic
systems are the rarefaction and shock curves. The rarefactions are obtained
from integral curves of the line fields, given by the eigenvectors of system
(49), where the matrices $\displaystyle A,B$ represent the Jacobian matrices
$\displaystyle A=D\hat{F}$ and $\displaystyle B=DG$ for the flux
$\displaystyle\hat{F}$ and the accumulation $\displaystyle G$. For system
$\displaystyle(\ref{eq1})$, the matrices $\displaystyle B$, $\displaystyle A$
are given by

$\displaystyle\displaystyle B_{i,1}=[\rho_{i}],\quad
B_{i,k+1}=\frac{\partial\rho_{wi}}{\partial
y_{k}}s+\frac{\partial\rho_{oi}}{\partial
y_{k}}s_{o}+\frac{\partial\rho_{ri}}{\partial y_{k}},\text{ and }B_{i,n+1}=0.$
(53) $\displaystyle\displaystyle A_{i,1}=(u/\phi)[\rho_{i}]\frac{\partial
f}{\partial s},\quad A_{i,k+1}=(u/\phi)\left(\frac{\partial\rho_{wi}}{\partial
y_{k}}f+\frac{\partial\rho_{oi}}{\partial y_{k}}f_{o}\right)\text{ and
}A_{i,n+1}=F_{i},$ (54)

in which $\displaystyle i=1,\cdots,4$, $\displaystyle k=1,2$ and

$[\rho_{i}]=\rho_{wi}-\rho_{oi}.$ (55)

From the Jacobian matrices, we obtain the eigenpairs which we have summarized
in

######  Proposition 1.

The eigenpairs of the eigenvalue problem $\displaystyle A\vec{r}=\lambda
B\vec{r}$, where the matrices $\displaystyle B$ and $\displaystyle A$
represent the Jacobian of the accumulation and flux terms of system
$\displaystyle(\ref{eq1})$ are the Buckley-Leverett eigenpairs
($\displaystyle\lambda_{s},\vec{r}_{s}$) given by

$\lambda_{s}=\frac{u}{\varphi}\frac{\partial f}{\partial s_{w}}\quad\text{ and
}\quad\vec{r}_{s}=(1,0,0,0)^{T},$ (56)

and two composition chemical with eigenpairs given by
$\displaystyle(\lambda_{\Lambda},\vec{r}_{\Lambda})$:

$\lambda_{\Lambda}=\frac{u}{\varphi}\frac{f-\Lambda}{s-\Lambda}.$ (57)

We obtain $\displaystyle\Lambda$ and $\displaystyle\vec{v}=(v_{1},v_{2})^{T}$
as the solutions of the compositional generalized eigenvalue problem

$((\mathcal{A}+\mathcal{A}_{1})-\Lambda(\mathcal{B}+\mathcal{B}_{1}))\vec{v}=0,$
(58)

where the matrices $\displaystyle\mathcal{A},\mathcal{A}_{1},\mathcal{B}$ and
$\displaystyle\mathcal{B}_{1}$ are given in (88) and (90).

Moreover, we obtain $\displaystyle\vec{r}_{\Lambda}=(r_{1},\cdots,r_{4})$ as

$\displaystyle\displaystyle
r_{4}=-\frac{\displaystyle\sum_{2}^{3}\mathcal{G}^{(4)}_{4,j}v_{j-1}}{\nu_{4}},~{}~{}r_{1}=\frac{\mathcal{G}_{1,4}r_{4}-\displaystyle\sum_{2}^{3}\mathcal{G}_{1,j}v_{j-1}}{\mathcal{G}_{1,1}},$
(59) $\displaystyle\displaystyle r_{j}=v_{j-1},\quad\text{for}\quad j=2,3.$
(60)

where the matrices $\displaystyle\mathcal{G}$ and
$\displaystyle\mathcal{G}^{(4)}_{4,j}$ evaluated at
$\displaystyle\lambda=\lambda_{\Lambda}$ are given in (61) and (77)-(79).

Proof of Proposition

The idea of the proof consist of reducing the matrix
$\displaystyle\mathcal{G}$ by continuous application of Gauss procedure
leaving the matrix as a reduced matrix where the eigenvalues and eigenvectors
are obtained easily.

To obtain the eigenvalues we solve $\displaystyle det(A-\lambda B)=0$, where
$\displaystyle A-\lambda B=\mathcal{G}$, and
$\displaystyle\mathcal{G}=(\mathcal{G}_{i,j})$ for $\displaystyle
i,j=1,\cdots,4$. Here we use an auxiliary variable $\displaystyle k$ that
ranges from $\displaystyle 1$ to $\displaystyle 2$ and we write
$\displaystyle\mathcal{G}_{i,j}$ as:

$\mathcal{G}_{i,1}=[\rho_{i}]\xi_{1},\quad\mathcal{G}_{i,k+1}=\frac{\partial\rho_{wi}}{\partial
y_{k}}\xi_{2}+\frac{\partial\rho_{oi}}{\partial
y_{k}}\xi_{3}-\lambda\frac{\partial\rho_{ri}}{\partial y_{k}}\quad\text{ and
}\quad\mathcal{G}_{i,4}=F_{i},$ (61)

where we define the auxiliary variables

$\xi_{1}=\left(u\frac{\partial f}{\partial
s}-\lambda\right),\quad\xi_{2}=(uf-\lambda s)\quad\text{ and
}\quad\xi_{3}=(uf_{o}-\lambda s_{o}).$ (62)

For calculations purposes, that we do in the next proofs, it is useful to
define the following functions $\displaystyle\gamma_{ij}$,
$\displaystyle\varrho_{ij}$, $\displaystyle\nu_{i}$,
$\displaystyle\vartheta_{ij}$ and $\displaystyle\varsigma_{ij}$ for
$\displaystyle i=1,\cdots,4$, $\displaystyle j=1,2$ as

$\displaystyle\displaystyle\gamma_{ij}=\frac{\partial\rho_{wi}}{\partial
y_{j}}[\rho_{1}]-\frac{\partial\rho_{w1}}{\partial
y_{j}}[\rho_{i}],\quad\varrho_{ij}=\frac{\partial\rho_{oi}}{\partial
y_{j}}[\rho_{1}]-\frac{\partial\rho_{o1}}{\partial y_{j}}[\rho_{i}],$ (63)
$\displaystyle\displaystyle\pi_{ij}=\frac{\partial\rho_{ri}}{\partial
y_{j}}[\rho_{1}]-\frac{\partial\rho_{r1}}{\partial
y_{j}}[\rho_{i}],\quad\nu_{i}=[\rho_{1}]F_{i}-[\rho_{i}]F_{1},$ (64)
$\displaystyle\displaystyle\vartheta_{ij}=\gamma_{ij}\nu_{4}-\gamma_{4,j}\nu_{i},\quad\varsigma_{ij}=\varrho_{ij}\nu_{4}-\varrho_{4,j}\nu_{i},$
(65) $\displaystyle\displaystyle\tau_{ij}=\pi_{ij}\nu_{4}-\pi_{4,j}\nu_{i}.$
(66)

All these coefficients ((63)- (66)) depend only of the variables
$\displaystyle y$. Now, we are able to perform the proof os results in the
text.

  * 1)

Substituting the $\displaystyle i$-th row of matrix
$\displaystyle(\ref{matrixnew})$, for $\displaystyle i=2,\cdots,4$, by the sum
of the first row of $\displaystyle(A-\lambda B)$ (the elements of which are
given by (61)) multiplied by $\displaystyle-[\rho_{i}]$ with its
$\displaystyle i$-th row multiplied by $\displaystyle[\rho_{1}]$, we obtain,
for $\displaystyle j=1,2$ and $\displaystyle i=2,3,4$

$\mathcal{G}\sim\begin{pmatrix}\mathcal{G}_{1,1}&\mathcal{G}_{1,j+1}&\mathcal{G}_{1,4}\\\
\mathbb{O}&\mathcal{G}^{(1)}_{i,j+1}&\mathcal{G}^{(1)}_{i,4}\\\
\end{pmatrix}.$ (67)

Here $\displaystyle\mathbb{O}$ is the column vector of three zeros and block
matrices $\displaystyle(\mathcal{G}^{(1)}_{i,j+1})$ and
$\displaystyle(\mathcal{G}^{(1)}_{i,4})$ for $\displaystyle i=2,3,4$ and
$\displaystyle j=1,2$ are the block matrices whose elements are given by

$\mathcal{G}^{(1)}_{i,j+1}=\left(\gamma_{ij}\xi_{2}+\varrho_{ij}\xi_{3}-\pi_{ij}\lambda\right)_{1\leq
j\leq 2}\quad\text{ and }\quad\mathcal{G}^{(1)}_{i,4}=\nu_{i},$ (68)

where $\displaystyle\gamma_{ij}$, $\displaystyle\varrho_{ij}$
$\displaystyle\pi_{ij}$ and $\displaystyle\nu_{i}$ are given by Eq. (63) and
(64).

Notice here that if $\displaystyle[\rho_{i}]$, given y
$\displaystyle(\ref{diferenca})$, is zero for some index $\displaystyle i$,
the corresponding position in the first column is zero and we do not need to
perform calculations to vanish this position. Moreover, if
$\displaystyle[\rho_{1}]$ is zero, we exchange the position of row
$\displaystyle 1$ with another row to obtain a non-zero pivot.

  * 2)

Now, we substitute $\displaystyle i$-th row of $\displaystyle(\ref{matrix3})$,
for $\displaystyle i=2,3$, by the sum of the four row of
$\displaystyle(\ref{matrix3})$ multiplied by $\displaystyle-\nu_{i}$ with the
$\displaystyle i$-th row of $\displaystyle(\ref{matrix3})$ multiplied by
$\displaystyle\nu_{4}$, and we obtain for $\displaystyle j=1,2$ and
$\displaystyle i=2,3$

$\mathcal{G}\sim\begin{pmatrix}\mathcal{G}_{1,1}&\mathcal{G}_{1,j+1}&\mathcal{G}_{1,4}\\\
\mathbb{O}&\mathcal{G}^{(2)}_{i,j+1}&\mathbb{O}\\\
0&\mathcal{G}^{(2)}_{4,j+1}&\nu_{4}\\\ \end{pmatrix}.$ (69)

Here $\displaystyle\mathbb{O}$ is the column vector of two zeros and
$\displaystyle(\mathcal{G}^{(2)}_{i,j+1})$ and
$\displaystyle(\mathcal{G}^{(1)}_{4,j+1})$ for $\displaystyle i=2,3$ and
$\displaystyle j=1,2$ are the block matrices given by

$\mathcal{G}^{(2)}_{i,j+1}=\vartheta_{ij}\xi_{2}+\varsigma_{ij}\xi_{3}-\tau_{ij}\lambda,~{}\mathcal{G}^{(2)}_{4,j+1}=\left(\gamma_{4j}\xi_{2}+\varrho_{4j}\xi_{3}-\pi_{4j}\lambda\right),$
(70)

where $\displaystyle\vartheta_{ij}$, $\displaystyle\varsigma_{ij}$ and
$\displaystyle\tau_{ij}$ are given by (65) and $\displaystyle(\ref{gamma})$,
for $\displaystyle i=1,\cdots,4$, $\displaystyle j=1,2$.

  * 3)

Now, we substitute the $\displaystyle i$-th row of $\displaystyle(\ref{pms})$,
for $\displaystyle i=2,3$, by the sum of the four row of
$\displaystyle(\ref{pms})$ multiplied by $\displaystyle-\tau_{i1}$ with the
$\displaystyle i$-th row of $\displaystyle(\ref{pms})$ multiplied by
$\displaystyle\pi_{41}$, and we obtain for $\displaystyle j=1,2$ and
$\displaystyle i=2,3$

$\mathcal{G}\sim\begin{pmatrix}\mathcal{G}_{1,1}&\mathcal{G}_{1,j+1}&\mathcal{G}_{1,4}\\\
\mathbb{O}&\mathcal{G}^{(3)}_{i,j+1}&\mathbb{O}\\\
0&\mathcal{G}^{(3)}_{4,j+1}&\nu_{4}\\\ \end{pmatrix}.$ (71)

Here $\displaystyle\mathbb{O}$ is the column vector of two zeros and
$\displaystyle(\mathcal{G}^{(3)}_{i,j+1})$ and
$\displaystyle(\mathcal{G}^{(3)}_{4,j+1})$ for $\displaystyle i=2,3$ and
$\displaystyle j=1,2$ are the block matrices given by

$\mathcal{G}^{(3)}_{i,2}=m_{i1}\xi_{2}+n_{i1}\xi_{3},~{}\mathcal{G}^{(3)}_{i,3}=m_{i2}\xi_{2}+n_{i2}\xi_{3}-\alpha_{i2}\lambda$
(72)
$\mathcal{G}^{(3)}_{4,j+1}=\left(\gamma_{4j}\xi_{2}+\varrho_{4j}\xi_{3}-\pi_{4j}\lambda\right),$
(73)

where

$\displaystyle\displaystyle
m_{ij}=\vartheta_{ij}\pi_{41}-\gamma_{4j}\tau_{i1},~{}n_{ij}=\varsigma_{ij}\pi_{41}-\rho_{4j}\tau_{i1},$
(74)
$\displaystyle\displaystyle\alpha_{i2}=\tau_{i2}\pi_{41}-\pi_{42}\tau_{i1},~{}\alpha_{i1}=0.$
(75)

and $\displaystyle\vartheta_{ij}$, $\displaystyle\varsigma_{ij}$ and
$\displaystyle\tau_{ij}$ are given by (64) and $\displaystyle(\ref{gamma})$,
for $\displaystyle i=1,\cdots,4$, $\displaystyle j=1,2$.

  * 4)

Now, we substitute third row of $\displaystyle(\ref{pmsa1})$, by the sum of
the second row of $\displaystyle(\ref{pmsa1})$ multiplied by
$\displaystyle-\alpha_{32}$ with the third row of $\displaystyle(\ref{pmsa1})$
multiplied by $\displaystyle\alpha_{21}$, and we obtain

$\mathcal{G}\sim\begin{pmatrix}\mathcal{G}_{1,1}&\mathcal{G}_{1,j+1}&\mathcal{G}_{1,4}\\\
\mathbb{O}&\mathcal{G}^{(4)}_{i,j+1}&\mathbb{O}\\\
0&\mathcal{G}^{(4)}_{4,j+1}&\nu_{4}\\\ \end{pmatrix}.$ (76)

Here $\displaystyle\mathbb{O}$ is the column vector of two zeros and
$\displaystyle(\mathcal{G}^{(4)}_{i,j+1})$ and
$\displaystyle(\mathcal{G}^{(4)}_{4,j+1})$ for $\displaystyle i=2,3$ and
$\displaystyle j=1,2$ are the block matrices given by:

$\mathcal{G}^{(4)}_{2,2}=\bar{m}_{21}\xi_{2}+\bar{n}_{21}\xi_{3},~{}\mathcal{G}^{(4)}_{2,3}=\bar{m}_{22}\xi_{2}+\bar{n}_{22}\xi_{3}+\alpha_{22}\lambda$
(77)
$\mathcal{G}^{(4)}_{3,2}=\bar{m}_{31}\xi_{2}+\bar{n}_{31}\xi_{3},~{}\mathcal{G}^{(4)}_{3,3}=\bar{m}_{32}\xi_{2}+\bar{n}_{32}\xi_{3}$
(78)
$\mathcal{G}^{(4)}_{4,j+1}=\left(\gamma_{4j}\xi_{2}+\varrho_{4j}\xi_{3}+\pi_{4j}\lambda\right),$
(79)

where

$\bar{m}_{21}=m_{21},\bar{n}_{21}=n_{21},\bar{m}_{22}=m_{22},\bar{n}_{22}=n_{22},$
(80)

and

$\displaystyle\displaystyle\bar{m}_{31}=m_{31}\alpha_{21}-m_{21}\alpha_{32},~{}\bar{n}_{31}=n_{31}\alpha_{21}-n_{21}\alpha_{32}.$
(81)
$\displaystyle\displaystyle\bar{m}_{32}=m_{32}\alpha_{21}-m_{22}\alpha_{32},~{}\bar{n}_{32}=n_{32}\alpha_{21}-n_{22}\alpha_{32}.$
(82)

  * 5)

Since $\displaystyle\mathcal{G}^{(4)}_{i,j+1}$ is block matrix, but this
matrix appears in the rows $\displaystyle 2$ to $\displaystyle 3$ and columns
$\displaystyle 2$ to $\displaystyle 3$ it is useful to define matrix
$\displaystyle\left(G^{(4)}_{l,r}\right)$ for $\displaystyle l,r=1,\cdots,2$
from $\displaystyle\left(\mathcal{G}^{(4)}_{i,j+1}\right)$

$G^{(4)}_{l,r}=\mathcal{G}^{(4)}_{l+1,r+1}\quad\text{ for }l,r=1,\cdots,2.$
(83)

From $\displaystyle(\ref{pmsa2})$ and $\displaystyle\det(A-\lambda B)=0$ a
solution is $\displaystyle\xi_{1}=0$. Since $\displaystyle\xi_{1}$ is given by
$\displaystyle(\ref{xis}.a)$, thus we obtain the eigenpair
$\displaystyle(\lambda_{s},\vec{r}_{s})$ given by
$\displaystyle(\ref{lambdasa})$. For this eigenpair, only saturation changes
and we identify this family wave as saturation wave, or Buckley-Leverett type
wave.

Using $\displaystyle(\ref{pmsa2})$ and $\displaystyle det(A-\lambda B)=0$, the
other eigenvalues are obtained by solving

$\det({G}^{(4)}_{l,r})=0\quad\text{ for }\quad l,r=1,2,$ (84)

with $\displaystyle{G}^{(4)}_{l,r}$ given by (77) and (78).

To obtain the corresponding eigenvalue, we substitute
$\displaystyle\lambda_{\Lambda}$ given by $\displaystyle(\ref{eig2a})$ into
$\displaystyle(\ref{eigeqs})$ and using that

$\xi_{2}=uf-\lambda_{\Lambda}s=u\frac{f(s-\Lambda)-(f-\Lambda)s}{s-\Lambda}=u\Lambda\frac{s-f}{s-\Lambda}.$
(85)
$\xi_{3}=uf_{o}-\lambda_{\Lambda}s_{o}=u(1-f)-\lambda_{\Lambda}(1-s)=u(1-\Lambda)\frac{s-f}{s-\Lambda}.$
(86)

Substituting $\displaystyle\xi_{2}$ and $\displaystyle\xi_{3}$ in (77) and
(78) the matrix in (84) can be rewrite as

$\left(\begin{array}[]{cc}\bar{n}_{21}+\Lambda(\,m_{21}-\,\bar{n}_{21})&\bar{n}_{22}+\Lambda(\,m_{22}-\,\bar{n}_{22})\\\
\bar{n}_{31}+\Lambda(\,m_{31}-\,\bar{n}_{31})&\bar{n}_{32}+\Lambda(\,m_{32}-\,\bar{n}_{32})\end{array}\right)+\left(\begin{array}[]{cc}0&\frac{\alpha_{22}(f-\Lambda)}{(s-f)}\\\
0&0\end{array}\right)$ (87)

Denoting by

$\mathcal{A}=\left(\begin{array}[]{cc}\bar{n}_{21}&\bar{n}_{22}\\\
\bar{n}_{31}&\bar{n}_{32}\end{array}\right)~{}~{}~{}~{}\text{and}~{}~{}~{}\mathcal{B}=-\left(\begin{array}[]{cc}\,m_{21}-\,\bar{n}_{21}&\,m_{22}-\,\bar{n}_{22}\\\
\,m_{31}-\,\bar{n}_{31}&\,m_{32}-\,\bar{n}_{32}\end{array}\right),$ (88)

we obtain

$\mathcal{A}-\Lambda\mathcal{B}+\left(\begin{array}[]{cc}0&\frac{\alpha_{22}(f-\Lambda)}{(s-f)}\\\
0&0\end{array}\right).$ (89)

Also, denoting by

$\mathcal{A}_{1}=\left(\begin{array}[]{cc}0&\frac{\alpha_{22}f}{(s-f)}\\\
0&0\end{array}\right)~{}~{}~{}~{}\text{and}~{}~{}~{}\mathcal{B}_{2}=\left(\begin{array}[]{cc}0&\frac{\alpha_{22}}{(s-f)}\\\
0&0\end{array}\right),$ (90)

we have

${G}^{(4)}_{l,r}=(\mathcal{A}-\Lambda\mathcal{B})+(\mathcal{A}_{1}-\Lambda\mathcal{B}_{2}),$
(91)

where the matrix $\displaystyle\mathcal{A}$ and $\displaystyle\mathcal{B}$
only depend the variables $\displaystyle y$, while the matrix
$\displaystyle\mathcal{A}_{1}$ and $\displaystyle\mathcal{B}_{1}$ depend on
$\displaystyle y$ and $\displaystyle s$.

In this notation is easy to verify that equation (84) is a quadratic equation
in the variable $\displaystyle\Lambda$. Thus equation (58) is a
$\displaystyle(2)\times(2)$ algebraic system. If we have two different roots,
we obtain two more eigenvalues of form $\displaystyle(\ref{eig2a})$. We call
each eigenvalues as $\displaystyle\lambda_{\Lambda_{i}}$ for $\displaystyle
i=1,2$.

Also, notice that if $\displaystyle\rho_{rj}=0$ $\displaystyle(j=1,\cdots,4)$
then $\displaystyle\alpha_{22}=0$ and the corresponding roots of equation (84)
only depend on the chemical variable $\displaystyle y$. In this case solution
can be found more easily (see [4]).

  * 6)

The eigenvectors related with $\displaystyle\lambda_{\Lambda}$ are obtained
solving

$\begin{pmatrix}\mathcal{G}_{1,1}&\mathcal{G}_{1,j+1}&\mathcal{G}_{1,4}\\\
\mathbb{O}&\mathcal{G}^{(4)}_{i,j+1}&\mathbb{O}\\\
0&\mathcal{G}^{(4)}_{4,j+1}&\nu_{4}\\\ \end{pmatrix}\vec{r}=0.$ (92)

We can split the calculation of $\displaystyle\vec{r}$. First we obtain the
coordinates $\displaystyle(r_{2},r_{3})$ of eigenvector
$\displaystyle\vec{r}_{\Lambda}=(r_{1},r_{2},r_{3},r_{4})$. To do our
calculations, we define the auxiliary vector $\displaystyle\vec{v}$ of two
coordinates, which is the solution of

$\left({G}^{(4)}_{l,r}\right)\vec{v}^{T}=0\quad\text{ for }\quad l,r=1,2.$
(93)

Using $\displaystyle(\ref{ueq1})$-$\displaystyle(\ref{ueq1n})$, we can see
that $\displaystyle(\ref{pmsnok})$, after simplifications is written again as
$\displaystyle(\ref{eigeqs2})$.

The coordinates $\displaystyle r_{n+1}$ and $\displaystyle r_{1}$ are obtained
by solving the first and the last equations of $\displaystyle(\ref{pmsno})$
and using $\displaystyle(\ref{ueq1})$-$\displaystyle(\ref{ueq1n})$, after some
tedious calculations, we obtain Eqs.
$\displaystyle(\ref{eiget1a})-(\ref{eiget1})$. $\displaystyle\square$

The integral curves $\displaystyle\mathcal{W}^{s}$ and
$\displaystyle\mathcal{W}^{\Lambda}$ associated to $\displaystyle\vec{r}_{s}$
and $\displaystyle\vec{r}_{\Lambda}$ are obtained by integrating the ODEs

$d(s,y,u)/d\xi=\vec{r}_{s}~{}\text{and}~{}\frac{d(s,y,u)}{d\xi}=\vec{r}_{\Lambda}.$
(94)

From (94)$\displaystyle a$, we obtain the Buckley-Leverett(B-L) wave, where
only the saturation changes, which is denotes by
$\displaystyle\mathcal{R}_{s}$. Moreover, the solution of (94)$\displaystyle
b$ determines the waves associated to $\displaystyle\lambda_{\Lambda}$ denotes
by $\displaystyle\mathcal{R}_{\Lambda}$. Here, we do the technical proofs of
results which appear in paper [5, 4]. These papers describe the numeric and
theory of the main wave interactions of the system of conservation laws
studied here.

###  5.1 Rankine-Hugoniot Locus

In this section we summarized how to calculated the discontinuous solution of
the system of conservation laws studied here. Details of proof can be found in
[5]. We denote $\displaystyle y=(pH,Cl)$ to gain in clearance.

The discontinuous solution of system of conservation laws
$\displaystyle(\ref{eq1})$ satisfy the Rankine-Hugoniot condition, i.e. for a
given left and right state $\displaystyle(s^{-},y^{-},u^{-})$ and
$\displaystyle(s^{+},y^{+},u^{+})$ respectively we have

$u^{+}F_{i}((s^{+},y^{+})-u^{-}F_{i}(s^{-},y^{-})=v^{s}(G_{i}(s^{+},y^{+})-G_{i}(s^{-},y^{-})),$
(95)

with $\displaystyle i=1,\dots,n+1$ and $\displaystyle v^{s}$ is the speed of
discontinuity. Equation (95) can be rewritten as

$\displaystyle\displaystyle\Phi_{i}.[v^{s},u^{+},u^{-}]=0,$ (96)

where $\displaystyle\Phi_{i}=(\Phi_{i1},\Phi_{i2},\Phi_{i2})$, with
$\displaystyle\Phi_{i1}=s^{+}(\rho_{wi}^{+}-\rho_{oi}^{+})+\rho_{oi}^{+}-(s^{-}(\rho_{wi}^{-}-\rho_{oi}^{-})+\rho_{oi}^{-})+(\rho_{ri}^{+}-\rho_{ri}^{+})$,
$\displaystyle\Phi_{i2}=-(f^{+}(\rho_{wi}^{+}-\rho_{oi}^{+})+\rho_{oi}^{+})$,
$\displaystyle\Phi_{i3}=(f^{+}(\rho_{wi}^{-}-\rho_{oi}^{-})+\rho_{oi}^{-})$
and $\displaystyle\rho^{+}=\rho(y^{+})$, $\displaystyle\rho^{-}=\rho(y^{-})$,
$\displaystyle f^{+}=f(s^{+},y^{+})$, $\displaystyle f^{-}=f(s^{-},y^{-})$.

For each fixed state $\displaystyle(s^{-},y^{-})$, the Hugoniot-locus
$\displaystyle\mathcal{HL}(s^{-},y^{-})$ consist of all states
$\displaystyle(s^{+},y^{+})$ satisfying $\displaystyle(\ref{rank1})$. In [43]
is prove that

$\mathcal{HL}(s^{-},y^{-})=\left\\{(s^{+},y^{+}):det(\Phi_{i}^{T},\Phi_{k}^{T},\Phi_{j}^{T})=0\right\\},$
(97)

for all the combination of distinct index
$\displaystyle\\{i,j,k\\}\in\\{1,2,\cdots,n+1\\}$. Also, we verify that
instead to consider all the combination it is enough to reduce $\displaystyle
n-1$ equations. This results can summarized as follows: let $\displaystyle
i_{1}$ and $\displaystyle i_{2}$ $\displaystyle\in\\{1,2,\cdots,n+1\\}$ two
indices such that $\displaystyle\Phi_{i_{1}}$ and $\displaystyle\Phi_{i2}$ are
independent linearly, then Equation (97) reduces to

$\mathcal{HL}(s^{-},y^{-})=\left\\{(s^{+},y^{+}):det(\Phi_{k}^{T},\Phi_{i_{1}}^{T},\Phi_{i_{2}}^{T})=0\right\\},$
(98)

for $\displaystyle k\in\\{1,2,\cdots,n+1\\}$ distinct of $\displaystyle i_{1}$
and $\displaystyle i_{2}$. Equation (98) represents a curve in the three
dimensional space $\displaystyle(s,y)$ consisting in the intersection of two
surfaces.

For fixed $\displaystyle(s^{-},y^{-})$, there exists a branch of
$\displaystyle\mathcal{HL}(s^{-},y^{-})$ consisting of the states of the form
$\displaystyle(s,y)\in\Omega$, with $\displaystyle s$ variable and
$\displaystyle y=y^{-}$. This branch is called Buckley-Leverett or saturation
branch, and it is denoted by $\displaystyle\mathcal{H}_{s}$. The other
branches are denoted by $\displaystyle\mathcal{H}_{\Lambda}$.

###  5.2 Numerical and semi analytical solution

In this work we find the solution for the system (24)-(27) by using the COMSOL
solver and the in house Riemann solver called RPNfilho. We are interested in
the Riemann-Goursat problem for these equations with piecewise constant
initial data

$\left\\{\begin{array}[c]{ll}L=(Sw_{L},pH_{L},Cl_{L},u_{L})&\textbf{if}\hskip
5.69046ptx<0,\\\ R=(Sw_{R},pH_{R},Cl_{R},\;\cdot\;)&\textbf{if}\hskip
5.69046ptx>0.\end{array}\right.$ (99)

The Riemann solution is obtained by the wave curve method. The solution is
constructed by means of a sequence of elementary waves $\displaystyle w_{k}$
(shocks and rarefactions) for $\displaystyle k=1,2,\cdots,m$ and constant
states $\displaystyle U_{k}$ for $\displaystyle k=1,2,\cdots,p$, in which
$\displaystyle p$ is not known a priori.

At any rate, this sequence of waves can be written as

$\mathcal{U}_{L}\equiv\mathcal{U}_{0}\overset{w_{1}}{\longrightarrow}\mathcal{U}_{1}\overset{w_{2}}{\longrightarrow}\cdots\overset{w_{m}}{\longrightarrow}\mathcal{U}_{m}\equiv\mathcal{U}_{R},$
(100)

where $\displaystyle\mathcal{U}=(s,pH,Cl,u)$. In the Riemann solution it is
necessary that the waves have increasing speed, satisfying the so called
geometrical compatibility. Sometimes, this geometrical compatibility is
sufficient to furnish existence and uniqueness of the solution. Moreover, this
condition is used to select the physical sequence of waves for the Riemann
solution. We do not go into details about the wave curve method in this paper.
Details can be found in [5].

The waves consist of the B-L saturation wave curve
$\displaystyle\mathcal{R}_{s}$ where only the saturation varies, the shock
wave curve $\displaystyle\mathcal{H}_{s}$, the locus where only the saturation
varies and the chemical saturation wave curves
$\displaystyle\mathcal{R}_{\Lambda_{i}}$ ($\displaystyle i=1,2$) associated to
the couple $\displaystyle(\lambda_{\Lambda_{i}},r_{\Lambda_{i}})$. Moreover
other compatible wave curves include the composite wave curve $\displaystyle
C_{s}=\mathcal{R}_{s}\cup\mathcal{H}_{s}$, formed by a characteristic B-L
shock curve, and the composite wave curves $\displaystyle
C_{\Lambda_{i}}=\mathcal{R}_{\Lambda_{i}}\cup\mathcal{H}_{\Lambda_{i}}$, where
$\displaystyle\mathcal{H}_{\Lambda_{i}}$ is a characteristic shock curve
associated to the branch of the family $\displaystyle\lambda_{\Lambda_{i}}$
with $\displaystyle i=1,2$ and the constant state $\displaystyle\mathcal{C}$.

##  6 Numerical solution

In this section we present the numerical solution for the system (24)-(27) by
using COMSOL and the Riemann solver RPNfilho. We solve the Riemann-Goursat
problem for some particular initial (right) and injection (left) conditions
for such system of equations. We do not impose any condition on the variable
$\displaystyle u_{R}$ in (99) because its value is obtained from the other
variables together with the solution of the system. Physically, the left state
represents the value of the variables $\displaystyle pH$, $\displaystyle[Cl]$,
$\displaystyle S_{w}$ and $\displaystyle u$ at the injection point. The right
state presents the initial condition of such variables in the reservoir. The
states studied here correspond to typical values of interest for petroleum
engineering.

We assume that the fluid is incompressible, but there is mass transfer between
phases of the carbon dioxide and the partial molar volume differs between
phases, thus a variable total Darcy velocity ensues.

###  6.1 RPNfilho and COMSOL solutions

We are interested in evaluating the salt concentrations reduction at the
injection point. To do so, first we cross-validated the method and second we
identify the main waves present in the simulation of relevant examples.
Besides calculating the speeds of the water saturation front and of the saline
front, we estimate the behavior of the $\displaystyle pH$.

In this simulation we take a reduction of 5$\displaystyle\%$ in the initial
chloride concentration in the reservoir, which we take as 0.3 mole per
kilogram of water. To the chloride concentration of 0.015 mole per kilogram of
water corresponds an oil saturation $\displaystyle S_{or}=0.228$ and an
initial water $\displaystyle S_{wi}=0.0398$. The initial condition for the
Riemann problem consists of $\displaystyle L=(0.7604,4.0,0.015,1.0\times
10^{-5})$ and $\displaystyle R=(0.0398,4.0,0.3)$.

Figures 9-10 show the water saturation, the chloride and the $\displaystyle
pH$ profiles corresponding to the Riemann problem for a given $\displaystyle
L$ and $\displaystyle R$. We verify that in this example the solution of
RPNfilho and COMSOL solvers match well. The saturation profile for both models
coincide. The pH and chloride profiles match well except in the shock front
solution.

The structure of the Riemann solution obtained by RPNfilho solver consists of
a chemical 2-rarefaction wave, a constant state, a B-L rarefaction, a constant
state, a 2-shock and a B-L shock. This sequence of waves can be represented
schematically as
$\displaystyle(L)\xrightarrow{\mathcal{R}_{\Lambda_{2}}}(A)\xrightarrow{\mathcal{R}_{s}}(B)\xrightarrow{\mathcal{H}_{\Lambda_{2}}}(C)\xrightarrow{\mathcal{H}_{s}}(R)$.
This scheme means that from state $\displaystyle L$ to the state
$\displaystyle A=(0.7604,4.12,0.015)$ there exists a rarefaction
$\displaystyle\mathcal{R}_{\Lambda_{2}}$ of the chemical family associated to
eigenvalues $\displaystyle\lambda_{\Lambda_{2}}$. This wave is followed by a
constant state from the state $\displaystyle A$ to the state $\displaystyle
B=(0.61,4.12,0.015)$, then we have a B-L rarefaction wave. Then we have
another constant state followed by a type 2 shock to the point $\displaystyle
C$. Immediately afterwards a constant state follows a B-L shock and reaches
the right state $\displaystyle R$.

Figure 9: a) Saturation profile from RPNfilho and COMSOL solvers, b)
$\displaystyle pH$ profile with RPNfilho and COMSOL solvers. Figure 10:
Chloride profile from RPNfilho and COMSOL solver

In the subsequent numerical examples the structure of the solution described
above does not change. In these simulations only the relative position of the
intermediate states vary when the initial conditions of the rock (right state)
and the injection parameters (formation water, left state) change.

###  6.2 Effect of brine reduction

In the following numerical experiments, we aim at investigating the effect of
reducing the salinity on the oil recovery. In the first case, we assume a
5$\displaystyle\%$ reduction in the initial salt concentration, and the
injection and initial states are defined by $\displaystyle
L=(0.76,4.0,0.015,1.0\times 10^{-5})$ and $\displaystyle R=(0.0398,4.0,0.3)$,
respectively. The corresponding residual oil saturation is $\displaystyle
S_{or}=0.228$ (as reported in [35]), and the initial water saturation is
$\displaystyle S_{wi}=0.0398$ (calculated using formula (47)).

In the second example, we assume a 20$\displaystyle\%$ reduction in the rock
salt concentration, and the initial and injection states are defined by
$\displaystyle R=(0.045,4.0,0.3)$ and $\displaystyle
L=(0.76,4.0,0.06,1.0\times 10^{-5})$, respectively. The corresponding residual
oil saturation is $\displaystyle S_{or}=0.239$, and the initial water
saturation is $\displaystyle S_{wi}=0.045$ (calculated using formula (47)).

Figure 11 displays the profiles of water saturation $\displaystyle S_{w}$,
$\displaystyle pH$, and chloride $\displaystyle[Cl]$ concentration for the two
examples. Figure 12 shows the corresponding oil recovery for a
5$\displaystyle\%$ and 20$\displaystyle\%$ reduction in the rock salt
concentration. We observed that a 5$\displaystyle\%$ reduction in salt
concentration resulted in an increase in the oil recovery fraction of 0.066
with respect to the case of a 20$\displaystyle\%$ of salt concentration
reduction. The oil recovery values obtained with our model are consistent with
those reported in the experiments, as shown in studies such as [35] and [39].

Figure 11: The figures show profiles of the saturation $\displaystyle S_{w}$,
scaled pH ($\displaystyle pH/4$), and scaled chloride concentration
($\displaystyle[Cl]*4$) at different velocities for a 5$\displaystyle\%$ and
20$\displaystyle\%$ reduction in salt concentration (left and right figures,
respectively). The velocity of the saline and pH fronts are $\displaystyle
4.05\times 10^{-5}$ and $\displaystyle 4.08\times 10^{-5}$, respectively, for
the 5$\displaystyle\%$ reduction case, and the velocity of the saturation
front is $\displaystyle 4.34\times 10^{-5}$. For the 20$\displaystyle\%$
reduction case, the velocity of the saline and pH fronts are $\displaystyle
4.08\times 10^{-5}$, while the velocity of the saturation front is
$\displaystyle 4.44\times 10^{-5}$. Figure 12: Evaluating oil recovery with a
decrease of 5$\displaystyle\%$ and 20$\displaystyle\%$ in salt concentration.

The Riemann solutions provide a numerical quantification and explanation of
the influence of low salinity injection on the system state. In case (a) (as
shown on the left side of Figure 11), we present the solution for a
5$\displaystyle\%$ reduction in the salt concentration, while in case (b) (as
shown on the right side of Figure 11), we present the solution for a
20$\displaystyle\%$ reduction in the salt concentration. The Riemann solutions
demonstrate the patterns of the system state, illustrating the effect of low
salinity injection on the system.

In both cases (a) and (b) (as shown in Figure 11), we observe the formation of
three fronts, i.e., saturation $\displaystyle S_{w}$, $\displaystyle pH$ and
clhoride concentration $\displaystyle[Cl]$. The saline and $\displaystyle pH$
fronts move at the same velocity, while the saturation front has a slightly
higher velocity. The velocity of the saline and $\displaystyle pH$ fronts is
$\displaystyle 4.05\times 10^{-5}$ for both cases, whereas the velocity of the
saturation front is $\displaystyle 4.34\times 10^{-5}$ for case (a) and
$\displaystyle 4.44\times 10^{-5}$ for case (b).

The increase in $\displaystyle pH$ and the location of the front in
$\displaystyle pH$ close to the saturation front can be explained by the fact
that $\displaystyle pH$ controls the number of surface species at the
interfaces of the oil/brine and the brine/carbonate (see e.g., [68]). This
behavior of $\displaystyle pH$ and of the saline front can also be explained
by the fact that $\displaystyle pH$ and salinity are among the most prominent
factors affecting the wettability state of a crude oil/brine/rock system
during waterflooding operations ([52]).

Another experiment was conducted to calculate the residual oil by changing the
types of oil 1 and 2 presented in Table 1. However, no significant changes
were observed for the type of oil studied in this experiment.

We confirmed that the coupling of SCM and compositional modeling accurately
reproduces the main effects observed in the experiments, including
wettability, residual oil, and connate water saturation as the main factors
responsible for oil recovery. The processes controlling wettability and
surface complexes are the change in $\displaystyle pH$ and the formation of a
$\displaystyle pH$ front.

##  7 Conclusion

This study quantifies the TBP of surface complexes and of wettability to
estimate the effect of reducing rock salt concentration on permeability, and
consequently, on oil recovery. We use a combination of SCM and compositional
modeling to analyze the influence of wettability, residual oil, and connate
water saturation on enhanced oil recovery. Our model reveals the presence of a
shock in the saturation profile, with higher oil recovery being responsible
for this phenomenon. We also observe a $\displaystyle pH$ wave in a constant
$\displaystyle pH$ flood, salt formation at the front, and a jump in water
saturation. The Riemann solution confirms the occurrence of $\displaystyle pH$
variations numerically using COMSOL and RPNfilho solvers, which is expected
when surface complexes are formed.

## Acknowledgments

The authors are gratefull to Ali A. Eftekhari for reviewing the calculations
carried out using the PHREEQC program. Additionally, they would like to thank
Sergio Pilotto for his support and acknowledge the funding received from CAPES
under grant 88881.156518/2017-01 and CAPES/NUFFIC grant 88887.156517/2017-00,
CNPq under grants 405366/2021-3 and 306566/2019-2, and FAPERJ under grants
E-26/210.738/2014, E-26/202.764/2017, and E-26/201.159/2021.

##  8 Declarations

There are no conflicts of competing interests.

## References

  * [1] Alvarez, A. C., Bruining, J., Lambert, W. J., & Marchesin, D.  Analytical and numerical solutions for carbonated waterflooding.  Computational Geosciences, 22(2):505–526, 2018.
  * [2] Al-Shalabi, E. W., Sepehrnoori, K., Delshad, M., & Pope, G.  A novel method to model low-salinity-water injection in carbonate oil reservoirs.  SPE Journal, 20(05):1154–1166, 2015.
  * [3] Al-Shalabi, E. W., Sepehrnoori, K., & Pope, G.  Geochemical interpretation of low-salinity-water injection in carbonate oil reservoirs.  SPE Journal, 20(06):1212–1226, 2015.
  * [4] Alvarez, A. C., Bruining, J., & Marchesin, D.  Nonlinear wave interactions in geochemical modeling.  Journal of Differential Equations, 359:1–22, 2023.
  * [5] Alvarez, A. C., Goedert, G. T., & Marchesin, D.  Resonance in rarefaction and shock curves: local analysis and numerics of the continuation method.  Journal of Hyperbolic Differential Equations, 17(04):639–676, 2020\.
  * [6] Alvarez, A. C., Blom, T., Lambert, W. J., Bruining, J., & Marchesin, D.  Analytical and numerical validation of a model for flooding by saline carbonated water.  Journal of Petroleum Science and Engineering, 167:900–917,2018.
  * [7] Appelo, C Anthony J and Postma, Dieke  Geochemistry, Groundwater and Pollution.  Taylor & Francis, 2005.
  * [8] Appelo, C. A. J.  Cation and proton exchange, $\displaystyle pH$ variations, and carbonate reactions in a freshening aquifer.  Water Resources Research, 30(10):2793–2805, 1994.
  * [9] Ayirala, S. C., & Yousef, A. A.  A state-of-the-art review to develop injection-water-chemistry requirement guidelines for ior/eor projects.  SPE Production & Operations, 30(01):26–42, 2015.
  * [10] Bonto, M., Eftekhari, A. A., & Nick, H. M.  An overview of the oil-brine interfacial behavior and a new surface complexation model.  Scientific reports, 9(1):1–16, 2019.
  * [11] Bordeaux-Rego, F., Mehrabi, M., Sanaei, A., & Sepehrnoori, K.  Improvements on modelling wettability alteration by engineered water injection: Surface complexation at the oil/brine/rock contact.  Fuel, 284:118991, 2021.
  * [12] Brady, P. V., Krumhansl, J. L., & Mariner, P. E.  Surface complexation modeling for improved oil recovery.  In SPE Improved Oil Recovery Symposium, volume SPE 153744,
  * [13] Hans Bruining.  Upscaling of Single-and Two-Phase Flow in Reservoir Engineering.  CRC Press, 2021.
  * [14] Bryant, S. L., Schechter, R. S., & Lake, L. W.  Interactions of precipitation/dissolution waves and ion exchange in flow through permeable media.  AIChE journal, 32(5):751–764, 1986.
  * [15] Bryant, S. L., Schechter, R. S., & Lake, L. W.  Mineral sequences in precipitation/dissolution waves.  AIChE journal, 33(8):1271–1287, 1987.
  * [16] Buckley, S. E., & Leverett, M.  Mechanism of fluid displacement in sands.  Transactions of the AIME, 146(01):107–116, 1942.
  * [17] Christensen, R. J.  Carbonated waterflood results–Texas and Oklahoma.  In Annual Meeting of Rocky Mountain Petroleum Engineers of AIME. OnePetro, 1961.
  * [18] De Nevers, N.  A calculation method for carbonated water flooding.  Society of Petroleum Engineers Journal, 4(01):9–20, 1964.
  * [19] Dong, Y., Dindoruk, B., Ishizawa, C., Lewis, E., & Kubicek, T.  An experimental investigation of carbonated water flooding.  In SPE Annual Technical Conference and Exhibition. OnePetro,
  * [20] Dubey, S. T., & Doe, P. H.  Base number and wetting properties of crude oils.  SPE Reservoir Engineering, 8(3):195–200, 1993.
  * [21] Dumoré, J. M., Hagoort, J., & Risseeuw, A. S.  An analytical model for one-dimensional, three-component condensing and vaporizing gas drives.  Society of Petroleum Engineers Journal, 24(02):169–179, 1984.
  * [22] Ebeltoft, E., Lomeland, F., Brautaset, A., & Haugen, Å  Parameter based scal-analysing relative permeability for full field application.  In International Symposium of the Society of Core Analysis, Avignon, France, pages 8–11, 2014.
  * [23] Elakneswaran, Y., Shimokawara, M., Nawa, T., & Takahashi, S.  Surface complexation and equilibrium modelling for low salinity waterflooding in sandstone reservoirs.  In Abu Dhabi International Petroleum Exhibition & Conference. OnePetro, 2017.
  * [24] Erzuah, S., Fjelde, I., & Omekeh, A. V.  Wettability estimation using surface-complexation simulations.  SPE Reservoir Evaluation & Engineering, 22(02):509–519, 2019.
  * [25] Erzuah, S., Fjelde, I., & Omekeh, A. V.  Modelling low-salinity waterflooding: Effect of divalent cations and capillary pressure.  Journal of Petroleum Science and Engineering, 149:1–8, 2017.
  * [26] Farajzadeh, R., Matsuura, T., van Batenburg, D., & Dijk, H.  Detailed modeling of the alkali/surfactant/polymer (ASP) process by coupling a multipurpose reservoir simulator to the chemistry package PHREEQC.  SPE Reservoir Evaluation & Engineering, 15(04):423–435, 2012.
  * [27] Foroozesh, J., Jamiolahmady, M., & Sohrabi, M.  Mathematical modeling of carbonated water injection for EOR and CO2 storage with a focus on mass transfer kinetics.  Fuel, 174:325–332, 2016.
  * [28] Ginn, Dennis.  Effects of potential determining ions and $\displaystyle pH$ on the wettability of intermediate wet outcrop limestone.  Master’s thesis, University of Stavanger, Norway, 2020.
  * [29] Glimm, J.  Solutions in the large for nonlinear hyperbolic systems of equations.  Communications on Pure and Applied Mathematics, 18(4):697–715, 1965\.
  * [30] Grogan, A. T., & Pinczewski, W. V.  The role of molecular diffusion processes in tertiary CO2 flooding.  Journal of Petroleum Technology, 39(05):591–602, 1987.
  * [31] Hassan, A. M., Ayoub, M., Eissa, M., Bruining, H., & Zitha, P.  Study of surface complexation modeling on a novel hybrid enhanced oil recovery (eor) method; smart-water assisted foam-flooding.  Journal of Petroleum Science and Engineering, 195:107563, 2020.
  * [32] Helfferich, F. G.  The theory of precipitation/dissolution waves.  AIChE journal, 35(1):75–87, 1989.
  * [33] Hirasaki, G., & Zhang, D. L.  Surface chemistry of oil recovery from fractured, oil-wet, carbonate formation.  In International Symposium on Oilfield Chemistry. OnePetro, 2003\.
  * [34] Honarpour, Mehdi Matt.  Relative Permeability of Petroleum Reservoirs.  CRC press, 1986.
  * [35] Jerauld, G. R., Webb, K. J., Lin, C. Y., & Seccombe, J. C.  Modeling low-salinity waterflooding.  SPE Reservoir Evaluation & Engineering, 11(06):1000–1012, 2008\.
  * [36] Jerauld, G. R., Lin, C. Y., Webb, K. J., & Seccombe, J. C.  Modeling low-salinity waterflooding.  SPE Reservoir Evaluation & Engineering, 11(06):1000–1012, 2008.
  * [37] Johns, R. T., Dindoruk, B., & Orr Jr, F. M.  Analytical theory of combined condensing/vaporizing gas drives.  SPE Advanced Technology Series, 1(02):7–16, 1993.
  * [38] Korrani, A. K., Jerauld, G. R., & Sepehrnoori, K.  Mechanistic modeling of low-salinity waterflooding through coupling a geochemical package with a compositional reservoir simulator.  SPE Reservoir Evaluation & Engineering, 19(01):142–162, 2016.
  * [39] Korrani, A. K., & Jerauld, G. R.  Modeling wettability change in sandstones and carbonates using a surface-complexation-based method.  Journal of Petroleum Science and Engineering, 174:1093–1112,2019.
  * [40] Kumar, S., & Mandal, A.  A comprehensive review on chemically enhanced water alternating gas/CO2 (CEWAG) injection for enhanced oil recovery.  Journal of Petroleum Science and Engineering, 157:696–715,2017.
  * [41] Lake, Larry W.  Enhanced Oil Recovery.  Prentice Hall Inc.  1989\.
  * [42] L.W. Lake, R.T. Johns, W.R. Rossen, and G.A. Pope.  Fundamentals of Enhanced Oil Recovery.  Society of Petroleum Engineers, Richardson, TX,  2014\.
  * [43] Lambert, W. J., Alvarez, A. C., Matos, V., Marchesin, D., & Bruining, J.  Nonlinear wave analysis of geochemical injection for multicomponent two phase flow in porous media.  Journal of Differential Equations, 266(1):406–454, 2019.
  * [44] Lax, P. D.  Hyperbolic systems of conservation laws II.  Communications on Pure and Applied Mathematics, 10(4):537–566,
  * [45] Leverett, M. C.  Flow of oil-water mixtures through unconsolidated sands.  Transactions of the AIME, 132(01):149–171, 1939.
  * [46] Liu, T. P.  The Riemann problem for general 2$\displaystyle\times$ 2 conservation laws.  Transactions of the American Mathematical Society, 199:89–112,
  * [47] Liu, T. P.  The Riemann problem for general systems of conservation laws.  Journal of Differential Equations, 18(1):218–234, 1975.
  * [48] Lomeland, F., Ebeltoft, E., & Hasanov, B.  A versatile representation of upscaled relative permeability for field applications (spe 154487).  In 74th EAGE Conference and Exhibition incorporating EUROPEC 2012, pages cp–293. European Association of Geoscientists $\displaystyle\&$ Engineers, 2012.
  * [49] Johannes Lutzenkirchen.  Surface Complexation Modelling.  Elsevier, 2006.
  * [50] Marmier, N., Dumonceau, J., & Fromage, F.  Surface Complexation Modeling of yb (iii) sorption and desorption on hematite and alumina.  Journal of Contaminant Hydrology, 26(1-4):159–167, 1997.
  * [51] Mehdiyev, F., Erzuah, S., Omekeh, A., & Fjelde, I.  Surface Complexation Modelling of wettability alteration during carbonated water flooding.  Energies, 15(9):3020, 2022.
  * [52] Mehraban, M. F., Ayatollahi, S., & Sharifi, M.  Experimental investigation on synergic effect of salinity and ph during low salinity water injection into carbonate oil reservoirs.  Journal of Petroleum Science and Engineering, 202:108555, 2021.
  * [53] Broder J. Merkel.  Groundwater Geochemistry.  Springer, 2005.
  * [54] Nowrouzi, I., Manshad, A. K., & Mohammadi, A. H.  Effects of dissolved binary ionic compounds and different densities of brine on interfacial tension (ift), wettability alteration, and contact angle in smart water and carbonated smart water injection processes in carbonate oil reservoirs.  Journal of Molecular Liquids, 254:83–92, 2018.
  * [55] Oleinik, O. A.  Discontinuous solutions of non-linear differential equations.  Uspekhi Matematicheskikh Nauk, 12(3):3–73, 1957.
  * [56] Omekeh, A., Friis, H. A., Fjelde, I., & Evje, S.  Modeling of ion-exchange and solubility in low salinity water flooding.  In SPE Improved Oil Recovery Symposium. OnePetro, 2012.
  * [57] David L Parkhurst and CAJ Appelo.  User’s guide to PHREEQC (version 2): A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations.  1999\.
  * [58] David L Parkhurst and CAJ Appelo.  Description of input and examples for PHREEQC version 3a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations.  2013\.
  * [59] Pope, G. A.  The application of fractional flow theory to enhanced oil recovery.  Society of Petroleum Engineers Journal, 20(03):191–205, 1980.
  * [60] Sanaei, A., Tavassoli, S., & Sepehrnoori, K.  Investigation of modified water chemistry for improved oil recovery: Application of DLVO theory and surface complexation model.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 574:131–145, 2019.
  * [61] Sanaei, A., Varavei, A., & Sepehrnoori, K.  Mechanistic modeling of carbonated waterflooding.  Journal of Petroleum Science and Engineering, 178:863–877,
  * [62] Sari, A., Chen, Y., Xie, Q., & Saeedi, A.  Low salinity water flooding in high acidic oil reservoirs: Impact of $\displaystyle pH$ on wettability of carbonate reservoirs.  Journal of Molecular Liquids, 281:444–450, 2019.
  * [63] James J. Sheng.  Enhanced oil recovery field case studies.  Gulf Professional Publishing, 2013.
  * [64] Sohrabi, M., Riazi, M., Jamiolahmady, M., Ireland, S., & Brown, C.  Mechanisms of oil recovery by carbonated water injection.  In SCA annual meeting, 2009.
  * [65] Werner Stumm and James J. Morgan.  Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters, volume 126.  John Wiley & Sons, 2012.
  * [66] Welge, H. J., Johnson, E. F., Ewing Jr, S. P., & Brinkman, F. H.  The linear displacement of oil from porous media by enriched gas.  Journal of Petroleum Technology, 13(08):787–796, 1961.
  * [67] Wolthers, M., Charlet, L., & Van Cappellen, P.  The surface chemistry of divalent metal carbonate minerals; a critical assessment of surface charge and potential data using the charge distribution multi-site ion complexation model.  American Journal of Science, 308(8):905–941, 2008.
  * [68] Xie, Q., Sari, A., Pu, W., Chen, Y., Brady, P. V., Al Maskari, N., & Saeedi, A.  pH effect on wettability of oil/brine/carbonate system: Implications for low salinity water flooding.  Journal of Petroleum Science and Engineering, 168:419–425, 2018.
  * [69] Yousef, A. A., Al-Saleh, S., & Al-Jawfi, M.  Smart waterflooding for carbonate reservoirs: Salinity and role of ions.  In SPE Middle East Oil and Gas Show and Conference. OnePetro, 2011.
  * [70] Yousef, A. A., Al-Saleh, S., & Al-Jawfi, M.  Improved/enhanced oil recovery from carbonate reservoirs by tuning injection water salinity and ionic content.  In SPE Improved Oil Recovery Symposium. OnePetro, 2012.